From be74892e11ec969facc91614ba2827cf6dd998c8 Mon Sep 17 00:00:00 2001 From: hrocho Date: Tue, 4 Aug 2020 17:35:33 +0200 Subject: [PATCH] backup --- CMakeLists.txt | 2 +- calculations.cpp | 3 +- catima.cpp | 33 ++- integrator.cpp | 2 + integrator.h | 132 +---------- libs/gkq_integrator.h | 440 ++++++++++++++++++++++++++++++++++++ libs/glq_integrator.h | 356 +++++++++++++++++++++++++++++ structures.h | 1 + tests/test_calculations.cpp | 13 +- 9 files changed, 846 insertions(+), 136 deletions(-) create mode 100644 libs/gkq_integrator.h create mode 100644 libs/glq_integrator.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 855b33e..4500831 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -72,7 +72,7 @@ if(GLOBAL) file(GLOB GLOBAL_SOURCES global/*.c) LIST (APPEND SOURCES ${GLOBAL_SOURCES}) endif(GLOBAL) -file(GLOB HEADERS *.h) +file(GLOB HEADERS *.h libs/*.h) add_library(catima ${SOURCES}) set_target_properties(catima PROPERTIES diff --git a/calculations.cpp b/calculations.cpp index a47c20f..bb71d0e 100644 --- a/calculations.cpp +++ b/calculations.cpp @@ -447,10 +447,11 @@ double sezi_dedx_e(const Projectile &p, const Material &mat, const Config &c){ double w; double sum=0.0; bool use95 = c.low_energy == low_energy_types::srim_95; + double T = p.T; for(int i=0;idouble{ double range = range_spline(T); - double tt = range - range_spline(x); - double t0 = std::min(range, t.thickness()); + double tt = range - range_spline(x); + double t0 = std::min(range, t.thickness()); return (tt-t0)*(tt-t0)*angular_variance_spline.derivative(x); }; - res.sigma_x = integrator.integrate(fx2,res.Eout, res.Ein)/t.density()/t.density(); + + + res.sigma_x = integrator_adaptive.integrate(fx2,res.Eout, res.Ein,1e-3, 1e-3,4)/t.density()/t.density(); + //res.sigma_x = integrator_adaptive.integrate(fx2p,0, t.thickness_cm(),1e-3, 1e-3,4)/t.density()/t.density(); res.sigma_x = sqrt(res.sigma_x); + + auto fx1 = [&](double x)->double{ + double range = range_spline(T); + double tt = range - range_spline(x); + double t0 = std::min(range, t.thickness()); + return (t0-tt)*angular_variance_spline.derivative(x); + }; + + res.cov = integrator_adaptive.integrate(fx1,res.Eout, res.Eout, 1e-6, 1e-3,4)/t.density(); + #ifdef REACTIONS res.sp = nonreaction_rate(p,t,c); #endif @@ -275,18 +288,20 @@ MultiResult calculate(Projectile &p, const Layers &layers, const Config &c){ res.total_result.Ein = e; res.results.reserve(layers.num()); double z = 0; + double cov = 0; for(auto&m:layers.get_materials()){ z += m.thickness_cm(); Result r = calculate(p,m,e,c); - e = r.Eout; - res.total_result.sigma_a += r.sigma_a*r.sigma_a; + e = r.Eout; res.total_result.Eloss += r.Eloss; res.total_result.sigma_E += r.sigma_E*r.sigma_E; res.total_result.tof += r.tof; - res.total_result.Eout = r.Eout; - double temp = r.sigma_a*(layers.thickness_cm() - z); - res.total_result.sigma_x += (r.sigma_x*r.sigma_x) + (temp*temp); - //res.total_result.sigma_x += r.sigma_x*r.sigma_x; + res.total_result.Eout = r.Eout; + double a2 = res.total_result.sigma_a; + res.total_result.sigma_x += (2*m.thickness_cm()*res.total_result.cov) + (a2*m.thickness_cm()*m.thickness_cm()) + r.sigma_x*r.sigma_x; + //res.total_result.sigma_x += (a2*m.thickness_cm()*m.thickness_cm()) + r.sigma_x*r.sigma_x; + res.total_result.cov += a2*m.thickness_cm() + r.cov; + res.total_result.sigma_a += r.sigma_a*r.sigma_a; #ifdef REACTIONS res.total_result.sp = (r.sp>=0.0)?res.total_result.sp*r.sp:-1; #endif diff --git a/integrator.cpp b/integrator.cpp index 20c7463..80386de 100644 --- a/integrator.cpp +++ b/integrator.cpp @@ -7,6 +7,8 @@ namespace catima{ integrator_type integrator; + integrator_adaptive_type integrator_adaptive; + #ifdef GSL_INTEGRATION double funcwrapper3(double x, void *_c){ std::function *f = (std::function *)_c; diff --git a/integrator.h b/integrator.h index 56bd398..a039055 100644 --- a/integrator.h +++ b/integrator.h @@ -23,11 +23,16 @@ //#ifdef USE_THREADS //#include //#endif - +#include "catima/glq_integrator.h" +#include "catima/gkq_integrator.h" #ifdef GSL_INTEGRATION #include "gsl/gsl_integration.h" #endif + +using integrators::GaussLegendreIntegration; +using integrators::GaussKronrodIntegration; + namespace catima{ #ifdef GSL_INTEGRATION @@ -52,136 +57,15 @@ class IntegratorGSL{ }; #endif -template -struct GL_data{ -}; - -template -class GaussLegendreIntegration{ -public: - template - double integrate(F f, double a, double b) const; - template - double operator()(F f, double a, double b) const {return integrate(f, a, b);} - double w(int i) const {return GL_data::w()[i];} - double x(int i) const {return GL_data::x()[i];} - int n() const {return order;} - std::array get_points(double a = -1.0, double b = 1.0)const; -}; - -template -template -double GaussLegendreIntegration::integrate(F f, double a, double b) const{ - double res=0.0; - double p = 0.5*(b-a); - double q = 0.5*(b+a); - - if(order%2){res+= w(0) * (f(p*x(0) + q));} // in case odd-order - for(int i=order%2;i -std::array GaussLegendreIntegration::get_points(double a, double b)const{ - std::array 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; -} - -// order = 8 -template<> -struct GL_data<8>{ - static const std::array& x(){ - static const std::array _x = - {0.1834346424956498049394761,0.5255324099163289858177390,0.7966664774136267395915539,0.9602898564975362316835609}; - return _x; - } - static const std::array& w(){ - static const std::array _w = - {0.3626837833783619829651504,0.3137066458778872873379622,0.2223810344533744705443560,0.1012285362903762591525314}; - return _w; - } -}; - -// order = 10 -template<> -struct GL_data<10>{ - static const std::array& x(){ - static const std::array _x = - {0.1488743389816312108848260,0.4333953941292471907992659,0.6794095682990244062343274,0.8650633666889845107320967,0.9739065285171717200779640}; - return _x; - } - static const std::array& w(){ - static const std::array _w = - {0.2955242247147528701738930,0.2692667193099963550912269,0.2190863625159820439955349,0.1494513491505805931457763,0.0666713443086881375935688}; - return _w; - } -}; - -// order = 12 -template<> -struct GL_data<12>{ - static const std::array& x(){ - static const std::array _x = - {0.1252334085114689154724414,0.3678314989981801937526915,0.5873179542866174472967024,0.7699026741943046870368938,0.9041172563704748566784659,0.9815606342467192506905491}; - return _x; - } - - static const std::array& w(){ - static const std::array _w = - {0.2491470458134027850005624,0.2334925365383548087608499,0.2031674267230659217490645,0.1600783285433462263346525,0.1069393259953184309602547,0.0471753363865118271946160}; - return _w; - } -}; - -// order = 14 -template<> -struct GL_data<14>{ - static const std::array& x(){ - static const std::array _x = - {0.1080549487073436620662447,0.3191123689278897604356718,0.5152486363581540919652907,0.6872929048116854701480198,0.8272013150697649931897947,0.9284348836635735173363911,0.9862838086968123388415973}; - return _x; - } - - static const std::array& w(){ - static const std::array _w = - {0.2152638534631577901958764,0.2051984637212956039659241,0.1855383974779378137417166,0.1572031671581935345696019,0.1215185706879031846894148,0.0801580871597602098056333,0.0351194603317518630318329}; - return _w; - } -}; - -// order = 16 -template<> -struct GL_data<16>{ - static const std::array& x(){ - static const std::array _x = - {0.0950125098376374401853193,0.2816035507792589132304605,0.4580167776572273863424194,0.6178762444026437484466718,0.7554044083550030338951012,0.8656312023878317438804679,0.9445750230732325760779884,0.9894009349916499325961542}; - return _x; - } - - static const std::array& w(){ - static const std::array _w = - {0.1894506104550684962853967,0.1826034150449235888667637,0.1691565193950025381893121,0.1495959888165767320815017,0.1246289712555338720524763,0.0951585116824927848099251,0.0622535239386478928628438,0.0271524594117540948517806}; - return _w; - } -}; - #ifdef GSL_INTEGRATION using integrator_type = IntegratorGSL; #else using integrator_type = GaussLegendreIntegration<8>; #endif +using integrator_adaptive_type = GaussKronrodIntegration<21>; extern integrator_type integrator; +extern integrator_adaptive_type integrator_adaptive; } #endif diff --git a/libs/gkq_integrator.h b/libs/gkq_integrator.h new file mode 100644 index 0000000..51db896 --- /dev/null +++ b/libs/gkq_integrator.h @@ -0,0 +1,440 @@ +/* + * Copyright(C) 2017, Andrej Prochazka + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Affero General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Affero General Public License for more details. + + * You should have received a copy of the GNU Affero General Public License + * along with this program. If not, see . + */ + +#ifndef GKQ_INTEGRATOR_H +#define GKQ_INTEGRATOR_H +#include +#include +#include +#include +#include + +namespace integrators{ + +template +struct GK_data{ +}; + +/** + * @brief adaptive integrator + * following orders are supported: 7,15,21,31,61 + * @tparam order Kronrod order + */ +template +class GaussKronrodIntegration{ +public: + template + static double integrate(F& f, double a, double b, double eps = 1e-3, double reps=1e-6, int N = 1); + + template + static double integrate_intervals(F& f, const std::vector>& intervals, double eps = 1e-3, double reps=1e-6); + + template + static std::pair integrate_nonadaptive(F& f, double a, double b); + + template + static double integrate_adaptive(F& f, double a, double b, double eps = 1e-3, double reps=1e-6, int level=49); + + static double w(int i) {return GK_data::w()[i];} + static double wg(int i) {return GK_data::wg()[i];} + static double x(int i) {return GK_data::x()[i];} + static int n() {return order;} + static std::array get_points(double a = -1.0, double b = 1.0); + +private: + static const size_t ngauss = (order-1)/2; + static const size_t xsize = order/2 + 1; + static const bool oddgauss = ((order-1)/2)&1; +}; + +template +template +std::pair GaussKronrodIntegration::integrate_nonadaptive(F& f, double a, double b){ + double res; + double gres = 0.0; + double val; + double p = 0.5*(b-a); + double q = 0.5*(b+a); + unsigned int gi = 1; //by default assume odd-order gauss + unsigned int ki = 2; + + // 1st kronrod point + val = f(p*x(0) + q); + res = w(0) * val; + if(oddgauss){ // 1st gaus point if odd order gauss + gres = wg(0) * val; + gi = 2; + ki = 1; + } + + for(unsigned int i=gi;i::epsilon()); + + return std::make_pair(p*res, p*err); +} + +template +template +double GaussKronrodIntegration::integrate_adaptive(F& f, double a, double b, double eps, double rprec, int level){ + double result = 0.0; + double err = 0.0; + const double numlimit = 10*std::numeric_limits::epsilon(); + + auto r = integrate_nonadaptive(f, a, b); + result = r.first; + err = r.second; + //printf("level = %d, I=%lf, e=%lf, %lf - %lf\n",level,result,err,a,b); + if( (std::abs(result) < numlimit) || ((b-a) aeps)){ + double mid = 0.5*(a+b); + result = integrate_adaptive(f, a, mid, 0.707*aeps, 0, level-1); + result += integrate_adaptive(f, mid, b, 0.707*aeps, 0, level-1); + } + + return result; + } + +template +template +double GaussKronrodIntegration::integrate(F& f, double a, double b, double eps, double reps, int N){ + double step = (b-a)/N; + double result = 0; + for(int i=0; i +template +double GaussKronrodIntegration::integrate_intervals(F& f, const std::vector>& intervals, double eps, double reps){ + double result = 0; + for(auto& i:intervals){ + result+=integrate_adaptive(f,i.first,i.second, eps/intervals.size(), reps); + } + + return result; +} + +template +std::array GaussKronrodIntegration::get_points(double a, double b){ + std::array 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; +} + + +/// weights and abscissas +// order = 7 +template<> +struct GK_data<7>{ + static std::array const & x(){ + static const std::array _x = + {0.0, 0.4342437493468025580021, 0.774596669241483377036,0.9604912687080202834235}; + return _x; + } + static std::array const & w(){ + static const std::array _w = + {0.4509165386584741423451, 0.4013974147759622229051, 0.2684880898683334407286, 0.1046562260264672651938}; + return _w; + } + static std::array const & wg(){ + static const std::array _wg = + {0.8888888888888888888889, 0.555555555555555555556}; + return _wg; + } +}; + +//order = 15 +template<> +struct GK_data<15>{ + static std::array const & x(){ + static const std::array _x = + {0.0, + 0.207784955007898467600689403773244913, + 0.405845151377397166906606412076961463, + 0.586087235467691130294144838258729598, + 0.741531185599394439863864773280788407, + 0.864864423359769072789712788640926201, + 0.949107912342758524526189684047851262, + 0.991455371120812639206854697526328517}; + return _x; + } + static std::array const & w(){ + static const std::array _w = + { + 0.209482141084727828012999174891714264, + 0.204432940075298892414161999234649085, + 0.190350578064785409913256402421013683, + 0.169004726639267902826583426598550284, + 0.140653259715525918745189590510237920, + 0.104790010322250183839876322541518017, + 0.0630920926299785532907006631892042867, + 0.0229353220105292249637320080589695920}; + return _w; + } + static std::array const & wg(){ + static const std::array _wg = + { 0.417959183673469387755102040816326531, + 0.381830050505118944950369775488975134, + 0.279705391489276667901467771423779582, + 0.129484966168869693270611432679082018}; + return _wg; + } +}; + +//order = 21 +template<> +struct GK_data<21>{ + static std::array const & x(){ + static const std::array _x = { + 0.00000000000000000e+00, + 1.48874338981631211e-01, + 2.94392862701460198e-01, + 4.33395394129247191e-01, + 5.62757134668604683e-01, + 6.79409568299024406e-01, + 7.80817726586416897e-01, + 8.65063366688984511e-01, + 9.30157491355708226e-01, + 9.73906528517171720e-01, + 9.95657163025808081e-01}; + return _x; + } + static std::array const & w(){ + static const std::array _w = + { + 1.49445554002916906e-01, + 1.47739104901338491e-01, + 1.42775938577060081e-01, + 1.34709217311473326e-01, + 1.23491976262065851e-01, + 1.09387158802297642e-01, + 9.31254545836976055e-02, + 7.50396748109199528e-02, + 5.47558965743519960e-02, + 3.25581623079647275e-02, + 1.16946388673718743e-02, + }; + return _w; + } + static std::array const & wg(){ + static const std::array _wg = + { + 2.95524224714752870e-01, + 2.69266719309996355e-01, + 2.19086362515982044e-01, + 1.49451349150580593e-01, + 6.66713443086881376e-02,}; + return _wg; + } +}; + + + +//order = 31 +template<> +struct GK_data<31>{ + static std::array const & x(){ + static const std::array _x = { + 0.00000000000000000e+00, + 1.01142066918717499e-01, + 2.01194093997434522e-01, + 2.99180007153168812e-01, + 3.94151347077563370e-01, + 4.85081863640239681e-01, + 5.70972172608538848e-01, + 6.50996741297416971e-01, + 7.24417731360170047e-01, + 7.90418501442465933e-01, + 8.48206583410427216e-01, + 8.97264532344081901e-01, + 9.37273392400705904e-01, + 9.67739075679139134e-01, + 9.87992518020485428e-01, + 9.98002298693397060e-01 + }; + return _x; + } + static std::array const & w(){ + static const std::array _w = + { + 1.01330007014791549e-01, + 1.00769845523875595e-01, + 9.91735987217919593e-02, + 9.66427269836236785e-02, + 9.31265981708253212e-02, + 8.85644430562117706e-02, + 8.30805028231330210e-02, + 7.68496807577203789e-02, + 6.98541213187282587e-02, + 6.20095678006706403e-02, + 5.34815246909280873e-02, + 4.45897513247648766e-02, + 3.53463607913758462e-02, + 2.54608473267153202e-02, + 1.50079473293161225e-02, + 5.37747987292334899e-03 + }; + return _w; + } + static std::array const & wg(){ + static const std::array _wg = + { + 2.02578241925561273e-01, + 1.98431485327111576e-01, + 1.86161000015562211e-01, + 1.66269205816993934e-01, + 1.39570677926154314e-01, + 1.07159220467171935e-01, + 7.03660474881081247e-02, + 3.07532419961172684e-02 + }; + return _wg; + } +}; + +//order = 61 +template<> +struct GK_data<61>{ + static std::array const & x(){ + static const std::array _x = { + 0.00000000000000000e+00, + 5.14718425553176958e-02, + 1.02806937966737030e-01, + 1.53869913608583547e-01, + 2.04525116682309891e-01, + 2.54636926167889846e-01, + 3.04073202273625077e-01, + 3.52704725530878113e-01, + 4.00401254830394393e-01, + 4.47033769538089177e-01, + 4.92480467861778575e-01, + 5.36624148142019899e-01, + 5.79345235826361692e-01, + 6.20526182989242861e-01, + 6.60061064126626961e-01, + 6.97850494793315797e-01, + 7.33790062453226805e-01, + 7.67777432104826195e-01, + 7.99727835821839083e-01, + 8.29565762382768397e-01, + 8.57205233546061099e-01, + 8.82560535792052682e-01, + 9.05573307699907799e-01, + 9.26200047429274326e-01, + 9.44374444748559979e-01, + 9.60021864968307512e-01, + 9.73116322501126268e-01, + 9.83668123279747210e-01, + 9.91630996870404595e-01, + 9.96893484074649540e-01, + 9.99484410050490638e-01 + }; + return _x; + } + static std::array const & w(){ + static const std::array _w = + { + 5.14947294294515676e-02, + 5.14261285374590259e-02, + 5.12215478492587722e-02, + 5.08817958987496065e-02, + 5.04059214027823468e-02, + 4.97956834270742064e-02, + 4.90554345550297789e-02, + 4.81858617570871291e-02, + 4.71855465692991539e-02, + 4.60592382710069881e-02, + 4.48148001331626632e-02, + 4.34525397013560693e-02, + 4.19698102151642461e-02, + 4.03745389515359591e-02, + 3.86789456247275930e-02, + 3.68823646518212292e-02, + 3.49793380280600241e-02, + 3.29814470574837260e-02, + 3.09072575623877625e-02, + 2.87540487650412928e-02, + 2.65099548823331016e-02, + 2.41911620780806014e-02, + 2.18280358216091923e-02, + 1.94141411939423812e-02, + 1.69208891890532726e-02, + 1.43697295070458048e-02, + 1.18230152534963417e-02, + 9.27327965951776343e-03, + 6.63070391593129217e-03, + 3.89046112709988405e-03, + 1.38901369867700762e-03 + }; + return _w; + } + static std::array const & wg(){ + static const std::array _wg = + { + 1.02852652893558840e-01, + 1.01762389748405505e-01, + 9.95934205867952671e-02, + 9.63687371746442596e-02, + 9.21225222377861287e-02, + 8.68997872010829798e-02, + 8.07558952294202154e-02, + 7.37559747377052063e-02, + 6.59742298821804951e-02, + 5.74931562176190665e-02, + 4.84026728305940529e-02, + 3.87991925696270496e-02, + 2.87847078833233693e-02, + 1.84664683110909591e-02, + 7.96819249616660562e-03 + }; + return _wg; + } +}; + + +} // end of namespace +#endif diff --git a/libs/glq_integrator.h b/libs/glq_integrator.h new file mode 100644 index 0000000..d3aaf3e --- /dev/null +++ b/libs/glq_integrator.h @@ -0,0 +1,356 @@ +/* + * Copyright(C) 2017, Andrej Prochazka + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Affero General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Affero General Public License for more details. + + * You should have received a copy of the GNU Affero General Public License + * along with this program. If not, see . + */ +#ifndef GLQ_INTEGRATOR_H +#define GLQ_INTEGRATOR_H +#include + +namespace integrators{ + +template +struct GL_data{ +}; + +template +class GaussLegendreIntegration{ +public: + template + double integrate(F& f, double a, double b) const; + template + double operator()(F& f, double a, double b) const {return integrate(f, a, b);} + double w(int i) const {return GL_data::w()[i];} + double x(int i) const {return GL_data::x()[i];} + int n() const {return order;} + std::array get_points(double a = -1.0, double b = 1.0)const; +}; + +template +template +double GaussLegendreIntegration::integrate(F& f, double a, double b) const{ + double res=0.0; + double p = 0.5*(b-a); + double q = 0.5*(b+a); + + if(order%2){res+= w(0) * (f(p*x(0) + q));} // in case odd-order + for(int i=order%2;i +std::array GaussLegendreIntegration::get_points(double a, double b)const{ + std::array 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; +} + + +template +class GaussLegendreIntegration2D{ +public: + template + double integrate(F& f, double a, double b, double c, double d) const; + template + double operator()(F& f, double a, double b, double c, double d) const {return integrate(f, a,b,c,d);} +private: + GaussLegendreIntegration integrator; +}; + +template +template +double GaussLegendreIntegration2D::integrate(F& f, double a, double b, double c, double d) const{ + double res=0.0; + double p = (b-a)/2.0; + double q = (b+a)/2.0; + double r = (d-c)/2.0; + double s = (d+c)/2.0; + + double xx, yy, sum; + for(int i=0;i +class GaussLegendreIntegration2DA{ +public: + template + double integrate(F& f, double a, double b, double c, double d) const; + template + double operator()(F& f, double a, double b, double c, double d) const {return integrate(f, a,b,c,d);} +private: + GaussLegendreIntegration integrator; + GaussLegendreIntegration integrator2; +}; + +template +template +double GaussLegendreIntegration2DA::integrate(F& f, double a, double b, double c, double d) const{ + double res=0.0; + double p = (b-a)/2.0; + double q = (b+a)/2.0; + double r = (d-c)/2.0; + double s = (d+c)/2.0; + + double xx, yy, sum; + for(int i=0;i +struct GL_data<3>{ + static const std::array& x(){ + static const std::array _x = {0.0,0.7745966692414833770359}; + return _x; + } + + static const std::array& w(){ + static const std::array _w = {0.8888888888888888888889,0.555555555555555555556}; + return _w; + } +}; + +// order = 4 +template<> +struct GL_data<4>{ + static const std::array& x(){ + static const std::array _x = + {0.3399810435848562648026658,0.8611363115940525752239465}; + return _x; + } + + static const std::array& w(){ + static const std::array _w = + {0.6521451548625461426269361,0.3478548451374538573730639}; + return _w; + } +}; + +// order = 5 +template<> +struct GL_data<5>{ + static const std::array& x(){ + static const std::array _x = + {0.0, 0.5384693101056830910363, 0.9061798459386639927976}; + return _x; + } + + static const std::array& w(){ + static const std::array _w = + {0.568888888888888888888, 0.4786286704993664680413, 0.2369268850561890875143}; + return _w; + } +}; + +// order = 6 +template<> +struct GL_data<6>{ + static const std::array& x(){ + static const std::array _x = + {0.2386191860831969086305017,0.6612093864662645136613996,0.9324695142031520278123016}; + return _x; + } + static const std::array& w(){ + static const std::array _w = + {0.4679139345726910473898703,0.3607615730481386075698335,0.1713244923791703450402961}; + return _w; + } +}; + +// order = 7 +template<> +struct GL_data<7>{ + static const std::array& x(){ + static const std::array _x = + {0.0, 0.4058451513773971669066, 0.7415311855993944398639,0.9491079123427585245262}; + return _x; + } + static const std::array& w(){ + static const std::array _w = + {0.417959183673469387755, 0.38183005050511894495, 0.279705391489276667901, 0.129484966168869693271}; + return _w; + } +}; + +// order = 8 +template<> +struct GL_data<8>{ + static const std::array& x(){ + static const std::array _x = + {0.1834346424956498049394761,0.5255324099163289858177390,0.7966664774136267395915539,0.9602898564975362316835609}; + return _x; + } + static const std::array& w(){ + static const std::array _w = + {0.3626837833783619829651504,0.3137066458778872873379622,0.2223810344533744705443560,0.1012285362903762591525314}; + return _w; + } +}; + +// order = 10 +template<> +struct GL_data<10>{ + static const std::array& x(){ + static const std::array _x = + {0.1488743389816312108848260,0.4333953941292471907992659,0.6794095682990244062343274,0.8650633666889845107320967,0.9739065285171717200779640}; + return _x; + } + static const std::array& w(){ + static const std::array _w = + {0.2955242247147528701738930,0.2692667193099963550912269,0.2190863625159820439955349,0.1494513491505805931457763,0.0666713443086881375935688}; + return _w; + } +}; + +// order = 12 +template<> +struct GL_data<12>{ + static const std::array& x(){ + static const std::array _x = + {0.1252334085114689154724414,0.3678314989981801937526915,0.5873179542866174472967024,0.7699026741943046870368938,0.9041172563704748566784659,0.9815606342467192506905491}; + return _x; + } + + static const std::array& w(){ + static const std::array _w = + {0.2491470458134027850005624,0.2334925365383548087608499,0.2031674267230659217490645,0.1600783285433462263346525,0.1069393259953184309602547,0.0471753363865118271946160}; + return _w; + } +}; + +// order = 14 +template<> +struct GL_data<14>{ + static const std::array& x(){ + static const std::array _x = + {0.1080549487073436620662447,0.3191123689278897604356718,0.5152486363581540919652907,0.6872929048116854701480198,0.8272013150697649931897947,0.9284348836635735173363911,0.9862838086968123388415973}; + return _x; + } + + static const std::array& w(){ + static const std::array _w = + {0.2152638534631577901958764,0.2051984637212956039659241,0.1855383974779378137417166,0.1572031671581935345696019,0.1215185706879031846894148,0.0801580871597602098056333,0.0351194603317518630318329}; + return _w; + } +}; + +// order = 16 +template<> +struct GL_data<16>{ + static const std::array& x(){ + static const std::array _x = + {0.0950125098376374401853193,0.2816035507792589132304605,0.4580167776572273863424194,0.6178762444026437484466718,0.7554044083550030338951012,0.8656312023878317438804679,0.9445750230732325760779884,0.9894009349916499325961542}; + return _x; + } + + static const std::array& w(){ + static const std::array _w = + {0.1894506104550684962853967,0.1826034150449235888667637,0.1691565193950025381893121,0.1495959888165767320815017,0.1246289712555338720524763,0.0951585116824927848099251,0.0622535239386478928628438,0.0271524594117540948517806}; + return _w; + } +}; + +// order = 20 +template<> +struct GL_data<20>{ + static const std::array& x(){ + static const std::array _x = + {0.0765265211334973337546404,0.2277858511416450780804962,0.3737060887154195606725482,0.5108670019508270980043641,0.6360536807265150254528367,0.7463319064601507926143051,0.8391169718222188233945291,0.9122344282513259058677524,0.9639719272779137912676661,0.9931285991850949247861224}; + return _x; + } + + static const std::array& w(){ + static const std::array _w = + {0.1527533871307258506980843,0.1491729864726037467878287,0.1420961093183820513292983,0.1316886384491766268984945,0.1181945319615184173123774,0.1019301198172404350367501,0.0832767415767047487247581,0.0626720483341090635695065,0.0406014298003869413310400,0.0176140071391521183118620}; + return _w; + } +}; + +// order = 32 +template<> +struct GL_data<32>{ + static const std::array& x(){ + static const std::array _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}; + return _x; + } + + static const std::array& w(){ + static const std::array _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}; + return _w; + } +}; + +//order = 64 +template<> +struct GL_data<64>{ + static const std::array& x(){ + static const std::array _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}; + return _x; + } + + static const std::array& w(){ + static const std::array _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}; + return _w; + } +}; + + +//order = 128 +template<> +struct GL_data<128>{ + static const std::array& x(){ + static const std::array _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}; + return _x; + } + + static const std::array& w(){ + static const std::array _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}; + return _w; + } +}; + +} // end of namespace +#endif diff --git a/structures.h b/structures.h index d63615d..8c46d17 100644 --- a/structures.h +++ b/structures.h @@ -239,6 +239,7 @@ namespace catima{ double sigma_a=0.0; double sigma_r=0.0; double sigma_x=0.0; + double cov = 0.0; double tof=0.0; #ifdef REACTIONS double sp = 1.0; diff --git a/tests/test_calculations.cpp b/tests/test_calculations.cpp index cb4ffa4..60cb286 100644 --- a/tests/test_calculations.cpp +++ b/tests/test_calculations.cpp @@ -241,6 +241,7 @@ using namespace std; water.thickness_cm(7.4); ll.add(water); res2 = catima::calculate(p,ll); + CHECK(ll.thickness_cm() == approx(17,0.00001)); CHECK(res2.total_result.sigma_x == approx(0.25,0.05)); CHECK(res2.total_result.sigma_x == approx(res17.sigma_x,0.001)); @@ -251,6 +252,16 @@ using namespace std; l.add(water); res2 = catima::calculate(p,l); CHECK(res2.total_result.sigma_x == approx(res.sigma_x).R(0.01)); + + catima::Layers lll; + water.thickness_cm(29.4/4); + lll.add(water); + lll.add(water); + lll.add(water); + lll.add(water); + res2 = catima::calculate(p,lll); + CHECK(res2.total_result.sigma_x == approx(res29.sigma_x).R(0.01)); + } TEST_CASE("result from stopped material"){ catima::Projectile p{12,6,6,1000}; @@ -461,4 +472,4 @@ using namespace std; CHECK(0.1*hbar*c_light/atomic_mass_unit == approx(0.21183,0.0001)); CHECK(16.0*dedx_constant*electron_mass*fine_structure/(atomic_mass_unit*3.0*4.0*M_PI) == approx(5.21721169334564e-7).R(1e-3)); } - \ No newline at end of file +