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