text:
#pragma once

/***************************************************
 *************** Auto-Generated File ***************
 ***************************************************/

#include <cstring>
#include <string>

#include <gsl/gsl_const_mksa.h>
#include <gsl/gsl_const_num.h>

#include "../../Utils/BinarySearch.h"

enum UnitType {
	NONE,
	DISTANCE,
	TIME,
	TIMESQUARED,
	SPEED,
	ACCELERATION,
	AREA,
	VOLUME,
	MASS,
	FORCE,
	MASSDISTANCE,
	MASSPERTIMESQUARED,
	ENERGY,
	MASSAREA
};

constexpr const int numUnits = 161;
constexpr const char* unitNames[numUnits] = {
	"acre", "astronomical_unit", "attogram", "attometer", "attometerscubed",
	"attometerssquared", "attosecond", "barn", "btu", "calorie", "canadian_gallon", "carat",
	"centimeter", "century", "cup", "decade", "decameter", "decimeter", "dyne",
	"electron_volt", "erg", "exagram", "exameter", "exameterscubed", "exameterssquared",
	"exasecond", "fathom", "femtogram", "femtometer", "femtometerscubed",
	"femtometerssquared", "femtosecond", "fluid_ounce", "foot", "gigagram", "gigameter",
	"gigameterscubed", "gigameterssquared", "gigasecond", "gram", "gram_force", "grammeter",
	"grammeterspersecondsquared", "grammeterssquared", "grammeterssquaredpersecondsquared",
	"grampersecondsquared", "grav_accel", "hectare", "hectometer", "hour", "inch", "joule",
	"kilogram", "kilometer", "kilometers_per_hour", "kilometerscubed", "kilometerssquared",
	"kilopound_force", "kilosecond", "knot", "light_year", "liter", "mass_electron",
	"mass_muon", "mass_neutron", "mass_proton", "megagram", "megameter", "megameterscubed",
	"megameterssquared", "megasecond", "meter", "meterscubed", "meterspersecond",
	"meterspersecondsquared", "meterssquared", "metric_ton", "microgram", "microliter",
	"micrometer", "micrometerscubed", "micrometerssquared", "microsecond", "mil", "mile",
	"miles_per_hour", "millenium", "milligram", "milliliter", "millimeter",
	"millimeterscubed", "millimeterssquared", "millisecond", "minute", "month", "nanogram",
	"nanometer", "nanometerscubed", "nanometerssquared", "nanosecond", "nautical_mile",
	"newton", "newton_meter", "ounce_mass", "petagram", "petameter", "petameterscubed",
	"petameterssquared", "petasecond", "picogram", "picometer", "picometerscubed",
	"picometerssquared", "picosecond", "pint", "pound_force", "pound_mass", "poundal",
	"quart", "rydberg", "second", "secondsquared", "solar_mass", "speed_of_light",
	"tablespoon", "teaspoon", "teragram", "terameter", "terameterscubed",
	"terameterssquared", "terasecond", "therm", "ton", "troy_ounce", "uk_gallon", "uk_ton",
	"unified_atomic_mass", "us_gallon", "week", "yard", "year", "yoctogram", "yoctometer",
	"yoctometerscubed", "yoctometerssquared", "yoctosecond", "yottagram", "yottameter",
	"yottameterscubed", "yottameterssquared", "yottasecond", "zeptogram", "zeptometer",
	"zeptometerscubed", "zeptometerssquared", "zeptosecond", "zettagram", "zettameter",
	"zettameterscubed", "zettameterssquared", "zettasecond"
};
constexpr const int unitIndices[numUnits] = {
	52, 0, 53, 54, 56, 55, 57, 58, 59, 61, 2, 1, 63, 62, 64, 66, 65, 67, 3, 4, 68, 5, 6, 8,
	7, 9, 76, 69, 71, 73, 72, 74, 70, 75, 11, 12, 14, 13, 15, 77, 83, 78, 79, 80, 81, 82,
	10, 85, 86, 84, 87, 16, 88, 90, 94, 92, 91, 89, 95, 93, 98, 17, 102, 103, 104, 105, 19,
	20, 22, 21, 23, 99, 107, 100, 101, 106, 118, 142, 143, 144, 146, 145, 147, 110, 109,
	116, 112, 108, 111, 113, 115, 114, 117, 18, 119, 120, 121, 123, 122, 125, 124, 24, 25,
	126, 26, 27, 29, 28, 30, 128, 129, 131, 130, 132, 133, 97, 96, 127, 134, 31, 135, 136,
	32, 60, 138, 141, 33, 34, 36, 35, 37, 139, 137, 140, 39, 40, 38, 41, 148, 149, 154, 150,
	151, 153, 152, 155, 42, 43, 45, 44, 46, 156, 157, 159, 158, 160, 47, 48, 50, 49, 51
};
constexpr const char * unitAbbreviations[numUnits] = {
	"AU", "CD", "CG", "Dy", "EV", "Eg", "Em", "Em^2", "Em^3", "Es", "G", "Gg", "Gm", "Gm^2",
	"Gm^3", "Gs", "J", "L", "M", "Mg", "Mm", "Mm^2", "Mm^3", "Ms", "N", "Nm", "Pg", "Pm",
	"Pm^2", "Pm^3", "Ps", "Ry", "SM", "Tg", "Tm", "Tm^2", "Tm^3", "Ts", "UAM", "UKG", "UKT",
	"USG", "Yg", "Ym", "Ym^2", "Ym^3", "Ys", "Zg", "Zm", "Zm^2", "Zm^3", "Zs", "ac", "ag",
	"am", "am^2", "am^3", "as", "bn", "btu", "c", "cal", "cen", "cm", "cp", "dcm", "dec",
	"dm", "erg", "fg", "floz", "fm", "fm^2", "fm^3", "fs", "ft", "fth", "g", "g*m",
	"g*m/s^2", "g*m^2", "g*m^2/s^2", "g/s^2", "gF", "h", "hec", "hm", "in", "kg", "klbF",
	"km", "km^2", "km^3", "knt", "kph", "ks", "lb", "lbF", "ly", "m", "m/s", "m/s^2", "mE",
	"mMu", "mNt", "mPt", "m^2", "m^3", "mg", "mi", "mil", "ml", "mln", "mm", "mm^2", "mm^3",
	"mph", "ms", "mt", "mth", "ng", "nm", "nm^2", "nm^3", "nmi", "ns", "oz", "pdl", "pg",
	"pm", "pm^2", "pm^3", "ps", "pt", "qt", "s", "s^2", "t", "tbsp", "thm", "toz", "tsp",
	"ug", "ul", "um", "um^2", "um^3", "us", "wk", "yd", "yg", "ym", "ym^2", "ym^3", "yr",
	"ys", "zg", "zm", "zm^2", "zm^3", "zs"
};
constexpr const UnitType unitTypes[numUnits] = {
	DISTANCE, MASS, VOLUME, FORCE, ENERGY, MASS, DISTANCE, AREA, VOLUME, TIME, ACCELERATION,
	MASS, DISTANCE, AREA, VOLUME, TIME, ENERGY, VOLUME, TIME, MASS, DISTANCE, AREA, VOLUME,
	TIME, FORCE, ENERGY, MASS, DISTANCE, AREA, VOLUME, TIME, ENERGY, MASS, MASS, DISTANCE,
	AREA, VOLUME, TIME, MASS, VOLUME, MASS, VOLUME, MASS, DISTANCE, AREA, VOLUME, TIME,
	MASS, DISTANCE, AREA, VOLUME, TIME, AREA, MASS, DISTANCE, AREA, VOLUME, TIME, AREA,
	ENERGY, SPEED, ENERGY, TIME, DISTANCE, VOLUME, DISTANCE, TIME, DISTANCE, ENERGY, MASS,
	VOLUME, DISTANCE, AREA, VOLUME, TIME, DISTANCE, DISTANCE, MASS, MASSDISTANCE, FORCE,
	MASSAREA, ENERGY, MASSPERTIMESQUARED, FORCE, TIME, AREA, DISTANCE, DISTANCE, MASS,
	FORCE, DISTANCE, AREA, VOLUME, SPEED, SPEED, TIME, MASS, FORCE, DISTANCE, DISTANCE,
	SPEED, ACCELERATION, MASS, MASS, MASS, MASS, AREA, VOLUME, MASS, DISTANCE, DISTANCE,
	VOLUME, TIME, DISTANCE, AREA, VOLUME, SPEED, TIME, MASS, TIME, MASS, DISTANCE, AREA,
	VOLUME, DISTANCE, TIME, MASS, FORCE, MASS, DISTANCE, AREA, VOLUME, TIME, VOLUME, VOLUME,
	TIME, TIMESQUARED, MASS, VOLUME, ENERGY, MASS, VOLUME, MASS, VOLUME, DISTANCE, AREA,
	VOLUME, TIME, TIME, DISTANCE, MASS, DISTANCE, AREA, VOLUME, TIME, TIME, MASS, DISTANCE,
	AREA, VOLUME, TIME
};
constexpr const double unitConversions[numUnits] = {
	GSL_CONST_MKSA_ASTRONOMICAL_UNIT, GSL_CONST_MKSA_CARAT*1000.0,
	GSL_CONST_MKSA_CANADIAN_GALLON, GSL_CONST_MKSA_DYNE*1000.0,
	GSL_CONST_MKSA_ELECTRON_VOLT*1000.0, GSL_CONST_NUM_EXA, GSL_CONST_NUM_EXA,
	GSL_CONST_NUM_EXA, GSL_CONST_NUM_EXA, GSL_CONST_NUM_EXA, GSL_CONST_MKSA_GRAV_ACCEL,
	GSL_CONST_NUM_GIGA, GSL_CONST_NUM_GIGA, GSL_CONST_NUM_GIGA, GSL_CONST_NUM_GIGA,
	GSL_CONST_NUM_GIGA, GSL_CONST_MKSA_JOULE*1000.0, GSL_CONST_MKSA_LITER,
	GSL_CONST_MKSA_MINUTE, GSL_CONST_NUM_MEGA, GSL_CONST_NUM_MEGA, GSL_CONST_NUM_MEGA,
	GSL_CONST_NUM_MEGA, GSL_CONST_NUM_MEGA, GSL_CONST_MKSA_NEWTON*1000.0, 1000,
	GSL_CONST_NUM_PETA, GSL_CONST_NUM_PETA, GSL_CONST_NUM_PETA, GSL_CONST_NUM_PETA,
	GSL_CONST_NUM_PETA, GSL_CONST_MKSA_RYDBERG*1000.0, GSL_CONST_MKSA_SOLAR_MASS*1000.0,
	GSL_CONST_NUM_TERA, GSL_CONST_NUM_TERA, GSL_CONST_NUM_TERA, GSL_CONST_NUM_TERA,
	GSL_CONST_NUM_TERA, GSL_CONST_MKSA_UNIFIED_ATOMIC_MASS*1000.0, GSL_CONST_MKSA_UK_GALLON,
	GSL_CONST_MKSA_UK_TON*1000.0, GSL_CONST_MKSA_US_GALLON, GSL_CONST_NUM_YOTTA,
	GSL_CONST_NUM_YOTTA, GSL_CONST_NUM_YOTTA, GSL_CONST_NUM_YOTTA, GSL_CONST_NUM_YOTTA,
	GSL_CONST_NUM_ZETTA, GSL_CONST_NUM_ZETTA, GSL_CONST_NUM_ZETTA, GSL_CONST_NUM_ZETTA,
	GSL_CONST_NUM_ZETTA, GSL_CONST_MKSA_ACRE, GSL_CONST_NUM_ATTO, GSL_CONST_NUM_ATTO,
	GSL_CONST_NUM_ATTO, GSL_CONST_NUM_ATTO, GSL_CONST_NUM_ATTO, GSL_CONST_MKSA_BARN,
	GSL_CONST_MKSA_BTU*1000.0, GSL_CONST_MKSA_SPEED_OF_LIGHT, GSL_CONST_MKSA_CALORIE*1000.0,
	3153600000, 1e-2, GSL_CONST_MKSA_CUP, 1e1, 315360000, 1e-1, GSL_CONST_MKSA_ERG*1000.0,
	GSL_CONST_NUM_FEMTO, GSL_CONST_MKSA_FLUID_OUNCE, GSL_CONST_NUM_FEMTO,
	GSL_CONST_NUM_FEMTO, GSL_CONST_NUM_FEMTO, GSL_CONST_NUM_FEMTO, GSL_CONST_MKSA_FOOT,
	GSL_CONST_MKSA_FATHOM, 1, 1, 1, 1, 1, 1, GSL_CONST_MKSA_GRAM_FORCE*1000.0,
	GSL_CONST_MKSA_HOUR, GSL_CONST_MKSA_HECTARE, 1e2, GSL_CONST_MKSA_INCH,
	GSL_CONST_NUM_KILO, GSL_CONST_MKSA_KILOPOUND_FORCE*1000.0, GSL_CONST_NUM_KILO,
	GSL_CONST_NUM_KILO, GSL_CONST_NUM_KILO, GSL_CONST_MKSA_KNOT,
	GSL_CONST_MKSA_KILOMETERS_PER_HOUR, GSL_CONST_NUM_KILO,
	GSL_CONST_MKSA_POUND_MASS*1000.0, GSL_CONST_MKSA_POUND_FORCE*1000.0,
	GSL_CONST_MKSA_LIGHT_YEAR, 1, 1, 1, GSL_CONST_MKSA_MASS_ELECTRON*1000.0,
	GSL_CONST_MKSA_MASS_MUON*1000.0, GSL_CONST_MKSA_MASS_NEUTRON*1000.0,
	GSL_CONST_MKSA_MASS_PROTON*1000.0, 1, 1, GSL_CONST_NUM_MILLI, GSL_CONST_MKSA_MILE,
	GSL_CONST_MKSA_MIL, 1e-6, 31536000000, GSL_CONST_NUM_MILLI, GSL_CONST_NUM_MILLI,
	GSL_CONST_NUM_MILLI, GSL_CONST_MKSA_MILES_PER_HOUR, GSL_CONST_NUM_MILLI,
	GSL_CONST_MKSA_METRIC_TON*1000.0, 2628288, GSL_CONST_NUM_NANO, GSL_CONST_NUM_NANO,
	GSL_CONST_NUM_NANO, GSL_CONST_NUM_NANO, GSL_CONST_MKSA_NAUTICAL_MILE,
	GSL_CONST_NUM_NANO, GSL_CONST_MKSA_OUNCE_MASS*1000.0, GSL_CONST_MKSA_POUNDAL*1000.0,
	GSL_CONST_NUM_PICO, GSL_CONST_NUM_PICO, GSL_CONST_NUM_PICO, GSL_CONST_NUM_PICO,
	GSL_CONST_NUM_PICO, GSL_CONST_MKSA_PINT, GSL_CONST_MKSA_QUART, 1, 1,
	GSL_CONST_MKSA_TON*1000.0, GSL_CONST_MKSA_TABLESPOON, GSL_CONST_MKSA_THERM*1000.0,
	GSL_CONST_MKSA_TROY_OUNCE*1000.0, GSL_CONST_MKSA_TEASPOON, GSL_CONST_NUM_MICRO, 1e-9,
	GSL_CONST_NUM_MICRO, GSL_CONST_NUM_MICRO, GSL_CONST_NUM_MICRO, GSL_CONST_NUM_MICRO,
	GSL_CONST_MKSA_WEEK, GSL_CONST_MKSA_YARD, GSL_CONST_NUM_YOCTO, GSL_CONST_NUM_YOCTO,
	GSL_CONST_NUM_YOCTO, GSL_CONST_NUM_YOCTO, 31536000, GSL_CONST_NUM_YOCTO,
	GSL_CONST_NUM_ZEPTO, GSL_CONST_NUM_ZEPTO, GSL_CONST_NUM_ZEPTO, GSL_CONST_NUM_ZEPTO,
	GSL_CONST_NUM_ZEPTO
};

