cloudy
trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 #include "cddefines.h" 00004 #include "physconst.h" 00005 #include "thirdparty.h" 00006 #include "continuum.h" 00007 00008 STATIC double RealF2_1( double alpha, double beta, double gamma, double chi ); 00009 STATIC complex<double> Hypergeometric2F1( complex<double> a, complex<double> b, complex<double> c, 00010 double chi, long *NumRenorms, long *NumTerms ); 00011 STATIC complex<double> F2_1( complex<double> alpha, complex<double> beta, complex<double> gamma, 00012 double chi, long *NumRenormalizations, long *NumTerms ); 00013 STATIC complex<double> HyperGeoInt( double v ); 00014 STATIC complex<double> qg32complex( double xl, double xu, complex<double> (*fct)(double) ); 00015 STATIC double GauntIntegrand( double y ); 00016 STATIC double FreeFreeGaunt( double x ); 00017 STATIC double DoBeckert_etal( double etai, double etaf, double chi ); 00018 STATIC double DoSutherland( double etai, double etaf, double chi ); 00019 00020 /* used to keep intermediate results from over- or underflowing. */ 00021 static const complex<double> Normalization(1e100, 1e100); 00022 static complex<double> CMinusBMinus1, BMinus1, MinusA; 00023 static double GlobalCHI; 00024 static double Zglobal, HNUglobal, TEglobal; 00025 00026 double cont_gaunt_calc( double temp, double z, double photon ) 00027 { 00028 double gaunt, u, gamma2; 00029 00030 Zglobal = z; 00031 HNUglobal = photon; 00032 TEglobal = temp; 00033 00034 u = TE1RYD*photon/temp; 00035 gamma2 = TE1RYD*z*z/temp; 00036 00037 if( log10(u)<-5. ) 00038 { 00039 /* this cutoff is where these two approximations are equal. */ 00040 if( log10( gamma2 ) < -0.75187 ) 00041 { 00042 /* given as eqn 3.2 in Hummer 88, original reference is 00043 * >>refer gaunt ff Elwert, G. 1954, Zs. Naturforshung, 9A, 637 */ 00044 gaunt = 0.551329 * ( 0.80888 - log(u) ); 00045 } 00046 else 00047 { 00048 /* given as eqn 3.1 in Hummer 88, original reference is 00049 * >>refer gaunt ff Scheuer, P.A.G. 1960, MNRAS, 120, 231 */ 00050 gaunt = -0.551329 * (0.5*log(gamma2) + log(u) + 0.056745); 00051 } 00052 } 00053 else 00054 { 00055 /* Perform integration. */ 00056 gaunt = qg32( 0.01, 1., GauntIntegrand ); 00057 gaunt += qg32( 1., 5., GauntIntegrand ); 00058 } 00059 ASSERT( gaunt>0. && gaunt<100. ); 00060 00061 return gaunt; 00062 } 00063 00064 STATIC double GauntIntegrand( double y ) 00065 { 00066 double value; 00067 value = FreeFreeGaunt( y ) * exp(-y); 00068 return value; 00069 } 00070 00071 STATIC double FreeFreeGaunt( double x ) 00072 { 00073 double Csum, zeta, etaf, etai, chi, gaunt, z, InitialElectronEnergy, FinalElectronEnergy, photon; 00074 bool lgSutherlandOn = false; 00075 long i; 00076 00077 z = Zglobal; 00078 photon = HNUglobal; 00079 ASSERT( z > 0. ); 00080 ASSERT( photon > 0. ); 00081 00082 /* The final electron energy should be xkT + hv and the initial, just xkT */ 00083 InitialElectronEnergy = sqrt(x) * TEglobal/TE1RYD; 00084 FinalElectronEnergy = photon + InitialElectronEnergy; 00085 ASSERT( InitialElectronEnergy > 0. ); 00086 00087 /* These are the free electron analog to a bound state principal quantum number.*/ 00088 etai = z/sqrt(InitialElectronEnergy); 00089 etaf = z/sqrt(FinalElectronEnergy); 00090 ASSERT( etai > 0. ); 00091 ASSERT( etaf > 0. ); 00092 chi = -4. * etai * etaf / POW2( etai - etaf ); 00093 zeta = etai-etaf; 00094 00095 if( etai>=130.) 00096 { 00097 if( etaf < 1.7 ) 00098 { 00099 /* Brussard and van de Hulst (1962), as given in Hummer 88, eqn 2.23b */ 00100 gaunt = 1.1027 * (1.-exp(-2.*PI*etaf)); 00101 } 00102 else if( etaf < 0.1*etai ) 00103 { 00104 /* Hummer 88, eqn 2.23a */ 00105 gaunt = 1. + 0.17282604*pow(etaf,-0.67) - 0.04959570*pow(etaf,-1.33) 00106 - 0.01714286*pow(etaf,-2.) + 0.00204498*pow(etaf,-2.67) 00107 - 0.00243945*pow(etaf,-3.33) - 0.00120387*pow(etaf,-4.) 00108 + 0.00071814*pow(etaf,-4.67) + 0.00026971*pow(etaf,-5.33); 00109 } 00110 else if( zeta > 0.5 ) 00111 { 00112 /* Grant 1958, as given in Hummer 88, eqn 2.25a */ 00113 gaunt = 1. + 0.21775*pow(zeta,-0.67) - 0.01312*pow(zeta, -1.33); 00114 } 00115 else 00116 { 00117 double a[10] = {1.47864486, -1.72329012, 0.14420320, 0.05744888, 0.01668957, 00118 0.01580779, 0.00464268, 0.00385156, 0.00116196, 0.00101967}; 00119 00120 Csum = 0.; 00121 for( i = 0; i <=9; i++ ) 00122 { 00123 /* The Chebyshev of the first kind is just a special kind of hypergeometric. */ 00124 Csum += a[i]*RealF2_1( (double)(-i), (double)i, 0.5, 0.5*(1.-zeta) ); 00125 } 00126 gaunt = fabs(0.551329*(0.57721 + log(zeta/2.))*exp(PI*zeta)*Csum); 00127 ASSERT( gaunt < 10. ); 00128 } 00129 } 00130 else if( lgSutherlandOn ) 00131 gaunt = DoSutherland( etai, etaf, chi ); 00132 else 00133 gaunt = DoBeckert_etal( etai, etaf, chi ); 00134 00135 /*if( gaunt*exp(-x) > 2. && TEglobal < 1000. ) 00136 fprintf( ioQQQ,"ni %.3e nf %.3e chi %.3e u %.3e gam2 %.3e gaunt %.3e x %.3e\n", 00137 etai, etaf, chi, TE1RYD*HNUglobal/TEglobal, 00138 TE1RYD*z*z/TEglobal, gaunt, x);*/ 00139 00142 ASSERT( gaunt > 0. && gaunt<BIGFLOAT ); 00143 00144 if( gaunt == 0. ) 00145 { 00146 fprintf( ioQQQ, "Uh-Oh! Gaunt is zero! Is this okay?\n"); 00147 /* assign some small value */ 00148 gaunt = 1e-5; 00149 } 00150 00151 return gaunt; 00152 } 00153 00154 00155 #if defined(__ICC) && defined(__i386) && __INTEL_COMPILER < 910 00156 #pragma optimize("", off) 00157 #endif 00158 /************************************ 00159 * This part calculates the gaunt factor as per Beckert et al 2000 */ 00161 STATIC double DoBeckert_etal( double etai, double etaf, double chi ) 00162 { 00163 double Delta, BeckertGaunt, MaxFReal, LnBeckertGaunt; 00164 long NumRenorms[2]={0,0}, NumTerms[2]={0,0}; 00165 int IndexMinNumRenorms, IndexMaxNumRenorms; 00166 complex<double> a,b,c,F[2]; 00167 00168 a = complex<double>( 1., -etai ); 00169 b = complex<double>( 0., -etaf ); 00170 c = 1.; 00171 00172 /* evaluate first hypergeometric function. */ 00173 F[0] = Hypergeometric2F1( a, b, c, chi, &NumRenorms[0], &NumTerms[0] ); 00174 00175 a = complex<double>( 1., -etaf ); 00176 b = complex<double>( 0., -etai ); 00177 00178 /* evaluate second hypergeometric function. */ 00179 F[1] = Hypergeometric2F1( a, b, c, chi, &NumRenorms[1], &NumTerms[1] ); 00180 00181 /* If there is a significant difference in the number of terms used, 00182 * they should be recalculated with the max number of terms in initial calculations */ 00183 /* If NumTerms[i]=-1, the hypergeometric was calculated by the use of an integral instead 00184 * of series summation...hence NumTerms has no meaning, and no need to recalculate. */ 00185 if( ( MAX2(NumTerms[1],NumTerms[0]) - MIN2(NumTerms[1],NumTerms[0]) >= 2 ) 00186 && NumTerms[1]!=-1 && NumTerms[0]!=-1) 00187 { 00188 a = complex<double>( 1., -etai ); 00189 b = complex<double>( 0., -etaf ); 00190 c = 1.; 00191 00192 NumTerms[0] = MAX2(NumTerms[1],NumTerms[0])+1; 00193 NumTerms[1] = NumTerms[0]; 00194 NumRenorms[0] = 0; 00195 NumRenorms[1] = 0; 00196 00197 /* evaluate first hypergeometric function. */ 00198 F[0] = Hypergeometric2F1( a, b, c, chi, &NumRenorms[0], &NumTerms[0] ); 00199 00200 a = complex<double>( 1., -etaf ); 00201 b = complex<double>( 0., -etai ); 00202 00203 /* evaluate second hypergeometric function. */ 00204 F[1] = Hypergeometric2F1( a, b, c, chi, &NumRenorms[1], &NumTerms[1] ); 00205 00206 ASSERT( NumTerms[0] == NumTerms[1] ); 00207 } 00208 00209 /* if magnitude of unNormalized F's are vastly different, zero out the lesser */ 00210 if( log10(abs(F[0])/abs(F[1])) + (NumRenorms[0]-NumRenorms[1])*log10(abs(Normalization)) > 10. ) 00211 { 00212 F[1] = 0.; 00213 /* no longer need to keep track of differences in NumRenorms */ 00214 NumRenorms[1] = NumRenorms[0]; 00215 } 00216 else if( log10(abs(F[1])/abs(F[0])) + (NumRenorms[1]-NumRenorms[0])*log10(abs(Normalization)) > 10. ) 00217 { 00218 F[0] = 0.; 00219 /* no longer need to keep track of differences in NumRenorms */ 00220 NumRenorms[0] = NumRenorms[1]; 00221 } 00222 00223 /* now must fix if NumRenorms[0] != NumRenorms[1], because the next calculation is the 00224 * difference of squares...information is lost and cannot be recovered if this calculation 00225 * is done with NumRenorms[0] != NumRenorms[1] */ 00226 MaxFReal = (fabs(F[1].real())>fabs(F[0].real())) ? fabs(F[1].real()):fabs(F[0].real()); 00227 while( NumRenorms[0] != NumRenorms[1] ) 00228 { 00229 /* but must be very careful to prevent both overflow and underflow. */ 00230 if( MaxFReal > 1e50 ) 00231 { 00232 IndexMinNumRenorms = ( NumRenorms[0] > NumRenorms[1] ) ? 1:0; 00233 F[IndexMinNumRenorms] /= Normalization; 00234 ++NumRenorms[IndexMinNumRenorms]; 00235 } 00236 else 00237 { 00238 IndexMaxNumRenorms = ( NumRenorms[0] > NumRenorms[1] ) ? 0:1; 00239 F[IndexMaxNumRenorms] = F[IndexMaxNumRenorms]*Normalization; 00240 --NumRenorms[IndexMaxNumRenorms]; 00241 } 00242 } 00243 00244 ASSERT( NumRenorms[0] == NumRenorms[1] ); 00245 00246 /* Okay, now we are guaranteed (?) a small dynamic range, but may still have to renormalize */ 00247 00248 /* Are we gonna have an overflow or underflow problem? */ 00249 ASSERT( (fabs(F[0].real())<1e+150) && (fabs(F[1].real())<1e+150) && 00250 (fabs(F[0].imag())<1e+150) && (fabs(F[1].real())<1e+150) ); 00251 ASSERT( (fabs(F[0].real())>1e-150) && ((fabs(F[0].imag())>1e-150) || (abs(F[0])==0.)) ); 00252 ASSERT( (fabs(F[1].real())>1e-150) && ((fabs(F[1].real())>1e-150) || (abs(F[1])==0.)) ); 00253 00254 /* guard against spurious underflow/overflow by braindead implementations of the complex class */ 00255 complex<double> CDelta = F[0]*F[0] - F[1]*F[1]; 00256 double renorm = MAX2(fabs(CDelta.real()),fabs(CDelta.imag())); 00257 ASSERT( renorm > 0. ); 00258 /* carefully avoid complex division here... */ 00259 complex<double> NCDelta( CDelta.real()/renorm, CDelta.imag()/renorm ); 00260 00261 Delta = renorm * abs( NCDelta ); 00262 00263 ASSERT( Delta > 0. ); 00264 00265 /* Now multiply by the coefficient in Beckert 2000, eqn 7 */ 00266 if( etaf > 100. ) 00267 { 00268 /* must compute logarithmically if etaf too big for linear computation. */ 00269 LnBeckertGaunt = 1.6940360 + log(Delta) + log(etaf) + log(etai) - log(fabs(etai-etaf)) - 6.2831853*etaf; 00270 LnBeckertGaunt += 2. * NumRenorms[0] * log(abs(Normalization)); 00271 BeckertGaunt = exp( LnBeckertGaunt ); 00272 NumRenorms[0] = 0; 00273 } 00274 else 00275 { 00276 BeckertGaunt = Delta*5.4413981*etaf*etai/fabs(etai - etaf) 00277 /(1.-exp(-6.2831853*etai) )/( exp(6.2831853*etaf) - 1.); 00278 00279 while( NumRenorms[0] > 0 ) 00280 { 00281 BeckertGaunt *= abs(Normalization); 00282 BeckertGaunt *= abs(Normalization); 00283 ASSERT( BeckertGaunt < BIGDOUBLE ); 00284 --NumRenorms[0]; 00285 } 00286 } 00287 00288 ASSERT( NumRenorms[0] == 0 ); 00289 00290 /*fprintf( ioQQQ,"etai %.3e etaf %.3e u %.3e B %.3e \n", 00291 etai, etaf, TE1RYD * HNUglobal / TEglobal, BeckertGaunt ); */ 00292 00293 return BeckertGaunt; 00294 } 00295 #if defined(__ICC) && defined(__i386) && __INTEL_COMPILER < 910 00296 #pragma optimize("", on) 00297 #endif 00298 00299 00300 /************************************ 00301 * This part calculates the gaunt factor as per Sutherland 98 */ 00303 STATIC double DoSutherland( double etai, double etaf, double chi ) 00304 { 00305 double Sgaunt, ICoef, weightI1, weightI0; 00306 long i, NumRenorms[2]={0,0}, NumTerms[2]={0,0}; 00307 complex<double> a,b,c,GCoef,kfac,etasum,G[2],I[2],ComplexFactors,GammaProduct; 00308 00309 kfac = complex<double>( fabs((etaf-etai)/(etaf+etai)), 0. ); 00310 etasum = complex<double>( 0., etai + etaf ); 00311 00312 GCoef = pow(kfac, etasum); 00313 /* GCoef is a complex vector that should be contained within the unit circle. 00314 * and have a non-zero magnitude. Or is it ON the unit circle? */ 00315 ASSERT( fabs(GCoef.real())<1.0 && fabs(GCoef.imag())<1.0 && ( GCoef.real()!=0. || GCoef.imag()!=0. ) ); 00316 00317 for( i = 0; i <= 1; i++ ) 00318 { 00319 a = complex<double>( i + 1., -etaf ); 00320 b = complex<double>( i + 1., -etai ); 00321 c = complex<double>( 2.*i + 2., 0. ); 00322 00323 /* First evaluate hypergeometric function. */ 00324 G[i] = Hypergeometric2F1( a, b, c, chi, &NumRenorms[i], &NumTerms[i] ); 00325 } 00326 00327 /* If there is a significant difference in the number of terms used, 00328 * they should be recalculated with the max number of terms in initial calculations */ 00329 /* If NumTerms[i]=-1, the hypergeometric was calculated by the use of an integral instead 00330 * of series summation...hence NumTerms has no meaning, and no need to recalculate. */ 00331 if( MAX2(NumTerms[1],NumTerms[0]) - MIN2(NumTerms[1],NumTerms[0]) > 2 00332 && NumTerms[1]!=-1 && NumTerms[0]!=-1 ) 00333 { 00334 NumTerms[0] = MAX2(NumTerms[1],NumTerms[0]); 00335 NumTerms[1] = NumTerms[0]; 00336 NumRenorms[0] = 0; 00337 NumRenorms[1] = 0; 00338 00339 for( i = 0; i <= 1; i++ ) 00340 { 00341 a = complex<double>( i + 1., -etaf ); 00342 b = complex<double>( i + 1., -etai ); 00343 c = complex<double>( 2.*i + 2., 0. ); 00344 00345 G[i] = Hypergeometric2F1( a, b, c, chi, &NumRenorms[i], &NumTerms[i] ); 00346 } 00347 00348 ASSERT( NumTerms[0] == NumTerms[1] ); 00349 } 00350 00351 for( i = 0; i <= 1; i++ ) 00352 { 00354 ASSERT( fabs(G[i].real())>0. && fabs(G[i].real())<1e100 && 00355 fabs(G[i].imag())>0. && fabs(G[i].imag())<1e100 ); 00356 00357 /* Now multiply by the coefficient in Sutherland 98, eqn 9 */ 00358 G[i] *= GCoef; 00359 00360 /* This is the coefficient in equation 8 in Sutherland */ 00361 /* Karzas and Latter give tgamma(2.*i+2.), Sutherland gives tgamma(2.*i+1.) */ 00362 ICoef = 0.25*pow(-chi, (double)i+1.)*exp( 1.5708*fabs(etai-etaf) )/factorial(2*i); 00363 GammaProduct = cdgamma(complex<double>(i+1.,etai))*cdgamma(complex<double>(i+1.,etaf)); 00364 ICoef *= abs(GammaProduct); 00365 00366 ASSERT( ICoef > 0. ); 00367 00368 I[i] = ICoef*G[i]; 00369 00370 while( NumRenorms[i] > 0 ) 00371 { 00372 I[i] *= Normalization; 00373 ASSERT( fabs(I[i].real()) < BIGDOUBLE && fabs(I[i].imag()) < BIGDOUBLE ); 00374 --NumRenorms[i]; 00375 } 00376 00377 ASSERT( NumRenorms[i] == 0 ); 00378 } 00379 00380 weightI0 = POW2(etaf+etai); 00381 weightI1 = 2.*etaf*etai*sqrt(1. + etai*etai)*sqrt(1. + etaf*etaf); 00382 00383 ComplexFactors = I[0] * ( weightI0*I[0] - weightI1*I[1] ); 00384 00385 /* This is Sutherland equation 13 */ 00386 Sgaunt = 1.10266 / etai / etaf * abs( ComplexFactors ); 00387 00388 return Sgaunt; 00389 } 00390 00391 #if defined(__ICC) && defined(__i386) && __INTEL_COMPILER < 910 00392 #pragma optimize("", off) 00393 #endif 00394 /* This routine is a wrapper for F2_1 */ 00395 STATIC complex<double> Hypergeometric2F1( complex<double> a, complex<double> b, complex<double> c, 00396 double chi, long *NumRenorms, long *NumTerms ) 00397 { 00398 complex<double> a1, b1, c1, a2, b2, c2, Result, Part[2], F[2]; 00399 complex<double> chifac, GammaProduct, Coef, FIntegral; 00401 double Interface1 = 0.4, Interface2 = 10.; 00402 long N_Renorms[2], N_Terms[2], IndexMaxNumRenorms, lgDoIntegral = false; 00403 00404 N_Renorms[0] = *NumRenorms; 00405 N_Renorms[1] = *NumRenorms; 00406 N_Terms[0] = *NumTerms; 00407 N_Terms[1] = *NumTerms; 00408 00409 /* positive and zero chi are not possible. */ 00410 ASSERT( chi < 0. ); 00411 00412 /* We want to be careful about evaluating the hypergeometric 00413 * in the vicinity of chi=1. So we employ three different methods...*/ 00414 00415 /* for small chi, we pass the parameters to the hypergeometric function as is. */ 00416 if( fabs(chi) < Interface1 ) 00417 { 00418 Result = F2_1( a, b, c, chi, &*NumRenorms, &*NumTerms ); 00419 } 00420 /* for large chi, we use a relation given as eqn 5 in Nicholas 89. */ 00421 else if( fabs(chi) > Interface2 ) 00422 { 00423 a1 = a; 00424 b1 = 1.-c+a; 00425 c1 = 1.-b+a; 00426 00427 a2 = b; 00428 b2 = 1.-c+b; 00429 c2 = 1.-a+b; 00430 00431 chifac = -chi; 00432 00433 F[0] = F2_1(a1,b1,c1,1./chi,&N_Renorms[0], &N_Terms[0]); 00434 F[1] = F2_1(a2,b2,c2,1./chi,&N_Renorms[1], &N_Terms[1]); 00435 00436 /* do it again if significant difference in number of terms. */ 00437 if( MAX2(N_Terms[1],N_Terms[0]) - MIN2(N_Terms[1],N_Terms[0]) >= 2 ) 00438 { 00439 N_Terms[0] = MAX2(N_Terms[1],N_Terms[0]); 00440 N_Terms[1] = N_Terms[0]; 00441 N_Renorms[0] = *NumRenorms; 00442 N_Renorms[1] = *NumRenorms; 00443 00444 F[0] = F2_1(a1,b1,c1,1./chi,&N_Renorms[0], &N_Terms[0]); 00445 F[1] = F2_1(a2,b2,c2,1./chi,&N_Renorms[1], &N_Terms[1]); 00446 ASSERT( N_Terms[0] == N_Terms[1] ); 00447 } 00448 00449 *NumTerms = MAX2(N_Terms[1],N_Terms[0]); 00450 00451 /************************************************************************/ 00452 /* Do the first part */ 00453 GammaProduct = (cdgamma(b-a)/cdgamma(b))*(cdgamma(c)/cdgamma(c-a)); 00454 00455 /* divide the hypergeometric by (-chi)^a and multiply by GammaProduct */ 00456 Part[0] = F[0]/pow(chifac,a)*GammaProduct; 00457 00458 /************************************************************************/ 00459 /* Do the second part */ 00460 GammaProduct = (cdgamma(a-b)/cdgamma(a))*(cdgamma(c)/cdgamma(c-b)); 00461 00462 /* divide the hypergeometric by (-chi)^b and multiply by GammaProduct */ 00463 Part[1] = F[1]/pow(chifac,b)*GammaProduct; 00464 00465 /************************************************************************/ 00466 /* Add the two parts to get the result. */ 00467 00468 /* First must fix it if N_Renorms[0] != N_Renorms[1] */ 00469 if( N_Renorms[0] != N_Renorms[1] ) 00470 { 00471 IndexMaxNumRenorms = ( N_Renorms[0] > N_Renorms[1] ) ? 0:1; 00472 Part[IndexMaxNumRenorms] *= Normalization; 00473 --N_Renorms[IndexMaxNumRenorms]; 00474 /* Only allow at most a difference of one in number of renormalizations... 00475 * otherwise something is really screwed up */ 00476 ASSERT( N_Renorms[0] == N_Renorms[1] ); 00477 } 00478 00479 *NumRenorms = N_Renorms[0]; 00480 00481 Result = Part[0] + Part[1]; 00482 } 00483 /* And for chi of order 1, we use Nicholas 89, eqn 27. */ 00484 else 00485 { 00486 /* the hypergeometric integral does not seem to work well. */ 00487 if( lgDoIntegral /* && fabs(chi+1.)>0.1 */) 00488 { 00489 /* a and b are always interchangeable, assign the lesser to b to 00490 * prevent Coef from blowing up */ 00491 if( abs(b) > abs(a) ) 00492 { 00493 complex<double> btemp = b; 00494 b = a; 00495 a = btemp; 00496 } 00497 Coef = cdgamma(c)/(cdgamma(b)*cdgamma(c-b)); 00498 CMinusBMinus1 = c-b-1.; 00499 BMinus1 = b-1.; 00500 MinusA = -a; 00501 GlobalCHI = chi; 00502 FIntegral = qg32complex( 0., 0.5, HyperGeoInt ); 00503 FIntegral += qg32complex( 0.5, 1., HyperGeoInt ); 00504 00505 Result = Coef + FIntegral; 00506 *NumTerms = -1; 00507 *NumRenorms = 0; 00508 } 00509 else 00510 { 00511 /* Near chi=1 solution */ 00512 a1 = a; 00513 b1 = c-b; 00514 c1 = c; 00515 chifac = 1.-chi; 00516 00517 Result = F2_1(a1,b1,c1,chi/(chi-1.),&*NumRenorms,&*NumTerms)/pow(chifac,a); 00518 } 00519 } 00520 00521 /* Limit the size of the returned value */ 00522 while( fabs(Result.real()) >= 1e50 ) 00523 { 00524 Result /= Normalization; 00525 ++*NumRenorms; 00526 } 00527 00528 return Result; 00529 } 00530 00531 /* This routine calculates hypergeometric functions */ 00532 STATIC complex<double> F2_1( 00533 complex<double> alpha, complex<double> beta, complex<double> gamma, 00534 double chi, long *NumRenormalizations, long *NumTerms ) 00535 { 00536 long i = 3, MinTerms; 00537 bool lgNotConverged = true; 00538 complex<double> LastTerm, Term, Sum; 00539 00540 MinTerms = MAX2( 3, *NumTerms ); 00541 00542 /* This is the first term of the hypergeometric series. */ 00543 Sum = 1./Normalization; 00544 ++*NumRenormalizations; 00545 00546 /* This is the second term */ 00547 LastTerm = Sum*alpha*beta/gamma*chi; 00548 00549 Sum += LastTerm; 00550 00551 /* Every successive term is easily found by multiplying the last term 00552 * by (alpha + i - 2)*(beta + i - 2)*chi/(gamma + i - 2)/(i-1.) */ 00553 do{ 00554 alpha += 1.; 00555 beta += 1.; 00556 gamma += 1.; 00557 00558 /* multiply old term by incremented alpha++*beta++/gamma++. Also multiply by chi/(i-1.) */ 00559 Term = LastTerm*alpha*beta/gamma*chi/(i-1.); 00560 00561 Sum += Term; 00562 00563 /* Renormalize if too big */ 00564 if( Sum.real() > 1e100 ) 00565 { 00566 Sum /= Normalization; 00567 LastTerm = Term/Normalization; 00568 ++*NumRenormalizations; 00569 /* notify of renormalization, and print the number of the term */ 00570 fprintf( ioQQQ,"Hypergeometric: Renormalized at term %li. Sum = %.3e %.3e\n", 00571 i, Sum.real(), Sum.imag()); 00572 } 00573 else 00574 LastTerm = Term; 00575 00576 /* Declare converged if this term does not affect Sum by much */ 00577 /* Must do this with abs because terms alternate sign. */ 00578 if( fabs(LastTerm.real()/Sum.real())<0.001 && fabs(LastTerm.imag()/Sum.imag())<0.001 ) 00579 lgNotConverged = false; 00580 00581 if( *NumRenormalizations >= 5 ) 00582 { 00583 fprintf( ioQQQ, "We've got too many (%li) renorms!\n",*NumRenormalizations ); 00584 } 00585 00586 ++i; 00587 00588 }while ( lgNotConverged || i<MinTerms ); 00589 00590 *NumTerms = i; 00591 00592 return Sum; 00593 } 00594 #if defined(__ICC) && defined(__i386) && __INTEL_COMPILER < 910 00595 #pragma optimize("", on) 00596 #endif 00597 00598 /* This routine calculates hypergeometric functions */ 00599 STATIC double RealF2_1( double alpha, double beta, double gamma, double chi ) 00600 { 00601 long i = 3; 00602 bool lgNotConverged = true; 00603 double LastTerm, Sum; 00604 00605 /* This is the first term of the hypergeometric series. */ 00606 Sum = 1.; 00607 00608 /* This is the second term */ 00609 LastTerm = alpha*beta*chi/gamma; 00610 00611 Sum += LastTerm; 00612 00613 /* Every successive term is easily found by multiplying the last term 00614 * by (alpha + i - 2)*(beta + i - 2)*chi/(gamma + i - 2)/(i-1.) */ 00615 do{ 00616 alpha++; 00617 beta++; 00618 gamma++; 00619 00620 /* multiply old term by incremented alpha++*beta++/gamma++. Also multiply by chi/(i-1.) */ 00621 LastTerm *= alpha*beta*chi/gamma/(i-1.); 00622 00623 Sum += LastTerm; 00624 00625 /* Declare converged if this term does not affect Sum by much */ 00626 /* Must do this with abs because terms alternate sign. */ 00627 if( fabs(LastTerm/Sum)<0.001 ) 00628 lgNotConverged = false; 00629 00630 ++i; 00631 00632 }while ( lgNotConverged ); 00633 00634 return Sum; 00635 } 00636 00637 STATIC complex<double> HyperGeoInt( double v ) 00638 { 00639 return pow(v,BMinus1)*pow(1.-v,CMinusBMinus1)*pow(1.-v*GlobalCHI,MinusA); 00640 } 00641 00642 /*complex 32 point Gaussian quadrature, originally given to Gary F by Jim Lattimer */ 00643 /* modified to handle complex numbers by Ryan Porter. */ 00644 STATIC complex<double> qg32complex( 00645 double xl, /*lower limit to integration range*/ 00646 double xu, /*upper limit to integration range*/ 00647 /*following is the pointer to the function that will be evaulated*/ 00648 complex<double> (*fct)(double) ) 00649 { 00650 double a, 00651 b, 00652 c; 00653 complex<double> y; 00654 00655 00656 /******************************************************************************** 00657 * * 00658 * 32-point Gaussian quadrature * 00659 * xl : the lower limit of integration * 00660 * xu : the upper limit * 00661 * fct : the (external) function * 00662 * returns the value of the integral * 00663 * * 00664 * simple call to integrate sine from 0 to pi * 00665 * double agn = qg32( 0., 3.141592654 , sin ); * 00666 * * 00667 *******************************************************************************/ 00668 00669 a = .5*(xu + xl); 00670 b = xu - xl; 00671 c = .498631930924740780*b; 00672 y = .35093050047350483e-2 * ( (*fct)(a+c) + (*fct)(a-c) ); 00673 c = b*.49280575577263417; 00674 y += .8137197365452835e-2 * ( (*fct)(a+c) + (*fct)(a-c) ); 00675 c = b*.48238112779375322; 00676 y += .1269603265463103e-1 * ( (*fct)(a+c) + (*fct)(a-c) ); 00677 c = b*.46745303796886984; 00678 y += .17136931456510717e-1* ( (*fct)(a+c) + (*fct)(a-c) ); 00679 c = b*.44816057788302606; 00680 y += .21417949011113340e-1* ( (*fct)(a+c) + (*fct)(a-c) ); 00681 c = b*.42468380686628499; 00682 y += .25499029631188088e-1* ( (*fct)(a+c) + (*fct)(a-c) ); 00683 c = b*.3972418979839712; 00684 y += .29342046739267774e-1* ( (*fct)(a+c) + (*fct)(a-c) ); 00685 c = b*.36609105937014484; 00686 y += .32911111388180923e-1* ( (*fct)(a+c) + (*fct)(a-c) ); 00687 c = b*.3315221334651076; 00688 y += .36172897054424253e-1* ( (*fct)(a+c) + (*fct)(a-c) ); 00689 c = b*.29385787862038116; 00690 y += .39096947893535153e-1* ( (*fct)(a+c) + (*fct)(a-c) ); 00691 c = b*.2534499544661147; 00692 y += .41655962113473378e-1* ( (*fct)(a+c) + (*fct)(a-c) ); 00693 c = b*.21067563806531767; 00694 y += .43826046502201906e-1* ( (*fct)(a+c) + (*fct)(a-c) ); 00695 c = b*.16593430114106382; 00696 y += .45586939347881942e-1* ( (*fct)(a+c) + (*fct)(a-c) ); 00697 c = b*.11964368112606854; 00698 y += .46922199540402283e-1* ( (*fct)(a+c) + (*fct)(a-c) ); 00699 c = b*.7223598079139825e-1; 00700 y += .47819360039637430e-1* ( (*fct)(a+c) + (*fct)(a-c) ); 00701 c = b*.24153832843869158e-1; 00702 y += .4827004425736390e-1 * ( (*fct)(a+c) + (*fct)(a-c) ); 00703 y *= b; 00704 00705 /* the answer */ 00706 00707 return( y ); 00708 }