| Geant4 Cross Reference |
1 /* 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.577350269189625764
20
21 #define n_3 3
22 static double weights_3[(n_3 + 1)/2] = { 8. /
23 static double xis_3[(n_3 + 1)/2] = { 0., 0.774
24
25 #define n_4 4
26 static double weights_4[(n_4 + 1)/2] = { 0.652
27 static double xis_4[(n_4 + 1)/2] = { 0.3399810
28
29 #define n_5 5
30 static double weights_5[(n_5 + 1)/2] = { 0.568
31 static double xis_5[(n_5 + 1)/2] = { 0.0, 0.53
32
33 #define n_10 10
34 static double weights_10[(n_10 + 1)/2] = { 0.2
35 static double xis_10[(n_10 + 1)/2] = { 0.14887
36
37 #define n_20 20
38 static double weights_20[(n_20 + 1)/2] = {
39 0.152753387130725850698, 0.149172986472603
40 0.101930119817240435037, 0.083276741576704
41 static double xis_20[(n_20 + 1)/2] = {
42 0.076526521133497333755, 0.227785851141645
43 0.746331906460150792614, 0.839116971822218
44
45 #define n_40 40
46 static double weights_40[(n_40 + 1)/2] = {
47 0.077505947978424811264, 0.077039818164247
48 0.070611647391286779696, 0.067912045815233
49 0.053227846983936824355, 0.048695807635072
50 0.027937006980023401099, 0.022245849194166
51 static double xis_40[(n_40 + 1)/2] = {
52 0.038772417506050821933, 0.116084070675255
53 0.413779204371605001525, 0.483075801686178
54 0.727318255189927103281, 0.778305651426519
55 0.932812808278676533361, 0.957916819213791
56
57 #define nSets 6
58 static struct nf_Legendre_GaussianQuadrature_d
59 { n_5, weights_5, xis_5 }, { n_10, weights
60 /*
61 **********************************************
62 */
63 nfu_status nf_Legendre_GaussianQuadrature( int
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 ), inte
72 *integral *= 2.; }
73 else if( degree < 4 ) {
74 x = 0.5 * ( -sqrt_inv3 * ( x2 - x1 ) +
75 if( ( status = func( x, integral, argL
76 x = 0.5 * ( sqrt_inv3 * ( x2 - x1
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].
83 }
84 n = ( GaussianQuadrature_degrees[i].n
85 weights = GaussianQuadrature_degrees[i
86 xis = GaussianQuadrature_degrees[i].xi
87 for( i = 0; i < n; i++ ) {
88 mu = xis[i];
89 x = 0.5 * ( x1 * ( 1 - mu ) + x2 *
90 if( ( status = func( x, &sum, argL
91 *integral += sum * weights[i];
92 if( mu == 0 ) continue;
93 x = x1 + x2 - x;
94 if( ( status = func( x, &sum, argL
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