cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cont_gaunt.cpp
Go to the documentation of this file.
1 /* This file is part of Cloudy and is copyright (C)1978-2008 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 
8 STATIC double RealF2_1( double alpha, double beta, double gamma, double chi );
9 STATIC complex<double> Hypergeometric2F1( complex<double> a, complex<double> b, complex<double> c,
10  double chi, long *NumRenorms, long *NumTerms );
11 STATIC complex<double> F2_1( complex<double> alpha, complex<double> beta, complex<double> gamma,
12  double chi, long *NumRenormalizations, long *NumTerms );
13 STATIC complex<double> HyperGeoInt( double v );
14 STATIC complex<double> qg32complex( double xl, double xu, complex<double> (*fct)(double) );
15 STATIC double GauntIntegrand( double y );
16 STATIC double FreeFreeGaunt( double x );
17 STATIC double DoBeckert_etal( double etai, double etaf, double chi );
18 STATIC double DoSutherland( double etai, double etaf, double chi );
19 
20 /* used to keep intermediate results from over- or underflowing. */
21 static const complex<double> Normalization(1e100, 1e100);
22 static complex<double> CMinusBMinus1, BMinus1, MinusA;
23 static double GlobalCHI;
24 static double Zglobal, HNUglobal, TEglobal;
25 
26 double 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 
64 STATIC double GauntIntegrand( double y )
65 {
66  double value;
67  value = FreeFreeGaunt( y ) * exp(-y);
68  return value;
69 }
70 
71 STATIC 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 */
161 STATIC 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 */
303 STATIC 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 */
395 STATIC 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 */
532 STATIC 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 */
599 STATIC 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 
637 STATIC 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. */
644 STATIC 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 }

Generated for cloudy by doxygen 1.8.1.2