/*
Returns index of the unit in the unitNames array.

Uses binary search under the hood to search for the index.

Parameters
----------
name: The name of the unit

Returns
-------
The index or -1 if the provided name is not a unit.
*/
CONSTEXPR_BINARY_SEARCH(getUnitIndex, unitNames, numUnits)

CONSTEXPR_BINARY_SEARCH(getAbbrIndex, unitAbbreviations, numUnits)

constexpr UnitType getUnitType(const char* name){
    int index = getUnitIndex(name);
    if (index != -1){
		return NONE;
	}
	return unitTypes[index];
}
inline UnitType getUnitType(const std::string& name){
    int index = getUnitIndex(name.c_str());
    if (index != -1){
		return NONE;
	}
	return unitTypes[index];
}
inline UnitType getUnitType(int index){
    if (index < 0 || index >= numUnits){
		return NONE;
	}
	return unitTypes[index];
}

constexpr double getUnitConversion(const char* name){
    int index = getUnitIndex(name);
    if (index != -1){
		return NONE;
	}
	return unitConversions[index];
}
inline double getUnitConversion(const std::string& name){
    int index = getUnitIndex(name.c_str());
    if (index != -1){
		return NONE;
	}
	return unitConversions[index];
}
inline double getUnitConversion(int index){
    if (index < 0 || index >= numUnits){
		return NONE;
	}
	return unitConversions[index];
}

ext: h

max_stars_repo_name: antoniojkim/CalcPlusPlus

max_issues_repo_name: antoniojkim/CalcPlusPlus

max_forks_repo_name: antoniojkim/CalcPlusPlus


text:
// @Copyright 2007-2017 Kristjan Haule
#include <iostream>
#include <fstream>
#include <vector>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>
#include "logging.h"
#include "smesh.h"
#include "sfunction.h"
#include "tanmesh.h"
#include "romberg.h"
#include <time.h>

#ifdef _MPI
#include <mpi.h>
#endif

using namespace std;

inline double fermi_kernel(double t, double w, double beta)
{
  double x = beta*w/2;
  double y = 2.*t/beta-1.;
  if (x>100.) return exp(-x*(y+1.));
  if (x<-100) return exp(x*(1.-y));
  return exp(-x*y)/(2*cosh(x));
}
inline double bose_kernel(double t, double w, double beta)
{
  double x = beta*w/2;
  double y = 2.*t/beta-1.;
  if (x>200.) return w*exp(-x*(y+1.));
  if (x<-200.) return -w*exp(x*(1.-y));
  return w*exp(-x*y)/(2*sinh(x));
}
inline void create_log_mesh(function1D<int>& ind_om, int nom_all, int nom, int ntail_)
{
  //Creates logarithmic mesh on Matsubara axis
  //     Takes first istart points from mesh om and the rest of om mesh is replaced by ntail poinst redistribued logarithmically.
  //     Input:
  //         om      -- original long mesh
  //         nom     -- number of points not in the tail
  //         ntail   -- tail replaced by ntail points only
  //     Output:
  //         ind_om  -- index array which conatins index to kept Matsubara points
  int istart = min(nom, nom_all);
  int ntail = min(ntail_, nom_all-istart);
  ind_om.resize(istart+ntail);
  double alpha = log((nom_all-1.)/istart)/(ntail-1.);
  for (int i=0; i< istart; i++) ind_om[i]=i;
  int ii=istart;
  for (int i=0; i < ntail; i++){
    int t = int(istart*exp(alpha*i)+0.5);
    if (t != ind_om[ii-1])
      ind_om[ii++] = t;
  }
  ind_om.resize(ii);
}

class SVDFunc{
public:
  vector<spline1D<double> > fU, fU_bose; // SVD functions in imaginary time
  mesh1D tau;                             // imaginary time optimized mesh
  int lmax, lmax_bose;                   // cutoff for the singular values
  int k_ovlp;                             // 2^{k_ovlp}+1 is the number of integration points in romberg method
  mesh1D om;                              // real axis mesh for SVD
  function1D<double> S, S_bose;          // eigenvalues of SVD decomposition
  function2D<double> Vt, Vt_bose;        // real axis basis functions
  function1D<double> dxi;                 // parts of fU splines stored in more efficient way for fast interpolation
  function3D<double> ff2;                 // fU splines stored in different way for fast interpolation
  double beta, L, x0;                     // parameters of the mesh should be remembered for reproducibility
  //
  SVDFunc() : lmax(0)  {  };
  //
  bool operator()() //  was it ever initialized?
  { return lmax>0; }
  double operator()(int l, int i){
    if (l<0) l=lmax-l;
    if (i<0) i=tau.size()+i;
    return fU[l][i];
  }
  double operator()(int l, int i) const{
    if (l<0) l=lmax-l;
    if (i<0) i=tau.size()+i;
    return fU[l][i];
  }
  spline1D<double>& operator[](int l){
    return fU[l];
  }

  inline vector<int> get_lll(int ii)
  { // conversion from a single combined index to the three l indices.
    //int ii = l0 + lmax*l1 + lmax*lmax*lb;
    int lb_ = ii / (lmax*lmax);
    int ir = ii % (lmax*lmax);
    int l1_ = ir / lmax;
    int l0_ = ir % lmax;
    int tmp[] = { l0_, l1_, lb_ };
    return vector<int>(tmp, tmp+3 );
  }
  void _cmp_(double beta, int& lmax, vector<spline1D<double> >& fU, function1D<double>& S, function2D<double>& Vt, function2D<double>& K, ostream& clog){
    clock_t t0 = clock();
    int Nl = min(om.size(),tau.size());
    S.resize(Nl);
    Vt.resize(om.size(),Nl);
    function2D<double> U(Nl,tau.size());

    int lwork;
    {
      int n = min(tau.size(),om.size());
      int mx = max(tau.size(),om.size());
      lwork =  4*n*n + 8*max(mx, n*(n+1));
    }
    function1D<double> work(lwork);
    function1D<int> iwork(8*min(om.size(),tau.size()));
    int info=0;
    info = xgesdd(true, tau.size(), om.size(), K.MemPt(), K.fullsize_Nd(), S.MemPt(), U.MemPt(), U.fullsize_Nd(), Vt.MemPt(), Vt.fullsize_Nd(),work.MemPt(), lwork, work.MemPt(), iwork.MemPt());
    if (info!=0){
      clog << "svd is returning "<<info<<" in svdfunc"<<endl;
    }
    double dt = static_cast<double>( clock() - t0 )/CLOCKS_PER_SEC;
    clog<<"svd time="<<dt<<endl;

    for (int l=0; l<lmax; l++)
      if (fabs(S[l])<1e-13){
	lmax=l;
	clog<<"lmax reduced to "<<lmax<<endl;
	break;
      }
    clog<<"last singular value="<<S[lmax-1]<<endl;

    for (int i=0; i<tau.size(); i++)
      for (int l=0; l<Nl; l++)
	U(l,i) *= 1./sqrt(tau.Dh(i));

    // Splines functions, which are singular values of the Kernel.
    vector<spline1D<double> > fu(lmax);
    for (int l=0; l<lmax; l++){
      fu[l].resize(tau.size());
      for (int i=0; i<tau.size(); i++) fu[l][i] = U(l,i);
      int n = tau.size();
      double x1 = 0.5*(tau[1]+tau[0]);
      double df1 = (U(l,1)-U(l,0))/(tau[1]-tau[0]);
      double x2 = 0.5*(tau[2]+tau[1]);
      double df2 = (U(l,2)-U(l,1))/(tau[2]-tau[1]);
      double df0 = df1 + (df2-df1)*(0.0-x1)/(x2-x1);
      x1 = 0.5*(tau[n-1]+tau[n-2]);
      df1 = (U(l,n-1)-U(l,n-2))/(tau[n-1]-tau[n-2]);
      x2 = 0.5*(tau[n-2]+tau[n-3]);
      df2 = (U(l,n-2)-U(l,n-3))/(tau[n-2]-tau[n-3]);
      double dfn = df1 + (df2-df1)*(beta-x1)/(x2-x1);
      fu[l].splineIt(tau, df0, dfn);
    }

    // Calculates overlap between spline interpolations
    function2D<double> overlap(lmax,lmax);
    cmpOverlap(overlap, fu, tau, beta, k_ovlp);
    // Calculates eigensystem of the overlap
    function1D<double> oE(lmax);
    SymEigensystem(overlap, oE);
    function2D<double> sou(overlap);// sou = 1/sqrt(overlap)
    // Computes  1/sqrt(overlap)
    for (int i=0; i<lmax; i++){
      for (int j=0; j<lmax; j++){
	double dsum =0.0;
	for (int k=0; k<lmax; k++) dsum += overlap(k,i)*1/sqrt(oE[k])*overlap(k,j);
	sou(i,j) = dsum;
      }
    }
    // Prepares new spline functions, which are more precisely orthonormal.
    fU.resize(lmax);
    for (int l=0; l<lmax; l++){
      fU[l].resize(tau.size());
      for (int i=0; i<tau.size(); i++){
	double ul=0;
	for (int lp=0; lp<lmax; lp++)  ul += sou(l,lp)*fu[lp][i];
	fU[l][i] = ul;
      }
      int n = tau.size();
      double x1 = 0.5*(tau[1]+tau[0]);
      double df1 = (fU[l][1]-fU[l][0])/(tau[1]-tau[0]);
      double x2 = 0.5*(tau[2]+tau[1]);
      double df2 = (fU[l][2]-fU[l][1])/(tau[2]-tau[1]);
      double df0 = df1 + (df2-df1)*(0.0-x1)/(x2-x1);
      x1 = 0.5*(tau[n-1]+tau[n-2]);
      df1 = (fU[l][n-1]-fU[l][n-2])/(tau[n-1]-tau[n-2]);
      x2 = 0.5*(tau[n-2]+tau[n-3]);
      df2 = (fU[l][n-2]-fU[l][n-3])/(tau[n-2]-tau[n-3]);
      double dfn = df1 + (df2-df1)*(beta-x1)/(x2-x1);
      fU[l].splineIt(tau, df0, dfn);
    }
  }
  //
  void cmp(const string& statistics, double beta_, int lmax_, int Ntau, double x0_, double L_, double Nw, ostream& clog, int k_ovlp_=5){
    k_ovlp=k_ovlp_;
    beta = beta_;
    L = L_;
    x0 = x0_;

    GiveTanMesh(om, x0, L, Nw);
    GiveDoubleExpMesh(tau, beta, Ntau);

    if (statistics=="fermi" || statistics=="both"){
      lmax = lmax_;
      clock_t t0 = clock();
      function2D<double> K(om.size(),tau.size());
      for (int j=0; j<om.size(); j++){
	for (int i=0; i<tau.size(); i++){
	  K(j,i) = fermi_kernel(tau[i],om[j],beta);
	  K(j,i) *= om.Dh(j)*sqrt(tau.Dh(i));
	}
      }
      double dt = static_cast<double>( clock() - t0 )/CLOCKS_PER_SEC;
      clog<<"setup time="<<dt<<endl;

      _cmp_(beta, lmax, fU, S, Vt, K, clog);
    }

    if (statistics=="bose" || statistics=="both"){
      lmax_bose = lmax_;
      clock_t t0 = clock();
      function2D<double> K(om.size(),tau.size());
      for (int j=0; j<om.size(); j++){
	for (int i=0; i<tau.size(); i++){
	  K(j,i) = bose_kernel(tau[i],om[j],beta);
	  K(j,i) *= om.Dh(j)*sqrt(tau.Dh(i));
	}
      }
      double dt = static_cast<double>( clock() - t0 )/CLOCKS_PER_SEC;
      clog<<"setup time="<<dt<<endl;

      _cmp_(beta, lmax_bose, fU_bose, S_bose, Vt_bose, K, clog);
      //cout<<"S_bose="<<S_bose<<endl;
    }
  }

  void Print(const string& filename, const string& statistics="fermi"){
    ofstream pul(filename.c_str());
    pul<<"# "<<lmax<<" "<<tau.size()<<" "<<L<<" "<<x0<<"  "<<(om.size()/2)<<endl;
    pul.precision(16);
    for (int i=0; i<tau.size(); i++){
      pul<<tau[i]<<"  ";
      if (statistics=="fermi")
	for (int l=0; l<lmax; l++) pul<<fU[l][i]<<"  ";
      else
	for (int l=0; l<lmax_bose; l++) pul<<fU_bose[l][i]<<"  ";
      pul<<endl;
    }
  }

  double CheckOverlap(/*double beta*/){
    // Now we check how well is the orthogonality obeyed after reorthogonalization
    function2D<double> overlap(lmax,lmax);
    cmpOverlap(overlap, fU, tau, beta, k_ovlp);
    double dsum=0;
    for (int l1=0; l1<lmax; l1++){
      for (int l2=0; l2<lmax; l2++){
	if (l1!=l2 && overlap(l1,l2)!=0) dsum += fabs(overlap(l1,l2));
	if (l1==l2) dsum += fabs(overlap(l1,l2)-1.0);
      }
    }
    return dsum;
  }

