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