mirror of
https://github.com/gwm17/catima.git
synced 2024-11-26 12:08:52 -05:00
integrator
This commit is contained in:
hrocho
parent
a2b3f4130d
commit
ac5a6d0048
|
|
@ -6,6 +6,7 @@ project(catima)
|
|||
option(CATIMA_PYTHON "compile the Catima python module(requires numpy and cython installed)" OFF)
|
||||
option(TESTS "build tests" OFF)
|
||||
option(EXAMPLES "build examples" ON)
|
||||
option(GSL_INTEGRATION "use GSL integration" ON)
|
||||
option(GENERATE_DATA "make data tables generator" OFF)
|
||||
option(THIN_TARGET_APPROXIMATION "thin target approximation" ON)
|
||||
option(DOCS "build documentation (requires doxygen)" OFF)
|
||||
|
|
|
|||
|
|
@ -2,5 +2,7 @@
|
|||
#define BUILD_CONFIG_H
|
||||
|
||||
#cmakedefine THIN_TARGET_APPROXIMATION
|
||||
#cmakedefine GSL_INTEGRATION
|
||||
|
||||
#endif
|
||||
|
||||
#endif
|
||||
|
|
|
|||
|
|
@ -631,7 +631,7 @@ double precalculated_lindhard_X(const Projectile &p){
|
|||
return v1+(dif*da/ls_coefficients::a_rel_increase);
|
||||
}
|
||||
|
||||
double dedx_rms(Projectile &p, Target &t, const Config &c){
|
||||
double dedx_variance(Projectile &p, Target &t, const Config &c){
|
||||
double zp_eff = z_effective(p,t,c);
|
||||
double gamma = gamma_from_T(p.T);
|
||||
double f = domega2dx_constant*pow(zp_eff,2)*t.Z*gamma*gamma;
|
||||
|
|
|
|||
|
|
@ -40,7 +40,7 @@ namespace catima{
|
|||
/**
|
||||
* returns energy loss straggling
|
||||
*/
|
||||
double dedx_rms(Projectile &p, Target &t, const Config &c=default_config);
|
||||
double dedx_variance(Projectile &p, Target &t, const Config &c=default_config);
|
||||
|
||||
/**
|
||||
* returns reduced energy loss unit for projectile-target combination
|
||||
|
|
|
|||
54
catima.cpp
54
catima.cpp
|
|
@ -253,12 +253,12 @@ std::vector<double> calculate_range(Projectile p, const Material &t, const Confi
|
|||
auto fdedx = [&](double x)->double{return 1.0/dedx(p,x,t);};
|
||||
|
||||
//calculate 1st point to have i-1 element ready for loop
|
||||
res = integratorGSL.Integrate(fdedx,Ezero,energy_table(0),int_eps_range,false);
|
||||
res = integrator.integrate(fdedx,Ezero,energy_table(0));
|
||||
res = p.A*res;
|
||||
values.push_back(res);
|
||||
|
||||
for(int i=1;i<max_datapoints;i++){
|
||||
res = integratorGSL.Integrate(fdedx,energy_table(i-1),energy_table(i),int_eps_range,false);
|
||||
res = integrator.integrate(fdedx,energy_table(i-1),energy_table(i));
|
||||
res = p.A*res;
|
||||
res += values[i-1];
|
||||
values.push_back(res);
|
||||
|
|
@ -278,11 +278,11 @@ std::vector<double> calculate_range_straggling(Projectile p, const Material &t,
|
|||
//return 1.0*domega2dx(p,x,t)/(de*de*de);
|
||||
//};
|
||||
//calculate 1st point to have i-1 element ready for loop
|
||||
res = integratorGSL.Integrate(function,Ezero,energy_table(0),int_eps_range_str,false);
|
||||
res = integrator.integrate(function,Ezero,energy_table(0));
|
||||
res = p.A*res;
|
||||
values.push_back(res);
|
||||
for(int i=1;i<max_datapoints;i++){
|
||||
res = integratorGSL.Integrate(function,energy_table(i-1),energy_table(i),int_eps_range_str,false);
|
||||
res = integrator.integrate(function,energy_table(i-1),energy_table(i));
|
||||
res = p.A*res;
|
||||
res += values[i-1];
|
||||
values.push_back(res);
|
||||
|
|
@ -297,11 +297,11 @@ std::vector<double> calculate_da2dx(Projectile p, const Material &t, const Confi
|
|||
values.reserve(max_datapoints);
|
||||
//auto function = [&](double x)->double{return p.A*da2dx(p,x,t)/dedx(p,x,t);};
|
||||
auto function = [&](double x)->double{return 1.0/dedx(p,x,t);};
|
||||
res = integratorGSL.Integrate(function,Ezero,energy_table(0),int_eps_ang_str,false);
|
||||
res = integrator.integrate(function,Ezero,energy_table(0));
|
||||
res = p.A*da2dx(p,energy_table(0),t)*res;
|
||||
values.push_back(res);
|
||||
for(int i=1;i<max_datapoints;i++){
|
||||
res = integratorGSL.Integrate(function,energy_table(i-1),energy_table(i),int_eps_ang_str,false);
|
||||
res = integrator.integrate(function,energy_table(i-1),energy_table(i));
|
||||
res = p.A*da2dx(p,energy_table(i),t)*res;
|
||||
res += values[i-1];
|
||||
values.push_back(res);
|
||||
|
|
@ -314,11 +314,11 @@ std::vector<double> calculate_tof(Projectile p, const Material &t, const Config
|
|||
std::vector<double> values;
|
||||
values.reserve(max_datapoints);
|
||||
auto function = [&](double x)->double{return 1.0/(dedx(p,x,t)*beta_from_T(x));};
|
||||
res = integratorGSL.Integrate(function,Ezero,energy_table(0),int_eps_tof,false);
|
||||
res = integrator.integrate(function,Ezero,energy_table(0));
|
||||
res = res*10.0*p.A/(c_light*t.density());
|
||||
values.push_back(res);
|
||||
for(int i=1;i<max_datapoints;i++){
|
||||
res = integratorGSL.Integrate(function,energy_table(i-1),energy_table(i),int_eps_tof,false);
|
||||
res = integrator.integrate(function,energy_table(i-1),energy_table(i));
|
||||
res = res*10.0*p.A/(c_light*t.density());
|
||||
res += values[i-1];
|
||||
values.push_back(res);
|
||||
|
|
@ -328,13 +328,9 @@ std::vector<double> calculate_tof(Projectile p, const Material &t, const Config
|
|||
|
||||
DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
|
||||
DataPoint dp(p,t,c);
|
||||
std::vector<double>range_values;
|
||||
range_values.reserve(max_datapoints);
|
||||
std::vector<double>range_straggling_values;
|
||||
range_straggling_values.reserve(max_datapoints);
|
||||
std::vector<double>angular_straggling_values;
|
||||
angular_straggling_values.reserve(max_datapoints);
|
||||
|
||||
dp.range.resize(max_datapoints);
|
||||
dp.range_straggling.resize(max_datapoints);
|
||||
dp.angular_variance.resize(max_datapoints);
|
||||
double dedxval;
|
||||
auto fdedx = [&](double x)->double{
|
||||
return 1.0/dedx(p,x,t);
|
||||
|
|
@ -346,32 +342,28 @@ DataPoint calculate_DataPoint(Projectile p, const Material &t, const Config &c){
|
|||
|
||||
double res;
|
||||
//calculate 1st point to have i-1 element ready for loop
|
||||
res = integratorGSL.Integrate(fdedx,Ezero,energy_table(0),int_eps_range,false);
|
||||
res = integrator.integrate(fdedx,Ezero,energy_table(0));
|
||||
res = p.A*res;
|
||||
range_values.push_back(res);
|
||||
dp.range[0] = res;
|
||||
res = da2dx(p,energy_table(0),t)*res;
|
||||
angular_straggling_values.push_back(res);
|
||||
dp.angular_variance[0] = res;
|
||||
|
||||
res = integratorGSL.Integrate(fomega,Ezero,energy_table(0),int_eps_range_str,false);
|
||||
res = integrator.integrate(fomega,Ezero,energy_table(0));
|
||||
res = p.A*res;
|
||||
range_straggling_values.push_back(res);
|
||||
dp.range_straggling[0]=res;
|
||||
|
||||
for(int i=1;i<max_datapoints;i++){
|
||||
res = p.A*integratorGSL.Integrate(fdedx,energy_table(i-1),energy_table(i),int_eps_range,false);
|
||||
range_values.push_back(res + range_values[i-1]);
|
||||
res = p.A*integrator.integrate(fdedx,energy_table(i-1),energy_table(i));
|
||||
dp.range[i] = res + dp.range[i-1];
|
||||
|
||||
res = da2dx(p,energy_table(i),t)*res;
|
||||
angular_straggling_values.push_back(res + angular_straggling_values[i-1]);
|
||||
dp.angular_variance[i] = res + dp.angular_variance[i-1];
|
||||
|
||||
res = integratorGSL.Integrate(fomega,energy_table(i-1),energy_table(i),int_eps_range_str,false);
|
||||
res = integrator.integrate(fomega,energy_table(i-1),energy_table(i));
|
||||
//res = integratorGSL.integrate(fomega,energy_table(i-1),energy_table(i));
|
||||
res = p.A*res;
|
||||
res += range_straggling_values[i-1];
|
||||
range_straggling_values.push_back(res);
|
||||
dp.range_straggling[i] = res + dp.range_straggling[i-1];
|
||||
}
|
||||
|
||||
dp.range = range_values;
|
||||
dp.range_straggling = range_straggling_values;
|
||||
dp.angular_variance = angular_straggling_values;
|
||||
return dp;
|
||||
}
|
||||
|
||||
|
|
@ -380,7 +372,7 @@ double calculate_tof_from_E(Projectile p, double Eout, const Material &t, const
|
|||
//double beta_in = beta_from_T(p.T);
|
||||
//double beta_out = beta_from_T(Eout);
|
||||
auto function = [&](double x)->double{return 1.0/(dedx(p,x,t)*beta_from_T(x));};
|
||||
res = integratorGSL.Integrate(function,Eout,p.T,int_eps_tof,false);
|
||||
res = integrator.integrate(function,Eout,p.T);
|
||||
res = res*10.0*p.A/(c_light*t.density());
|
||||
return res;
|
||||
}
|
||||
|
|
|
|||
126
integrator.cpp
126
integrator.cpp
|
|
@ -4,23 +4,28 @@
|
|||
|
||||
namespace catima{
|
||||
|
||||
IntegratorGSL integratorGSL;
|
||||
integrator_type integrator;
|
||||
IntegratorGSL integratorGSL(true);
|
||||
|
||||
double funcwrapper3(double x, void *_c){
|
||||
std::function<double(double)> *f = (std::function<double(double)> *)_c;
|
||||
return (*f)(x);
|
||||
}
|
||||
|
||||
IntegratorGSL::IntegratorGSL(){
|
||||
gsl_set_error_handler_off();
|
||||
w=gsl_integration_workspace_alloc(500);
|
||||
IntegratorGSL::IntegratorGSL(bool adapt):adaptive(adapt){
|
||||
//gsl_set_error_handler_off();
|
||||
if(adaptive){
|
||||
w=gsl_integration_workspace_alloc(100);
|
||||
}
|
||||
};
|
||||
|
||||
|
||||
IntegratorGSL::~IntegratorGSL(){
|
||||
gsl_integration_workspace_free(w);
|
||||
if(adaptive){
|
||||
gsl_integration_workspace_free(w);
|
||||
}
|
||||
};
|
||||
|
||||
double IntegratorGSL::Integrate(std::function<double(double)> f, double _min, double _max, double prec, bool adaptive){
|
||||
double IntegratorGSL::integrate(std::function<double(double)> f, double _min, double _max, double prec){
|
||||
gsl_function F;
|
||||
F.function = funcwrapper3;
|
||||
F.params = &f;
|
||||
|
|
@ -32,10 +37,111 @@ namespace catima{
|
|||
#ifdef USE_THREADS
|
||||
std::lock_guard<std::mutex> lock(integration_mutex);
|
||||
#endif
|
||||
gsl_integration_qag(&F,_min,_max,prec,prec,500,6,w,&result,&error);
|
||||
gsl_integration_qag(&F,_min,_max,1e-6,prec,100,6,w,&result,&error);
|
||||
}
|
||||
else
|
||||
gsl_integration_qng(&F,_min,_max,0,prec,&result,&error,&num);
|
||||
else{
|
||||
gsl_integration_qng(&F,_min,_max,1e-6,prec,&result,&error,&num);
|
||||
}
|
||||
return result;
|
||||
};
|
||||
|
||||
//////////////////// GaussLegendre weights /////////////////////
|
||||
// order = 4
|
||||
template<>
|
||||
std::array<double,2> GaussLegendreIntegration<4>::
|
||||
x = {0.3399810435848562648026658,0.8611363115940525752239465};
|
||||
|
||||
template<>
|
||||
std::array<double,2> GaussLegendreIntegration<4>::
|
||||
w = {0.6521451548625461426269361,0.3478548451374538573730639};
|
||||
|
||||
// order = 6
|
||||
template<>
|
||||
std::array<double,3> GaussLegendreIntegration<6>::
|
||||
x = {0.2386191860831969086305017,0.6612093864662645136613996,0.9324695142031520278123016};
|
||||
|
||||
template<>
|
||||
std::array<double,3> GaussLegendreIntegration<6>::
|
||||
w = {0.4679139345726910473898703,0.3607615730481386075698335,0.1713244923791703450402961};
|
||||
|
||||
// order = 8
|
||||
template<>
|
||||
std::array<double,4> GaussLegendreIntegration<8>::
|
||||
x ={0.1834346424956498049394761,0.5255324099163289858177390,0.7966664774136267395915539,0.9602898564975362316835609};
|
||||
|
||||
template<>
|
||||
std::array<double,4> GaussLegendreIntegration<8>::
|
||||
w = {0.3626837833783619829651504,0.3137066458778872873379622,0.2223810344533744705443560,0.1012285362903762591525314};
|
||||
|
||||
// order = 10
|
||||
template<>
|
||||
std::array<double,5> GaussLegendreIntegration<10>::
|
||||
x = {0.1488743389816312108848260,0.4333953941292471907992659,0.6794095682990244062343274,0.8650633666889845107320967,0.9739065285171717200779640};
|
||||
|
||||
template<>
|
||||
std::array<double,5> GaussLegendreIntegration<10>::
|
||||
w = {0.2955242247147528701738930,0.2692667193099963550912269,0.2190863625159820439955349,0.1494513491505805931457763,0.0666713443086881375935688};
|
||||
|
||||
// order = 12
|
||||
template<>
|
||||
std::array<double,6> GaussLegendreIntegration<12>::
|
||||
x = {0.1252334085114689154724414,0.3678314989981801937526915,0.5873179542866174472967024,0.7699026741943046870368938,0.9041172563704748566784659,0.9815606342467192506905491};
|
||||
|
||||
template<>
|
||||
std::array<double,6> GaussLegendreIntegration<12>::
|
||||
w = {0.2491470458134027850005624,0.2334925365383548087608499,0.2031674267230659217490645,0.1600783285433462263346525,0.1069393259953184309602547,0.0471753363865118271946160};
|
||||
|
||||
// order = 14
|
||||
template<>
|
||||
std::array<double,7> GaussLegendreIntegration<14>::
|
||||
x = {0.1080549487073436620662447,0.3191123689278897604356718,0.5152486363581540919652907,0.6872929048116854701480198,0.8272013150697649931897947,0.9284348836635735173363911,0.9862838086968123388415973};
|
||||
template<>
|
||||
std::array<double,7> GaussLegendreIntegration<14>::
|
||||
w = {0.2152638534631577901958764,0.2051984637212956039659241,0.1855383974779378137417166,0.1572031671581935345696019,0.1215185706879031846894148,0.0801580871597602098056333,0.0351194603317518630318329};
|
||||
|
||||
|
||||
// order = 16
|
||||
template<>
|
||||
std::array<double,8> GaussLegendreIntegration<16>::
|
||||
x = {0.0950125098376374401853193,0.2816035507792589132304605,0.4580167776572273863424194,0.6178762444026437484466718,0.7554044083550030338951012,0.8656312023878317438804679,0.9445750230732325760779884,0.9894009349916499325961542};
|
||||
|
||||
template<>
|
||||
std::array<double,8> GaussLegendreIntegration<16>::
|
||||
w = {0.1894506104550684962853967,0.1826034150449235888667637,0.1691565193950025381893121,0.1495959888165767320815017,0.1246289712555338720524763,0.0951585116824927848099251,0.0622535239386478928628438,0.0271524594117540948517806};
|
||||
|
||||
// order = 20
|
||||
template<>
|
||||
std::array<double,10> GaussLegendreIntegration<20>::
|
||||
x = {0.0765265211334973337546404,0.2277858511416450780804962,0.3737060887154195606725482,0.5108670019508270980043641,0.6360536807265150254528367,0.7463319064601507926143051,0.8391169718222188233945291,0.9122344282513259058677524,0.9639719272779137912676661,0.9931285991850949247861224};
|
||||
template<>
|
||||
std::array<double,10> GaussLegendreIntegration<20>::
|
||||
w = {0.1527533871307258506980843,0.1491729864726037467878287,0.1420961093183820513292983,0.1316886384491766268984945,0.1181945319615184173123774,0.1019301198172404350367501,0.0832767415767047487247581,0.0626720483341090635695065,0.0406014298003869413310400,0.0176140071391521183118620};
|
||||
|
||||
// order = 32
|
||||
template<>
|
||||
std::array<double,16> GaussLegendreIntegration<32>::
|
||||
x = {0.0483076656877383162348126,0.1444719615827964934851864,0.2392873622521370745446032,0.3318686022821276497799168,0.4213512761306353453641194,0.5068999089322293900237475,0.5877157572407623290407455,0.6630442669302152009751152,0.7321821187402896803874267,0.7944837959679424069630973,0.8493676137325699701336930,0.8963211557660521239653072,0.9349060759377396891709191,0.9647622555875064307738119,0.9856115115452683354001750,0.9972638618494815635449811};
|
||||
|
||||
template<>
|
||||
std::array<double,16> GaussLegendreIntegration<32>::
|
||||
w = {0.0965400885147278005667648,0.0956387200792748594190820,0.0938443990808045656391802,0.0911738786957638847128686,0.0876520930044038111427715,0.0833119242269467552221991,0.0781938957870703064717409,0.0723457941088485062253994,0.0658222227763618468376501,0.0586840934785355471452836,0.0509980592623761761961632,0.0428358980222266806568786,0.0342738629130214331026877,0.0253920653092620594557526,0.0162743947309056706051706,0.0070186100094700966004071};
|
||||
|
||||
//order = 64
|
||||
template<>
|
||||
std::array<double,32> GaussLegendreIntegration<64>::
|
||||
x = {0.0243502926634244325089558,0.0729931217877990394495429,0.1214628192961205544703765,0.1696444204239928180373136,0.2174236437400070841496487,0.2646871622087674163739642,0.3113228719902109561575127,0.3572201583376681159504426,0.4022701579639916036957668,0.4463660172534640879849477,0.4894031457070529574785263,0.5312794640198945456580139,0.5718956462026340342838781,0.6111553551723932502488530,0.6489654712546573398577612,0.6852363130542332425635584,0.7198818501716108268489402,0.7528199072605318966118638,0.7839723589433414076102205,0.8132653151227975597419233,0.8406292962525803627516915,0.8659993981540928197607834,0.8893154459951141058534040,0.9105221370785028057563807,0.9295691721319395758214902,0.9464113748584028160624815,0.9610087996520537189186141,0.9733268277899109637418535,0.9833362538846259569312993,0.9910133714767443207393824,0.9963401167719552793469245,0.9993050417357721394569056};
|
||||
|
||||
template<>
|
||||
std::array<double,32> GaussLegendreIntegration<64>::
|
||||
w = {0.0486909570091397203833654,0.0485754674415034269347991,0.0483447622348029571697695,0.0479993885964583077281262,0.0475401657148303086622822,0.0469681828162100173253263,0.0462847965813144172959532,0.0454916279274181444797710,0.0445905581637565630601347,0.0435837245293234533768279,0.0424735151236535890073398,0.0412625632426235286101563,0.0399537411327203413866569,0.0385501531786156291289625,0.0370551285402400460404151,0.0354722132568823838106931,0.0338051618371416093915655,0.0320579283548515535854675,0.0302346570724024788679741,0.0283396726142594832275113,0.0263774697150546586716918,0.0243527025687108733381776,0.0222701738083832541592983,0.0201348231535302093723403,0.0179517157756973430850453,0.0157260304760247193219660,0.0134630478967186425980608,0.0111681394601311288185905,0.0088467598263639477230309,0.0065044579689783628561174,0.0041470332605624676352875,0.0017832807216964329472961};
|
||||
|
||||
//order = 128
|
||||
template<>
|
||||
std::array<double,64> GaussLegendreIntegration<128>::
|
||||
x = {0.0122236989606157641980521,0.0366637909687334933302153,0.0610819696041395681037870,0.0854636405045154986364980,0.1097942311276437466729747,0.1340591994611877851175753,0.1582440427142249339974755,0.1823343059853371824103826,0.2063155909020792171540580,0.2301735642266599864109866,0.2538939664226943208556180,0.2774626201779044028062316,0.3008654388776772026671541,0.3240884350244133751832523,0.3471177285976355084261628,0.3699395553498590266165917,0.3925402750332674427356482,0.4149063795522750154922739,0.4370245010371041629370429,0.4588814198335521954490891,0.4804640724041720258582757,0.5017595591361444642896063,0.5227551520511754784539479,0.5434383024128103634441936,0.5637966482266180839144308,0.5838180216287630895500389,0.6034904561585486242035732,0.6228021939105849107615396,0.6417416925623075571535249,0.6602976322726460521059468,0.6784589224477192593677557,0.6962147083695143323850866,0.7135543776835874133438599,0.7304675667419088064717369,0.7469441667970619811698824,0.7629743300440947227797691,0.7785484755064119668504941,0.7936572947621932902433329,0.8082917575079136601196422,0.8224431169556438424645942,0.8361029150609068471168753,0.8492629875779689691636001,0.8619154689395484605906323,0.8740527969580317986954180,0.8856677173453972174082924,0.8967532880491581843864474,0.9073028834017568139214859,0.9173101980809605370364836,0.9267692508789478433346245,0.9356743882779163757831268,0.9440202878302201821211114,0.9518019613412643862177963,0.9590147578536999280989185,0.9656543664319652686458290,0.9717168187471365809043384,0.9771984914639073871653744,0.9820961084357185360247656,0.9864067427245862088712355,0.9901278184917343833379303,0.9932571129002129353034372,0.9957927585349811868641612,0.9977332486255140198821574,0.9990774599773758950119878,0.9998248879471319144736081};
|
||||
|
||||
template<>
|
||||
std::array<double,64> GaussLegendreIntegration<128>::
|
||||
w = {0.0244461801962625182113259,0.0244315690978500450548486,0.0244023556338495820932980,0.0243585572646906258532685,0.0243002001679718653234426,0.0242273192228152481200933,0.0241399579890192849977167,0.0240381686810240526375873,0.0239220121367034556724504,0.0237915577810034006387807,0.0236468835844476151436514,0.0234880760165359131530253,0.0233152299940627601224157,0.0231284488243870278792979,0.0229278441436868469204110,0.0227135358502364613097126,0.0224856520327449668718246,0.0222443288937997651046291,0.0219897106684604914341221,0.0217219495380520753752610,0.0214412055392084601371119,0.0211476464682213485370195,0.0208414477807511491135839,0.0205227924869600694322850,0.0201918710421300411806732,0.0198488812328308622199444,0.0194940280587066028230219,0.0191275236099509454865185,0.0187495869405447086509195,0.0183604439373313432212893,0.0179603271850086859401969,0.0175494758271177046487069,0.0171281354231113768306810,0.0166965578015892045890915,0.0162550009097851870516575,0.0158037286593993468589656,0.0153430107688651440859909,0.0148731226021473142523855,0.0143943450041668461768239,0.0139069641329519852442880,0.0134112712886163323144890,0.0129075627392673472204428,0.0123961395439509229688217,0.0118773073727402795758911,0.0113513763240804166932817,0.0108186607395030762476596,0.0102794790158321571332153,0.0097341534150068058635483,0.0091830098716608743344787,0.0086263777986167497049788,0.0080645898904860579729286,0.0074979819256347286876720,0.0069268925668988135634267,0.0063516631617071887872143,0.0057726375428656985893346,0.0051901618326763302050708,0.0046045842567029551182905,0.0040162549837386423131943,0.0034255260409102157743378,0.0028327514714579910952857,0.0022382884309626187436221,0.0016425030186690295387909,0.0010458126793403487793129,0.0004493809602920903763943};
|
||||
|
||||
}
|
||||
|
|
|
|||
63
integrator.h
63
integrator.h
|
|
@ -16,8 +16,10 @@
|
|||
|
||||
#ifndef INTEGRATOR_H
|
||||
#define INTEGRATOR_H
|
||||
#include "catima/build_config.h"
|
||||
#include "gsl/gsl_integration.h"
|
||||
#include <functional>
|
||||
#include <array>
|
||||
#ifdef USE_THREADS
|
||||
#include <mutex>
|
||||
#endif
|
||||
|
|
@ -27,15 +29,15 @@ namespace catima{
|
|||
/// helper class to integrate functions using the GSL library
|
||||
class IntegratorGSL{
|
||||
public:
|
||||
IntegratorGSL();
|
||||
IntegratorGSL(bool adapt=true);
|
||||
~IntegratorGSL();
|
||||
|
||||
double Integrate(std::function<double(double)> f, double min, double max, double precision, bool adaptive=true);
|
||||
double integrate(std::function<double(double)> f, double min, double max, double precision=0.001);
|
||||
|
||||
private:
|
||||
gsl_integration_workspace *w;
|
||||
gsl_integration_workspace *w;
|
||||
bool adaptive;
|
||||
double error;
|
||||
double precision;
|
||||
double result;
|
||||
double min;
|
||||
double max;
|
||||
|
|
@ -44,8 +46,59 @@ class IntegratorGSL{
|
|||
#endif
|
||||
};
|
||||
|
||||
extern IntegratorGSL integratorGSL;
|
||||
// built in integrator
|
||||
template<int order>
|
||||
class GaussLegendreIntegration{
|
||||
public:
|
||||
template<typename F>
|
||||
double integrate(F f, double a, double b) const;
|
||||
template<typename F>
|
||||
double operator()(F f, double a, double b) const {return integrate(f, a, b);}
|
||||
double get_w(int i) const {return w[i];}
|
||||
double get_x(int i) const {return x[i];}
|
||||
int n() const {return order;}
|
||||
std::array<double,order> get_points(double a = -1.0, double b = 1.0)const;
|
||||
|
||||
public:
|
||||
static std::array<double,order/2> w;
|
||||
static std::array<double,order/2> x;
|
||||
};
|
||||
|
||||
template<int order>
|
||||
template<typename F>
|
||||
double GaussLegendreIntegration<order>::integrate(F f, double a, double b) const{
|
||||
double res=0.0;
|
||||
double p = 0.5*(b-a);
|
||||
double q = 0.5*(b+a);
|
||||
|
||||
for(int i=0;i<order/2;i++){
|
||||
res += w[i] * (f(p*x[i] + q) + f(-p*x[i] + q));
|
||||
}
|
||||
return p*res;
|
||||
}
|
||||
|
||||
template<int order>
|
||||
std::array<double,order> GaussLegendreIntegration<order>::get_points(double a, double b)const{
|
||||
std::array<double,order> points;
|
||||
double p = 0.5*(b-a);
|
||||
double q = 0.5*(b+a);
|
||||
|
||||
int num = (order/2);
|
||||
for(int i=0;i< num;i++){
|
||||
points[num-i-1] = -p*x[i] + q;
|
||||
points[num+i] = p*x[i] + q;
|
||||
}
|
||||
return points;
|
||||
}
|
||||
|
||||
#ifdef GSL_INTEGRATION
|
||||
using integrator_type = IntegratorGSL;
|
||||
#else
|
||||
using integrator_type = GaussLegendreIntegration<8>;
|
||||
#endif
|
||||
|
||||
extern integrator_type integrator;
|
||||
extern IntegratorGSL integratorGSL;
|
||||
}
|
||||
|
||||
#endif
|
||||
|
|
|
|||
|
|
@ -12,27 +12,28 @@ namespace catima{
|
|||
|
||||
Material get_compound(material m){
|
||||
switch(m){
|
||||
case material::Plastics: return Material({{0,1,10},{0,6,9}},1.032);
|
||||
case material::Air: return Material({{0,7,0.781},{0,8,0.2095},{0,18,0.0095}},0.0012);
|
||||
case material::Plastics: return Material({{0,1,0.085},{0,6,0.915}},1.032);
|
||||
case material::Air: return Material({{0,7,0.755267},{0,8,0.231781},{0,18,0.012827},{0,6,0.000124}},0.001205);
|
||||
case material::CH2: return Material({{0,6,1},{0,1,2}},0.94);
|
||||
case material::LH2: return Material({{0,1,1}},0.0708);
|
||||
case material::LD2: return Material({{2.01355,1,1}},0.162);
|
||||
case material::LD2: return Material({{2.0141,1,1}},0.162);
|
||||
case material::Water: return Material({{0,1,2},{0,8,1}},1);
|
||||
case material::Diamond: return Material({{0,6,1}},3.52);
|
||||
case material::Glass: return Material({{0,14,1},{0,8,2}},2.4);
|
||||
case material::ALMG3: return Material({{0,13,1}},2.67);
|
||||
case material::ArCO2_30: return Material({{0,18,7},{0,8,6},{0,6,3}},0.00171);
|
||||
case material::Glass: return Material({{0,14,0.37722},{0,8,0.539562},{0,5,0.040064},{0,11,0.028191},{0,13,0.011644},{0,19,0.003321}},2.4);
|
||||
case material::ALMG3: return Material({{0,13,0.97},{0,12,0.03}},2.67);
|
||||
case material::ArCO2_30: return Material({{0,18,0.679},{0,8,0.2321},{0,6,0.0889}},0.00171);
|
||||
case material::CF4: return Material({{0,6,1},{0,9,4}},0.00372);
|
||||
case material::Isobutane: return Material({{0,6,4},{0,1,10}},0.00251);
|
||||
case material::Kapton: return Material({{0,1,10},{0,6,22},{0,7,2},{0,8,5}},1.42);
|
||||
case material::Mylar: return Material({{0,6,5},{0,1,4},{0,8,2}},1.38);
|
||||
case material::Kapton: return Material({{0,1,0.026362},{0,6,0.691133},{0,7,0.07327},{0,8,0.209235}},1.42);
|
||||
case material::Mylar: return Material({{0,6,10},{0,1,8},{0,8,4}},1.38);
|
||||
case material::NaF: return Material({{0,11,1},{0,9,1}},2.56);
|
||||
case material::P10: return Material({{0,18,9},{0,6,1},{0,1,4}},0.00166);
|
||||
case material::PolyOlefin: return Material({{0,8,10},{0,1,16}},0.9);
|
||||
case material::P10: return Material({{0,18,0.64171},{0,6,0.268247},{0,1,0.09004}},0.00166);
|
||||
case material::Polyolefin: return Material({{0,8,8},{0,1,16}},0.9);
|
||||
case material::CmO2: return Material({{0,96,1},{0,8,2}},12);
|
||||
case material::Suprasil: return Material({{0,14,1},{0,8,2}},2.2);
|
||||
case material::HAVAR: return Material({{0,27,42},{0,24,20},{0,28,13},{0,74,3},{0,42,2},{0,26,20}},8.3);
|
||||
case material::Steel: return Material({{0,26,74},{0,24,18},{0,28,8}},8);
|
||||
case material::Suprasil: return Material({{0,14,0.37722},{0,8,0.539562},{0,5,0.040064},{0,11,0.028191},{0,13,0.011644},{0,19,0.003321}},2.2);
|
||||
case material::HAVAR: return Material({{0,27,0.403228},{0,24,0.169412},{0,28,0.124301},{0,74,0.089847},{0,42,0.031259},{0,26,0.181952}},8.3);
|
||||
case material::Steel: return Material({{0,26,0.74621},{0,24,0.169},{0,28,0.08479}},8);
|
||||
case material::CH4: return Material({{0,1,4},{0,6,1}},0.0006);
|
||||
default:break;
|
||||
}
|
||||
return Material();
|
||||
|
|
|
|||
|
|
@ -21,11 +21,12 @@ namespace catima{
|
|||
Mylar = 214,
|
||||
NaF = 215,
|
||||
P10 = 216,
|
||||
PolyOlefin = 217,
|
||||
Polyolefin = 217,
|
||||
CmO2 = 218,
|
||||
Suprasil = 219,
|
||||
HAVAR = 220,
|
||||
Steel = 221
|
||||
Steel = 221,
|
||||
CH4 = 222
|
||||
};
|
||||
|
||||
Material get_compound(material m);
|
||||
|
|
|
|||
19
storage.h
19
storage.h
|
|
@ -29,25 +29,6 @@ namespace catima{
|
|||
|
||||
enum interpolation_t {cspline, linear};
|
||||
|
||||
//inline double energy_function( int i ) { return exp(M_LN10*(logEmin + ((double)i)*(logEmax-logEmin)/(max_datapoints - 1.0))); }
|
||||
|
||||
//enum DataType{TYPE_RANGE,TYPE_LS};
|
||||
/*
|
||||
template<int N>
|
||||
struct EnergyTable{
|
||||
constexpr EnergyTable(double logmin, double logmax):values(),step(0.0),num(N){
|
||||
step = (logmax-logmin)/(N - 1.0);
|
||||
for(auto i=0;i<N;i++){
|
||||
values[i]=exp(M_LN10*(logmin + ((double)i)*step));
|
||||
}
|
||||
}
|
||||
double operator()(int i)const{return values[i];}
|
||||
double values[N];
|
||||
double step;
|
||||
std::size_t num;
|
||||
};
|
||||
*/
|
||||
|
||||
template<int N>
|
||||
struct EnergyTable{
|
||||
EnergyTable(double logmin, double logmax):values(),step(0.0),num(N){
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user