  void cmpOverlap(function2D<double>& overlap, const vector<spline1D<double> >& fu, const mesh1D& tau, double beta, int k=5){
    int lmax = fu.size();

    int M = static_cast<int>(pow(2.0,k))+1;    // Number of points in Romberg integration routine
    vector<double> utu(M); // storage for Romberg

    overlap.resize(lmax,lmax);
    overlap=0.0;
    int nt = tau.size();
    for (int l1=0; l1<lmax; l1++){
      double odd_1 = fu[l1][0]*fu[l1][nt-1];
      for (int l2=0; l2<=l1; l2++){
	double odd_2 = fu[l2][0]*fu[l2][nt-1];
	if (odd_1*odd_2 > 0){// both are even or odd
	  double oij=0; // overlap between u_{l1} and u_{l2} functions
	  for (int i=0; i<tau.size()-1; i++){
	    // here we integrate with a fine mesh of M points using romberg routine
	    double a = tau[i];  // integral only between t_i and t_{i+1} points
	    double b = tau[i+1];
	    double fa = fu[l1][i]   * fu[l2][i];
	    double fb = fu[l1][i+1] * fu[l2][i+1];
	    utu[0]   = fa;
	    utu[M-1] = fb;
	    tint ia=0;
	    for (int j=1; j<M-1; j++){
	      intpar p = tau.InterpLeft( a + (b-a)*j/(M-1.0), ia);
	      utu[j] = fu[l1](p) * fu[l2](p);
	    }
	    oij += romberg(utu, b-a);
	  }
	  overlap(l1,l2) = oij;
	  overlap(l2,l1) = oij;
	}
      }
    }
  }

  template <class container>
  void CreateGfromCoeff(container& gf, const functionb<double>& gl, const string& statistics="fermi"){

    vector<spline1D<double> >* pfU = (statistics=="fermi") ? &fU : &fU_bose;

    gf.resize(tau.size());
    for (int i=0; i<tau.size(); i++) gf[i] = 0.0;

    for (int l=0; l<lmax; l++)
      for (int i=0; i<tau.size(); i++)
	gf[i] += gl[l]*(*pfU)[l][i];
  }

  void CreateSplineFromCoeff(spline1D<double>& gf, const functionb<double>& gl, const string& statistics="fermi"){
    CreateGfromCoeff(gf,gl, statistics);
    {
      double x1 = 0.5*(tau[1]+tau[0]);
      double df1 = (gf[1]-gf[0])/(tau[1]-tau[0]);
      double x2 = 0.5*(tau[2]+tau[1]);
      double df2 = (gf[2]-gf[1])/(tau[2]-tau[1]);
      double df0 = df1 + (df2-df1)*(0.0-x1)/(x2-x1);
      int n = tau.size();
      x1 = 0.5*(tau[n-1]+tau[n-2]);
      df1 = (gf[n-1]-gf[n-2])/(tau[n-1]-tau[n-2]);
      x2 = 0.5*(tau[n-2]+tau[n-3]);
      df2 = (gf[n-2]-gf[n-3])/(tau[n-2]-tau[n-3]);
      double dfn = df1 + (df2-df1)*(tau[n-1]-x1)/(x2-x1);
      gf.splineIt(tau, df0, dfn);
    }
  }

  template <class container>
  void MatsubaraFrequency(container& giom, const spline1D<double>& gf, int nmax, const string& statistics="fermi"){
    giom.resize(nmax);
    //double beta=tau[tau.size()-1];
    double one = (statistics=="fermi") ? 1.0 : 0.0;
    for (int in=0; in<nmax; in++){
      double iom = (2*in+one)*M_PI/beta;
      giom[in] = gf.Fourier_(iom, tau);
    }
  }

  template <class functor>
  void ExpandAnalyticFunction(const functor& fG, function1D<double>& gl, const string& statistics="fermi") const
  {
    const vector<spline1D<double> >* pfU = (statistics=="fermi") ? &fU : &fU_bose;

    gl.resize(lmax);
    int M = static_cast<int>(pow(2.0,k_ovlp))+1;    // Number of points in Romberg integration routine
    vector<double> utu(M); // storage for Romberg
    //double beta = tau[tau.size()-1];
    for (int l=0; l<lmax; l++){
      gl[l]=0;
      for (int i=0; i<tau.size()-1; i++){
	double a = tau[i];  // integral only between t_i and t_{i+1} points
	double b = tau[i+1];
	double fa = (*pfU)[l][i]   * fG(a, beta);
	double fb = (*pfU)[l][i+1] * fG(b, beta);
	utu[0]   = fa;
	utu[M-1] = fb;
	tint ia=0;
	for (int j=1; j<M-1; j++){
	  double t = a + (b-a)*j/(M-1.0);
	  intpar p = tau.InterpLeft( t, ia);
	  utu[j] = (*pfU)[l](p) * fG(t, beta);
	}
	gl[l] += romberg(utu, b-a);
      }
    }
  }

  int l_critical(const functionb<double>& gl, double max_ratio=3., const string& statistics="fermi"){
    int nsa=5;

    int _lmax_ = statistics=="fermi" ? lmax : lmax_bose;

    if (_lmax_<=nsa) return _lmax_;

    function1D<double> gos(_lmax_);

    function1D<double>* pS = (statistics=="fermi") ? &S : &S_bose;

    for (int i=0; i<_lmax_; i++) gos[i] = fabs(gl[i]/(*pS)[i]);

    //cout<<"gos="<<gos<<endl;

    double* p = gos.begin() + nsa+1; // pointer to the fifth entry

    for (int l=nsa+1; l<_lmax_; ++l && ++p){
      double sa = *max_element(p-nsa,p);  // maximum in the past five coefficients
      //cout<<l<<" g[l]="<<gl[l]<<" gos[l]="<<gos[l]<<" gmax="<<sa<<" g/gmax="<<fabs(gos[l]/sa)<<endl;
      if (fabs(gos[l]/sa)>max_ratio){
	//cout<<"lcritical should be "<<l<<endl;
	return l;
      }
    }
    return _lmax_;
  }

  void GiveRealAxis(function2D<double>& Aw, const functionb<double>& gl, double max_ratio=3., const string& statistics="fermi")
  {
    int lcritical = l_critical(gl,max_ratio,statistics);
    Aw.resize(lcritical,om.size());

    function1D<double> gos(gl.size());

    function1D<double>* pS = (statistics=="fermi") ? &S : &S_bose;

    for (int i=0; i<gl.size(); i++) gos[i] = gl[i]/(*pS)[i];

    function2D<double>* pVt = (statistics=="fermi") ? &Vt : &Vt_bose;

    for(int i=0; i<om.size(); i++)
      Aw(0,i) = -gos[0]*(*pVt)(i,0);

    for (int l=1; l<lcritical; l++){
      for (int i=0; i<om.size(); i++)
	Aw(l,i) = Aw(l-1,i) - gos[l]*(*pVt)(i,l);
    }

    if (statistics=="bose"){
      for (int l=0; l<lcritical; l++)
	for (int iw=0; iw<om.size(); iw++)
	  Aw(l,iw) *= M_PI*om[iw];
    }
  }

