/*gffsub compute gaunt factors for any charge, drive Hummer routine, mod by J Ferguson */ #include "cddefines.h" #include "gffsub.h" #include "logte.h" void gffsub(double z, float xlf[], float g[], long int low, long int n, long int *iflag) { long int i, ir, j, _r; double b[11], c[8], con, gamma2, slope, txg, txu, u, xlrkt; static int _aini = 1; static double _es0[88]; double (*const d)[8] = (double(*)[8])_es0; double *const dd = (double*)_es0; if( _aini ){ /* Do 1 TIME INITIALIZATIONS! */ { static double _itmp0[] = {8.986940175e00,-4.009515855e00, 8.808871266e-01,2.640245111e-02,-4.580645915e-02,-3.568055702e-03, 2.827798067e-03,3.365860195e-04,-8.006936989e-01,9.466021705e-01, 9.043402532e-02,-9.608451450e-02,-1.885629865e-02,1.050313890e-02, 2.800889961e-03,-1.078209202e-03,-3.781305103e-01,1.102726332e-01, -1.543619180e-02,8.310561114e-03,2.179620525e-02,4.259726289e-03, -4.181588794e-03,-1.770208330e-03,1.877213132e-02,-1.004885705e-01, -5.483366378e-02,-4.520154409e-03,8.366530426e-03,3.700273930e-03, 6.889320423e-04,9.460313195e-05,7.300158392e-02,3.576785497e-03, -4.545307025e-03,-1.017965604e-02,-9.530211924e-03,-3.450186162e-03, 1.040482914e-03,1.407073544e-03,-1.744671550e-03,2.864013856e-02, 1.903394837e-02,7.091074494e-03}; for( i=1, _r = 0; i <= 44; i++ ) { dd[i-1] = _itmp0[_r++]; } } { static double _itmp1[] = {-9.668371391e-04,-2.999107465e-03, -1.820642230e-03,-3.874082085e-04,-1.707268366e-02,-4.694254776e-03, 1.311691517e-03,5.316703136e-03,5.178193095e-03,2.451228935e-03, -2.277321615e-05,-8.182359057e-04,2.567331664e-04,-9.155339970e-03, -6.997479192e-03,-3.571518641e-03,-2.096101038e-04,1.553822487e-03, 1.509584686e-03,6.212627837e-04,4.098322531e-03,1.635218463e-03, -5.918883504e-04,-2.333091048e-03,-2.484138313e-03,-1.359996060e-03, -5.371426147e-05,5.553549563e-04,3.837562402e-05,2.938325230e-03, 2.393747064e-03,1.328839809e-03,9.135013312e-05,-7.137252303e-04, -7.656848158e-04,-3.504683798e-04,-8.491991820e-04,-3.615327726e-04, 3.148015257e-04,8.909207650e-04,9.869737522e-04,6.134671184e-04, 1.068883394e-04,-2.046080100e-04}; for( i=45, _r = 0; i <= 88; i++ ) { dd[i-1] = _itmp1[_r++]; } } _aini = 0; } # ifdef DEBUG_FUN fputs( "<+>gffsub()\n", debug_fp ); # endif /*compute gaunt factors for any charge, drive Hummer routine, mod by J Ferguson * * GENERATES THERMALLY AVERAGED FREE-FREE NON-RELATIVISTIC GAUNT * FACTOR FOR A HYDROGENIC ION OF CHARGE Z, WITH A MAXIMUM RELATIVE * ERROR OF .OO7, (RMS FITTING ERROR = .001) FOR TEMPERATURE AND * FREQUENCY IN INTERVALS: * 10**-4 LE U LE 10**1.5, * 10**-3 LE GAMS LE 10**3, * WHERE U = H*NU/K*T AND GAMS = Z**2*RYD/K*T. TO OBTAIN THE * STATED ACCURACY, THE FULL NUMBER OF SIGNIFICANT FIGURES IN THE * COEFFICIENTS MUST BE RETAINED. * * THIS SUBROUTINE USES A TWO-DIMENSIONAL CHEBYSHEV EXPANSION * COMPUTED FROM EXPRESSIONS GIVEN BY KARZAS AND LATTER (AP.J. * SUPPL., V.6, P.167, 1961) AUGMENTED BY VARIOUS LIMITING FORMS * OF ENERGY-SPECIFIC GAUNT FACTORS. * D.G.HUMMER, JILA, MAY 1987. Ap.J. 327, 477 * modified with correct limits, J Ferguson, July 94 * * SUBROUTINE ARGUMENTS ARE: * Z = base 10 log of nuclear charge * => T = TEMPERATURE IN DEGREES KELVIN (now in common) * must set alogte in common logte * XLF = ARRAY CONTAINING INPUT VALUES OF LOG10(H*NU/RYD) * G = ARRAY CONTAINING OUTPUT VALUES OF GFF * LOW = FIRST ELEMENT OF ARRAY TO BE EVALUATED * N = LENGTH OF XLF AND G ARRAYS * NDIM = DIMENSION OF XLF AND G ARRAYS IN CALLING PROGRAM * IFLAG= explanation given in each area. * * */ /* COMPUTE TEMPERATURE-DEPENDENT COEFFICIENTS FOR U EXPANSION */ *iflag = 0; /* XLRXT IS LOG(RYD/KT), code not valid for TE> 10**8 * XLRXT IS LOG(RYD/KT), code not valid for TE< 10**2.5 */ if( logte.alogte < 2.5 ) { xlrkt = 5.1983649 - 2.5; *iflag = 1; } else if( logte.alogte > 8. ) { xlrkt = 5.1983649 - 8.; *iflag = 8; } else { xlrkt = 5.1983649 - logte.alogte; } /* Z is the base 10 log of the nuclear charge * TXG (a Hummer variable) is related to gamma^2 * TXU is related to U */ txg = 0.66666667*(2.0*z + xlrkt); gamma2 = pow((double)10,txg*1.5); con = 0.72727273*xlrkt + 0.90909091; for( j=0; j < 8; j++ ) { ir = 9; b[10] = d[10][j]; b[9] = txg*b[10] + d[9][j]; for( i=0; i < 9; i++ ) { b[ir-1] = txg*b[ir] - b[ir+1] + d[ir-1][j]; ir -= 1; } c[j] = 0.25*(b[0] - b[2]); } /* sum U expansion * loop through energy at fixed temperature */ for( i=low-1; i < n; i++ ) { /* TXU is related to U */ txu = 0.72727273*xlf[i] + con; u = pow((double)10,(txu - .90909091)/.72727273); /* criteria set by Hummer limits. It is a wedge from * (LOG(hnu),log(T)) = * (-5,2.5) to (-5,4) to (-1,8) to (4,8) to (-1.5,2.5). * These limits correspond to the gamma^2 and U limits given above */ if( fabs(txu) <= 2.0 ) { ir = 6; b[7] = c[7]; b[6] = txu*b[7] + c[6]; for( j=0; j < 6; j++ ) { b[ir-1] = txu*b[ir] - b[ir+1] + c[ir-1]; ir -= 1; } g[i] = (float)(b[0] - b[2]); if( logte.alogte >= 2.5 && logte.alogte <= 8.0 ) *iflag = 0; /* On the bottom side of the hummer box,u<-4 and gamma2>.33 */ } else if( log10(u) < -4.0 && gamma2 < 0.3 ) { g[i] = (float)(0.551329*log(2.24593/u)); if( logte.alogte >= 2.5 && logte.alogte <= 8.0 ) *iflag = 2; /* On the bottom side of the box,u<-4 and gamma2<.33 */ } else if( log10(u) < -4.0 && gamma2 >= 0.3 ) { g[i] = (float)(0.551329*log(0.944931/u/sqrt(gamma2))); if( logte.alogte >= 2.5 && logte.alogte <= 8.0 ) *iflag = 3; /* Now on the bottom side of the box. */ } else if( txu > 2.0 ) { /* Top of the box high T first */ if( log10(gamma2) < -3. ) { g[i] = (float)sqrt(0.9549297/u); if( logte.alogte >= 2.5 && logte.alogte <= 8.0 ) *iflag = 4; /* Must interpolate between two asymptotes */ } else if( log10(gamma2) < 0.0 ) { slope = (sqrt(12*gamma2/u) - sqrt(0.9549297/u))/3.0; g[i] = (float)(sqrt(12*gamma2/u) + slope*log10(gamma2)); if( logte.alogte >= 2.5 && logte.alogte <= 8.0 ) *iflag = 7; } else { /* The top side with wedge of 1.05 where region 6 fails */ g[i] = (float)(sqrt(12.*gamma2/u)); g[i] = (float)MIN2(1.05,g[i]); if( (g[i] == 1.05 && logte.alogte >= 2.5) && logte.alogte <= 8.0 ) *iflag = 5; if( (g[i] < 1.05 && logte.alogte >= 2.5) && logte.alogte <= 8.0 ) *iflag = 6; } } u = 0.0; } # ifdef DEBUG_FUN fputs( " <->gffsub()\n", debug_fp ); # endif return; }