#include <iostream> #include "catima/catima.h" #include "catima/nucdata.h" #include <functional> #include <fstream> #include <sstream> using namespace catima; using std::cout; using std::endl; //// configuration //// constexpr double logEmin_forLS = 0.0; constexpr double logEmax_forLS = catima::logEmax; constexpr double num_LS_datapoints = 200; ///////////////////////////////////////// inline double ls_energy_table( int i ) { return exp(M_LN10*(logEmin_forLS + ((double)i)*(logEmax_forLS-logEmin_forLS)/(num_LS_datapoints - 1.0))); } int main(){ double r; double r2; double dx; double x1,x2; #define LS_COEFFICENTS #ifdef LS_COEFFICENTS std::ofstream file("generated_LS_coeff.h",std::ofstream::out); catima::Projectile p; std::vector<double> ve; std::stringstream sv1; std::stringstream sv2; std::stringstream sx1; std::stringstream sx2; file<<"//LS precalculated LS coefficients \n "; file<<"#include \"catima/storage.h\""; file<<"\n\n"; file<<"#define LS_NUM_ENERGY_POINTS "<<num_LS_datapoints<<"\n"; file<<"#define LS_MAX_Z "<<ATOMIC_WEIGHT_MAXZ<<"\n"; file<<"namespace catima{\n\n"; file<<"namespace ls_coefficients{\n\n"; file<<"// relative difference of mass number for 2nd mass point ls_coefficients_ahi \n\n"; file<<"constexpr double a_rel_increase=0.05;\n"; file<<"constexpr double logEmin = "<<logEmin_forLS<<";\n"; file<<"constexpr double logEmax = "<<logEmax_forLS<<";\n"; file<<"// energy points array in MeV/u \n"; //file<<"constexpr EnergyTable<"<<num_LS_datapoints<<"> ls_energy_table("<<logEmin_forLS<<","<<logEmax_forLS<<");\n"; file<<"EnergyTable<"<<num_LS_datapoints<<"> ls_energy_table("<<logEmin_forLS<<","<<logEmax_forLS<<");\n"; /* file<<"double ls_energy_points["<<num_LS_datapoints<<"]={"; for(int i=0;i<num_LS_datapoints;i++){ if(i!=0)file<<","; file<<ls_energy_table(i); } file<<"};\n"; */ sv1<<"{"; sv2<<"{"; sx1<<"{"; sx2<<"{"; for(int z=1;z<ATOMIC_WEIGHT_MAXZ-1;z++){ cout<<z<<"/"<<ATOMIC_WEIGHT_MAXZ-2<<endl; p.Z = z; if(z!=1)sv1<<","; if(z!=1)sv2<<","; if(z!=1)sx1<<","; if(z!=1)sx2<<","; sv1<<"\n\t{"; sv2<<"\n\t{"; sx1<<"\n\t{"; sx2<<"\n\t{"; for(int i=0;i<num_LS_datapoints;i++){ p.T = ls_energy_table(i); p.A = element_atomic_weight((int)p.Z); r = bethek_lindhard(p); x1 = bethek_lindhard_X(p); dx = p.A*0.05; p.A +=dx; r2 = bethek_lindhard(p); x2 = bethek_lindhard_X(p); if(i!=0)sv1<<","; sv1<<r; if(i!=0)sv2<<","; sv2<<r2; if(i!=0)sx1<<","; sx1<<x1; if(i!=0)sx2<<","; sx2<<x2; } sv1<<"}"; sv2<<"}"; sx1<<"}"; sx2<<"}"; } sv1<<"\n};"; sv2<<"\n};"; sx1<<"\n};"; sx2<<"\n};"; file<<"\n"; file<<"//arrays dimensions are [z][energy], z=1 starts from index=0\n"; file<<"\n//LS coefficient for A=atomic weight\n"; file<<"double ls_coefficients_a[]["<<num_LS_datapoints<<"]=\n"; file<<sv1.str(); file<<"\n"; file<<"\n//LS coefficient for A=atomic weight * 1.05\n"; file<<"double ls_coefficients_ahi[]["<<num_LS_datapoints<<"]=\n"; file<<sv2.str(); file<<"\n"; file<<"\n"; file<<"//arrays dimensions are [z][energy], z=1 starts from index=0\n"; file<<"\n//LS X coefficient (dE straggling) for A=atomic weight\n"; file<<"double ls_X_coefficients_a[]["<<num_LS_datapoints<<"]=\n"; file<<sx1.str(); file<<"\n"; file<<"\n//LS X coefficient for A=atomic weight * 1.05\n"; file<<"double ls_X_coefficients_ahi[]["<<num_LS_datapoints<<"]=\n"; file<<sx2.str(); file<<"\n"; file<<"\n}\n"; // end of ls_coefficient namespace file<<"\n}\n"; // end of catima namespace file.close(); #endif return 0; }