  void SetUpFastInterp()
  {
    ff2.resize(tau.size(),lmax,2);
    dxi.resize(tau.size());
    for (int i=0; i<tau.size(); i++) dxi[i] = fU[0].dxi[i];
    for (int l=0; l<lmax; l++){
      for (int i=0; i<tau.size(); i++){
	ff2(i,l,0) = fU[l][i];
	ff2(i,l,1) = fU[l].f2[i];
      }
    }
  }
  void FastInterp(functionb<double>& res, const intpar& ip, double cst) const
  {
    int i= ip.i;
    double p = ip.p, q=1-ip.p;
    double dx26 = dxi[i]*dxi[i]/6.;
    double dq = dx26*q*(q*q-1);
    double dp = dx26*p*(p*p-1);
    q  *= cst;
    dq *= cst;
    p  *= cst;
    dp *= cst;
    double* __restrict__ _res = res.MemPt();
    const double* __restrict__ _ff2 = ff2[i];
    // this is fast equivalent of
    // res[l] += q * fU[l].f[i] + dq * fU[l].f2[i];
    for (int l=0; l<lmax; l++) _res[l] += q * _ff2[2*l] + dq * _ff2[2*l+1];
    _ff2 += 2*lmax;
    // this is fast equivalent of
    // res[l] += p * fu[l].f[i+1] + dp * fU[l].f2[i+1];
    for (int l=0; l<lmax; l++) _res[l] += p * _ff2[2*l] + dp * _ff2[2*l+1];
  }
  void FastInterp_(functionb<double>& res, const intpar& ip, double cst) const
  {// This is the same as FastInterp, except res is zeroth first.
    int i= ip.i;
    double p = ip.p, q=1-ip.p;
    double dx26 = dxi[i]*dxi[i]/6.;
    double dq = dx26*q*(q*q-1);
    double dp = dx26*p*(p*p-1);
    q  *= cst;
    dq *= cst;
    p  *= cst;
    dp *= cst;
    double* __restrict__ _res = res.MemPt();
    const double* __restrict__ _ff2 = ff2[i];
    for (int l=0; l<lmax; l++) _res[l] = q * _ff2[2*l] + dq * _ff2[2*l+1];
    _ff2 += 2*lmax;
    for (int l=0; l<lmax; l++) _res[l] += p * _ff2[2*l] + dp * _ff2[2*l+1];
  }

#ifdef _MPI
  void ComputeClCoefficients(function2D<double>& Cl, int my_rank, int mpi_size, int Master, int nw=0, int ntail=0, int cutoff_iom=0){
#else
  void ComputeClCoefficients(function2D<double>& Cl, int nw=0, int ntail=0, int cutoff_iom=0){
#endif
    bool BRISI = false;

    int Nt = tau.size();
    //double L = om[om.size()-1];
    //double beta = tau[tau.size()-1];
    if (nw==0)    nw = L*beta/3;      // number of Matsubara points treated exactly
    if (ntail==0) ntail=150;          // number of points used for the tail
    if (cutoff_iom==0) cutoff_iom=100;// how far the tail will extend. We will go to Matsubara index nw*cutoff_iom

    function1D<int> ind_om;
    create_log_mesh(ind_om, nw*cutoff_iom, nw, ntail); // log mesh with nw exact points in [0,...nw-1], and ntail points distributed between [nw,nw*cutoff_iom]
    mesh1D iom(ind_om.size());
    for (int in=0; in<ind_om.size(); in++)  iom[in] = ind_om[in];
    iom.SetUp(0);

    if (BRISI){
      cout.precision(11);
      cout<<"# Starting single loop for u(iw) "<<endl;
    }
    // This creates u_l(iom) from u_l(tau)
    function2D<complex<double> > u_om(lmax,iom.size());
    for (int lf=0; lf<lmax; lf++){
      for (int in=0; in<iom.size(); in++){
	double om = (2*iom[in]+1.)*M_PI/beta;
	u_om(lf,in) = fU[lf].Fourier(om, tau);
      }
    }
    if (BRISI) cout<<"# Starting double loop for w(iw) "<<endl;
    // This creates w_l(iom) from u_l(tau).
    // w_{lb,lf}(iom) = 1/beta \sum_{iOm} u^b_{lb}(iOm) * u^f_{lf}(iom+iOm)
    // which is equivalent to
    // w_{lb,lf}(iom) = Integrate[ u^b_{lb}(beta-t) * u^f_{lf}(t) e^{i*om*t}, {t,0,beta}]
    int Nl2 = lmax_bose*lmax;
#ifdef _MPI
    int ipr_proc = (Nl2 % mpi_size==0) ? Nl2/mpi_size : Nl2/mpi_size+1;
    int iistart = ipr_proc*my_rank;
    int iiend   = ipr_proc*(my_rank+1);
    if (iistart>Nl2) iistart=Nl2;
    if (iiend  >Nl2) iiend=Nl2;
    int _Nl2_ = max(Nl2, mpi_size*ipr_proc);
    if (BRISI) cout<<"istart="<<iistart<<" iend="<<iiend<<" pr_proc="<<ipr_proc<<" pr_proc*mpi_size="<<ipr_proc*mpi_size<<" Nl2="<<Nl2<<endl;
#else
    int iistart = 0;
    int iiend = Nl2;
    int ipr_proc = Nl2;
    int _Nl2_ = Nl2;
#endif
    function2D<complex<double> > w_iom(_Nl2_, iom.size());
    w_iom = 0.0;
    for (int ii=iistart; ii<iiend; ii++){
      // ii = lb * lmax + lf
      int lb = ii / lmax;
      int lf = ii % lmax;
      spline1D<double> wt(Nt);
      for (int it=0; it<tau.size(); it++)
	wt[it] = fU_bose[lb][Nt-1-it]*fU[lf][it]; // wt[t] = u^b_{lb}(beta-t) * u^f_{lf}(t)

      int n = tau.size();
      double x1 = 0.5*(tau[1]+tau[0]);
      double df1 = (wt[1]-wt[0])/(tau[1]-tau[0]);
      double x2 = 0.5*(tau[2]+tau[1]);
      double df2 = (wt[2]-wt[1])/(tau[2]-tau[1]);
      double df0 = df1 + (df2-df1)*(0.0-x1)/(x2-x1);
      x1 = 0.5*(tau[n-1]+tau[n-2]);
      df1 = (wt[n-1]-wt[n-2])/(tau[n-1]-tau[n-2]);
      x2 = 0.5*(tau[n-2]+tau[n-3]);
      df2 = (wt[n-2]-wt[n-3])/(tau[n-2]-tau[n-3]);
      double dfn = df1 + (df2-df1)*(beta-x1)/(x2-x1);
      wt.splineIt(tau, df0, dfn);
      // w(iom) = Fourier[ wt[t] ]
      for (int in=0; in<iom.size(); in++){
	double om = (2*iom[in]+1.)*M_PI/beta;
	//w_om(lb,lf,in) = wt.Fourier(om, tau);
	w_iom(ii,in) = wt.Fourier(om, tau);
      }
    }
#ifdef _MPI
    if (mpi_size>1)
      MPI_Allgather(MPI_IN_PLACE, ipr_proc*iom.size()*2, MPI_DOUBLE, w_iom.MemPt(), ipr_proc*iom.size()*2, MPI_DOUBLE, MPI_COMM_WORLD);
#endif
    if (BRISI) cout<<"# Starting last part"<<endl;
    int Nl = (lmax*lmax*lmax_bose);
    int offset=2;
    spline1D<double> ftmp(ntail+offset);
    mesh1D wom(ntail+offset);
    int ifirst = iom.size()-ntail-offset;
    for (int in=ifirst; in<iom.size(); in++) wom[in-ifirst] = iom[in];
    wom.SetUp(0);

    deque<pair<int,int> > icase;
    for (int i=0; i<Nl; i++){
      for (int j=i; j<Nl; j++){
	vector<int> li = get_lll(i);
	vector<int> lj = get_lll(j);
	if ( (li[0] + li[1] + li[2] + lj[0] + lj[1] + lj[2])%2 ) continue; // this vanishes due to symmetry
	icase.push_back( make_pair(i,j) );
      }
    }

#ifdef _MPI
    int pr_proc = (icase.size() % mpi_size==0) ? icase.size()/mpi_size : icase.size()/mpi_size+1;
    int istart = pr_proc*my_rank;
    int iend   = pr_proc*(my_rank+1);
    if (istart>icase.size()) istart=icase.size();
    if (iend  >icase.size()) iend=icase.size();
    if (BRISI) cout<<"istart="<<istart<<" iend="<<iend<<" pr_proc="<<pr_proc<<" pr_proc*mpi_size="<<pr_proc*mpi_size<<" icase.size="<<icase.size()<<endl;
#else
    int istart = 0;
    int iend = icase.size();
    int pr_proc = icase.size();
#endif
    function1D<double> Cl0(pr_proc);
    Cl0=0;
    // This loop is parallelized
    for (int ii=istart; ii<iend; ii++){
      int i = icase[ii].first;
      int j = icase[ii].second;
      vector<int> li = get_lll(i);
      vector<int> lj = get_lll(j);
      if ( (li[0] + li[1] + li[2] + lj[0] + lj[1] + lj[2])%2 ) continue; // this vanishes due to symmetry

      // lj[2],lj[1],lj[0]  :  bose,fermi,fermi
      //                       lb * lmax + lf
      //complex<double> * w_1 = w_iom[ lj[2]*lmax + li[1] ].MemPt();
      //complex<double> * w_2 = w_iom[ li[2]*lmax + lj[1] ].MemPt();
      //complex<double> * u_1 = u_om[ li[0] ].MemPt();
      //complex<double> * u_2 = u_om[ lj[0] ].MemPt();
      complex<double> * w_1 = w_iom[ lj[2]*lmax + li[0] ].MemPt();
      complex<double> * w_2 = w_iom[ li[2]*lmax + lj[0] ].MemPt();
      complex<double> * u_1 = u_om[ li[1] ].MemPt();
      complex<double> * u_2 = u_om[ lj[1] ].MemPt();

      double wsum = 0.0;
      for (int in=0; in<iom.size()-ntail; in++)
	wsum += (w_1[in] * conj(w_2[in]) * u_1[in] * conj(u_2[in])).real();
      // fill int the 1D-interpolation for the tail
      for (int in=ifirst; in<iom.size(); in++)
	ftmp[in-ifirst] = (w_1[in] * conj(w_2[in]) * u_1[in] * conj(u_2[in])).real();
      // intepolate the tail
      int n=ftmp.size();
      ftmp.splineIt(wom, (ftmp[1]-ftmp[0])/(wom[1]-wom[0]), (ftmp[n-1]-ftmp[n-2])/(wom[n-1]-wom[n-2]) );
      // Now anlytically evaluate the sum over all points, but loop only over a few points in the tail.
      // The sum of the spline in the interval [n1,n2] can be analytically evaluated:
      //  1/2*(f[n2]+f[n1])*(n2-n1)  - 1/24*(n2-n1)*((n2-n1)^2-1) * (f2[n2]+f2[n1])
      //  where f2[n1], f2[n2] is the second derivative at n1 and n2 points.
      //  Note that the first and the last point come in with the weight 1/2, hence we need to correct for that.
      wsum += 0.5*ftmp[iom.size()-ntail-ifirst] + 0.5*ftmp[iom.size()-1-ifirst];
      for (int in=iom.size()-ntail; in<iom.size()-1; in++){
	double dh = iom[in+1]-iom[in];
	wsum += 0.5*(ftmp[in-ifirst]+ftmp[in+1-ifirst])*dh - dh/24.0*(dh*dh - 1.0)*(ftmp.f2[in+1-ifirst]+ftmp.f2[in-ifirst]);
      }
      /* //The above few lines are equivalent to
	 int iom_last = iom[iom.size()-1]; tint ia=0;
	 for (int iw=iom.size()-ntail; iw<iom_last; iw++){
	 double ww = ftmp( wom.InterpLeft(iw, ia) );
	 wsum += ww;
	 }
      */
      /*
      Cl(i,j) = 2*(wsum)/beta;
      Cl(j,i) = Cl(i,j);
      */
      Cl0[ii-istart] = 2*(wsum)/beta;
      //cout<<setw(4)<<i<<" "<<setw(4)<<j<<"     "<<setw(3)<<li[2]<<" "<<setw(3)<<li[1]<<" "<<setw(3)<<li[0]<<"; "<<setw(3)<<lj[2]<<" "<<setw(3)<<lj[1]<<" "<<setw(3)<<lj[0]<<"    "<<setw(10)<<left<<Cl(i,j)<<right<<endl;
    }
#ifdef _MPI
    if (mpi_size>1){
      if (BRISI){
	cout<<"MPI_Gather"<<endl;
	cout.flush();
      }
      function1D<double> cCl0;
      cCl0.resize(pr_proc*mpi_size);
      MPI_Gather(Cl0.MemPt(), pr_proc, MPI_DOUBLE, cCl0.MemPt(), pr_proc, MPI_DOUBLE, Master, MPI_COMM_WORLD);
      if (my_rank==Master){
	Cl0.resize(icase.size());
	for (int i=0; i<icase.size(); i++) Cl0[i] = cCl0[i];
      }
    }
    if (my_rank == Master){
#endif
      Cl.resize(Nl,Nl);
      Cl=0.0;
      for (int ii=0; ii<icase.size(); ii++){
	int i = icase[ii].first;
	int j = icase[ii].second;
	Cl(i,j) = Cl0[ii];
	Cl(j,i) = Cl0[ii];
      }
#ifdef _MPI
    }
#endif
  }
};

/*
double fGs(double tau, double beta){
  // On real axis, this corresponds to
  // A(w) = 1/2*(delta(w-x0)+delta(w+x0))
  // with x0=1.
  return -exp(-beta/2.)*cosh(beta/2.-tau);
}

double chis(double tau, double beta){
// On real axis this corresponds to
// chi''(x) = pi*[delta(w-x0)-delta(w+x0)]
// with x0=1/beta
  const double bx2 = 0.5;
  return cosh(bx2*(1-2.*tau/beta))/sinh(bx2);
}

class fGspl{
public:
  spline1D<double>& fg;
  mesh1D t;
  fGspl(spline1D<double>& fg_, mesh1D& t_) : fg(fg_), t(t_) {}
  double operator()(double x, double beta) const{
    return fg(t.Interp(x));
  }
};

int main(){
  double x0=0.005;
  double L=10.;
  int Nw = 500;
  //double beta=9.35887;
  double beta=100;
  int lmax = 40;
  int k_ovlp=5;
  int Ntau = static_cast<int>(5*beta+100);
  bool readGFile=false;//true;

  spline1D<double> ginp;
  mesh1D tc;
  if (readGFile){
    ifstream  inp("../Gtau.dat_4");
    deque<double> ct, gt;
    while (inp){
      double t, g;
      inp >> t >> g;
      //cout<<t<<" "<<g<<endl;
      if (!inp) break;
      ct.push_back(t);
      gt.push_back(g);
      inp.ignore(numeric_limits<streamsize>::max(), '\n');
    }

    beta = ct.back()+ct[0];
    cout<<"beta="<<beta<<endl;
    ginp.resize(gt.size());
    tc.resize(ct.size());
    {
      for (int i=0; i<ct.size(); i++){
	tc[i] = ct[i];
	ginp[i] = gt[i];
      }
      tc.SetUp(0);
      double df0 = (ginp[1]-ginp[0])/(tc[1]-tc[0]);
      int n = tc.size();
      double dfn = (ginp[n-1]-ginp[n-2])/(tc[n-1]-tc[n-2]);
      cout<<"df0="<<df0<<" dfn="<<dfn<<endl;
      ginp.splineIt(tc, df0, dfn);
    }
  }
  fGspl Fg(ginp, tc);

  clog<<"beta="<<beta<<" Ntau="<<Ntau<<endl;

  SVDFunc svdf;

  svdf.cmp("both", beta, lmax, Ntau, x0, L, Nw, clog);

  // Print SVD functions
  svdf.Print("uls.dat",beta);
  svdf.Print("ulb.dat",beta, "bose");

  // Now we check how well is the orthogonality obeyed after reorthogonalization
  clog<<"Final overlap is smaller than "<<svdf.CheckOverlap(beta)<<endl;

  ofstream fo("cc.dat");
  mesh1D& tau = svdf.tau;
  function1D<double> Gts(tau.size());
  for (int i=0; i<tau.size(); i++){
    if (readGFile) Gts[i] = Fg(tau[i],beta);
    else Gts[i] = fGs(tau[i], beta);
    fo<<tau[i]<<" "<<Gts[i]<<" "<<svdf[0][i]<<" "<<svdf[2][i]<<endl;
  }
  fo.close();

  function1D<double> gl(lmax);
  if (readGFile) svdf.ExpandAnalyticFunction(Fg, gl);
  else svdf.ExpandAnalyticFunction(fGs, gl);

  function1D<double> gl2(lmax);
  svdf.SetUpFastInterp();
  {
    k_ovlp=5;
    gl2.resize(lmax);
    gl2 = 0.0;
    int M = pow(2,k_ovlp)+1;    // Number of points in Romberg integration routine
    function2D<double> utu_T(M,lmax); // storage for Romberg
    double beta = tau[tau.size()-1];
    for (int i=0; i<tau.size()-1; i++){
      utu_T=0;
      double a = tau[i];  // integral only between t_i and t_{i+1} points
      double b = tau[i+1];
      tint ia=0;
      for (int j=0; j<M; j++){
	double t = a + (b-a)*j/(M-1.0);
	intpar p = tau.InterpLeft( t, ia);
	double cst;
	if (readGFile) cst = Fg(t,beta);
	else cst = fGs(t, beta);
	///// Slow equivalent ////
	//for (int l=0; l<svdf.lmax; l++) utu_T(j,l) += cst * svdf.fU[l](p);
	svdf.FastInterp(utu_T[j], p, cst );
      }
      function1D<double> utu(M);
      for (int l=0; l<lmax; l++){
	for (int j=0; j<M; j++) utu[j] = utu_T(j,l);
	gl2[l] += romberg(utu, b-a);
      }
    }
  }
  cout<<"Should be equal:"<<endl;
  for (int l=0; l<lmax; l++)
    cout<<l<<" "<<gl[l]<<" "<<gl2[l]<<" "<<gl[l]-gl2[l]<<endl;

  spline1D<double> gf;
  svdf.CreateSplineFromCoeff(gf, gl);
  int nmax=5000;
  function1D<complex<double> > giom;
  svdf.MatsubaraFrequency(giom, gf, nmax);

  ofstream fig("giom.dat");
  fig.precision(16);
  for (int in=0; in<nmax; in++){
    double iom = (2*in+1)*M_PI/beta;
    double gi_exact = -iom/(iom*iom+1.);
    fig<<iom<<" "<<giom[in]<<" "<<gi_exact<<endl;
  }

  ofstream fog("gl.dat");
  fog.precision(16);
  for (int l=0; l<svdf.lmax; l++)
    if (fabs(gl[l])>1e-10)
      fog<<l<<" "<<gl[l]<<endl;
  fog.close();

  ofstream foG("Gapp.dat");
  foG.precision(16);
  for (int i=0; i<tau.size(); i++){
    double dval=0;
    for (int l=0; l<svdf.lmax; l++) dval += svdf[l][i]*gl[l];
    foG<<tau[i]<<" "<<dval<<" "<<Gts[i]<<" "<<fabs(dval-Gts[i])<<endl;
  }
  foG.close();

  //int lc = svdf.l_critical(gl);
  //cout<<lc<<endl;

  function2D<double> Aw;
  ofstream fAw("Aw.out");
  svdf.GiveRealAxis(Aw, gl, 2.0);
  for (int i=0; i<svdf.om.size(); i++){
    fAw << svdf.om[i]<<" ";
    for (int l=0; l<Aw.size_N(); l++){
      fAw << Aw(l,i) <<" ";
    }
    fAw << endl;
  }

  svdf.ExpandAnalyticFunction(chis, gl, "bose");
  svdf.CreateSplineFromCoeff(gf, gl, "bose");
  ofstream foc("Chid.dat");
  foc.precision(16);
  for (int i=0; i<svdf.tau.size(); i++){
    double val = chis(svdf.tau[i],beta);
    foc<<svdf.tau[i]<<" "<<gf[i]<<" "<<val<<" "<<fabs(val-gf[i])<<endl;
  }
  foc.close();

  svdf.MatsubaraFrequency(giom, gf, nmax, "bose");
  ofstream sig("siom.dat");
  sig.precision(16);
  for (int in=0; in<nmax; in++){
    double iom = (2*in+0.0)*M_PI/beta;
    double x0 = 1./beta;
    double gi_exact = 2*x0/(x0*x0+iom*iom);
    sig<<iom<<" "<<giom[in]<<" "<<gi_exact<<" "<<giom[in]-gi_exact<<endl;
  }

  cout<<"bose coefficients are"<<gl<<endl;

  svdf.GiveRealAxis(Aw, gl, 10.0, "bose");
  ofstream sAw("sw.out");
  for (int i=0; i<svdf.om.size(); i++){
    sAw << svdf.om[i]<<" ";
    for (int l=0; l<Aw.size_N(); l++) sAw << Aw(l,i) <<" ";
    sAw << endl;
  }

  svd.Print("usb.dat",common::beta,"bose");
  for (int l=0; l<svdf.lmax_bose; l++) svdf.MatsubaraFrequency(bUw[l], svdf.fU_bose[l], N_W);

  return 0;
}
*/


text:
/* Author:  G. Jungman */

/* Convenience header */
#ifndef __GSL_SPECFUNC_H__
#define __GSL_SPECFUNC_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <gsl/gsl_sf.h>

#endif /* __GSL_SPECFUNC_H__ */

max_stars_repo_name: zzpwahaha/Chimera-Control-Trim

max_issues_repo_name: zzpwahaha/Chimera-Control-Trim

max_forks_repo_name: zzpwahaha/Chimera-Control-Trim


text:
/* -----------------------------------------------------------------------------
 * Copyright 2021 Jonathan Haigh
 * SPDX-License-Identifier: MIT
 * ---------------------------------------------------------------------------*/

#ifndef SQ_INCLUDE_GUARD_system_linux_SqFieldSchemaImpl_h_
#define SQ_INCLUDE_GUARD_system_linux_SqFieldSchemaImpl_h_

#include "core/typeutil.h"
#include "system/SqFieldSchema.gen.h"
#include "system/schema.h"

#include <gsl/gsl>

namespace sq::system::linux {

class SqFieldSchemaImpl : public SqFieldSchema<SqFieldSchemaImpl> {
public:
  explicit SqFieldSchemaImpl(const FieldSchema &field_schema);

  SQ_ND Result get_name() const;
  SQ_ND Result get_doc() const;
  SQ_ND Result get_params() const;
  SQ_ND Result get_return_type() const;
  SQ_ND Result get_return_list() const;
  SQ_ND Result get_null() const;

