cloudy trunk
Loading...
Searching...
No Matches
cont_gaunt.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3#include "cddefines.h"
4#include "physconst.h"
5#include "thirdparty.h"
6#include "continuum.h"
7
8STATIC double RealF2_1( double alpha, double beta, double gamma, double chi );
9STATIC complex<double> Hypergeometric2F1( complex<double> a, complex<double> b, complex<double> c,
10 double chi, long *NumRenorms, long *NumTerms );
11STATIC complex<double> F2_1( complex<double> alpha, complex<double> beta, complex<double> gamma,
12 double chi, long *NumRenormalizations, long *NumTerms );
13STATIC complex<double> HyperGeoInt( double v );
14STATIC complex<double> qg32complex( double xl, double xu, complex<double> (*fct)(double) );
15STATIC double GauntIntegrand( double y );
16STATIC double FreeFreeGaunt( double x );
17STATIC double DoBeckert_etal( double etai, double etaf, double chi );
18STATIC double DoSutherland( double etai, double etaf, double chi );
19
20/* used to keep intermediate results from over- or underflowing. */
21static const complex<double> Normalization(1e100, 1e100);
22static complex<double> CMinusBMinus1, BMinus1, MinusA;
23static double GlobalCHI;
24static double Zglobal, HNUglobal, TEglobal;
25
26double cont_gaunt_calc( double temp, double z, double photon )
27{
28 double gaunt, u, gamma2;
29
30 Zglobal = z;
31 HNUglobal = photon;
32 TEglobal = temp;
33
34 u = TE1RYD*photon/temp;
35 gamma2 = TE1RYD*z*z/temp;
36
37 if( log10(u)<-5. )
38 {
39 /* this cutoff is where these two approximations are equal. */
40 if( log10( gamma2 ) < -0.75187 )
41 {
42 /* given as eqn 3.2 in Hummer 88, original reference is
43 * >>refer gaunt ff Elwert, G. 1954, Zs. Naturforshung, 9A, 637 */
44 gaunt = 0.551329 * ( 0.80888 - log(u) );
45 }
46 else
47 {
48 /* given as eqn 3.1 in Hummer 88, original reference is
49 * >>refer gaunt ff Scheuer, P. A. G. 1960, MNRAS, 120, 231 */
50 gaunt = -0.551329 * (0.5*log(gamma2) + log(u) + 0.056745);
51 }
52 }
53 else
54 {
55 /* Perform integration. */
56 gaunt = qg32( 0.01, 1., GauntIntegrand );
57 gaunt += qg32( 1., 5., GauntIntegrand );
58 }
59 ASSERT( gaunt>0. && gaunt<100. );
60
61 return gaunt;
62}
63
64STATIC double GauntIntegrand( double y )
65{
66 double value;
67 value = FreeFreeGaunt( y ) * exp(-y);
68 return value;
69}
70
71STATIC double FreeFreeGaunt( double x )
72{
73 double Csum, zeta, etaf, etai, chi, gaunt, z, InitialElectronEnergy, FinalElectronEnergy, photon;
74 bool lgSutherlandOn = false;
75 long i;
76
77 z = Zglobal;
78 photon = HNUglobal;
79 ASSERT( z > 0. );
80 ASSERT( photon > 0. );
81
82 /* The final electron energy should be xkT + hv and the initial, just xkT */
83 InitialElectronEnergy = sqrt(x) * TEglobal/TE1RYD;
84 FinalElectronEnergy = photon + InitialElectronEnergy;
85 ASSERT( InitialElectronEnergy > 0. );
86
87 /* These are the free electron analog to a bound state principal quantum number.*/
88 etai = z/sqrt(InitialElectronEnergy);
89 etaf = z/sqrt(FinalElectronEnergy);
90 ASSERT( etai > 0. );
91 ASSERT( etaf > 0. );
92 chi = -4. * etai * etaf / POW2( etai - etaf );
93 zeta = etai-etaf;
94
95 if( etai>=130.)
96 {
97 if( etaf < 1.7 )
98 {
99 /* Brussard and van de Hulst (1962), as given in Hummer 88, eqn 2.23b */
100 gaunt = 1.1027 * (1.-exp(-2.*PI*etaf));
101 }
102 else if( etaf < 0.1*etai )
103 {
104 /* Hummer 88, eqn 2.23a */
105 gaunt = 1. + 0.17282604*pow(etaf,-0.67) - 0.04959570*pow(etaf,-1.33)
106 - 0.01714286*pow(etaf,-2.) + 0.00204498*pow(etaf,-2.67)
107 - 0.00243945*pow(etaf,-3.33) - 0.00120387*pow(etaf,-4.)
108 + 0.00071814*pow(etaf,-4.67) + 0.00026971*pow(etaf,-5.33);
109 }
110 else if( zeta > 0.5 )
111 {
112 /* Grant 1958, as given in Hummer 88, eqn 2.25a */
113 gaunt = 1. + 0.21775*pow(zeta,-0.67) - 0.01312*pow(zeta, -1.33);
114 }
115 else
116 {
117 double a[10] = {1.47864486, -1.72329012, 0.14420320, 0.05744888, 0.01668957,
118 0.01580779, 0.00464268, 0.00385156, 0.00116196, 0.00101967};
119
120 Csum = 0.;
121 for( i = 0; i <=9; i++ )
122 {
123 /* The Chebyshev of the first kind is just a special kind of hypergeometric. */
124 Csum += a[i]*RealF2_1( (double)(-i), (double)i, 0.5, 0.5*(1.-zeta) );
125 }
126 gaunt = fabs(0.551329*(0.57721 + log(zeta/2.))*exp(PI*zeta)*Csum);
127 ASSERT( gaunt < 10. );
128 }
129 }
130 else if( lgSutherlandOn )
131 gaunt = DoSutherland( etai, etaf, chi );
132 else
133 gaunt = DoBeckert_etal( etai, etaf, chi );
134
135 /*if( gaunt*exp(-x) > 2. && TEglobal < 1000. )
136 fprintf( ioQQQ,"ni %.3e nf %.3e chi %.3e u %.3e gam2 %.3e gaunt %.3e x %.3e\n",
137 etai, etaf, chi, TE1RYD*HNUglobal/TEglobal,
138 TE1RYD*z*z/TEglobal, gaunt, x);*/
139
142 ASSERT( gaunt > 0. && gaunt<BIGFLOAT );
143
144 if( gaunt == 0. )
145 {
146 fprintf( ioQQQ, "Uh-Oh! Gaunt is zero! Is this okay?\n");
147 /* assign some small value */
148 gaunt = 1e-5;
149 }
150
151 return gaunt;
152}
153
154
155#if defined(__ICC) && defined(__i386) && __INTEL_COMPILER < 910
156#pragma optimize("", off)
157#endif
158/************************************
159 * This part calculates the gaunt factor as per Beckert et al 2000 */
161STATIC double DoBeckert_etal( double etai, double etaf, double chi )
162{
163 double Delta, BeckertGaunt, MaxFReal, LnBeckertGaunt;
164 long NumRenorms[2]={0,0}, NumTerms[2]={0,0};
165 int IndexMinNumRenorms, IndexMaxNumRenorms;
166 complex<double> a,b,c,F[2];
167
168 a = complex<double>( 1., -etai );
169 b = complex<double>( 0., -etaf );
170 c = 1.;
171
172 /* evaluate first hypergeometric function. */
173 F[0] = Hypergeometric2F1( a, b, c, chi, &NumRenorms[0], &NumTerms[0] );
174
175 a = complex<double>( 1., -etaf );
176 b = complex<double>( 0., -etai );
177
178 /* evaluate second hypergeometric function. */
179 F[1] = Hypergeometric2F1( a, b, c, chi, &NumRenorms[1], &NumTerms[1] );
180
181 /* If there is a significant difference in the number of terms used,
182 * they should be recalculated with the max number of terms in initial calculations */
183 /* If NumTerms[i]=-1, the hypergeometric was calculated by the use of an integral instead
184 * of series summation...hence NumTerms has no meaning, and no need to recalculate. */
185 if( ( MAX2(NumTerms[1],NumTerms[0]) - MIN2(NumTerms[1],NumTerms[0]) >= 2 )
186 && NumTerms[1]!=-1 && NumTerms[0]!=-1)
187 {
188 a = complex<double>( 1., -etai );
189 b = complex<double>( 0., -etaf );
190 c = 1.;
191
192 NumTerms[0] = MAX2(NumTerms[1],NumTerms[0])+1;
193 NumTerms[1] = NumTerms[0];
194 NumRenorms[0] = 0;
195 NumRenorms[1] = 0;
196
197 /* evaluate first hypergeometric function. */
198 F[0] = Hypergeometric2F1( a, b, c, chi, &NumRenorms[0], &NumTerms[0] );
199
200 a = complex<double>( 1., -etaf );
201 b = complex<double>( 0., -etai );
202
203 /* evaluate second hypergeometric function. */
204 F[1] = Hypergeometric2F1( a, b, c, chi, &NumRenorms[1], &NumTerms[1] );
205
206 ASSERT( NumTerms[0] == NumTerms[1] );
207 }
208
209 /* if magnitude of unNormalized F's are vastly different, zero out the lesser */
210 if( log10(abs(F[0])/abs(F[1])) + (NumRenorms[0]-NumRenorms[1])*log10(abs(Normalization)) > 10. )
211 {
212 F[1] = 0.;
213 /* no longer need to keep track of differences in NumRenorms */
214 NumRenorms[1] = NumRenorms[0];
215 }
216 else if( log10(abs(F[1])/abs(F[0])) + (NumRenorms[1]-NumRenorms[0])*log10(abs(Normalization)) > 10. )
217 {
218 F[0] = 0.;
219 /* no longer need to keep track of differences in NumRenorms */
220 NumRenorms[0] = NumRenorms[1];
221 }
222
223 /* now must fix if NumRenorms[0] != NumRenorms[1], because the next calculation is the
224 * difference of squares...information is lost and cannot be recovered if this calculation
225 * is done with NumRenorms[0] != NumRenorms[1] */
226 MaxFReal = (fabs(F[1].real())>fabs(F[0].real())) ? fabs(F[1].real()):fabs(F[0].real());
227 while( NumRenorms[0] != NumRenorms[1] )
228 {
229 /* but must be very careful to prevent both overflow and underflow. */
230 if( MaxFReal > 1e50 )
231 {
232 IndexMinNumRenorms = ( NumRenorms[0] > NumRenorms[1] ) ? 1:0;
233 F[IndexMinNumRenorms] /= Normalization;
234 ++NumRenorms[IndexMinNumRenorms];
235 }
236 else
237 {
238 IndexMaxNumRenorms = ( NumRenorms[0] > NumRenorms[1] ) ? 0:1;
239 F[IndexMaxNumRenorms] = F[IndexMaxNumRenorms]*Normalization;
240 --NumRenorms[IndexMaxNumRenorms];
241 }
242 }
243
244 ASSERT( NumRenorms[0] == NumRenorms[1] );
245
246 /* Okay, now we are guaranteed (?) a small dynamic range, but may still have to renormalize */
247
248 /* Are we gonna have an overflow or underflow problem? */
249 ASSERT( (fabs(F[0].real())<1e+150) && (fabs(F[1].real())<1e+150) &&
250 (fabs(F[0].imag())<1e+150) && (fabs(F[1].real())<1e+150) );
251 ASSERT( (fabs(F[0].real())>1e-150) && ((fabs(F[0].imag())>1e-150) || (abs(F[0])==0.)) );
252 ASSERT( (fabs(F[1].real())>1e-150) && ((fabs(F[1].real())>1e-150) || (abs(F[1])==0.)) );
253
254 /* guard against spurious underflow/overflow by braindead implementations of the complex class */
255 complex<double> CDelta = F[0]*F[0] - F[1]*F[1];
256 double renorm = MAX2(fabs(CDelta.real()),fabs(CDelta.imag()));
257 ASSERT( renorm > 0. );
258 /* carefully avoid complex division here... */
259 complex<double> NCDelta( CDelta.real()/renorm, CDelta.imag()/renorm );
260
261 Delta = renorm * abs( NCDelta );
262
263 ASSERT( Delta > 0. );
264
265 /* Now multiply by the coefficient in Beckert 2000, eqn 7 */
266 if( etaf > 100. )
267 {
268 /* must compute logarithmically if etaf too big for linear computation. */
269 LnBeckertGaunt = 1.6940360 + log(Delta) + log(etaf) + log(etai) - log(fabs(etai-etaf)) - 6.2831853*etaf;
270 LnBeckertGaunt += 2. * NumRenorms[0] * log(abs(Normalization));
271 BeckertGaunt = exp( LnBeckertGaunt );
272 NumRenorms[0] = 0;
273 }
274 else
275 {
276 BeckertGaunt = Delta*5.4413981*etaf*etai/fabs(etai - etaf)
277 /(1.-exp(-6.2831853*etai) )/( exp(6.2831853*etaf) - 1.);
278
279 while( NumRenorms[0] > 0 )
280 {
281 BeckertGaunt *= abs(Normalization);
282 BeckertGaunt *= abs(Normalization);
283 ASSERT( BeckertGaunt < BIGDOUBLE );
284 --NumRenorms[0];
285 }
286 }
287
288 ASSERT( NumRenorms[0] == 0 );
289
290 /*fprintf( ioQQQ,"etai %.3e etaf %.3e u %.3e B %.3e \n",
291 etai, etaf, TE1RYD * HNUglobal / TEglobal, BeckertGaunt ); */
292
293 return BeckertGaunt;
294}
295#if defined(__ICC) && defined(__i386) && __INTEL_COMPILER < 910
296#pragma optimize("", on)
297#endif
298
299
300/************************************
301 * This part calculates the gaunt factor as per Sutherland 98 */
303STATIC double DoSutherland( double etai, double etaf, double chi )
304{
305 double Sgaunt, ICoef, weightI1, weightI0;
306 long i, NumRenorms[2]={0,0}, NumTerms[2]={0,0};
307 complex<double> a,b,c,GCoef,kfac,etasum,G[2],I[2],ComplexFactors,GammaProduct;
308
309 kfac = complex<double>( fabs((etaf-etai)/(etaf+etai)), 0. );
310 etasum = complex<double>( 0., etai + etaf );
311
312 GCoef = pow(kfac, etasum);
313 /* GCoef is a complex vector that should be contained within the unit circle.
314 * and have a non-zero magnitude. Or is it ON the unit circle? */
315 ASSERT( fabs(GCoef.real())<1.0 && fabs(GCoef.imag())<1.0 && ( GCoef.real()!=0. || GCoef.imag()!=0. ) );
316
317 for( i = 0; i <= 1; i++ )
318 {
319 a = complex<double>( i + 1., -etaf );
320 b = complex<double>( i + 1., -etai );
321 c = complex<double>( 2.*i + 2., 0. );
322
323 /* First evaluate hypergeometric function. */
324 G[i] = Hypergeometric2F1( a, b, c, chi, &NumRenorms[i], &NumTerms[i] );
325 }
326
327 /* If there is a significant difference in the number of terms used,
328 * they should be recalculated with the max number of terms in initial calculations */
329 /* If NumTerms[i]=-1, the hypergeometric was calculated by the use of an integral instead
330 * of series summation...hence NumTerms has no meaning, and no need to recalculate. */
331 if( MAX2(NumTerms[1],NumTerms[0]) - MIN2(NumTerms[1],NumTerms[0]) > 2
332 && NumTerms[1]!=-1 && NumTerms[0]!=-1 )
333 {
334 NumTerms[0] = MAX2(NumTerms[1],NumTerms[0]);
335 NumTerms[1] = NumTerms[0];
336 NumRenorms[0] = 0;
337 NumRenorms[1] = 0;
338
339 for( i = 0; i <= 1; i++ )
340 {
341 a = complex<double>( i + 1., -etaf );
342 b = complex<double>( i + 1., -etai );
343 c = complex<double>( 2.*i + 2., 0. );
344
345 G[i] = Hypergeometric2F1( a, b, c, chi, &NumRenorms[i], &NumTerms[i] );
346 }
347
348 ASSERT( NumTerms[0] == NumTerms[1] );
349 }
350
351 for( i = 0; i <= 1; i++ )
352 {
354 ASSERT( fabs(G[i].real())>0. && fabs(G[i].real())<1e100 &&
355 fabs(G[i].imag())>0. && fabs(G[i].imag())<1e100 );
356
357 /* Now multiply by the coefficient in Sutherland 98, eqn 9 */
358 G[i] *= GCoef;
359
360 /* This is the coefficient in equation 8 in Sutherland */
361 /* Karzas and Latter give tgamma(2.*i+2.), Sutherland gives tgamma(2.*i+1.) */
362 ICoef = 0.25*pow(-chi, (double)i+1.)*exp( 1.5708*fabs(etai-etaf) )/factorial(2*i);
363 GammaProduct = cdgamma(complex<double>(i+1.,etai))*cdgamma(complex<double>(i+1.,etaf));
364 ICoef *= abs(GammaProduct);
365
366 ASSERT( ICoef > 0. );
367
368 I[i] = ICoef*G[i];
369
370 while( NumRenorms[i] > 0 )
371 {
372 I[i] *= Normalization;
373 ASSERT( fabs(I[i].real()) < BIGDOUBLE && fabs(I[i].imag()) < BIGDOUBLE );
374 --NumRenorms[i];
375 }
376
377 ASSERT( NumRenorms[i] == 0 );
378 }
379
380 weightI0 = POW2(etaf+etai);
381 weightI1 = 2.*etaf*etai*sqrt(1. + etai*etai)*sqrt(1. + etaf*etaf);
382
383 ComplexFactors = I[0] * ( weightI0*I[0] - weightI1*I[1] );
384
385 /* This is Sutherland equation 13 */
386 Sgaunt = 1.10266 / etai / etaf * abs( ComplexFactors );
387
388 return Sgaunt;
389}
390
391#if defined(__ICC) && defined(__i386) && __INTEL_COMPILER < 910
392#pragma optimize("", off)
393#endif
394/* This routine is a wrapper for F2_1 */
395STATIC complex<double> Hypergeometric2F1( complex<double> a, complex<double> b, complex<double> c,
396 double chi, long *NumRenorms, long *NumTerms )
397{
398 complex<double> a1, b1, c1, a2, b2, c2, Result, Part[2], F[2];
399 complex<double> chifac, GammaProduct, Coef, FIntegral;
401 double Interface1 = 0.4, Interface2 = 10.;
402 long N_Renorms[2], N_Terms[2], IndexMaxNumRenorms, lgDoIntegral = false;
403
404 N_Renorms[0] = *NumRenorms;
405 N_Renorms[1] = *NumRenorms;
406 N_Terms[0] = *NumTerms;
407 N_Terms[1] = *NumTerms;
408
409 /* positive and zero chi are not possible. */
410 ASSERT( chi < 0. );
411
412 /* We want to be careful about evaluating the hypergeometric
413 * in the vicinity of chi=1. So we employ three different methods...*/
414
415 /* for small chi, we pass the parameters to the hypergeometric function as is. */
416 if( fabs(chi) < Interface1 )
417 {
418 Result = F2_1( a, b, c, chi, &*NumRenorms, &*NumTerms );
419 }
420 /* for large chi, we use a relation given as eqn 5 in Nicholas 89. */
421 else if( fabs(chi) > Interface2 )
422 {
423 a1 = a;
424 b1 = 1.-c+a;
425 c1 = 1.-b+a;
426
427 a2 = b;
428 b2 = 1.-c+b;
429 c2 = 1.-a+b;
430
431 chifac = -chi;
432
433 F[0] = F2_1(a1,b1,c1,1./chi,&N_Renorms[0], &N_Terms[0]);
434 F[1] = F2_1(a2,b2,c2,1./chi,&N_Renorms[1], &N_Terms[1]);
435
436 /* do it again if significant difference in number of terms. */
437 if( MAX2(N_Terms[1],N_Terms[0]) - MIN2(N_Terms[1],N_Terms[0]) >= 2 )
438 {
439 N_Terms[0] = MAX2(N_Terms[1],N_Terms[0]);
440 N_Terms[1] = N_Terms[0];
441 N_Renorms[0] = *NumRenorms;
442 N_Renorms[1] = *NumRenorms;
443
444 F[0] = F2_1(a1,b1,c1,1./chi,&N_Renorms[0], &N_Terms[0]);
445 F[1] = F2_1(a2,b2,c2,1./chi,&N_Renorms[1], &N_Terms[1]);
446 ASSERT( N_Terms[0] == N_Terms[1] );
447 }
448
449 *NumTerms = MAX2(N_Terms[1],N_Terms[0]);
450
451 /************************************************************************/
452 /* Do the first part */
453 GammaProduct = (cdgamma(b-a)/cdgamma(b))*(cdgamma(c)/cdgamma(c-a));
454
455 /* divide the hypergeometric by (-chi)^a and multiply by GammaProduct */
456 Part[0] = F[0]/pow(chifac,a)*GammaProduct;
457
458 /************************************************************************/
459 /* Do the second part */
460 GammaProduct = (cdgamma(a-b)/cdgamma(a))*(cdgamma(c)/cdgamma(c-b));
461
462 /* divide the hypergeometric by (-chi)^b and multiply by GammaProduct */
463 Part[1] = F[1]/pow(chifac,b)*GammaProduct;
464
465 /************************************************************************/
466 /* Add the two parts to get the result. */
467
468 /* First must fix it if N_Renorms[0] != N_Renorms[1] */
469 if( N_Renorms[0] != N_Renorms[1] )
470 {
471 IndexMaxNumRenorms = ( N_Renorms[0] > N_Renorms[1] ) ? 0:1;
472 Part[IndexMaxNumRenorms] *= Normalization;
473 --N_Renorms[IndexMaxNumRenorms];
474 /* Only allow at most a difference of one in number of renormalizations...
475 * otherwise something is really screwed up */
476 ASSERT( N_Renorms[0] == N_Renorms[1] );
477 }
478
479 *NumRenorms = N_Renorms[0];
480
481 Result = Part[0] + Part[1];
482 }
483 /* And for chi of order 1, we use Nicholas 89, eqn 27. */
484 else
485 {
486 /* the hypergeometric integral does not seem to work well. */
487 if( lgDoIntegral /* && fabs(chi+1.)>0.1 */)
488 {
489 /* a and b are always interchangeable, assign the lesser to b to
490 * prevent Coef from blowing up */
491 if( abs(b) > abs(a) )
492 {
493 complex<double> btemp = b;
494 b = a;
495 a = btemp;
496 }
497 Coef = cdgamma(c)/(cdgamma(b)*cdgamma(c-b));
498 CMinusBMinus1 = c-b-1.;
499 BMinus1 = b-1.;
500 MinusA = -a;
501 GlobalCHI = chi;
502 FIntegral = qg32complex( 0., 0.5, HyperGeoInt );
503 FIntegral += qg32complex( 0.5, 1., HyperGeoInt );
504
505 Result = Coef + FIntegral;
506 *NumTerms = -1;
507 *NumRenorms = 0;
508 }
509 else
510 {
511 /* Near chi=1 solution */
512 a1 = a;
513 b1 = c-b;
514 c1 = c;
515 chifac = 1.-chi;
516
517 Result = F2_1(a1,b1,c1,chi/(chi-1.),&*NumRenorms,&*NumTerms)/pow(chifac,a);
518 }
519 }
520
521 /* Limit the size of the returned value */
522 while( fabs(Result.real()) >= 1e50 )
523 {
524 Result /= Normalization;
525 ++*NumRenorms;
526 }
527
528 return Result;
529}
530
531/* This routine calculates hypergeometric functions */
532STATIC complex<double> F2_1(
533 complex<double> alpha, complex<double> beta, complex<double> gamma,
534 double chi, long *NumRenormalizations, long *NumTerms )
535{
536 long i = 3, MinTerms;
537 bool lgNotConverged = true;
538 complex<double> LastTerm, Term, Sum;
539
540 MinTerms = MAX2( 3, *NumTerms );
541
542 /* This is the first term of the hypergeometric series. */
543 Sum = 1./Normalization;
544 ++*NumRenormalizations;
545
546 /* This is the second term */
547 LastTerm = Sum*alpha*beta/gamma*chi;
548
549 Sum += LastTerm;
550
551 /* Every successive term is easily found by multiplying the last term
552 * by (alpha + i - 2)*(beta + i - 2)*chi/(gamma + i - 2)/(i-1.) */
553 do{
554 alpha += 1.;
555 beta += 1.;
556 gamma += 1.;
557
558 /* multiply old term by incremented alpha++*beta++/gamma++. Also multiply by chi/(i-1.) */
559 Term = LastTerm*alpha*beta/gamma*chi/(i-1.);
560
561 Sum += Term;
562
563 /* Renormalize if too big */
564 if( Sum.real() > 1e100 )
565 {
566 Sum /= Normalization;
567 LastTerm = Term/Normalization;
568 ++*NumRenormalizations;
569 /* notify of renormalization, and print the number of the term */
570 fprintf( ioQQQ,"Hypergeometric: Renormalized at term %li. Sum = %.3e %.3e\n",
571 i, Sum.real(), Sum.imag());
572 }
573 else
574 LastTerm = Term;
575
576 /* Declare converged if this term does not affect Sum by much */
577 /* Must do this with abs because terms alternate sign. */
578 if( fabs(LastTerm.real()/Sum.real())<0.001 && fabs(LastTerm.imag()/Sum.imag())<0.001 )
579 lgNotConverged = false;
580
581 if( *NumRenormalizations >= 5 )
582 {
583 fprintf( ioQQQ, "We've got too many (%li) renorms!\n",*NumRenormalizations );
584 }
585
586 ++i;
587
588 }while ( lgNotConverged || i<MinTerms );
589
590 *NumTerms = i;
591
592 return Sum;
593}
594#if defined(__ICC) && defined(__i386) && __INTEL_COMPILER < 910
595#pragma optimize("", on)
596#endif
597
598/* This routine calculates hypergeometric functions */
599STATIC double RealF2_1( double alpha, double beta, double gamma, double chi )
600{
601 long i = 3;
602 bool lgNotConverged = true;
603 double LastTerm, Sum;
604
605 /* This is the first term of the hypergeometric series. */
606 Sum = 1.;
607
608 /* This is the second term */
609 LastTerm = alpha*beta*chi/gamma;
610
611 Sum += LastTerm;
612
613 /* Every successive term is easily found by multiplying the last term
614 * by (alpha + i - 2)*(beta + i - 2)*chi/(gamma + i - 2)/(i-1.) */
615 do{
616 alpha++;
617 beta++;
618 gamma++;
619
620 /* multiply old term by incremented alpha++*beta++/gamma++. Also multiply by chi/(i-1.) */
621 LastTerm *= alpha*beta*chi/gamma/(i-1.);
622
623 Sum += LastTerm;
624
625 /* Declare converged if this term does not affect Sum by much */
626 /* Must do this with abs because terms alternate sign. */
627 if( fabs(LastTerm/Sum)<0.001 )
628 lgNotConverged = false;
629
630 ++i;
631
632 }while ( lgNotConverged );
633
634 return Sum;
635}
636
637STATIC complex<double> HyperGeoInt( double v )
638{
639 return pow(v,BMinus1)*pow(1.-v,CMinusBMinus1)*pow(1.-v*GlobalCHI,MinusA);
640}
641
642/*complex 32 point Gaussian quadrature, originally given to Gary F by Jim Lattimer */
643/* modified to handle complex numbers by Ryan Porter. */
644STATIC complex<double> qg32complex(
645 double xl, /*lower limit to integration range*/
646 double xu, /*upper limit to integration range*/
647 /*following is the pointer to the function that will be evaulated*/
648 complex<double> (*fct)(double) )
649{
650 double a,
651 b,
652 c;
653 complex<double> y;
654
655
656 /********************************************************************************
657 * *
658 * 32-point Gaussian quadrature *
659 * xl : the lower limit of integration *
660 * xu : the upper limit *
661 * fct : the (external) function *
662 * returns the value of the integral *
663 * *
664 * simple call to integrate sine from 0 to pi *
665 * double agn = qg32( 0., 3.141592654 , sin ); *
666 * *
667 *******************************************************************************/
668
669 a = .5*(xu + xl);
670 b = xu - xl;
671 c = .498631930924740780*b;
672 y = .35093050047350483e-2 * ( (*fct)(a+c) + (*fct)(a-c) );
673 c = b*.49280575577263417;
674 y += .8137197365452835e-2 * ( (*fct)(a+c) + (*fct)(a-c) );
675 c = b*.48238112779375322;
676 y += .1269603265463103e-1 * ( (*fct)(a+c) + (*fct)(a-c) );
677 c = b*.46745303796886984;
678 y += .17136931456510717e-1* ( (*fct)(a+c) + (*fct)(a-c) );
679 c = b*.44816057788302606;
680 y += .21417949011113340e-1* ( (*fct)(a+c) + (*fct)(a-c) );
681 c = b*.42468380686628499;
682 y += .25499029631188088e-1* ( (*fct)(a+c) + (*fct)(a-c) );
683 c = b*.3972418979839712;
684 y += .29342046739267774e-1* ( (*fct)(a+c) + (*fct)(a-c) );
685 c = b*.36609105937014484;
686 y += .32911111388180923e-1* ( (*fct)(a+c) + (*fct)(a-c) );
687 c = b*.3315221334651076;
688 y += .36172897054424253e-1* ( (*fct)(a+c) + (*fct)(a-c) );
689 c = b*.29385787862038116;
690 y += .39096947893535153e-1* ( (*fct)(a+c) + (*fct)(a-c) );
691 c = b*.2534499544661147;
692 y += .41655962113473378e-1* ( (*fct)(a+c) + (*fct)(a-c) );
693 c = b*.21067563806531767;
694 y += .43826046502201906e-1* ( (*fct)(a+c) + (*fct)(a-c) );
695 c = b*.16593430114106382;
696 y += .45586939347881942e-1* ( (*fct)(a+c) + (*fct)(a-c) );
697 c = b*.11964368112606854;
698 y += .46922199540402283e-1* ( (*fct)(a+c) + (*fct)(a-c) );
699 c = b*.7223598079139825e-1;
700 y += .47819360039637430e-1* ( (*fct)(a+c) + (*fct)(a-c) );
701 c = b*.24153832843869158e-1;
702 y += .4827004425736390e-1 * ( (*fct)(a+c) + (*fct)(a-c) );
703 y *= b;
704
705 /* the answer */
706
707 return( y );
708}
static double a2[63]
static realnum b1[83]
static double b2[63]
static double a1[83]
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
#define MIN2
Definition cddefines.h:761
double qg32(double, double, double(*)(double))
Definition service.cpp:1053
#define POW2
Definition cddefines.h:929
#define MAX2
Definition cddefines.h:782
STATIC double DoSutherland(double etai, double etaf, double chi)
STATIC complex< double > qg32complex(double xl, double xu, complex< double >(*fct)(double))
double cont_gaunt_calc(double temp, double z, double photon)
static double GlobalCHI
static complex< double > BMinus1
STATIC complex< double > Hypergeometric2F1(complex< double > a, complex< double > b, complex< double > c, double chi, long *NumRenorms, long *NumTerms)
static const complex< double > Normalization(1e100, 1e100)
static double HNUglobal
STATIC double FreeFreeGaunt(double x)
static complex< double > MinusA
STATIC complex< double > HyperGeoInt(double v)
STATIC complex< double > F2_1(complex< double > alpha, complex< double > beta, complex< double > gamma, double chi, long *NumRenormalizations, long *NumTerms)
STATIC double RealF2_1(double alpha, double beta, double gamma, double chi)
STATIC double GauntIntegrand(double y)
static double Zglobal
static double TEglobal
STATIC double DoBeckert_etal(double etai, double etaf, double chi)
static complex< double > CMinusBMinus1
UNUSED const realnum BIGFLOAT
Definition cpu.h:189
const double BIGDOUBLE
Definition cpu.h:194
UNUSED const double PI
Definition physconst.h:29
UNUSED const double TE1RYD
Definition physconst.h:183
double factorial(long n)
complex< double > cdgamma(complex< double > x)