| 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_from_ptwXY_callback_s { 13 struct nf_Legendre_from_ptwXY_callback_s {
14 int l; 14 int l;
15 double mu1, mu2, f1, f2; 15 double mu1, mu2, f1, f2;
16 }; 16 };
17 17
18 static nfu_status nf_Legendre_to_ptwXY2( doubl 18 static nfu_status nf_Legendre_to_ptwXY2( double mu, double *P, void *argList );
19 static nfu_status nf_Legendre_from_ptwXY_callb 19 static nfu_status nf_Legendre_from_ptwXY_callback( double mu, double *f, void *argList );
20 /* 20 /*
21 ********************************************** 21 ************************************************************
22 */ 22 */
23 nf_Legendre *nf_Legendre_new( int initialSize, 23 nf_Legendre *nf_Legendre_new( int initialSize, int maxOrder, double *Cls, nfu_status *status ) {
24 24
25 int l; 25 int l;
26 nf_Legendre *Legendre = (nf_Legendre *) nf 26 nf_Legendre *Legendre = (nf_Legendre *) nfu_malloc( sizeof( nf_Legendre ) );
27 27
28 *status = nfu_mallocError; 28 *status = nfu_mallocError;
29 if( Legendre == NULL ) return( NULL ); 29 if( Legendre == NULL ) return( NULL );
30 if( ( *status = nf_Legendre_setup( Legendr 30 if( ( *status = nf_Legendre_setup( Legendre, initialSize, maxOrder ) ) != nfu_Okay ) {
31 nfu_free( Legendre ); 31 nfu_free( Legendre );
32 return( NULL ); 32 return( NULL );
33 } 33 }
34 for( l = 0; l <= Legendre->maxOrder; l++ ) 34 for( l = 0; l <= Legendre->maxOrder; l++ ) Legendre->Cls[l] = Cls[l];
35 return( Legendre ); 35 return( Legendre );
36 } 36 }
37 /* 37 /*
38 ********************************************** 38 ************************************************************
39 */ 39 */
40 nfu_status nf_Legendre_setup( nf_Legendre *Leg 40 nfu_status nf_Legendre_setup( nf_Legendre *Legendre, int initialSize, int maxOrder ) {
41 41
42 memset( Legendre, 0, sizeof( nf_Legendre ) 42 memset( Legendre, 0, sizeof( nf_Legendre ) );
43 if( maxOrder < 0 ) maxOrder = -1; 43 if( maxOrder < 0 ) maxOrder = -1;
44 if( maxOrder > nf_Legendre_maxMaxOrder ) m 44 if( maxOrder > nf_Legendre_maxMaxOrder ) maxOrder = nf_Legendre_maxMaxOrder;
45 Legendre->maxOrder = maxOrder; 45 Legendre->maxOrder = maxOrder;
46 if( initialSize < ( maxOrder + 1 ) ) initi 46 if( initialSize < ( maxOrder + 1 ) ) initialSize = maxOrder + 1;
47 return( nf_Legendre_reallocateCls( Legendr 47 return( nf_Legendre_reallocateCls( Legendre, initialSize, 0 ) );
48 } 48 }
49 /* 49 /*
50 ********************************************** 50 ************************************************************
51 */ 51 */
52 nfu_status nf_Legendre_release( nf_Legendre *L 52 nfu_status nf_Legendre_release( nf_Legendre *Legendre ) {
53 53
54 if( Legendre->allocated > 0 ) nfu_free( Le 54 if( Legendre->allocated > 0 ) nfu_free( Legendre->Cls );
55 memset( Legendre, 0, sizeof( nf_Legendre ) 55 memset( Legendre, 0, sizeof( nf_Legendre ) );
56 return( nfu_Okay ); 56 return( nfu_Okay );
57 } 57 }
58 /* 58 /*
59 ********************************************** 59 ************************************************************
60 */ 60 */
61 nf_Legendre *nf_Legendre_free( nf_Legendre *Le 61 nf_Legendre *nf_Legendre_free( nf_Legendre *Legendre ) {
62 62
63 nf_Legendre_release( Legendre ); 63 nf_Legendre_release( Legendre );
64 nfu_free( Legendre ); 64 nfu_free( Legendre );
65 return( NULL ); 65 return( NULL );
66 } 66 }
67 /* 67 /*
68 ********************************************** 68 ************************************************************
69 */ 69 */
70 nf_Legendre *nf_Legendre_clone( nf_Legendre *n 70 nf_Legendre *nf_Legendre_clone( nf_Legendre *nfL, nfu_status *status ) {
71 71
72 return( nf_Legendre_new( 0, nfL->maxOrder, 72 return( nf_Legendre_new( 0, nfL->maxOrder, nfL->Cls, status ) );
73 } 73 }
74 /* 74 /*
75 ********************************************** 75 ************************************************************
76 */ 76 */
77 nfu_status nf_Legendre_reallocateCls( nf_Legen 77 nfu_status nf_Legendre_reallocateCls( nf_Legendre *Legendre, int size, int forceSmallerResize ) {
78 78
79 nfu_status status = nfu_Okay; 79 nfu_status status = nfu_Okay;
80 80
81 if( size < nf_Legendre_minMaxOrder ) size 81 if( size < nf_Legendre_minMaxOrder ) size = nf_Legendre_minMaxOrder;
82 if( size > ( nf_Legendre_maxMaxOrder + 1 ) 82 if( size > ( nf_Legendre_maxMaxOrder + 1 ) ) size = nf_Legendre_maxMaxOrder + 1;
83 if( size != Legendre->allocated ) { 83 if( size != Legendre->allocated ) {
84 if( size > Legendre->allocated ) { 84 if( size > Legendre->allocated ) {
85 Legendre->Cls = (double *) nfu_rea 85 Legendre->Cls = (double *) nfu_realloc( size * sizeof( double ), Legendre->Cls ); }
86 else { 86 else {
87 if( size < ( Legendre->maxOrder + 87 if( size < ( Legendre->maxOrder + 1 ) ) size = Legendre->maxOrder + 1;
88 if( ( Legendre->allocated > 2 * si 88 if( ( Legendre->allocated > 2 * size ) || forceSmallerResize ) {
89 Legendre->Cls = (double *) 89 Legendre->Cls = (double *) nfu_realloc( size * sizeof( double ), Legendre->Cls ); }
90 else { 90 else {
91 size = Legendre->allocated; 91 size = Legendre->allocated;
92 } 92 }
93 } 93 }
94 if( Legendre->Cls == NULL ) { 94 if( Legendre->Cls == NULL ) {
95 size = 0; 95 size = 0;
96 status = nfu_mallocError; 96 status = nfu_mallocError;
97 } 97 }
98 Legendre->allocated = size; 98 Legendre->allocated = size;
99 } 99 }
100 return( status ); 100 return( status );
101 } 101 }
102 /* 102 /*
103 ********************************************** 103 ************************************************************
104 */ 104 */
105 int nf_Legendre_maxOrder( nf_Legendre *Legendr 105 int nf_Legendre_maxOrder( nf_Legendre *Legendre ) {
106 106
107 return( Legendre->maxOrder ); 107 return( Legendre->maxOrder );
108 } 108 }
109 /* 109 /*
110 ********************************************** 110 ************************************************************
111 */ 111 */
112 int nf_Legendre_allocated( nf_Legendre *Legend 112 int nf_Legendre_allocated( nf_Legendre *Legendre ) {
113 113
114 return( Legendre->allocated ); 114 return( Legendre->allocated );
115 } 115 }
116 /* 116 /*
117 ********************************************** 117 ************************************************************
118 */ 118 */
119 double nf_Legendre_getCl( nf_Legendre *Legendr 119 double nf_Legendre_getCl( nf_Legendre *Legendre, int l, nfu_status *status ) {
120 120
121 *status = nfu_Okay; 121 *status = nfu_Okay;
122 if( ( l < 0 ) || ( l > Legendre->maxOrder 122 if( ( l < 0 ) || ( l > Legendre->maxOrder ) ) {
123 *status = nfu_badIndex; 123 *status = nfu_badIndex;
124 return( 0. ); 124 return( 0. );
125 } 125 }
126 return( Legendre->Cls[l] ); 126 return( Legendre->Cls[l] );
127 } 127 }
128 /* 128 /*
129 ********************************************** 129 ************************************************************
130 */ 130 */
131 nfu_status nf_Legendre_setCl( nf_Legendre *Leg 131 nfu_status nf_Legendre_setCl( nf_Legendre *Legendre, int l, double Cl ) {
132 132
133 nfu_status status; 133 nfu_status status;
134 134
135 if( ( l < 0 ) || ( l > ( Legendre->maxOrde 135 if( ( l < 0 ) || ( l > ( Legendre->maxOrder + 1 ) ) ) return( nfu_badIndex );
136 if( Legendre->allocated <= l ) { 136 if( Legendre->allocated <= l ) {
137 if( ( status = nf_Legendre_reallocateC 137 if( ( status = nf_Legendre_reallocateCls( Legendre, l + nf_Legendre_sizeIncrement, 0 ) ) != nfu_Okay ) return( status );
138 } 138 }
139 if( l > Legendre->maxOrder ) Legendre->max 139 if( l > Legendre->maxOrder ) Legendre->maxOrder = l;
140 Legendre->Cls[l] = Cl; 140 Legendre->Cls[l] = Cl;
141 return( nfu_Okay ); 141 return( nfu_Okay );
142 } 142 }
143 /* 143 /*
144 ********************************************** 144 ************************************************************
145 */ 145 */
146 nfu_status nf_Legendre_normalize( nf_Legendre 146 nfu_status nf_Legendre_normalize( nf_Legendre *Legendre ) {
147 147
148 int l; 148 int l;
149 double norm; 149 double norm;
150 150
151 if( Legendre->maxOrder >= 0 ) { 151 if( Legendre->maxOrder >= 0 ) {
152 if( ( norm = Legendre->Cls[0] ) == 0 ) 152 if( ( norm = Legendre->Cls[0] ) == 0 ) return( nfu_divByZero );
153 for( l = 0; l <= Legendre->maxOrder; l 153 for( l = 0; l <= Legendre->maxOrder; l++ ) Legendre->Cls[l] /= norm;
154 } 154 }
155 return( nfu_Okay ); 155 return( nfu_Okay );
156 } 156 }
157 /* 157 /*
158 ********************************************** 158 ************************************************************
159 */ 159 */
160 double nf_Legendre_evauluateAtMu( nf_Legendre 160 double nf_Legendre_evauluateAtMu( nf_Legendre *Legendre, double mu, nfu_status *status ) {
161 161
162 int l; 162 int l;
163 double P = 0.; 163 double P = 0.;
164 164
165 *status = nfu_XOutsideDomain; 165 *status = nfu_XOutsideDomain;
166 if( ( mu >= -1. ) && ( mu <= 1. ) ) { 166 if( ( mu >= -1. ) && ( mu <= 1. ) ) {
167 *status = nfu_Okay; 167 *status = nfu_Okay;
168 for( l = 0; l <= Legendre->maxOrder; l 168 for( l = 0; l <= Legendre->maxOrder; l++ ) P += ( l + 0.5 ) * Legendre->Cls[l] * nf_Legendre_PofL_atMu( l, mu );
169 } 169 }
170 return( P ); 170 return( P );
171 } 171 }
172 /* 172 /*
173 ********************************************** 173 ************************************************************
174 */ 174 */
175 double nf_Legendre_PofL_atMu( int l, double mu 175 double nf_Legendre_PofL_atMu( int l, double mu ) {
176 176
177 int l_, twoL_plus1; 177 int l_, twoL_plus1;
178 double Pl_minus1, Pl, Pl_plus1; 178 double Pl_minus1, Pl, Pl_plus1;
179 179
180 if( l == 0 ) { 180 if( l == 0 ) {
181 return( 1. ); } 181 return( 1. ); }
182 else if( l == 1 ) { 182 else if( l == 1 ) {
183 return( mu ); } 183 return( mu ); }
184 /* 184 /*
185 else if( l <= 9 ) { 185 else if( l <= 9 ) {
186 double mu2 = mu * mu; 186 double mu2 = mu * mu;
187 if ( l == 2 ) { 187 if ( l == 2 ) {
188 return( 1.5 * mu2 - 188 return( 1.5 * mu2 - 0.5 ); }
189 else if( l == 3 ) { 189 else if( l == 3 ) {
190 return( 2.5 * mu2 - 190 return( 2.5 * mu2 - 1.5 ) * mu; }
191 else if( l == 4 ) { 191 else if( l == 4 ) {
192 return( 4.375 * mu2 - 192 return( 4.375 * mu2 - 3.75 ) * mu2 + 0.375; }
193 else if( l == 5 ) { 193 else if( l == 5 ) {
194 return( ( 7.875 * mu2 - 194 return( ( 7.875 * mu2 - 8.75 ) * mu2 + 1.875 ) * mu; }
195 else if( l == 6 ) { 195 else if( l == 6 ) {
196 return( ( 14.4375 * mu2 - 196 return( ( 14.4375 * mu2 - 19.6875 ) * mu2 + 6.5625 ) * mu2 - 0.3125; }
197 else if( l == 7 ) { 197 else if( l == 7 ) {
198 return( ( ( 26.8125 * mu2 - 198 return( ( ( 26.8125 * mu2 - 43.3125 ) * mu2 + 19.6875 ) * mu2 - 2.1875 ) * mu; }
199 else if( l == 8 ) { 199 else if( l == 8 ) {
200 return( ( ( 50.2734375 * mu2 - 200 return( ( ( 50.2734375 * mu2 - 93.84375 ) * mu2 + 54.140625 ) * mu2 - 9.84375 ) * mu2 + 0.2734375; }
201 else { 201 else {
202 return( ( ( ( 94.9609375 * mu2 - 202 return( ( ( ( 94.9609375 * mu2 - 201.09375 ) * mu2 + 140.765625 ) * mu2 - 36.09375 ) * mu2 + 2.4609375 ) * mu;
203 } 203 }
204 } 204 }
205 */ 205 */
206 206
207 Pl = 0.; 207 Pl = 0.;
208 Pl_plus1 = 1.; 208 Pl_plus1 = 1.;
209 for( l_ = 0, twoL_plus1 = 1; l_ < l; l_++, 209 for( l_ = 0, twoL_plus1 = 1; l_ < l; l_++, twoL_plus1 += 2 ) {
210 Pl_minus1 = Pl; 210 Pl_minus1 = Pl;
211 Pl = Pl_plus1; 211 Pl = Pl_plus1;
212 Pl_plus1 = ( twoL_plus1 * mu * Pl - l_ 212 Pl_plus1 = ( twoL_plus1 * mu * Pl - l_ * Pl_minus1 ) / ( l_ + 1 );
213 } 213 }
214 return( Pl_plus1 ); 214 return( Pl_plus1 );
215 } 215 }
216 /* 216 /*
217 ********************************************** 217 ************************************************************
218 */ 218 */
219 ptwXYPoints *nf_Legendre_to_ptwXY( nf_Legendre 219 ptwXYPoints *nf_Legendre_to_ptwXY( nf_Legendre *Legendre, double accuracy, int biSectionMax, int checkForRoots, nfu_status *status ) {
220 220
221 int i, n = 1; 221 int i, n = 1;
222 double dx, xs[1000]; 222 double dx, xs[1000];
223 void *argList = (void *) Legendre; 223 void *argList = (void *) Legendre;
224 224
225 *status = nfu_Okay; 225 *status = nfu_Okay;
226 xs[0] = -1; 226 xs[0] = -1;
227 if( Legendre->maxOrder > 1 ) { 227 if( Legendre->maxOrder > 1 ) {
228 n = Legendre->maxOrder - 1; 228 n = Legendre->maxOrder - 1;
229 if( n > 249 ) n = 249; 229 if( n > 249 ) n = 249;
230 n = 4 * n + 1; 230 n = 4 * n + 1;
231 dx = 2. / n; 231 dx = 2. / n;
232 for( i = 1; i < n; i++ ) xs[i] = xs[i- 232 for( i = 1; i < n; i++ ) xs[i] = xs[i-1] + dx;
233 } 233 }
234 xs[n] = 1.; 234 xs[n] = 1.;
235 return( ptwXY_createFromFunction( n + 1, x 235 return( ptwXY_createFromFunction( n + 1, xs, nf_Legendre_to_ptwXY2, (void *) argList, accuracy, checkForRoots, biSectionMax, status ) );
236 } 236 }
237 /* 237 /*
238 ********************************************** 238 ************************************************************
239 */ 239 */
240 static nfu_status nf_Legendre_to_ptwXY2( doubl 240 static nfu_status nf_Legendre_to_ptwXY2( double mu, double *P, void *argList ) {
241 241
242 nfu_status status; /* Set by nf_Legen 242 nfu_status status; /* Set by nf_Legendre_evauluateAtMu. */
243 243
244 *P = nf_Legendre_evauluateAtMu( (nf_Legend 244 *P = nf_Legendre_evauluateAtMu( (nf_Legendre *) argList, mu, &status );
245 return( status ); 245 return( status );
246 } 246 }
247 /* 247 /*
248 ********************************************** 248 ************************************************************
249 */ 249 */
250 nf_Legendre *nf_Legendre_from_ptwXY( ptwXYPoin 250 nf_Legendre *nf_Legendre_from_ptwXY( ptwXYPoints *ptwXY, int maxOrder, nfu_status *status ) {
251 251
252 int l, i, n = (int) ptwXY_length( ptwXY ); 252 int l, i, n = (int) ptwXY_length( ptwXY );
253 nf_Legendre *Legendre; 253 nf_Legendre *Legendre;
254 double mu1, mu2, f1, f2, Cl, Cls[1] = { 0 254 double mu1, mu2, f1, f2, Cl, Cls[1] = { 0 }, integral;
255 struct nf_Legendre_from_ptwXY_callback_s a 255 struct nf_Legendre_from_ptwXY_callback_s argList;
256 256
257 if( ( *status = ptwXY_getStatus( ptwXY ) ) 257 if( ( *status = ptwXY_getStatus( ptwXY ) ) != nfu_Okay ) return( NULL );
258 258
259 ptwXY_getXYPairAtIndex( ptwXY, 0, &mu1, &f 259 ptwXY_getXYPairAtIndex( ptwXY, 0, &mu1, &f1 );
260 if( mu1 < -1 ) { 260 if( mu1 < -1 ) {
261 *status = nfu_XOutsideDomain; 261 *status = nfu_XOutsideDomain;
262 return( NULL ); 262 return( NULL );
263 } 263 }
264 264
265 ptwXY_getXYPairAtIndex( ptwXY, 0, &mu2, &f 265 ptwXY_getXYPairAtIndex( ptwXY, 0, &mu2, &f2 );
266 if( mu2 > 1 ) { 266 if( mu2 > 1 ) {
267 *status = nfu_XOutsideDomain; 267 *status = nfu_XOutsideDomain;
268 return( NULL ); 268 return( NULL );
269 } 269 }
270 270
271 if( ( Legendre = nf_Legendre_new( maxOrder 271 if( ( Legendre = nf_Legendre_new( maxOrder + 1, -1, Cls, status ) ) == NULL ) return( NULL );
272 272
273 if( maxOrder > nf_Legendre_maxMaxOrder ) m 273 if( maxOrder > nf_Legendre_maxMaxOrder ) maxOrder = nf_Legendre_maxMaxOrder;
274 for( l = 0; l <= maxOrder; l++ ) { 274 for( l = 0; l <= maxOrder; l++ ) {
275 ptwXY_getXYPairAtIndex( ptwXY, 0, &mu1 275 ptwXY_getXYPairAtIndex( ptwXY, 0, &mu1, &f1 );
276 argList.l = l; 276 argList.l = l;
277 for( i = 1, Cl = 0; i < n; i++ ) { 277 for( i = 1, Cl = 0; i < n; i++ ) {
278 ptwXY_getXYPairAtIndex( ptwXY, i, 278 ptwXY_getXYPairAtIndex( ptwXY, i, &mu2, &f2 );
279 argList.mu1 = mu1; 279 argList.mu1 = mu1;
280 argList.f1 = f1; 280 argList.f1 = f1;
281 argList.mu2 = mu2; 281 argList.mu2 = mu2;
282 argList.f2 = f2; 282 argList.f2 = f2;
283 if( ( *status = nf_Legendre_Gaussi 283 if( ( *status = nf_Legendre_GaussianQuadrature( l + 1, mu1, mu2, nf_Legendre_from_ptwXY_callback, (void *) &argList, &integral ) ) != nfu_Okay )
284 goto err; 284 goto err;
285 Cl += integral; 285 Cl += integral;
286 mu1 = mu2; 286 mu1 = mu2;
287 f1 = f2; 287 f1 = f2;
288 } 288 }
289 if( ( *status = nf_Legendre_setCl( Leg 289 if( ( *status = nf_Legendre_setCl( Legendre, l, Cl ) ) != nfu_Okay ) goto err;
290 } 290 }
291 return( Legendre ); 291 return( Legendre );
292 292
293 err: 293 err:
294 nf_Legendre_free( Legendre ); 294 nf_Legendre_free( Legendre );
295 return( NULL ); 295 return( NULL );
296 } 296 }
297 /* 297 /*
298 ********************************************** 298 ************************************************************
299 */ 299 */
300 static nfu_status nf_Legendre_from_ptwXY_callb 300 static nfu_status nf_Legendre_from_ptwXY_callback( double mu, double *f, void *argList ) {
301 301
302 struct nf_Legendre_from_ptwXY_callback_s * 302 struct nf_Legendre_from_ptwXY_callback_s *args = (struct nf_Legendre_from_ptwXY_callback_s *) argList;
303 303
304 *f = ( args->f1 * ( args->mu2 - mu ) + arg 304 *f = ( args->f1 * ( args->mu2 - mu ) + args->f2 * ( mu - args->mu1 ) ) / ( args->mu2 - args->mu1 );
305 *f *= nf_Legendre_PofL_atMu( args->l, mu ) 305 *f *= nf_Legendre_PofL_atMu( args->l, mu );
306 return( nfu_Okay ); 306 return( nfu_Okay );
307 } 307 }
308 308
309 #if defined __cplusplus 309 #if defined __cplusplus
310 } 310 }
311 #endif 311 #endif
312 312