  SQ_ND Primitive to_primitive() const override;

private:
  gsl::not_null<const FieldSchema *> field_schema_;
};

} // namespace sq::system::linux

#endif // SQ_INCLUDE_GUARD_system_linux_SqFieldSchemaImpl_h_

max_stars_repo_name: jonathanhaigh/sq

max_issues_repo_name: jonathanhaigh/sq

max_forks_repo_name: jonathanhaigh/sq


text:
/*--------------------------------------------------------------------------------------------------
  FUNCTION: ssblock_lu

  DESCRIPTION:
--------------------------------------------------------------------------------------------------
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>

#ifdef OPENBLAS
#include <cblas.h>
#elif MKL
#include <mkl.h>
#endif

#define max(a, b) ( a < b ? b : a )

extern void getlvl2minors_lu_(int*, const int*, double*, int*, double*, int*);

void ssblock_lu(double *csc, int no, const int nv1, const int nv2, const int ns1, const int ns2, double *wf1, double *wf2, double *l1minor, double *ss, const int print_level)
{
	int ld_csc;
        int ld_msum, ld_tmp;
        int ld_l2minors;
	int tmp_cols;
        int i, o1, o2;
        int pos_tmp, pos_msum;
        int st1, st2;
	int num_threads;
	double *msum = NULL;  //Auxiliary matrix to store computed determinants
	double *tmp = NULL;  // Auxiliary matrix for intermediate results
	double *l2minors = NULL;
	double zero = 0.0, one = 1.0, neg_one = -1.0;
	int    one_i = 1;
	double coef;
	double ss_local;
	int v1;

        /* Timing routines */
        double start, finish;
        double time00, time0;
        double time_l2m, time_det, time_sum, time_tot;

	/* Auxiliary variables */
	int ndevices;
	int devices[10];

	if (print_level >= 2) {
	    time00 = omp_get_wtime();
            printf("\n         ---- start ssblock subroutine ----\n");
	}

	/* Set times to 0 */
	time_l2m = 0.0;
	time_det = 0.0;
	time_sum = 0.0;

        /* Set leading dimensions */
        ld_csc =  no + nv1;
        ld_l2minors = no*no*no;
	ld_tmp = max(no, nv1);
	tmp_cols = nv2;
	ld_msum = nv1;

        /* Allocate arrays */
        msum = (double*) malloc (ld_msum * nv2 * sizeof(double));
        tmp = (double*) malloc(ld_tmp * tmp_cols * ns2 * sizeof(double));
	l2minors = (double*) malloc(ld_l2minors * no * sizeof(double));

#ifdef OPENBLAS
	num_threads = openblas_get_num_threads();
#elif MKL
	num_threads = mkl_get_max_threads();
#endif

	/* Compute the determinants of all possible l2minors of the reference matrix. Pass through all combinations of rows and columns (r1, c1, r2, c2) to removed in order to form a l2minor */
	time0 = omp_get_wtime();
	getlvl2minors_lu_(&no, &nv2, csc, &ld_csc, l2minors, &ld_l2minors);
	time_l2m = omp_get_wtime() - time0;

        /* Scale lower csc matrix block with alternating array (1,-1,1,-1...) i.e. starting from the second
         * column, scale every second column with -1.
         * Arrays wf1 and wf2 are used as auxiliary arrays and are not referenced (unchaged on exit).
	 */
	cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, nv1, no/2, nv2, zero, wf1, nv1, wf2, nv2, neg_one, &csc[no+nv1 + no], 2*ld_csc);

	/* Scale ss matrix to zero */
        cblas_dscal(ns1*ns2, zero, ss, 1);

	/* Start iterate through the blocks */
	for( o1 = 0; o1 < no; o1++ ){
	    for( o2 = 0; o2 < no; o2++ ){

	/*	printf("[ssblock_lu] Computing block (%d, %d)\n", o1, o2); */

		/* Apply rows from the right of the ref matrix, cscmat positions (1:no, no+1:np+nv2)
		 * Compute first product l2minors * csc(1:no, no+1:no+nv2) */
        	time0 = omp_get_wtime();
		pos_tmp = o1*no*no + o2*no;
		ld_tmp = no;

        	cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, no, nv2, no, one, &l2minors[pos_tmp], ld_l2minors, &csc[no*ld_csc], ld_csc, zero, tmp, ld_tmp);

                pos_msum = o1*no*nv1 + o2*nv1;

		cblas_dgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, nv1, nv2, no, pow(neg_one, o1+2), &csc[no], ld_csc, tmp, ld_tmp, zero, msum, ld_msum);

		/* Get l1minor determinants and apply permutation sign */
		coef = pow(neg_one, o1+o2) * l1minor[o2*no + o1];

		/* Get coefficients (no+1:no+nv1, no+1:no+nv2) into auxiliary storage msum */
		for( i = 0; i < nv2; i++ ){
			cblas_daxpy(nv1, coef, &csc[(no+i)*ld_csc+no], 1, &msum[i*ld_msum], 1);
		}
        	time_det += omp_get_wtime() - time0;

	        time0 = omp_get_wtime();
		/* Compute ss matrix  */
		for( v1 = 0; v1 < nv1; v1++ ){
			cblas_dgemv( CblasColMajor, CblasTrans, nv2, ns2, one, &wf2[o2*nv2], no*nv2, &msum[v1], nv1, zero, tmp, 1 );
			cblas_dger( CblasColMajor, ns1, ns2, one, &wf1[o1*nv1+v1], no*nv1, tmp, one, ss, ns1 );
		}
        	finish = omp_get_wtime();
		time_sum +=  omp_get_wtime() - time0;
	    }
	}

	/* Free allocated spaces */
	free(msum);
	free(tmp);
	free(l2minors);

	/* Print timings */
	if (print_level >= 2) {
            printf("\n         ssblock time:\n");
	    printf("             L2 minors                  - time (sec): %.4lf \n",  time_l2m);
	    printf("             Determinants from minors   - time (sec): %.4lf \n",  time_det);
	    printf("             Final sum                  - time (sec): %.4lf \n",  time_sum);
	    printf("                                                     --------------\n");
	    printf("             Total                      - time (sec): %.4lf \n",  time_tot);
            printf("         ---- end ssblock subroutine ----\n");
	}

}

ext: c


text:
#include <stdio.h>
#include <stdlib.h>
//#include <gsl_matrix.h>

#define IKFAST_HAS_LIBRARY // Build IKFast with API functions
#define IKFAST_NO_MAIN // Don't include main() from IKFast

#include "ikfast.cpp"
#include "kin.h"

int KinCalcFowardKinematics(jointAngles_t* angles, eeCoord_t* ee) {
	int i;
	double eetrans[3];
	double eerot[9];

	for(i=0;i<9;i++){
		eerot[i] = -48;
	}

	ComputeFk((double*)angles,eetrans,eerot);

	printf("Found fk solution for end frame: \n\n");

	printf("  Translation:  x: %f  y: %f  z: %f  \n", eetrans[0], eetrans[1], eetrans[2] );
	printf("\n");
	printf("     Rotation     %f   %f   %f  \n", eerot[0], eerot[1], eerot[2] );
	printf("       Matrix:    %f   %f   %f  \n", eerot[3], eerot[4], eerot[5] );
	printf("                  %f   %f   %f  \n", eerot[6], eerot[7], eerot[8] );
	printf("\n");

	return 0;
}

#ifdef STANDALONE

int main(int argc, char **argv){
	int i;
	if(argc != 9) {
		printf("usage: %s x1 y1 z1 w1 x2 y2 z2 w2\n", argv[0]);
		return 1;
	}

	const int points = 100;

	waypoint_t waypoint[points];

	waypoint[0].x = atof(argv[1]);
	waypoint[0].y = atof(argv[2]);
	waypoint[0].z = atof(argv[3]);
	waypoint[0].w = atof(argv[4]);
	waypoint[points-1].x = atof(argv[5]);
	waypoint[points-1].y = atof(argv[6]);
	waypoint[points-1].z = atof(argv[7]);
	waypoint[points-1].w = atof(argv[8]);

	KinCalcFowardKinematics(&angles, &ee);

}

#endif

max_stars_repo_name: bcmetz/arm

max_issues_repo_name: bcmetz/arm

max_forks_repo_name: bcmetz/arm


text:
/*
# Program to run a Monte Carlo radiation transfer through the 2D
# simulations of GRB jets.
#
# Python code written by D. Lazzati at Oregonstate, C code written by Tyler Parsotan @ Oregon State
# ver 0.1 July 8, 2015
# ver 1.1 July 20, 2015: added record of number of scatterings, included
# 	all terms in weight. Should now give correct light curves.
# ver 1.2 July 21, 2015: added parameter file to keep track of input
# 	params of each simulation

# ver 2.0 July 22, 2015: corrected the problem that arises when there is
# 	no scattering in the time span of one frame. Fixed output arrays dimension.

# ver 2.1 July 25, 2015: fixed bug that did not make the number of
# 	scattering grow with the number of photons.

# ver 3.0 July 28, 2015: using scipy nearest neighbor interpolation to
# 	speed things up. Gained about factor 2

# ver 3.1 July 29, 2015: added radial spread of photon injection points
# ver 3.2 July 31, 2015: added Gamma to the weight of photons!!!

# ver 4.0 Aug 5, 2015: try to speed up by inverting cycle
# ver 4.1 Aug 8, 2015: add spherical test as an option
# ver 4.2 Aug 9, 2015: saving files appending rather than re-writing
# ver 4.3 Aug 11, 2015: corrected error in the calculation of the local temperature
# ver 4.4 Aug 13, 2015: added cylindrical test
# ver 4.5 Aug 18, 2015: fixd various problems pointed by the cylindrical test
# ver 4.6 Aug 21, 2015: corrected mean free path for large radii

# ver 5.0 Aug 25, 2015: corrected problem with high-T electrons and excess scatterings
# ver 5.1 Aug 25, 2015: cleaned-up coding
# ver 5.2 Sept 3, 2015: fixed problem with number of scatterings for multiple injections
 *
 * ver 6.0 Dec 28, 2016: rewrote the code in C, added checkpoint file so if the code is interrupted all the progress wont be lost, made the code only need to be compiled once for a given MC_XXX directory path
                                            so you just need to supply the sub directory of MC_XXX as a command line argument
*  version 7.0 used OpenMP to parallelize the code by angle and the function findminmfp()

 version 8.0 added 3D capabilities for RIKEN hydro data  and 2D capablities for RIKEN 2D hydro data and made it more efficient with grid selection to speed it up
*/

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <sys/types.h>
#include <sys/stat.h>
#include <unistd.h>
#include <dirent.h>
#include <math.h>
#include <gsl/gsl_rng.h>
#include "mclib_3d.h"
#include <omp.h>
#include "mpi.h"

/*
#define THISRUN "Science"
#define FILEPATH "/home/physics/parsotat/16OI/"
#define FILEROOT "rhd_jet_big_16OI_hdf5_plt_cnt_"
#define MC_PATH "MPI_CMC_16OI_SPHERICAL/"

#define THISRUN "Science"
#define FILEPATH "/Users/Tylerparsotan/Documents/Box Sync/1spike/"
#define FILEROOT "m0_rhop0.1big_hdf5_plt_cnt_"
#define MC_PATH "CMC_1spike/"
//#define MC_PATH "MC_16OI/Single_Photon_Cy_mc_total/"
 * */
 /*
 #define THISRUN "Science"
#define FILEPATH "/home/physics/parsotat/16OM/"
#define FILEROOT "rhd_jet_big_16OM_hdf5_plt_cnt_"
#define MC_PATH "DIR_TEST/"

 #define THISRUN "Science"
#define FILEPATH "/Volumes/DATA6TB/Collapsars/2D/HUGE_BOXES/VARY/40spikes/"
#define FILEROOT "m0_rhop0.1big_hdf5_plt_cnt_"
#define MC_PATH "CMC_40spikes_TEST/"
 * */

 #define THISRUN "Spherical"
#define FILEPATH "/Volumes/DATA6TB/Collapsars/2D/HUGE_BOXES/CONSTANT/16OI/"
//#define FILEPATH "/Users/Tylerparsotan//Documents/16OI_TEST/"
#define FILEROOT "rhd_jet_big_16OI_hdf5_plt_cnt_"
#define MC_PATH "TEST/"

#define MCPAR "mc.par"
#define RIKEN_SWITCH 0

