Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
nf_Legendre.h File Reference
#include <nf_utilities.h>
#include <ptwXY.h>

Go to the source code of this file.

Classes

struct  nf_Legendre_s
 

Macros

#define nf_Legendre_minMaxOrder   4
 
#define nf_Legendre_maxMaxOrder   64
 
#define nf_Legendre_sizeIncrement   8
 

Typedefs

typedef struct nf_Legendre_s nf_Legendre
 
typedef nfu_status(* nf_Legendre_GaussianQuadrature_callback) (double x, double *y, void *argList)
 

Functions

nf_Legendrenf_Legendre_new (int initialSize, int maxOrder, double *Cls, nfu_status *status)
 
nfu_status nf_Legendre_setup (nf_Legendre *nfL, int initialSize, int maxOrder)
 
nfu_status nf_Legendre_release (nf_Legendre *nfL)
 
nf_Legendrenf_Legendre_free (nf_Legendre *nfL)
 
nf_Legendrenf_Legendre_clone (nf_Legendre *nfL, nfu_status *status)
 
nfu_status nf_Legendre_reallocateCls (nf_Legendre *Legendre, int size, int forceSmallerResize)
 
int nf_Legendre_maxOrder (nf_Legendre *Legendre)
 
int nf_Legendre_allocated (nf_Legendre *Legendre)
 
double nf_Legendre_getCl (nf_Legendre *Legendre, int l, nfu_status *status)
 
nfu_status nf_Legendre_setCl (nf_Legendre *Legendre, int l, double Cl)
 
nfu_status nf_Legendre_normalize (nf_Legendre *Legendre)
 
double nf_Legendre_evauluateAtMu (nf_Legendre *nfL, double mu, nfu_status *status)
 
double nf_Legendre_PofL_atMu (int l, double mu)
 
ptwXYPointsnf_Legendre_to_ptwXY (nf_Legendre *nfL, double accuracy, int biSectionMax, int checkForRoots, nfu_status *status)
 
nf_Legendrenf_Legendre_from_ptwXY (ptwXYPoints *ptwXY, int maxOrder, nfu_status *status)
 
nfu_status nf_Legendre_GaussianQuadrature (int degree, double x1, double x2, nf_Legendre_GaussianQuadrature_callback func, void *argList, double *integral)
 

Macro Definition Documentation

◆ nf_Legendre_maxMaxOrder

#define nf_Legendre_maxMaxOrder   64

Definition at line 18 of file nf_Legendre.h.

◆ nf_Legendre_minMaxOrder

#define nf_Legendre_minMaxOrder   4

Definition at line 17 of file nf_Legendre.h.

◆ nf_Legendre_sizeIncrement

#define nf_Legendre_sizeIncrement   8

Definition at line 19 of file nf_Legendre.h.

Typedef Documentation

◆ nf_Legendre

typedef struct nf_Legendre_s nf_Legendre

Definition at line 21 of file nf_Legendre.h.

◆ nf_Legendre_GaussianQuadrature_callback

typedef nfu_status(* nf_Legendre_GaussianQuadrature_callback) (double x, double *y, void *argList)

Definition at line 29 of file nf_Legendre.h.

Function Documentation

◆ nf_Legendre_allocated()

int nf_Legendre_allocated ( nf_Legendre Legendre)

Definition at line 112 of file nf_Legendre.cc.

112 {
113
114 return( Legendre->allocated );
115}

◆ nf_Legendre_clone()

nf_Legendre * nf_Legendre_clone ( nf_Legendre nfL,
nfu_status status 
)

Definition at line 70 of file nf_Legendre.cc.

70 {
71
72 return( nf_Legendre_new( 0, nfL->maxOrder, nfL->Cls, status ) );
73}
nf_Legendre * nf_Legendre_new(int initialSize, int maxOrder, double *Cls, nfu_status *status)
Definition: nf_Legendre.cc:23
double * Cls
Definition: nf_Legendre.h:26

