ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
nf_Legendre.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file nf_Legendre.cc
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 
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 */
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 */
62 
63  nf_Legendre_release( Legendre );
64  nfu_free( Legendre );
65  return( NULL );
66 }
67 /*
68 ************************************************************
69 */
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 
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 */
106 
107  return( Legendre->maxOrder );
108 }
109 /*
110 ************************************************************
111 */
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 */
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 
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