Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/lend/src/nf_Legendre.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_from_ptwXY_callback_s {
 14     int l;
 15     double mu1, mu2, f1, f2;
 16 };
 17 
 18 static nfu_status nf_Legendre_to_ptwXY2( double mu, double *P, void *argList );
 19 static nfu_status nf_Legendre_from_ptwXY_callback( double mu, double *f, void *argList );
 20 /*
 21 ************************************************************
 22 */
 23 nf_Legendre *nf_Legendre_new( int initialSize, int maxOrder, double *Cls, nfu_status *status ) {
 24 
 25     int l;
 26     nf_Legendre *Legendre = (nf_Legendre *) nfu_malloc( sizeof( nf_Legendre ) );
 27 
 28     *status = nfu_mallocError;
 29     if( Legendre == NULL ) return( NULL );
 30     if( ( *status = nf_Legendre_setup( Legendre, initialSize, maxOrder ) ) != nfu_Okay ) {
 31         nfu_free( Legendre );
 32         return( NULL );
 33     }
 34     for( l = 0; l <= Legendre->maxOrder; l++ ) Legendre->Cls[l] = Cls[l];
 35     return( Legendre );
 36 }
 37 /*
 38 ************************************************************
 39 */
 40 nfu_status nf_Legendre_setup( nf_Legendre *Legendre, int initialSize, int maxOrder ) {
 41 
 42     memset( Legendre, 0, sizeof( nf_Legendre ) );
 43     if( maxOrder < 0 ) maxOrder = -1;
 44     if( maxOrder > nf_Legendre_maxMaxOrder ) maxOrder = nf_Legendre_maxMaxOrder;
 45     Legendre->maxOrder = maxOrder;
 46     if( initialSize < ( maxOrder + 1 ) ) initialSize = maxOrder + 1;
 47     return( nf_Legendre_reallocateCls( Legendre, initialSize, 0 ) );
 48 }
 49 /*
 50 ************************************************************
 51 */
 52 nfu_status nf_Legendre_release( nf_Legendre *Legendre ) {
 53 
 54     if( Legendre->allocated > 0 ) nfu_free( Legendre->Cls );
 55     memset( Legendre, 0, sizeof( nf_Legendre ) );
 56     return( nfu_Okay );
 57 }
 58 /*
 59 ************************************************************
 60 */
 61 nf_Legendre *nf_Legendre_free( nf_Legendre *Legendre ) {
 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 *nfL, nfu_status *status ) {
 71 
 72     return( nf_Legendre_new( 0, nfL->maxOrder, nfL->Cls, status ) );
 73 }
 74 /*
 75 ************************************************************
 76 */
 77 nfu_status nf_Legendre_reallocateCls( nf_Legendre *Legendre, int size, int forceSmallerResize ) {
 78 
 79     nfu_status status = nfu_Okay;
 80 
 81     if( size < nf_Legendre_minMaxOrder ) size = nf_Legendre_minMaxOrder;
 82     if( size > ( nf_Legendre_maxMaxOrder + 1 ) ) size = nf_Legendre_maxMaxOrder + 1;
 83     if( size != Legendre->allocated ) {
 84         if( size > Legendre->allocated ) {
 85             Legendre->Cls = (double *) nfu_realloc( size * sizeof( double ), Legendre->Cls ); }
 86         else {
 87             if( size < ( Legendre->maxOrder + 1 ) ) size = Legendre->maxOrder + 1;
 88             if( ( Legendre->allocated > 2 * size ) || forceSmallerResize ) {
 89                     Legendre->Cls = (double *) nfu_realloc( size * sizeof( double ), Legendre->Cls ); } 
 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 *Legendre ) {
106 
107     return( Legendre->maxOrder );
108 }
109 /*
110 ************************************************************
111 */
112 int nf_Legendre_allocated( nf_Legendre *Legendre ) {
113 
114     return( Legendre->allocated );
115 }
116 /*
117 ************************************************************
118 */
119 double nf_Legendre_getCl( nf_Legendre *Legendre, int l, nfu_status *status ) {
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 *Legendre, int l, double Cl ) {
132 
133     nfu_status status;
134 
135     if( ( l < 0 ) || ( l > ( Legendre->maxOrder + 1 ) ) ) return( nfu_badIndex );
136     if( Legendre->allocated <= l ) {
137         if( ( status = nf_Legendre_reallocateCls( Legendre, l + nf_Legendre_sizeIncrement, 0 ) ) != nfu_Okay ) return( status );
138     }
139     if( l > Legendre->maxOrder ) Legendre->maxOrder = l;
140     Legendre->Cls[l] = Cl;
141     return( nfu_Okay );
142 }
143 /*
144 ************************************************************
145 */
146 nfu_status nf_Legendre_normalize( nf_Legendre *Legendre ) {
147 
148     int l;
149     double norm;
150 
151     if( Legendre->maxOrder >= 0 ) {
152         if( ( norm = Legendre->Cls[0] ) == 0 ) return( nfu_divByZero );
153         for( l = 0; l <= Legendre->maxOrder; l++ ) Legendre->Cls[l] /= norm;
154     }
155     return( nfu_Okay );
156 }
157 /*
158 ************************************************************
159 */
160 double nf_Legendre_evauluateAtMu( nf_Legendre *Legendre, double mu, nfu_status *status ) {
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++ ) P += ( l + 0.5 ) * Legendre->Cls[l] * nf_Legendre_PofL_atMu( l, mu );
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 - 0.5 ); }
189         else if( l == 3 ) {
190             return(                2.5 * mu2 -          1.5 ) * mu; }
191         else if( l == 4 ) {
192             return(              4.375 * mu2 -         3.75 ) * mu2 +       0.375; }
193         else if( l == 5 ) {
194             return( (            7.875 * mu2 -         8.75 ) * mu2 +       1.875 ) * mu; }
195         else if( l == 6 ) {
196             return( (          14.4375 * mu2 -      19.6875 ) * mu2 +      6.5625 ) * mu2 -      0.3125; }
197         else if( l == 7 ) {
198             return( ( (        26.8125 * mu2 -      43.3125 ) * mu2 +     19.6875 ) * mu2 -      2.1875 ) * mu; }
199         else if( l == 8 ) {
200             return( ( (     50.2734375 * mu2 -     93.84375 ) * mu2 +   54.140625 ) * mu2 -     9.84375 ) * mu2 + 0.2734375; }
201         else {
202             return( ( ( (   94.9609375 * mu2 -    201.09375 ) * mu2 +  140.765625 ) * mu2 -    36.09375 ) * mu2 + 2.4609375 ) * mu;
203         }
204     }
205 */
206 
207     Pl = 0.;
208     Pl_plus1 = 1.;
209     for( l_ = 0, twoL_plus1 = 1; l_ < l; l_++, twoL_plus1 += 2 ) {
210         Pl_minus1 = Pl;
211         Pl = Pl_plus1;
212         Pl_plus1 = ( twoL_plus1 * mu * Pl - l_ * Pl_minus1 ) / ( l_ + 1 );
213     }
214     return( Pl_plus1 );
215 }
216 /*
217 ************************************************************
218 */
219 ptwXYPoints *nf_Legendre_to_ptwXY( nf_Legendre *Legendre, double accuracy, int biSectionMax, int checkForRoots, nfu_status *status ) {
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-1] + dx;
233     }
234     xs[n] = 1.;
235     return( ptwXY_createFromFunction( n + 1, xs, nf_Legendre_to_ptwXY2, (void *) argList, accuracy, checkForRoots, biSectionMax, status ) );
236 }
237 /*
238 ************************************************************
239 */
240 static nfu_status nf_Legendre_to_ptwXY2( double mu, double *P, void *argList ) {
241 
242     nfu_status status;      /* Set by nf_Legendre_evauluateAtMu. */
243 
244     *P = nf_Legendre_evauluateAtMu( (nf_Legendre *) argList, mu, &status );
245     return( status );
246 }
247 /*
248 ************************************************************
249 */
250 nf_Legendre *nf_Legendre_from_ptwXY( ptwXYPoints *ptwXY, int maxOrder, nfu_status *status ) {
251 
252     int l, i, n = (int) ptwXY_length( ptwXY );
253     nf_Legendre *Legendre;
254     double mu1, mu2, f1, f2, Cl, Cls[1] = { 0 }, integral;
255     struct nf_Legendre_from_ptwXY_callback_s argList;
256 
257     if( ( *status = ptwXY_getStatus( ptwXY ) ) != nfu_Okay ) return( NULL );
258 
259     ptwXY_getXYPairAtIndex( ptwXY, 0, &mu1, &f1 );
260     if( mu1 < -1 ) {
261         *status = nfu_XOutsideDomain;
262         return( NULL );
263     }
264 
265     ptwXY_getXYPairAtIndex( ptwXY, 0, &mu2, &f2 );
266     if( mu2 > 1 ) {
267         *status = nfu_XOutsideDomain;
268         return( NULL );
269     }
270 
271     if( ( Legendre = nf_Legendre_new( maxOrder + 1, -1, Cls, status ) ) == NULL ) return( NULL );
272 
273     if( maxOrder > nf_Legendre_maxMaxOrder ) maxOrder = nf_Legendre_maxMaxOrder;
274     for( l = 0; l <= maxOrder; l++ ) {
275         ptwXY_getXYPairAtIndex( ptwXY, 0, &mu1, &f1 );
276         argList.l = l;
277         for( i = 1, Cl = 0; i < n; i++ ) {
278             ptwXY_getXYPairAtIndex( ptwXY, i, &mu2, &f2 );
279             argList.mu1 = mu1;
280             argList.f1 = f1;
281             argList.mu2 = mu2;
282             argList.f2 = f2;
283             if( ( *status = nf_Legendre_GaussianQuadrature( l + 1, mu1, mu2, nf_Legendre_from_ptwXY_callback, (void *) &argList, &integral ) ) != nfu_Okay )
284                 goto err;
285             Cl += integral;
286             mu1 = mu2;
287             f1 = f2;
288         }
289         if( ( *status = nf_Legendre_setCl( Legendre, l, Cl ) ) != nfu_Okay ) goto err;
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_callback( double mu, double *f, void *argList ) {
301 
302     struct nf_Legendre_from_ptwXY_callback_s *args = (struct nf_Legendre_from_ptwXY_callback_s *) argList;
303 
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 );
306     return( nfu_Okay );
307 }
308 
309 #if defined __cplusplus
310 }
311 #endif
312