◆ nf_Legendre_evauluateAtMu()

double nf_Legendre_evauluateAtMu ( nf_Legendre nfL,
double  mu,
nfu_status status 
)

Definition at line 160 of file nf_Legendre.cc.

160 {
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}
double nf_Legendre_PofL_atMu(int l, double mu)
Definition: nf_Legendre.cc:175
@ nfu_Okay
Definition: nf_utilities.h:25
@ nfu_XOutsideDomain
Definition: nf_utilities.h:26

◆ nf_Legendre_free()

nf_Legendre * nf_Legendre_free ( nf_Legendre nfL)

Definition at line 61 of file nf_Legendre.cc.

61 {
62
63 nf_Legendre_release( Legendre );
64 nfu_free( Legendre );
65 return( NULL );
66}
nfu_status nf_Legendre_release(nf_Legendre *Legendre)
Definition: nf_Legendre.cc:52
void * nfu_free(void *p)

Referenced by nf_Legendre_from_ptwXY().

◆ nf_Legendre_from_ptwXY()

nf_Legendre * nf_Legendre_from_ptwXY ( ptwXYPoints ptwXY,
int  maxOrder,
nfu_status status 
)

Definition at line 250 of file nf_Legendre.cc.

250 {
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;
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
293err:
294 nf_Legendre_free( Legendre );
295 return( NULL );
296}
nf_Legendre * nf_Legendre_free(nf_Legendre *Legendre)
Definition: nf_Legendre.cc:61
nfu_status nf_Legendre_setCl(nf_Legendre *Legendre, int l, double Cl)
Definition: nf_Legendre.cc:131
#define nf_Legendre_maxMaxOrder
Definition: nf_Legendre.h:18
nfu_status nf_Legendre_GaussianQuadrature(int degree, double x1, double x2, nf_Legendre_GaussianQuadrature_callback func, void *argList, double *integral)
int64_t ptwXY_length(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:583
nfu_status ptwXY_getStatus(ptwXYPoints *ptwXY)
Definition: ptwXY_core.cc:351
nfu_status ptwXY_getXYPairAtIndex(ptwXYPoints *ptwXY, int64_t index, double *x, double *y)
Definition: ptwXY_core.cc:698

◆ nf_Legendre_GaussianQuadrature()

nfu_status nf_Legendre_GaussianQuadrature ( int  degree,
double  x1,
double  x2,
nf_Legendre_GaussianQuadrature_callback  func,
void *  argList,
double *  integral 
)

Definition at line 63 of file nf_Legendre_GaussianQuadrature.cc.

63 {
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 ), integral, argList );
72 *integral *= 2.; }
73 else if( degree < 4 ) {
74 x = 0.5 * ( -sqrt_inv3 * ( x2 - x1 ) + x1 + x2 );
75 if( ( status = func( x, integral, argList ) ) == nfu_Okay ) {
76 x = 0.5 * ( sqrt_inv3 * ( x2 - x1 ) + x1 + x2 );
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].n > ( degree + 1 ) / 2 ) break;
83 }
84 n = ( GaussianQuadrature_degrees[i].n + 1 ) / 2;
85 weights = GaussianQuadrature_degrees[i].weights;
86 xis = GaussianQuadrature_degrees[i].xis;
87 for( i = 0; i < n; i++ ) {
88 mu = xis[i];
89 x = 0.5 * ( x1 * ( 1 - mu ) + x2 * ( mu + 1 ) );
90 if( ( status = func( x, &sum, argList ) ) != nfu_Okay ) break;
91 *integral += sum * weights[i];
92 if( mu == 0 ) continue;
93 x = x1 + x2 - x;
94 if( ( status = func( x, &sum, argList ) ) != nfu_Okay ) break;
95 *integral += sum * weights[i];
96 }
97 }
98 *integral *= 0.5 * ( x2 - x1 );
99 return( status );
100}
enum nfu_status_e nfu_status

Referenced by nf_Legendre_from_ptwXY().

