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