8 STATIC double RealF2_1(
double alpha,
double beta,
double gamma,
double chi );
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 );
14 STATIC complex<double>
qg32complex(
double xl,
double xu, complex<double> (*fct)(
double) );
28 double gaunt, u, gamma2;
40 if( log10( gamma2 ) < -0.75187 )
44 gaunt = 0.551329 * ( 0.80888 - log(u) );
50 gaunt = -0.551329 * (0.5*log(gamma2) + log(u) + 0.056745);
59 ASSERT( gaunt>0. && gaunt<100. );
73 double Csum, zeta, etaf, etai, chi, gaunt, z, InitialElectronEnergy, FinalElectronEnergy, photon;
74 bool lgSutherlandOn =
false;
83 InitialElectronEnergy = sqrt(x) * TEglobal/
TE1RYD;
84 FinalElectronEnergy = photon + InitialElectronEnergy;
85 ASSERT( InitialElectronEnergy > 0. );
88 etai = z/sqrt(InitialElectronEnergy);
89 etaf = z/sqrt(FinalElectronEnergy);
92 chi = -4. * etai * etaf /
POW2( etai - etaf );
100 gaunt = 1.1027 * (1.-exp(-2.*
PI*etaf));
102 else if( etaf < 0.1*etai )
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);
110 else if( zeta > 0.5 )
113 gaunt = 1. + 0.21775*pow(zeta,-0.67) - 0.01312*pow(zeta, -1.33);
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};
121 for( i = 0; i <=9; i++ )
124 Csum += a[i]*
RealF2_1( (
double)(-i), (
double)i, 0.5, 0.5*(1.-zeta) );
126 gaunt = fabs(0.551329*(0.57721 + log(zeta/2.))*exp(
PI*zeta)*Csum);
130 else if( lgSutherlandOn )
146 fprintf(
ioQQQ,
"Uh-Oh! Gaunt is zero! Is this okay?\n");
155 #if defined(__ICC) && defined(__i386) && __INTEL_COMPILER < 910
156 #pragma optimize("", off)
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];
168 a = complex<double>( 1., -etai );
169 b = complex<double>( 0., -etaf );
175 a = complex<double>( 1., -etaf );
176 b = complex<double>( 0., -etai );
185 if( (
MAX2(NumTerms[1],NumTerms[0]) -
MIN2(NumTerms[1],NumTerms[0]) >= 2 )
186 && NumTerms[1]!=-1 && NumTerms[0]!=-1)
188 a = complex<double>( 1., -etai );
189 b = complex<double>( 0., -etaf );
192 NumTerms[0] =
MAX2(NumTerms[1],NumTerms[0])+1;
193 NumTerms[1] = NumTerms[0];
200 a = complex<double>( 1., -etaf );
201 b = complex<double>( 0., -etai );
206 ASSERT( NumTerms[0] == NumTerms[1] );
210 if( log10(abs(F[0])/abs(F[1])) + (NumRenorms[0]-NumRenorms[1])*log10(abs(
Normalization)) > 10. )
214 NumRenorms[1] = NumRenorms[0];
216 else if( log10(abs(F[1])/abs(F[0])) + (NumRenorms[1]-NumRenorms[0])*log10(abs(
Normalization)) > 10. )
220 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] )
230 if( MaxFReal > 1e50 )
232 IndexMinNumRenorms = ( NumRenorms[0] > NumRenorms[1] ) ? 1:0;
234 ++NumRenorms[IndexMinNumRenorms];
238 IndexMaxNumRenorms = ( NumRenorms[0] > NumRenorms[1] ) ? 0:1;
240 --NumRenorms[IndexMaxNumRenorms];
244 ASSERT( NumRenorms[0] == NumRenorms[1] );
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.)) );
255 complex<double> CDelta = F[0]*F[0] - F[1]*F[1];
256 double renorm =
MAX2(fabs(CDelta.real()),fabs(CDelta.imag()));
259 complex<double> NCDelta( CDelta.real()/renorm, CDelta.imag()/renorm );
261 Delta = renorm * abs( NCDelta );
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 );
276 BeckertGaunt = Delta*5.4413981*etaf*etai/fabs(etai - etaf)
277 /(1.-exp(-6.2831853*etai) )/( exp(6.2831853*etaf) - 1.);
279 while( NumRenorms[0] > 0 )
288 ASSERT( NumRenorms[0] == 0 );
295 #if defined(__ICC) && defined(__i386) && __INTEL_COMPILER < 910
296 #pragma optimize("", on)
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;
309 kfac = complex<double>( fabs((etaf-etai)/(etaf+etai)), 0. );
310 etasum = complex<double>( 0., etai + etaf );
312 GCoef = pow(kfac, etasum);
315 ASSERT( fabs(GCoef.real())<1.0 && fabs(GCoef.imag())<1.0 && ( GCoef.real()!=0. || GCoef.imag()!=0. ) );
317 for( i = 0; i <= 1; i++ )
319 a = complex<double>( i + 1., -etaf );
320 b = complex<double>( i + 1., -etai );
321 c = complex<double>( 2.*i + 2., 0. );
331 if(
MAX2(NumTerms[1],NumTerms[0]) -
MIN2(NumTerms[1],NumTerms[0]) > 2
332 && NumTerms[1]!=-1 && NumTerms[0]!=-1 )
334 NumTerms[0] =
MAX2(NumTerms[1],NumTerms[0]);
335 NumTerms[1] = NumTerms[0];
339 for( i = 0; i <= 1; i++ )
341 a = complex<double>( i + 1., -etaf );
342 b = complex<double>( i + 1., -etai );
343 c = complex<double>( 2.*i + 2., 0. );
348 ASSERT( NumTerms[0] == NumTerms[1] );
351 for( i = 0; i <= 1; i++ )
354 ASSERT( fabs(G[i].real())>0. && fabs(G[i].real())<1e100 &&
355 fabs(G[i].imag())>0. && fabs(G[i].imag())<1e100 );
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);
370 while( NumRenorms[i] > 0 )
377 ASSERT( NumRenorms[i] == 0 );
380 weightI0 =
POW2(etaf+etai);
381 weightI1 = 2.*etaf*etai*sqrt(1. + etai*etai)*sqrt(1. + etaf*etaf);
383 ComplexFactors = I[0] * ( weightI0*I[0] - weightI1*I[1] );
386 Sgaunt = 1.10266 / etai / etaf * abs( ComplexFactors );
391 #if defined(__ICC) && defined(__i386) && __INTEL_COMPILER < 910
392 #pragma optimize("", off)
396 double chi,
long *NumRenorms,
long *NumTerms )
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;
404 N_Renorms[0] = *NumRenorms;
405 N_Renorms[1] = *NumRenorms;
406 N_Terms[0] = *NumTerms;
407 N_Terms[1] = *NumTerms;
416 if( fabs(chi) < Interface1 )
418 Result =
F2_1( a, b, c, chi, &*NumRenorms, &*NumTerms );
421 else if( fabs(chi) > Interface2 )
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]);
437 if(
MAX2(N_Terms[1],N_Terms[0]) -
MIN2(N_Terms[1],N_Terms[0]) >= 2 )
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;
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] );
449 *NumTerms =
MAX2(N_Terms[1],N_Terms[0]);
456 Part[0] = F[0]/pow(chifac,a)*GammaProduct;
463 Part[1] = F[1]/pow(chifac,b)*GammaProduct;
469 if( N_Renorms[0] != N_Renorms[1] )
471 IndexMaxNumRenorms = ( N_Renorms[0] > N_Renorms[1] ) ? 0:1;
473 --N_Renorms[IndexMaxNumRenorms];
476 ASSERT( N_Renorms[0] == N_Renorms[1] );
479 *NumRenorms = N_Renorms[0];
481 Result = Part[0] + Part[1];
491 if( abs(b) > abs(a) )
493 complex<double> btemp = b;
498 CMinusBMinus1 = c-b-1.;
505 Result = Coef + FIntegral;
517 Result =
F2_1(a1,b1,c1,chi/(chi-1.),&*NumRenorms,&*NumTerms)/pow(chifac,a);
522 while( fabs(Result.real()) >= 1e50 )
533 complex<double> alpha, complex<double> beta, complex<double> gamma,
534 double chi,
long *NumRenormalizations,
long *NumTerms )
536 long i = 3, MinTerms;
537 bool lgNotConverged =
true;
538 complex<double> LastTerm, Term, Sum;
540 MinTerms =
MAX2( 3, *NumTerms );
544 ++*NumRenormalizations;
547 LastTerm = Sum*alpha*beta/gamma*chi;
559 Term = LastTerm*alpha*beta/gamma*chi/(i-1.);
564 if( Sum.real() > 1e100 )
568 ++*NumRenormalizations;
570 fprintf(
ioQQQ,
"Hypergeometric: Renormalized at term %li. Sum = %.3e %.3e\n",
571 i, Sum.real(), Sum.imag());
578 if( fabs(LastTerm.real()/Sum.real())<0.001 && fabs(LastTerm.imag()/Sum.imag())<0.001 )
579 lgNotConverged =
false;
581 if( *NumRenormalizations >= 5 )
583 fprintf(
ioQQQ,
"We've got too many (%li) renorms!\n",*NumRenormalizations );
588 }
while ( lgNotConverged || i<MinTerms );
594 #if defined(__ICC) && defined(__i386) && __INTEL_COMPILER < 910
595 #pragma optimize("", on)
602 bool lgNotConverged =
true;
603 double LastTerm, Sum;
609 LastTerm = alpha*beta*chi/gamma;
621 LastTerm *= alpha*beta*chi/gamma/(i-1.);
627 if( fabs(LastTerm/Sum)<0.001 )
628 lgNotConverged =
false;
632 }
while ( lgNotConverged );
639 return pow(v,BMinus1)*pow(1.-v,CMinusBMinus1)*pow(1.-v*GlobalCHI,MinusA);
648 complex<double> (*fct)(
double) )
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) );