◆ nf_Legendre_getCl()

double nf_Legendre_getCl ( nf_Legendre Legendre,
int  l,
nfu_status status 
)

Definition at line 119 of file nf_Legendre.cc.

119 {
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}
@ nfu_badIndex
Definition: nf_utilities.h:26

◆ nf_Legendre_maxOrder()

int nf_Legendre_maxOrder ( nf_Legendre Legendre)

Definition at line 105 of file nf_Legendre.cc.

105 {
106
107 return( Legendre->maxOrder );
108}

◆ nf_Legendre_new()

nf_Legendre * nf_Legendre_new ( int  initialSize,
int  maxOrder,
double *  Cls,
nfu_status status 
)

Definition at line 23 of file nf_Legendre.cc.

23 {
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}
nfu_status nf_Legendre_setup(nf_Legendre *Legendre, int initialSize, int maxOrder)
Definition: nf_Legendre.cc:40
@ nfu_mallocError
Definition: nf_utilities.h:25
void * nfu_malloc(size_t size)

Referenced by nf_Legendre_clone(), and nf_Legendre_from_ptwXY().

◆ nf_Legendre_normalize()

nfu_status nf_Legendre_normalize ( nf_Legendre Legendre)

Definition at line 146 of file nf_Legendre.cc.

146 {
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}
@ nfu_divByZero
Definition: nf_utilities.h:27

◆ nf_Legendre_PofL_atMu()

double nf_Legendre_PofL_atMu ( int  l,
double  mu 
)

Definition at line 175 of file nf_Legendre.cc.

175 {
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}

Referenced by nf_Legendre_evauluateAtMu().

◆ nf_Legendre_reallocateCls()

nfu_status nf_Legendre_reallocateCls ( nf_Legendre Legendre,
int  size,
int  forceSmallerResize 
)

Definition at line 77 of file nf_Legendre.cc.

77 {
78
79 nfu_status status = nfu_Okay;
80
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}
#define nf_Legendre_minMaxOrder
Definition: nf_Legendre.h:17
void * nfu_realloc(size_t size, void *old)

Referenced by nf_Legendre_setCl(), and nf_Legendre_setup().

◆ nf_Legendre_release()

nfu_status nf_Legendre_release ( nf_Legendre nfL)

Definition at line 52 of file nf_Legendre.cc.

52 {
53
54 if( Legendre->allocated > 0 ) nfu_free( Legendre->Cls );
55 memset( Legendre, 0, sizeof( nf_Legendre ) );
56 return( nfu_Okay );
57}

Referenced by nf_Legendre_free().

◆ nf_Legendre_setCl()

nfu_status nf_Legendre_setCl ( nf_Legendre Legendre,
int  l,
double  Cl 
)

Definition at line 131 of file nf_Legendre.cc.

131 {
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}
nfu_status nf_Legendre_reallocateCls(nf_Legendre *Legendre, int size, int forceSmallerResize)
Definition: nf_Legendre.cc:77
#define nf_Legendre_sizeIncrement
Definition: nf_Legendre.h:19

Referenced by nf_Legendre_from_ptwXY().

◆ nf_Legendre_setup()

nfu_status nf_Legendre_setup ( nf_Legendre nfL,
int  initialSize,
int  maxOrder 
)

Definition at line 40 of file nf_Legendre.cc.

40 {
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}

Referenced by nf_Legendre_new().

◆ nf_Legendre_to_ptwXY()

ptwXYPoints * nf_Legendre_to_ptwXY ( nf_Legendre nfL,
double  accuracy,
int  biSectionMax,
int  checkForRoots,
nfu_status status 
)

Definition at line 219 of file nf_Legendre.cc.

219 {
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}
ptwXYPoints * ptwXY_createFromFunction(int n, double *xs, ptwXY_createFromFunction_callback func, void *argList, double accuracy, int checkForRoots, int biSectionMax, nfu_status *status)
Definition: ptwXY_misc.cc:40