Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/lend/src/nf_Legendre_GaussianQuadrature.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ] 

  1 /*
  2 # <<BEGIN-copyright>>
  3 # <<END-copyright>>
  4 */
  5 
  6 #include "nf_Legendre.h"
  7 
  8 #if defined __cplusplus
  9 namespace GIDI {
 10 using namespace GIDI;
 11 #endif
 12 
 13 struct nf_Legendre_GaussianQuadrature_degree {
 14     int n;
 15     double *weights;
 16     double *xis;
 17 };
 18 
 19 static double sqrt_inv3 = 0.57735026918962576451;        /* sqrt( 1. / 3. ); */
 20 
 21 #define n_3 3
 22 static double weights_3[(n_3 + 1)/2] = { 8. / 9., 5. / 9. };
 23 static double xis_3[(n_3 + 1)/2] = { 0., 0.77459666924148337704 };
 24 
 25 #define n_4 4
 26 static double weights_4[(n_4 + 1)/2] = { 0.65214515486254614263, 0.34785484513745385737 };
 27 static double xis_4[(n_4 + 1)/2] = { 0.33998104358485626480, 0.86113631159405257522 };
 28 
 29 #define n_5 5
 30 static double weights_5[(n_5 + 1)/2] = { 0.568888888888889, 0.478628670499366, 0.236926885056189 };
 31 static double xis_5[(n_5 + 1)/2] = { 0.0, 0.538469310105683, 0.906179845938664 };
 32 
 33 #define n_10 10
 34 static double weights_10[(n_10 + 1)/2] = { 0.295524224714752870, 0.269266719309996355, 0.219086362515982044, 0.149451349150580593, 0.066671344308688138 };
 35 static double xis_10[(n_10 + 1)/2] = { 0.148874338981631211, 0.433395394129247191, 0.679409568299024406, 0.865063366688984511, 0.973906528517171720 };
 36 
 37 #define n_20 20
 38 static double weights_20[(n_20 + 1)/2] = { 
 39     0.152753387130725850698, 0.149172986472603746788, 0.142096109318382051329, 0.131688638449176626898, 0.118194531961518417312,
 40     0.101930119817240435037, 0.083276741576704748725, 0.062672048334109063570, 0.040601429800386941331, 0.017614007139152118312 };
 41 static double xis_20[(n_20 + 1)/2] = { 
 42     0.076526521133497333755, 0.227785851141645078080, 0.373706088715419560673, 0.510867001950827098004, 0.636053680726515025453,
 43     0.746331906460150792614, 0.839116971822218823395, 0.912234428251325905868, 0.963971927277913791268, 0.993128599185094924786 };
 44 
 45 #define n_40 40
 46 static double weights_40[(n_40 + 1)/2] = {
 47     0.077505947978424811264, 0.077039818164247965588, 0.076110361900626242372, 0.074723169057968264200, 0.072886582395804059061, 
 48     0.070611647391286779696, 0.067912045815233903826, 0.064804013456601038075, 0.061306242492928939167, 0.057439769099391551367, 
 49     0.053227846983936824355, 0.048695807635072232061, 0.043870908185673271992, 0.038782167974472017640, 0.033460195282547847393, 
 50     0.027937006980023401099, 0.022245849194166957262, 0.016421058381907888713, 0.010498284531152813615, 0.004521277098533191258 };
 51 static double xis_40[(n_40 + 1)/2] = { 
 52     0.038772417506050821933, 0.116084070675255208483, 0.192697580701371099716, 0.268152185007253681141, 0.341994090825758473007, 
 53     0.413779204371605001525, 0.483075801686178712909, 0.549467125095128202076, 0.612553889667980237953, 0.671956684614179548379, 
 54     0.727318255189927103281, 0.778305651426519387695, 0.824612230833311663196, 0.865959503212259503821, 0.902098806968874296728, 
 55     0.932812808278676533361, 0.957916819213791655805, 0.977259949983774262663, 0.990726238699457006453, 0.998237709710559200350 };
 56 
 57 #define nSets 6
 58 static struct nf_Legendre_GaussianQuadrature_degree GaussianQuadrature_degrees[nSets] = { { n_3, weights_3, xis_3 }, { n_4, weights_4, xis_4 }, 
 59     { n_5, weights_5, xis_5 }, { n_10, weights_10, xis_10 }, { n_20, weights_20, xis_20 }, { n_40, weights_40, xis_40 } };
 60 /*
 61 ************************************************************
 62 */
 63 nfu_status nf_Legendre_GaussianQuadrature( int degree, double x1, double x2, nf_Legendre_GaussianQuadrature_callback func, void *argList, double *integral ) {
 64 
 65     int i, n;
 66     double x, mu, sum, *weights, *xis;
 67     nfu_status status = nfu_Okay;
 68 
 69     *integral = 0;
 70     if( degree < 2 ) {
 71         status = func( 0.5 * ( x1 + x2 ), integral, argList );
 72         *integral *= 2.; }
 73     else if( degree < 4 ) {
 74         x = 0.5 * ( -sqrt_inv3 * ( x2 - x1 ) + x1 + x2 );
 75         if( ( status = func( x, integral, argList ) ) == nfu_Okay ) {
 76             x = 0.5 * (  sqrt_inv3 * ( x2 - x1 ) + x1 + x2 );
 77             status = func( x, &sum, argList );
 78             *integral += sum;
 79         } }
 80     else {
 81         for( i = 0; i < nSets - 1; i++ ) {
 82             if( GaussianQuadrature_degrees[i].n > ( degree + 1 ) / 2 ) break;
 83         }
 84         n = ( GaussianQuadrature_degrees[i].n + 1 ) / 2;
 85         weights = GaussianQuadrature_degrees[i].weights;
 86         xis = GaussianQuadrature_degrees[i].xis;
 87         for( i = 0; i < n; i++ ) {
 88             mu = xis[i];
 89             x = 0.5 * ( x1 * ( 1 - mu ) + x2 * ( mu + 1 ) );
 90             if( ( status = func( x, &sum, argList ) ) != nfu_Okay ) break;
 91             *integral += sum * weights[i];
 92             if( mu == 0 ) continue;
 93             x = x1 + x2 - x;
 94             if( ( status = func( x, &sum, argList ) ) != nfu_Okay ) break;
 95             *integral += sum * weights[i];
 96         }
 97     }
 98     *integral *= 0.5 * ( x2 - x1 );
 99     return( status );
100 }
101 
102 #if defined __cplusplus
103 }
104 #endif
105