int main(int argc, char **argv)
{
    //compile each time a macro is changed, have to supply the subfolder within the MC_PATH directory as a command line argument to the C program eg. MCRAT 1/

	// Define variables
	char flash_prefix[200]="";
	char mc_file[200]="" ;
    char this_run[200]=THISRUN;
    char *cyl="Cylindrical";
    char *sph="Spherical";
    char spect;//type of spectrum
    char restrt;//restart or not
    double fps, fps_modified, theta_jmin, theta_jmax, hydro_domain_y,hydro_domain_x ;//frames per second of sim, min opening angle of jet, max opening angle of jet in radians, max y value of fluid simulation domain
    double inj_radius_small, inj_radius_large,  ph_weight_suggest, ph_weight_small, ph_weight_large ;//radius at chich photons are injected into sim
    int frm0,last_frm, frm2_small, frm2_large, j=0, min_photons, max_photons, frm0_small, frm0_large ;//frame starting from, last frame of sim, frame of last injection
    int dim_switch=0;
    int find_nearest_grid_switch=0;
    int increment_inj=1, increment_scatt=1; //increments for injection loop and scattering loop, outer and inner loops respectively, the increment can change for RIKEN 3D hydro files

    double inj_radius;
    int frm2;
    char mc_filename[200]="";
    char mc_filename_2[200]="";
    char mc_operation[200]="";
    char mc_dir[200]="" ;
    int file_count = 0;
    DIR * dirp;
    struct dirent * entry;
    struct stat st = {0};
    double theta_jmin_thread=0, theta_jmax_thread=0;

    char flash_file[200]="";
    char log_file[200]="";
    FILE *fPtr=NULL; //pointer to log file for each thread
    double *xPtr=NULL,  *yPtr=NULL,  *rPtr=NULL,  *thetaPtr=NULL,  *velxPtr=NULL,  *velyPtr=NULL,  *densPtr=NULL,  *presPtr=NULL,  *gammaPtr=NULL,  *dens_labPtr=NULL;
    double *szxPtr=NULL,*szyPtr=NULL, *tempPtr=NULL; //pointers to hold data from FLASH files
    double *phiPtr=NULL, *velzPtr=NULL, *zPtr=NULL;
    int num_ph=0, array_num=0, ph_scatt_index=0, max_scatt=0, min_scatt=0,i=0; //number of photons produced in injection algorithm, number of array elleemnts from reading FLASH file, index of photon whch does scattering, generic counter
    double dt_max=0, thescatt=0, accum_time=0;
    double  gamma_infinity=0, time_now=0, time_step=0, avg_scatt=0, avg_r=0; //gamma_infinity not used?
    double ph_dens_labPtr=0, ph_vxPtr=0, ph_vyPtr=0, ph_tempPtr=0, ph_vzPtr=0;;// *ph_cosanglePtr=NULL ;
    double min_r=0, max_r=0, min_theta=0, max_theta=0;
    int frame=0, scatt_frame=0, frame_scatt_cnt=0, scatt_framestart=0, framestart=0;
    struct photon *phPtr=NULL; //pointer to array of photons

    int num_thread=0, angle_count=0;
    int num_angles=0, old_num_angle_procs=0; //old_num_angle_procs is to hold the old number of procs in each angle when cont sims, if  restarting sims this gets set to angle_procs
    int *frame_array=NULL, *proc_frame_array=NULL, *element_num=NULL, proc_frame_size=0;
    double *thread_theta=NULL; //saves ranges of thetas for each thread to go through
    double delta_theta=1;

    int myid, numprocs, angle_procs, angle_id, procs_per_angle;

   //new OpenMPI stuff
   MPI_Init(NULL,NULL);
   MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
   MPI_Comm_rank(MPI_COMM_WORLD, &myid);

    //new muliple threads injecting and propagating photons
    const gsl_rng_type *rng_t;
    gsl_rng *rng;
    gsl_rng_env_setup();
    rng_t = gsl_rng_ranlxs0;
    rng = gsl_rng_alloc (rng_t); //initalize random number generator to seed the others with random numbers

    //want to break up simulation by angle and injection frame & have each thread save data in its own folder
    //have each thread check if its directory is made and if its restarting (delete evrything) or if its continuing with a previous simulation
    //the angle and the injection frames will be the names of mc_dir, therefore read mc.par first in MC_XXX directory

    //make strings of proper directories etc.
	snprintf(flash_prefix,sizeof(flash_prefix),"%s%s",FILEPATH,FILEROOT );
    snprintf(mc_file,sizeof(flash_prefix),"%s%s%s",FILEPATH, MC_PATH,MCPAR);

    printf(">> mc.py:  Reading mc.par: %s\n", mc_file);

    readMcPar(mc_file, &hydro_domain_x, &hydro_domain_y, &fps, &theta_jmin, &theta_jmax, &delta_theta, &inj_radius_small,&inj_radius_large, &frm0_small,&frm0_large, &last_frm ,&frm2_small, &frm2_large, &ph_weight_small, &ph_weight_large, &min_photons, &max_photons, &spect, &restrt, &num_thread,&dim_switch); //thetas that comes out is in degrees
    //printf("%c\n", restrt);

    //divide up angles and frame injections among threads DONT WANT NUMBER OF THREADS TO BE ODD
    //assign ranges to array that hold them

    //leave angles in degrees here
    num_angles=(int) (((theta_jmax-theta_jmin)/delta_theta)) ;//*(180/M_PI));
    thread_theta=malloc( num_angles *sizeof(double) );
    *(thread_theta+0)=theta_jmin;//*(180/M_PI);
    //printf("%e\n", *(thread_theta+0));

    for (j=1;j<(num_angles); j++)
    {
        *(thread_theta+j)=*(thread_theta+(j-1))+delta_theta;
        //printf("%e\n", *(thread_theta+j));
    }

    //make comm without the procs that deal with angle
     //comm for angles

     procs_per_angle= numprocs/num_angles;
     //printf("%d\n", procs_per_angle);

     MPI_Comm angle_comm;
     if (restrt=='r') //uncomment this when I run MCRAT for sims that didnt originally save angle_procs
     {

        MPI_Comm_split(MPI_COMM_WORLD, myid/procs_per_angle , myid, &angle_comm);
        MPI_Comm_rank(angle_comm, &angle_id);
        MPI_Comm_size(angle_comm, &angle_procs);

        //printf("WORLD RANK/SIZE: %d/%d \t ROW RANK/SIZE: %d/%d\n", myid, numprocs, angle_id, angle_procs);

        theta_jmin_thread= (*(thread_theta+  (myid/procs_per_angle))) *(M_PI/180);
        theta_jmax_thread= theta_jmin_thread+(delta_theta*(M_PI/180));

        snprintf(mc_dir,sizeof(flash_prefix),"%s%s%0.1lf-%0.1lf/",FILEPATH,MC_PATH, theta_jmin_thread*180/M_PI, theta_jmax_thread*180/M_PI ); //have to add angle into this

        old_num_angle_procs=angle_procs;
     }
     else
     {
         MPI_Group sub_world_group;
         MPI_Comm sub_world_comm;
         int incl_procs[procs_per_angle*num_angles], count, sub_world_id;
         int total_num_to_restart=0;
         int color=1;
         int  *all_cont_process_idPtr=NULL, *each_num_to_restart_per_anglePtr=NULL, *tmp=NULL;
        //for restart='c' case if the number of processes isnt a multiple of procs_per_angle*num_angles make a comm out of those that are in order to analyze files and count number of processes for each angle range need to con't
        count=0;
        for (j=0;j<numprocs;j++)
        {
            if (j<procs_per_angle*num_angles)
            {
                incl_procs[count]=j;
                count++;
            }
        }

        if (myid<procs_per_angle*num_angles)
        {
            int myid_2=0;
            // Get the group of processes in MPI_COMM_WORLD and make a sub group to go through checkpoint files
            MPI_Group world_group;
            MPI_Comm root_angle_comm;

            MPI_Comm_group(MPI_COMM_WORLD, &world_group);
            MPI_Group_incl(world_group, procs_per_angle*num_angles, incl_procs, &sub_world_group);
            MPI_Comm_create_group(MPI_COMM_WORLD, sub_world_group, 0, &sub_world_comm);
            MPI_Comm_rank(sub_world_comm, &myid_2);

            MPI_Comm_split(sub_world_comm, myid_2/procs_per_angle , myid_2, &angle_comm);
            MPI_Comm_rank(angle_comm, &angle_id);
            MPI_Comm_size(angle_comm, &angle_procs);

            //create group of all the processes that have angle_id==0
            if (angle_id==0)
            {
                color=0; //set different color for root processes in each group of angle_comm
            }
            MPI_Comm_split(sub_world_comm, color , myid_2, &root_angle_comm); //create comm to exchange info about number of processes to restart for each angle range

            printf("WORLD RANK/SIZE: %d/%d \t ROW RANK/SIZE: %d/%d\n", myid, numprocs, angle_id, angle_procs);

            theta_jmin_thread= (*(thread_theta+  (myid_2/procs_per_angle))) *(M_PI/180);
            theta_jmax_thread= theta_jmin_thread+(delta_theta*(M_PI/180));

            snprintf(mc_dir,sizeof(flash_prefix),"%s%s%0.1lf-%0.1lf/",FILEPATH,MC_PATH, theta_jmin_thread*180/M_PI, theta_jmax_thread*180/M_PI ); //have to add angle into this

            //call the function to count the num of processes for each angle range that need to be con't
            int  count_cont_procs=0, total_cont_procs_angle=0, global_cont_procs=0;
            int *cont_proc_idsPtr=NULL, *total_cont_procs_angle_Ptr=NULL, *displPtr=NULL; //becomes the size of the number of old procceses
            int *cont_proc_ids_anglePtr=NULL;

            old_num_angle_procs=getOrigNumProcesses(&count_cont_procs,  &cont_proc_idsPtr, mc_dir, angle_id,  angle_procs,  last_frm, dim_switch, RIKEN_SWITCH);

            if (old_num_angle_procs==-1)
            {
                printf("MCRAT wasnt able to get a value of old_num_angle_procs to continue the simulation. Now exiting to prevent data corruption.\n" );
                MPI_Abort(MPI_COMM_WORLD, 1);
            }

            total_cont_procs_angle_Ptr=malloc(angle_procs*sizeof(int));
            displPtr=malloc(angle_procs*sizeof(int));
            MPI_Gather(&count_cont_procs,1,MPI_INT, total_cont_procs_angle_Ptr, 1, MPI_INT, 0,angle_comm );//hold the number of elements that each process will send the root process

            MPI_Barrier(angle_comm);
            MPI_Barrier(sub_world_comm);
            if (angle_id==0)
            {
                printf("Angle_procs: %d 1st gather: %d, %d, %d\n", angle_procs, *(total_cont_procs_angle_Ptr), *(total_cont_procs_angle_Ptr+1), *(total_cont_procs_angle_Ptr+2));
            }

            MPI_Reduce(&count_cont_procs, &total_cont_procs_angle, 1, MPI_INT, MPI_SUM, 0, angle_comm); //for each angle sum the number of procs to continue and pass it to the root for angle_comm

            cont_proc_ids_anglePtr=malloc(total_cont_procs_angle*sizeof(int)); //each root proc in angle comm has to hold the id's of the old set of processes to cont

            *(displPtr+0)=0;
            if (angle_id==0)
            {
                for (j=1;j<angle_procs;j++)
                {
                    *(displPtr+j)=(*(displPtr+j-1))+(*(total_cont_procs_angle_Ptr+j-1 )); //set the displacement for each proces to put its vector of pprocess IDs that need to be continued
                    printf("Displacement: %d\n", *(displPtr+j));
                }

            }

            MPI_Gatherv(cont_proc_idsPtr,count_cont_procs,MPI_INT, cont_proc_ids_anglePtr, total_cont_procs_angle_Ptr, displPtr , MPI_INT, 0,angle_comm ); //send the vectors with the ids of the old processes that need to be cont to root in angle_comm

            MPI_Barrier(angle_comm);
            MPI_Barrier(sub_world_comm);

            if (angle_id==0)
            {
                printf("Total Cont Procs: %d\n", total_cont_procs_angle);
                for (j=0;j<total_cont_procs_angle;j++)
                {
                    {
                        printf("ID: %d\n", *(cont_proc_ids_anglePtr+j));
                    }
                }
            }

            //each root for angle_comm has the number of processes each angle range needs to restart and the array of what the IDs of those processes used to be
            //now have to combine all that info for rank 0 in MPI_COMM_WORLD and then end it to all processes in MPI_COMM_WORLD
            //if (myid==0)
            {
                free(displPtr);
                displPtr=NULL;
                //initalize variables to hold all data

                each_num_to_restart_per_anglePtr=malloc(num_angles*sizeof(int));
                displPtr=malloc(num_angles*sizeof(int));
                *(displPtr+0)=0;
            }

            MPI_Barrier(angle_comm);
            MPI_Barrier(sub_world_comm);

            if (angle_id==0)
            {
                //this is the part where all the root processes of angle_comm transfer thier info to the root proc of MPI_WORLD
                MPI_Reduce(&total_cont_procs_angle, &total_num_to_restart, 1, MPI_INT, MPI_SUM, 0, root_angle_comm); //for each angle sum the number of procs to continue and pass it to the root for MPI_COMM_WORLD

                MPI_Gather(&total_cont_procs_angle,1,MPI_INT, each_num_to_restart_per_anglePtr, 1, MPI_INT, 0,root_angle_comm );//hold the number of elements that each process sent the root  for MPI_COMM_WORLD

                if (myid==0)
                {
                    for (j=1;j<num_angles;j++)
                    {
                        *(displPtr+j)=(*(displPtr+j-1))+(*(each_num_to_restart_per_anglePtr+j-1 )); //set the displacement for each proces to put its vector of pprocess IDs that need to be continued
                    }
                }

                all_cont_process_idPtr=malloc(total_num_to_restart*sizeof(int));

                MPI_Gatherv(cont_proc_ids_anglePtr, total_cont_procs_angle, MPI_INT, all_cont_process_idPtr, each_num_to_restart_per_anglePtr,   displPtr, MPI_INT, 0, root_angle_comm);
            }

            MPI_Barrier(angle_comm);
            MPI_Barrier(sub_world_comm);

            if (myid==0)
            {
                printf("Global Cont Procs: %d\n", total_num_to_restart);
                for (j=0;j<total_num_to_restart;j++)
                {
                    {
                        printf("Global ID: %d\n", *(all_cont_process_idPtr+j));
                    }
                }
            }

            //destroy the old comms
            MPI_Barrier(angle_comm);
            MPI_Barrier(sub_world_comm);
            //destroy current angle comm and recreate a new one
            MPI_Comm_free(&root_angle_comm);
            MPI_Comm_free(&angle_comm);
            MPI_Comm_free(&sub_world_comm);
            MPI_Group_free(&sub_world_group);
            MPI_Group_free(&world_group);
            free(cont_proc_idsPtr);
            free(cont_proc_ids_anglePtr);
            free(total_cont_procs_angle_Ptr);
            free(displPtr);
            //free(each_num_to_restart_per_anglePtr);
            //free(all_cont_process_idPtr);
        }

        //send all of myid==0 data to all processes in MPI_COMM_WORLD
        MPI_Bcast( &total_num_to_restart, 1, MPI_INT, 0, MPI_COMM_WORLD );

        if (total_num_to_restart>0)
        {
            if (myid != 0 )
            {
                //allocate data of appropriate size for all processes to hold the data from MPI_Bcast
                tmp=realloc(all_cont_process_idPtr,total_num_to_restart *sizeof(int));
                if (tmp!=NULL)
                {
                    all_cont_process_idPtr=tmp;
                }
                else
                {
                    printf("Error with reserving space to hold data about restarting process ID's\n");
                }
                //free(tmp);
                tmp=realloc(each_num_to_restart_per_anglePtr, num_angles*sizeof(int));
                if (tmp!=NULL)
                {
                    each_num_to_restart_per_anglePtr=tmp;
                }
                else
                {
                    printf("Error with reserving space to hold data about restarting process numbers for each angle range\n");
                }
                //free(tmp);
            }

            MPI_Bcast( all_cont_process_idPtr, total_num_to_restart, MPI_INT, 0, MPI_COMM_WORLD );
            MPI_Bcast( each_num_to_restart_per_anglePtr, num_angles, MPI_INT, 0, MPI_COMM_WORLD );
            MPI_Bcast( &old_num_angle_procs, 1, MPI_INT, 0, MPI_COMM_WORLD );

            MPI_Barrier(MPI_COMM_WORLD);
            if (myid==numprocs-1)
            {
                printf("Number of processes: %d\n", old_num_angle_procs);
                printf("restarting process numbers for each angle range: %d, %d, %d\n", *(each_num_to_restart_per_anglePtr), *(each_num_to_restart_per_anglePtr+1), *(each_num_to_restart_per_anglePtr+2));
            }

            //assign proper number of processes to each angle range to con't sims and then reset angle_id to original value from when simulation was first started
            color=0; //by default all processes have this value

            count=0;
            for (j=0;j<num_angles;j++)
            {
                if (myid>=count   &&   myid<count+(*(each_num_to_restart_per_anglePtr+j)) )
                {
                    color=j;
                }
                count+=(*(each_num_to_restart_per_anglePtr+j));
                printf("Myid: %d, Color: %d, Count %d, Num To Start Per Angle: %d\n", myid, color, count, (*(each_num_to_restart_per_anglePtr+j)));
            }

            MPI_Comm_split(MPI_COMM_WORLD, color , myid, &angle_comm);
            MPI_Comm_rank(angle_comm, &angle_id);
            MPI_Comm_size(angle_comm, &angle_procs);

            printf("WORLD RANK/SIZE: %d/%d \t ROW RANK/SIZE: %d/%d\n", myid, numprocs, angle_id, angle_procs);

            angle_procs=old_num_angle_procs;

            //reset the angle for each process
            theta_jmin_thread= (*(thread_theta+  color)) *(M_PI/180);
            theta_jmax_thread= theta_jmin_thread+(delta_theta*(M_PI/180));

            //reset the angle_id for each process
            count=0;
            for (j=0;j<color;j++)
            {
                count+=(*(each_num_to_restart_per_anglePtr+j));
            }

            angle_id=(*(all_cont_process_idPtr+count+angle_id));

            snprintf(mc_dir,sizeof(flash_prefix),"%s%s%0.1lf-%0.1lf/",FILEPATH,MC_PATH, theta_jmin_thread*180/M_PI, theta_jmax_thread*180/M_PI ); //have to add angle into this
        }
        else
        {
            //if there are no more processes to continue just break up processes normally  so they read in checkpoint files of completed processes and jump to merging files
            MPI_Comm_split(MPI_COMM_WORLD, myid/procs_per_angle , myid, &angle_comm);
            MPI_Comm_rank(angle_comm, &angle_id);
            MPI_Comm_size(angle_comm, &angle_procs);
        }
        free(all_cont_process_idPtr);
        free(each_num_to_restart_per_anglePtr);
    }

         MPI_Barrier(MPI_COMM_WORLD);

        if ((theta_jmin_thread >= 0) &&  (theta_jmax_thread <= (2*M_PI/180) )) //if within small angle (0-2 degrees) use _small inj_radius and frm2 have to think about this for larger domains
        {
            inj_radius=inj_radius_small;
            frm2=frm2_small;
            frm0=frm0_small;
            ph_weight_suggest=ph_weight_small;
        }
        else
        {
            inj_radius=inj_radius_large;
            frm2=frm2_large;
            frm0=frm0_large;
            ph_weight_suggest=ph_weight_large;
        }

        //make vector to hold the frames we are injecting in, vector should have (frm2-frm0)/angle_procs slots, if fps is const
            proc_frame_size=ceil((frm2-frm0)/ (float) angle_procs);
            frame_array=malloc(((frm2-frm0)+1)*sizeof(int));

            for (j=0;j<((frm2-frm0)+1); j++)
            {
                *(frame_array+j)=frm0+j ;
            	//printf("proc: %d frame: %d\n", angle_id, *(frame_array+j));
            }

        {
            //set this now incase there is no checkpoint file, then this wont be overwritten and the corretc values will be passed even if the user decides to restart
            framestart=(*(frame_array +(angle_id*proc_frame_size)));
            scatt_framestart=framestart;

            if (angle_id != (angle_procs-1))
            {
                frm2=(*(frame_array +((angle_id*proc_frame_size) + proc_frame_size-1) )); //section off blocks of the frame_array to give to each angle_id
            }
            else
            {
                frm2=(*(frame_array + (frm2-frm0) )); //if angle_id is last give it the last set, even if its uneven
            }

            if (restrt=='c')
            {
                printf(">> mc.py:  Reading checkpoint\n");
                //#pragma omp critical

                    readCheckpoint(mc_dir, &phPtr, &frm2, &framestart, &scatt_framestart, &num_ph, &restrt, &time_now, angle_id, &angle_procs, dim_switch, RIKEN_SWITCH);

                /*
                for (i=0;i<num_ph;i++)
                {
                    printf("%e,%e,%e, %e,%e,%e, %e, %e\n",(phPtr+i)->p0, (phPtr+i)->p1, (phPtr+i)->p2, (phPtr+i)->p3, (phPtr+i)->r0, (phPtr+i)->r1, (phPtr+i)->r2, (phPtr+i)->num_scatt );
                }
                */
                if (restrt=='c')
                {
                    printf(">> Rank %d: Starting from photons injected at frame: %d out of %d\n", angle_id,framestart, frm2);
                    printf(">> Rank %d with angles %0.1lf-%0.1lf: Continuing scattering %d photons from frame: %d\n", angle_id, theta_jmin_thread*180/M_PI, theta_jmax_thread*180/M_PI,num_ph, scatt_framestart);
                    printf(">> Rank %d with angles %0.1lf-%0.1lf: The time now is: %e\n", angle_id, theta_jmin_thread*180/M_PI, theta_jmax_thread*180/M_PI,time_now);
                }
                else
                {
                    printf(">> Rank %d with angles %0.1lf-%0.1lf: Continuing simulation by injecting photons at frame: %d out of %d\n", angle_id, theta_jmin_thread*180/M_PI, theta_jmax_thread*180/M_PI,framestart, frm2); //starting with new photon injection is same as restarting sim
                }

            }
            else if ((stat(mc_dir, &st) == -1) && (restrt=='r'))
            {
                mkdir(mc_dir, 0777); //make the directory with full permissions

            }
            else
            {
                if (angle_id==0)
                {
                    printf(">> proc %d with angles %0.1lf-%0.1lf:  Cleaning directory \n",angle_id, theta_jmin_thread*180/M_PI, theta_jmax_thread*180/M_PI);
                    dirp = opendir(mc_dir);
                    while ((entry = readdir(dirp)) != NULL)
                    {
                        if (entry->d_type == DT_REG) { /* If the entry is a regular file */
                             file_count++; //count how many files are in dorectory
                        }
                    }
                    printf("File count %d\n", file_count);

                    if (file_count>0)
                    {
                        for (i=0;i<=last_frm;i++)
                        {

                            snprintf(mc_filename,sizeof(mc_filename),"%s%s%d%s", mc_dir,"mcdata_",i,"_P0.dat");
                            //snprintf(mc_filename_2,sizeof(mc_filename),"%s%s%d%s", mc_dir,"mcdata_",i,"_P0_0.dat");
                            for  (j=0;j<angle_procs;j++)
                            {
                                snprintf(mc_filename_2,sizeof(mc_filename),"%s%s%d%s%d%s", mc_dir,"mcdata_",i,"_P0_",j, ".dat");

                                if(( access( mc_filename, F_OK ) != -1 )  || ( access( mc_filename_2, F_OK ) != -1 ) )
                                {
                                    snprintf(mc_operation,sizeof(flash_prefix),"%s%s%s%d%s","exec rm ", mc_dir,"mcdata_",i,"_*.dat"); //prepares string to remove *.dat in mc_dir
                                    //printf("%s\n",mc_operation);
                                    system(mc_operation);

                                    snprintf(mc_operation,sizeof(flash_prefix),"%s%s%s%d%s","exec rm ", mc_dir,"mcdata_",i,"_*");
                                    system(mc_operation);
                                }

                            }
                        }
                        snprintf(mc_operation,sizeof(flash_prefix),"%s%s%s","exec rm ", mc_dir,"mcdata_PW_*.dat"); //prepares string to remove *.dat in mc_dir
                        system(mc_operation);

                        snprintf(mc_operation,sizeof(flash_prefix),"%s%s%s","exec rm ", mc_dir,"mcdata_PW.dat"); //prepares string to remove *.dat in mc_dir
                        system(mc_operation);

                        snprintf(mc_operation,sizeof(flash_prefix),"%s%s%s","exec rm ", mc_dir,"mc_chkpt_*.dat"); //prepares string to remove *.dat in mc_dir
                        system(mc_operation);

                        snprintf(mc_operation,sizeof(flash_prefix),"%s%s%s","exec rm ", mc_dir,"mc_output_*.log"); //prepares string to remove *.log in mc_dir
                        system(mc_operation);

                    }
                }
            }

            if ((RIKEN_SWITCH==1) && (dim_switch==1) && (framestart>=3000))
            {
                increment_inj=10; //when the frame ==3000 for RIKEN 3D hydro files, increment file numbers by 10 instead of by 1
                fps_modified=1;
            }
            else
            {
                increment_inj=1;
                fps_modified=fps;
            }

            dt_max=1.0/fps_modified;

            MPI_Barrier(angle_comm);
            snprintf(log_file,sizeof(log_file),"%s%s%d%s",mc_dir,"mc_output_", angle_id,".log" );
            printf("%s\n",log_file);
            fPtr=fopen(log_file, "a");

            printf( "Im Proc %d with angles %0.1lf-%0.1lf proc_frame_size is %d Starting on Frame: %d Injecting until %d scatt_framestart: %d\n", angle_id, theta_jmin_thread*180/M_PI, theta_jmax_thread*180/M_PI, proc_frame_size, framestart, frm2, scatt_framestart);

            fprintf(fPtr, "Im Proc %d with angles %0.1lf-%0.1lf  Starting on Frame: %d scatt_framestart: %d\n", angle_id, theta_jmin_thread*180/M_PI, theta_jmax_thread*180/M_PI, framestart, scatt_framestart);
            fflush(fPtr);

            free(frame_array);

            //for a checkpoint implementation, start from the last saved "frame" value and go to the saved "frm2" value

            //#pragma omp for

            for (frame=framestart;frame<=frm2;frame=frame+increment_inj)
            {

                if ((RIKEN_SWITCH==1) && (dim_switch==1) && (frame>=3000))
                {
                    increment_inj=10; //when the frame ==3000 for RIKEN 3D hydro files, increment file numbers by 10 instead of by 1
                    fps_modified=1;
                }
                else
                {
                    increment_inj=1;
                    fps_modified=fps;
                }

                 if (restrt=='r')
                 {
                    time_now=frame/fps; //for a checkpoint implmentation, load the saved "time_now" value when reading the ckeckpoint file otherwise calculate it normally
                 }

                //printf(">> mc.py: Working on Frame %d\n", frame);
                fprintf(fPtr,"Im Proc: %d with angles %0.1lf - %0.1lf Working on Frame: %d\n", angle_id, theta_jmin_thread*180/M_PI, theta_jmax_thread*180/M_PI, frame);
                fflush(fPtr);

                if (restrt=='r')
                {

                    if (dim_switch==0)
                    {
                        if (RIKEN_SWITCH==0)
                        {
                            //if using FLASH data for 2D
                        //put proper number at the end of the flash file
                        modifyFlashName(flash_file, flash_prefix, frame, dim_switch);

                        fprintf(fPtr,">> Im Proc: %d with angles %0.1lf-%0.1lf: Opening FLASH file %s\n",angle_id, theta_jmin_thread*180/M_PI, theta_jmax_thread*180/M_PI, flash_file);
                        fflush(fPtr);

                        readAndDecimate(flash_file, inj_radius, fps_modified, &xPtr,  &yPtr,  &szxPtr, &szyPtr, &rPtr,\
                                &thetaPtr, &velxPtr,  &velyPtr,  &densPtr,  &presPtr,  &gammaPtr,  &dens_labPtr, &tempPtr, &array_num, 1, min_r, max_r, min_theta, max_theta, fPtr);
                        }
                        else
                        {
                            //if using RIKEN hydro data for 2D szx becomes delta r szy becomes delta theta
                            readHydro2D(FILEPATH, frame, inj_radius, fps_modified, &xPtr,  &yPtr,  &szxPtr, &szyPtr, &rPtr,\
                                        &thetaPtr, &velxPtr,  &velyPtr,  &densPtr,  &presPtr,  &gammaPtr,  &dens_labPtr, &tempPtr, &array_num, 1, min_r, max_r, fPtr);
                            //fprintf(fPtr, "%d\n\n", array_num);
                        }
                    }
                    else
                    {
                        fprintf(fPtr,">> Im Proc: %d with angles %0.1lf-%0.1lf\n",angle_id, theta_jmin_thread*180/M_PI, theta_jmax_thread*180/M_PI);
                        fflush(fPtr);

                        read_hydro(FILEPATH, frame, inj_radius, &xPtr,  &yPtr, &zPtr,  &szxPtr, &szyPtr, &rPtr,\
                                   &thetaPtr, &phiPtr, &velxPtr,  &velyPtr, &velzPtr,  &densPtr,  &presPtr,  &gammaPtr,  &dens_labPtr, &tempPtr, &array_num, 1, min_r, max_r, fps_modified, fPtr);
                    }

                    //check for run type
                    if(strcmp(cyl, this_run)==0)
                    {
                        //printf("In cylindrical prep\n");
                        cylindricalPrep(gammaPtr, velxPtr, velyPtr, densPtr, dens_labPtr, presPtr, tempPtr, array_num);
                    }
                    else if (strcmp(sph, this_run)==0)
                    {
                        printf("In Spherical\n");
                        sphericalPrep(rPtr, xPtr, yPtr,gammaPtr, velxPtr, velyPtr, densPtr, dens_labPtr, presPtr, tempPtr, array_num , fPtr);
                    }

                    //determine where to place photons and how many should go in a given place
                    //for a checkpoint implmentation, dont need to inject photons, need to load photons' last saved data
                    fprintf(fPtr,">>  Proc: %d with angles %0.1lf-%0.1lf: Injecting photons\n",angle_id, theta_jmin_thread*180/M_PI, theta_jmax_thread*180/M_PI);
                    fflush(fPtr);

                    if (dim_switch==0)
                    {
                        photonInjection(&phPtr, &num_ph, inj_radius, ph_weight_suggest, min_photons, max_photons,spect, array_num, fps_modified, theta_jmin_thread, theta_jmax_thread, xPtr, yPtr, szxPtr, szyPtr,rPtr,thetaPtr, tempPtr, velxPtr, velyPtr,rng, RIKEN_SWITCH, fPtr );
                    }
                    else
                    {
                        photonInjection3D(&phPtr, &num_ph, inj_radius, ph_weight_suggest, min_photons, max_photons,spect, array_num, fps_modified, theta_jmin_thread, theta_jmax_thread, xPtr, yPtr, zPtr, szxPtr, szyPtr,rPtr,thetaPtr, phiPtr, tempPtr, velxPtr, velyPtr, velzPtr, rng, fPtr);

                    }

                    //printf("This many Photons: %d\n",num_ph); //num_ph is one more photon than i actually have

                    //for (i=0;i<num_ph;i++)
                    //    printf("%e,%e,%e \n",(phPtr+i)->r0, (phPtr+i)->r1, (phPtr+i)->r2 );

                }

                //scatter photons all the way thoughout the jet
                //for a checkpoint implmentation, start from the last saved "scatt_frame" value eh start_frame=frame or start_frame=cont_frame
                if (restrt=='r')
                {
                    scatt_framestart=frame; //have to make sure that once the inner loop is done and the outer loop is incrememnted by one the inner loop starts at that new value and not the one read by readCheckpoint()
                }

                for (scatt_frame=scatt_framestart;scatt_frame<=last_frm;scatt_frame=scatt_frame+increment_scatt)
                {
                    if ((RIKEN_SWITCH==1) && (dim_switch==1) && (scatt_frame>=3000))
                    {
                        increment_scatt=10; //when the frame ==3000 for RIKEN 3D hydro files, increment file numbers by 10 instead of by 1
                        fps_modified=1; //therefore dt between files become 1 second

                    }
                    else
                    {
                        increment_scatt=1;
                        fps_modified=fps;
                    }

                    dt_max=1.0/fps_modified; //if working with RIKEN files and scatt_frame>=3000 dt  is 1 second between each subsequent frame

                    fprintf(fPtr,">>\n");
                    fprintf(fPtr,">> Proc %d with angles %0.1lf-%0.1lf: Working on photons injected at frame: %d out of %d\n", angle_id, theta_jmin_thread*180/M_PI, theta_jmax_thread*180/M_PI,frame, frm2);
                    fprintf(fPtr,">> Proc %d with angles %0.1lf-%0.1lf: %s - Working on frame %d\n",angle_id, theta_jmin_thread*180/M_PI, theta_jmax_thread*180/M_PI, THISRUN, scatt_frame);
                    fprintf(fPtr,">> Proc %d with angles %0.1lf-%0.1lf: Opening file...\n", angle_id, theta_jmin_thread*180/M_PI, theta_jmax_thread*180/M_PI);
                    fflush(fPtr);

                    //set new seed to increase randomness?
                    gsl_rng_set(rng, gsl_rng_get(rng));

                    if (dim_switch==0)
                    {
                        if (RIKEN_SWITCH==0)
                        {
                            //put proper number at the end of the flash file
                            modifyFlashName(flash_file, flash_prefix, scatt_frame, dim_switch);
                            phMinMax(phPtr, num_ph, &min_r, &max_r, &min_theta, &max_theta, fPtr);
                            readAndDecimate(flash_file, inj_radius, fps_modified, &xPtr,  &yPtr,  &szxPtr, &szyPtr, &rPtr,\
                                    &thetaPtr, &velxPtr,  &velyPtr,  &densPtr,  &presPtr,  &gammaPtr,  &dens_labPtr, &tempPtr, &array_num, 0, min_r, max_r, min_theta, max_theta, fPtr);
                        }
                        else
                        {
                            phMinMax(phPtr, num_ph, &min_r, &max_r, &min_theta, &max_theta, fPtr);
                            //if using RIKEN hydro data for 2D szx becomes delta r szy becomes delta theta
                            readHydro2D(FILEPATH, scatt_frame, inj_radius, fps_modified, &xPtr,  &yPtr,  &szxPtr, &szyPtr, &rPtr,\
                                        &thetaPtr, &velxPtr,  &velyPtr,  &densPtr,  &presPtr,  &gammaPtr,  &dens_labPtr, &tempPtr, &array_num, 0, min_r, max_r, fPtr);

                        }
                    }
                    else
                    {
                        phMinMax(phPtr, num_ph, &min_r, &max_r, &min_theta, &max_theta, fPtr);
                        read_hydro(FILEPATH, scatt_frame, inj_radius, &xPtr,  &yPtr, &zPtr,  &szxPtr, &szyPtr, &rPtr,\
                                   &thetaPtr, &phiPtr, &velxPtr,  &velyPtr, &velzPtr,  &densPtr,  &presPtr,  &gammaPtr,  &dens_labPtr, &tempPtr, &array_num, 0, min_r, max_r, fps_modified, fPtr);
                    }
                    fprintf(fPtr, "Number of Flash Elements %d\n", array_num);

                    //check for run type
                    if(strcmp(cyl, this_run)==0)
                    {
                        //printf("In cylindrical prep\n");
                        cylindricalPrep(gammaPtr, velxPtr, velyPtr, densPtr, dens_labPtr, presPtr, tempPtr, array_num);
                    }
                    else if (strcmp(sph, this_run)==0)
                    {
                        sphericalPrep(rPtr, xPtr, yPtr,gammaPtr, velxPtr, velyPtr, densPtr, dens_labPtr, presPtr, tempPtr, array_num, fPtr );
                    }
                        //printf("The result of read and decimate are arrays with %d elements\n", array_num);

                    fprintf(fPtr,">> Proc %d with angles %0.1lf-%0.1lf: propagating and scattering %d photons\n",angle_id, theta_jmin_thread*180/M_PI, theta_jmax_thread*180/M_PI,num_ph);
                    fflush(fPtr);

                    frame_scatt_cnt=0;
                    find_nearest_grid_switch=1; // set to true so the function findNearestPropertiesAndMinMFP by default finds the index of the grid block closest to each photon since we just read in a file and the prior index is invalid

                    while (time_now<((scatt_frame+increment_scatt)/fps))
                    {
                        //if simulation time is less than the simulation time of the next frame, keep scattering in this frame
                        //for RIKEN hydro data, theres still 10 fps but after frame 3000, file increment is 10 not 1, therefore modify dt_max not fps

                        //go through each photon and find blocks closest to each photon and properties of those blocks to calulate mean free path
                        //and choose the photon with the smallest mfp and calculate the timestep

                        ph_scatt_index=findNearestPropertiesAndMinMFP(phPtr, num_ph, array_num, hydro_domain_x, hydro_domain_y, &time_step, xPtr,  yPtr, zPtr, szxPtr, szyPtr, velxPtr,  velyPtr,  velzPtr, dens_labPtr, tempPtr,\
                                                                      &ph_dens_labPtr, &ph_vxPtr, &ph_vyPtr, &ph_vzPtr, &ph_tempPtr, rng, dim_switch, find_nearest_grid_switch, RIKEN_SWITCH, fPtr);

                        find_nearest_grid_switch=0; //set to zero (false) since we do not absolutely need to refind the index, this makes the function findNearestPropertiesAndMinMFP just check if the photon is w/in the given grid box still

                        //fprintf(fPtr, "In main: %e, %d, %e, %e\n",((phPtr+ph_scatt_index)->num_scatt), ph_scatt_index, time_step, time_now);
                        //fflush(fPtr);

                         if (time_step<dt_max)
                        {

                            //update number of scatterings and time
                            ((phPtr+ph_scatt_index)->num_scatt)+=1;
                            frame_scatt_cnt+=1;
                            time_now+=time_step;

                            updatePhotonPosition(phPtr, num_ph, time_step, fPtr);

                            //scatter the photon
                            //fprintf(fPtr, "Passed Parameters: %e, %e, %e\n", (ph_vxPtr), (ph_vyPtr), (ph_tempPtr));

                            photonScatter( (phPtr+ph_scatt_index), (ph_vxPtr), (ph_vyPtr),ph_vzPtr, (ph_tempPtr), rng, dim_switch, fPtr );

                            if (frame_scatt_cnt%1000 == 0)
                            {
                                fprintf(fPtr,"Scattering Number: %d\n", frame_scatt_cnt);
                                fprintf(fPtr,"The local temp is: %e\n", (ph_tempPtr));
                                fprintf(fPtr,"Average photon energy is: %e\n", averagePhotonEnergy(phPtr, num_ph)); //write function to average over the photons p0 and then do (*3e10/1.6e-9)
                                fprintf(fPtr,"The last time step was: %e.\nThe time now is: %e\n", time_step,time_now);
                                fflush(fPtr);
                            }

                        }
                        else
                        {
                            time_now+=dt_max;

                            //for each photon update its position based on its momentum

                            updatePhotonPosition(phPtr, num_ph, dt_max, fPtr);
                        }

                        //printf("In main 2: %e, %d, %e, %e\n", ((phPtr+ph_scatt_index)->num_scatt), ph_scatt_index, time_step, time_now);

                    }

                    //get scattering statistics
                    phScattStats(phPtr, num_ph, &max_scatt, &min_scatt, &avg_scatt, &avg_r);

                    fprintf(fPtr,"The number of scatterings in this frame is: %d\n", frame_scatt_cnt);
                    fprintf(fPtr,"The last time step was: %e.\nThe time now is: %e\n", time_step,time_now);
                    fprintf(fPtr,"The maximum number of scatterings for a photon is: %d\nThe minimum number of scattering for a photon is: %d\n", max_scatt, min_scatt);
                    fprintf(fPtr,"The average number of scatterings thus far is: %lf\nThe average position of photons is %e\n", avg_scatt, avg_r);
                    fflush(fPtr);

                    printPhotons(phPtr, num_ph,  scatt_frame , frame, mc_dir, angle_id);
                    //exit(0);
                    fprintf(fPtr, ">> Proc %d with angles %0.1lf-%0.1lf: Making checkpoint file\n", angle_id, theta_jmin_thread*180/M_PI, theta_jmax_thread*180/M_PI);
                    fflush(fPtr);

                    fprintf(fPtr, " mc_dir: %s\nframe %d\nfrm2: %d\nscatt_frame: %d\n num_photon: %d\ntime_now: %e\nlast_frame: %d\n", mc_dir, frame, frm2, scatt_frame, num_ph, time_now, last_frm  );
                    fflush(fPtr);

                    saveCheckpoint(mc_dir, frame, frm2, scatt_frame, num_ph, time_now, phPtr, last_frm, angle_id, old_num_angle_procs);

                     if (dim_switch==1)
                    {
                        if (RIKEN_SWITCH==1)
                        {
                            free(zPtr);free(phiPtr);free(velzPtr);
                            zPtr=NULL; phiPtr=NULL; velzPtr=NULL;
                        }
                    }

                    free(xPtr);free(yPtr);free(szxPtr);free(szyPtr);free(rPtr);free(thetaPtr);free(velxPtr);free(velyPtr);free(densPtr);free(presPtr);
                    free(gammaPtr);free(dens_labPtr);free(tempPtr);
                    xPtr=NULL; yPtr=NULL;  rPtr=NULL;thetaPtr=NULL;velxPtr=NULL;velyPtr=NULL;densPtr=NULL;presPtr=NULL;gammaPtr=NULL;dens_labPtr=NULL;
                    szxPtr=NULL; szyPtr=NULL; tempPtr=NULL;
                }

                restrt='r';//set this to make sure that the next iteration of propogating photons doesnt use the values from the last reading of the checkpoint file
                free(phPtr);
                phPtr=NULL;

            }
            saveCheckpoint(mc_dir, frame, frm2, scatt_frame, 0, time_now, phPtr, last_frm, angle_id, old_num_angle_procs); //this is for processes using the old code that didnt restart efficiently
            fprintf(fPtr, "Process %d has completed the MC calculation.\n", angle_id);
            fflush(fPtr);

        }//end omp parallel inner section

        MPI_Barrier(angle_comm);

        //merge files from each worker thread within a directory
        {

             increment_scatt=1;
             file_count=0;

             //count number of files
             for (i=frm0;i<=last_frm;i=i+increment_scatt)
             {
                if ((RIKEN_SWITCH==1) && (dim_switch==1) && (i>=3000))
                {
                    increment_scatt=10; //when the frame ==3000 for RIKEN 3D hydro files, increment file numbers by 10 instead of by 1
                }
                file_count++;
             }

             //holds number of files for each process to merge
             MPI_Comm_size(angle_comm, &angle_procs); //to get the proper number of processes within the group
             MPI_Comm_rank(angle_comm, &angle_id); //reset the value of angle_id to what it should actualy be to properly distribute files to merge

             proc_frame_size=floor(file_count/ (float) angle_procs);
             frame_array=malloc(file_count*sizeof(int));
             proc_frame_array=malloc(angle_procs*sizeof(int)); //sets index of each proceesed acquired value
             element_num=malloc(angle_procs*sizeof(int));

             for (i=0;i<angle_procs;i++)
             {
                *(proc_frame_array+i)=i*proc_frame_size;
                *(element_num+i)=1;
             }

             //make vector with the files in order to pass them to each of the processes
             increment_scatt=1;
             file_count=0;
             for (i=frm0;i<=last_frm;i=i+increment_scatt)
             {
                if ((RIKEN_SWITCH==1) && (dim_switch==1) && (i>=3000))
                {
                    increment_scatt=10; //when the frame ==3000 for RIKEN 3D hydro files, increment file numbers by 10 instead of by 1
                }

                *(frame_array+file_count)=i ;
                file_count++;
                //printf("file_count: %d frame: %d\n",  file_count-1, *(frame_array+file_count-1));
             }
             //pass  first frame number that each rpocess should start to merge, can calulate the file it should merge until
             MPI_Scatterv(frame_array, element_num, proc_frame_array, MPI_INT, &frm0, 1, MPI_INT, 0, angle_comm);

             //fprintf(fPtr, "Value: last_frm: ,%d\n", file_count);
             //fflush(fPtr);

             //make sure all files get merged by giving the rest to the last process
             if (angle_id==angle_procs-1)
             {
                proc_frame_size=file_count-proc_frame_size*(angle_procs-1); //for last process take over the remaining number of files
             }
             //calculate what the last file the preocess should merge up to
             i=0;
             last_frm=frm0;
             while(i<proc_frame_size)
             {
                if ((RIKEN_SWITCH==1) && (dim_switch==1) && (last_frm>=3000))
                {
                    increment_scatt=10; //when the frame ==3000 for RIKEN 3D hydro files, increment file numbers by 10 instead of by 1
                }
                else
                {
                    increment_scatt=1;
                }

                last_frm+=increment_scatt;
                i++;
             }

            //if (angle_id==0)
            {
                //fprintf(fPtr, ">> Proc %d with angles %0.1lf-%0.1lf: Merging Files from %d to %d\n", angle_id, theta_jmin_thread*180/M_PI, theta_jmax_thread*180/M_PI, frm0, last_frm);
                fprintf(fPtr, ">> Proc %d with angles %0.1lf-%0.1lf: Merging Files from %d to %d\n", angle_id, theta_jmin_thread*180/M_PI, theta_jmax_thread*180/M_PI, frm0, last_frm);
                fflush(fPtr);

                dirFileMerge(mc_dir, frm0, last_frm, old_num_angle_procs, angle_id, dim_switch, RIKEN_SWITCH, fPtr);
            }
        }

        fprintf(fPtr, "Process %d has completed merging files.\n", angle_id);
        fflush(fPtr);

    fclose(fPtr);
    gsl_rng_free (rng);

    MPI_Finalize();
    //free(rng);
    //free(thread_theta);

	return 0;
}

max_stars_repo_name: outflows/MCRaT

max_issues_repo_name: outflows/MCRaT

max_forks_repo_name: outflows/MCRaT