/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and * others. For conditions of distribution and use see copyright notice in license.txt */ /*HCSAR_interp interpolate on collision strengths */ /*C6cs123 line collision rates for lower levels of hydrogenic carbon, n=1,2,3 */ /*Ca20cs123 line collision rates for lower levels of hydrogenic calcium, n=1,2,3 */ /*Hydcs123 Hydrogenic de-excitation collision rates n=1,2,3 */ /*H1cs123 hydrogen collision data levels involving 1s,2s,2p,3. */ /*Ne10cs123 line collision rates for lower levels of hydrogenic neon, n=1,2,3 */ /*He2cs123 line collision strengths for lower levels of helium ion, n=1,2,3, by K Korista */ /*Fe26cs123 line collision rates for lower levels of hydrogenic iron, n=1,2,3 */ #include "cddefines.h" #include "atmdat.h" #include "dense.h" #include "helike_cs.h" #include "hydro_vs_rates.h" #include "iso.h" #include "opacity.h" #include "phycon.h" #include "physconst.h" #include "taulines.h" STATIC double Fe26cs123(long int i, long int j); STATIC double He2cs123(long int i, long int j); STATIC double Hydcs123(long int ilow, long int ihigh, long int iz, long int chType); STATIC double C6cs123(long int i, long int j); STATIC double Ca20cs123(long int i, long int j); STATIC double Ne10cs123(long int i, long int j); STATIC realnum HCSAR_interp( int ipLo , int ipHi ); STATIC double CS_ThermAve_PR78(long ipISO, long nelem, long nHi, long nLo, double deltaE, double temp ); STATIC double Therm_ave_coll_str_int_PR78( double EOverKT ); STATIC double CS_PercivalRichards78( double Ebar ); static long global_ipISO, global_nelem, global_nHi, global_nLo; static double kTRyd, global_deltaE; static const realnum HCSTE[NHCSTE] = {5802.f,11604.f,34812.f,58020.f,116040.f,174060.f,232080.f,290100.f}; /*HCSAR_interp interpolate on collision strengths */ STATIC realnum HCSAR_interp( int ipLo , int ipHi ) { static int ip=1; realnum cs; DEBUG_ENTRY( "HCSAR_interp()" ); if( ipLo==1 && ipHi==2 ) { fprintf(ioQQQ,"HCSAR_interp was called for the 2s-2p transition, which it cannot do\n"); cdEXIT(EXIT_FAILURE); } if( phycon.te <= HCSTE[0] ) { cs = t_ADfA::Inst().h_coll_str( ipLo , ipHi , 0 ); } else if( phycon.te >= HCSTE[NHCSTE-1] ) { cs = t_ADfA::Inst().h_coll_str( ipLo , ipHi , NHCSTE-1 ); } else { /* the ip index is most likely correct since it points to the last temperature */ if( (HCSTE[ip-1] >= phycon.te ) || ( phycon.te > HCSTE[ip]) ) { /* we must find the temperature in the array */ for( ip=1; ip>refer h1 cs Callaway, J. 1983, Phys Let A, 96, 83 * >>refer h1 cs Zygelman, B., & Dalgarno, A. 1987, Phys Rev A, 35, 4085 * for 2p-2s only. * The fit from Callaway is in nuclear charge for 1s - 2s,2p only. * For transtions involving level 3, interpolation in Z involving * the functions He2cs123,C6cs123,Ne10cs123,Ca20cs123, Fe26cs123. * * The fits from ZD are for 2p-2s for Z=2,6,12,16,18 only other charges are * interpolated, both electron and proton rates are included, * the variable chType is either 'e' or 'p'. * * ipLow is the lower level and runs from 1 to 3 (1s, 2s, 2p) * ipHi is the upper level and runs from 2 to 6 (2s, 2p, 3s, 3p, 3d) */ /* for Callaway fit: */ /* for Zygelman and Dalgarno: */ /* first entry is 2p, then 2s */ /* fit in nuclear charge Z for 2p-2s collisions in hydrogenic species * equation is a+bx+cx^2ln(x)+dexp(x)+eln(x)/x^2, where x=te/Z^2 in a.u. * first are the proton rates: */ /* these are electron rates: */ /* following is for charged species */ /* charge is on scale with He+=1, Li++=2, etc */ ASSERT( nelem > ipHYDROGEN ); ASSERT( nelem < LIMELM ); /* these are the pointers to upper and lower levels. 1=1s, 2=2s, 3=2p, 4=3 */ ASSERT( ipLow > 0); ASSERT( ipLow <= 3); ASSERT( ipHi > 1 ); ASSERT( ipHi <=6 ); /* set quantum numbers and stat. weights of the transitions: */ if( ipHi == 6 ) { /* upper is n=3 then set level, stat. weight */ QuanNUp = 3.; gHi = 10.; /* following will be set here even though it is not used in this case, * to prevent good compilers from falsing on i not set, * there is assert when used to make sure it is ok */ i = -1; } else if( ipHi == 5 ) { /* upper is n=3 then set level, stat. weight */ QuanNUp = 3.; gHi = 6.; /* following will be set here even though it is not used in this case, * to prevent good compilers from falsing on i not set, * there is assert when used to make sure it is ok */ i = -1; } else if( ipHi == 4 ) { /* upper is n=3 then set level, stat. weight */ QuanNUp = 3.; gHi = 2.; /* following will be set here even though it is not used in this case, * to prevent good compilers from falsing on i not set, * there is assert when used to make sure it is ok */ i = -1; } else if( ipHi == 3 ) { /* upper is nl=2p then set level, stat. weight */ QuanNUp = 2.; gHi = 6.; /* used to point within vectors defined above */ i = 0; } else if( ipHi == 2 ) { /* upper is nl=2s then set level, stat. weight */ QuanNUp = 2.; gHi = 2.; /* used to point within vectors defined above */ i = 1; } else { fprintf( ioQQQ, " Insane levels in Hydcs123\n" ); cdEXIT(EXIT_FAILURE); } /* which lower level? */ if( ipLow == 1 ) { /* lower is n=1 then set level, stat. weight */ QuanNLo = 1.; gLo = 2.; } else if( ipLow == 2 ) { /* lower is nl=2s then set level, stat. weight */ QuanNLo = 2.; gLo = 2.; } else if( ipLow == 3 ) { QuanNLo = 2.; gLo = 6.; } else { fprintf( ioQQQ, " Insane levels in Hydcs123\n" ); cdEXIT(EXIT_FAILURE); } /* following is the physical charge */ Charge = (double)(nelem + 1); /* square of charge */ ChargeSquared = Charge*Charge; if( ipLow == 2 && ipHi == 3 ) { /*************************************** this is 2p-2s: * series of if statements determines which entries Charge is between. */ if( nelem < 1 ) { /* this can't happen since routine returned above when ip=1 and * special atomic hydrogen routine called */ fprintf( ioQQQ, " insane charge given to Hydcs123\n" ); cdEXIT(EXIT_FAILURE); } else if( nelem == 1 ) { zlow = 2.; j = 1; zhigh = 2.; k = 1; } /* Li through C */ else if( nelem <= 5 ) { zlow = 2.; j = 1; zhigh = 6.; k = 2; } else if( nelem <= 11 ) { zlow = 6.; j = 2; zhigh = 12.; k = 3; } else if( nelem <= 15 ) { zlow = 12.; j = 3; zhigh = 16.; k = 4; } else if( nelem <= 17 ) { zlow = 16.; j = 4; zhigh = 18.; k = 5; } /* following changed to else from else if, * to prevent false comment in good compilers */ /*else if( nelem > 18 )*/ else { zlow = 18.; j = 5; zhigh = 18.; k = 5; } /* convert Te to a.u./Z^2 * determine rate at the low Charge */ x = EVRYD/TE1RYD*phycon.te/(27.211396*pow2(zlow)); TeUse = MIN2(x,0.80); x = MAX2(0.025,TeUse); /* what type of collision are we dealing with? */ if( chType == 'e' ) { /* electron collisions */ Ratelow = ae[j-1] + be[j-1]*x + ce[j-1]*pow2(x)*log(x) + de[j-1]* exp(x) + ee[j-1]*log(x)/pow2(x); } else if( chType == 'p' ) { Ratelow = ap[j-1] + bp[j-1]*x + cp[j-1]*pow2(x)*log(x) + dp[j-1]* exp(x) + ep[j-1]*log(x)/pow2(x); } else { /* this can't happen */ fprintf( ioQQQ, " insane collision species given to Hydcs123\n" ); cdEXIT(EXIT_FAILURE); } /* determine rate at the high Charge */ x = EVRYD/TE1RYD*phycon.te/(27.211396*pow2(zhigh)); TeUse = MIN2(x,0.80); x = MAX2(0.025,TeUse); if( chType == 'e' ) { Ratehigh = ae[k-1] + be[k-1]*x + ce[k-1]*pow2(x)*log(x) + de[k-1]*exp(x) + ee[k-1]*log(x)/pow2(x); } else { Ratehigh = ap[k-1] + bp[k-1]*x + cp[k-1]*pow2(x)*log(x) + dp[k-1]*exp(x) + ep[k-1]*log(x)/pow2(x); } /* linearly interpolate in charge */ if( fp_equal( zlow, zhigh ) ) { rate = Ratelow; } else { slope = (Ratehigh - Ratelow)/(zhigh - zlow); rate = slope*(Charge - zlow) + Ratelow; } rate = rate/ChargeSquared/Charge*1.0e-7; /* must convert to cs and need to know the valid temp range */ templow = 0.025*27.211396*TE1RYD/EVRYD*ChargeSquared; temphigh = 0.80*27.211396*TE1RYD/EVRYD*ChargeSquared; TeUse = MIN2((double)phycon.te,temphigh); temp = MAX2(TeUse,templow); Hydcs123_v = rate*gHi*sqrt(temp)/COLL_CONST; if( chType == 'p' ) { /* COLL_CONST is incorrect for protons, correct here */ Hydcs123_v *= pow( PROTON_MASS/ELECTRON_MASS, 1.5 ); } } else if( ipHi == 4 || ipHi == 5 || ipHi == 6 ) { /* n = 3 * for the rates involving n = 3, must do something different. */ if( nelem < 1 ) { fprintf( ioQQQ, " insane charge given to Hydcs123\n" ); cdEXIT(EXIT_FAILURE); } else if( nelem == 1 ) { zlow = 2.; Ratelow = He2cs123(ipLow,ipHi); zhigh = 2.; Ratehigh = Ratelow; } else if( nelem > 1 && nelem <= 5 ) { zlow = 2.; Ratelow = He2cs123(ipLow,ipHi); zhigh = 6.; Ratehigh = C6cs123(ipLow,ipHi); } else if( nelem > 5 && nelem <= 9 ) { zlow = 6.; Ratelow = C6cs123(ipLow,ipHi); zhigh = 10.; Ratehigh = Ne10cs123(ipLow,ipHi); } else if( nelem > 9 && nelem <= 19 ) { zlow = 10.; Ratelow = Ne10cs123(ipLow,ipHi); zhigh = 20.; Ratehigh = Ca20cs123(ipLow,ipHi); } else if( nelem > 19 && nelem <= 25 ) { zlow = 20.; Ratelow = Ca20cs123(ipLow,ipHi); zhigh = 26.; Ratehigh = Fe26cs123(ipLow,ipHi); } /*>>>chng 98 dec 17, to else to stop comment from good compilers*/ /*else if( nelem > 26 )*/ else { Charge = 26.; zlow = 26.; Ratelow = Fe26cs123(ipLow,ipHi); zhigh = 26.; Ratehigh = Ratelow; } /* linearly interpolate */ if( fp_equal( zlow, zhigh ) ) { rate = Ratelow; } else { slope = (Ratehigh - Ratelow)/(zhigh - zlow); rate = slope*(Charge - zlow) + Ratelow; } /** NB - all of these actually calculate EIE collision strengths */ Hydcs123_v = rate; /* the routines C6cs123, Ne10cs123, etc... do not resolve L for n>2 */ /* dividing by N should roughly recover the original l-resolved data */ Hydcs123_v /= 3.; } else { /* this branch 1-2s, 1-2p */ if( nelem == 1 ) { /* this brance for helium, then return */ Hydcs123_v = He2cs123(ipLow,ipHi); return( Hydcs123_v ); } /* electron temperature in eV */ tev = phycon.te / EVDEGK; /* energy in eV for hydrogenic species and these quantum numbers */ EE = ChargeSquared*EVRYD*(1./QuanNLo/QuanNLo - 1./QuanNUp/QuanNUp); /* EE/kT for this transion */ q = EE/tev; TeUse = MIN2(q,10.); /* q is now EE/kT but between 1 and 10 */ q = MAX2(1.,TeUse); expq = exp(q); /* i must be 0 or 1 */ ASSERT( i==0 || i==1 ); C = C0[i] + C1[i]/Charge + C2[i]/ChargeSquared; D = D0[i] + D1[i]/Charge + D2[i]/ChargeSquared; /* following code changed so that ee1 always returns e1, * orifinal version only returned e1 for x < 1 */ /* use disabled e1: */ /*if( q < 1. )*/ /*{*/ /*rate = (B[i-1] + D*q)*exp(-q) + (A[i-1] + C*q - D*q*q)**/ /*ee1(q);*/ /*}*/ /*else*/ /*{*/ /*rate = (B[i-1] + D*q) + (A[i-1] + C*q - D*q*q)*ee1(q);*/ /*}*/ /*rate *= 8.010e-8/2./ChargeSquared/tev*sqrt(tev);*/ /* convert to de-excitation */ /*if( q < 1. )*/ /*{*/ /*rate = rate*exp(q)*gLo/gHi;*/ /*}*/ /*else*/ /*{*/ /*rate = rate*gLo/gHi;*/ /*}*/ /*>>>chng 98 dec 17, ee1 always returns e1 */ rate = (B[i] + D*q)/expq + (A[i] + C*q - D*q*q)* ee1(q); rate *= 8.010e-8/2./ChargeSquared/tev*sqrt(tev); /* convert to de-excitation */ rate *= expq*gLo/gHi; /* convert to cs */ Hydcs123_v = rate*gHi*phycon.sqrte/COLL_CONST; } return( Hydcs123_v ); } /*C6cs123 line collision rates for lower levels of hydrogenic carbon, n=1,2,3 */ STATIC double C6cs123(long int i, long int j) { double C6cs123_v, TeUse, t, x; static const double a[3] = {-92.23774,-1631.3878,-6326.4947}; static const double b[3] = {-11.93818,-218.3341,-849.8927}; static const double c[3] = {0.07762914,1.50127,5.847452}; static const double d[3] = {78.401154,1404.8475,5457.9291}; static const double e[3] = {332.9531,5887.4263,22815.211}; DEBUG_ENTRY( "C6cs123()" ); /* These are fits to Table 5 of * >>refer c6 cs Aggarwal, K.M., & Kingston, A.E. 1991, J Phys B, 24, 4583 * C VI collision rates for 1s-3l, 2s-3l, and 2p-3l, * principal quantum numbers n and l. * * i is the lower level and runs from 1 to 3 (1s, 2s, 2p) * j is the upper level and runs from 2 to 6 (2s, 2p, 3s, 3p, 3d) * 1s-2s,2p is not done here. * check temperature: fits only good between 3.8 < log Te < 6.2 */ /* arrays for fits of 3 transitions see the code below for key: */ TeUse = MAX2(phycon.te,6310.); t = MIN2(TeUse,1.6e6); x = log10(t); if( i == 1 && j == 2 ) { /* 1s - 2s (first entry) */ fprintf( ioQQQ, " Carbon VI 2s-1s not done in C6cs123\n" ); cdEXIT(EXIT_FAILURE); } else if( i == 1 && j == 3 ) { /* 1s - 2p (second entry) */ fprintf( ioQQQ, " Carbon VI 2p-1s not done in C6cs123\n" ); cdEXIT(EXIT_FAILURE); } else if( i == 1 && ( j == 4 || j == 5 || j == 6 ) ) { /* 1s - 3 (first entry) */ C6cs123_v = a[0] + b[0]*x + c[0]*pow2(x)*sqrt(x) + d[0]*log(x) + e[0]*log(x)/pow2(x); } else if( i == 2 && ( j == 4 || j == 5 || j == 6 ) ) { /* 2s - 3 (second entry) */ C6cs123_v = a[1] + b[1]*x + c[1]*pow2(x)*sqrt(x) + d[1]*log(x) + e[1]*log(x)/pow2(x); } else if( i == 3 && ( j == 4 || j == 5 || j == 6 ) ) { /* 2p - 3s (third entry) */ C6cs123_v = a[2] + b[2]*x + c[2]*pow2(x)*sqrt(x) + d[2]*log(x) + e[2]*log(x)/pow2(x); } else { fprintf( ioQQQ, " insane levels for C VI n=1,2,3 !!!\n" ); cdEXIT(EXIT_FAILURE); } return( C6cs123_v ); } /*Ca20cs123 line collision rates for lower levels of hydrogenic calcium, n=1,2,3 */ STATIC double Ca20cs123(long int i, long int j) { double Ca20cs123_v, TeUse, t, x; static const double a[3] = {-12.5007,-187.2303,-880.18896}; static const double b[3] = {-1.438749,-22.17467,-103.1259}; static const double c[3] = {0.008219688,0.1318711,0.6043752}; static const double d[3] = {10.116516,153.2650,717.4036}; static const double e[3] = {45.905343,685.7049,3227.2836}; DEBUG_ENTRY( "Ca20cs123()" ); /* * These are fits to Table 5 of * >>refer ca20 cs Aggarwal, K.M., & Kingston, A.E. 1992, J Phys B, 25, 751 * Ca XX collision rates for 1s-3l, 2s-3l, and 2p-3l, * principal quantum numbers n and l. * * i is the lower level and runs from 1 to 3 (1s, 2s, 2p) * j is the upper level and runs from 2 to 6 (2s, 2p, 3s, 3p, 3d) * 1s-2s,2p is not done here. * check temperature: fits only good between 5.0 < log Te < 7.2 */ /* arrays for fits of 3 transitions see the code below for key: */ TeUse = MAX2(phycon.te,1.0e5); t = MIN2(TeUse,1.585e7); x = log10(t); if( i == 1 && j == 2 ) { /* 1s - 2s (first entry) */ fprintf( ioQQQ, " Ca XX 2s-1s not done in Ca20cs123\n" ); cdEXIT(EXIT_FAILURE); } else if( i == 1 && j == 3 ) { /* 1s - 2p (second entry) */ fprintf( ioQQQ, " Ca XX 2p-1s not done in Ca20cs123\n" ); cdEXIT(EXIT_FAILURE); } else if( i == 1 && ( j == 4 || j == 5 || j == 6 )) { /* 1s - 3 (first entry) */ Ca20cs123_v = a[0] + b[0]*x + c[0]*pow2(x)*sqrt(x) + d[0]*log(x) + e[0]*log(x)/pow2(x); } else if( i == 2 && ( j == 4 || j == 5 || j == 6 )) { /* 2s - 3 (second entry) */ Ca20cs123_v = a[1] + b[1]*x + c[1]*pow2(x)*sqrt(x) + d[1]*log(x) + e[1]*log(x)/pow2(x); } else if( i == 3 && ( j == 4 || j == 5 || j == 6 )) { /* 2p - 3s (third entry) */ Ca20cs123_v = a[2] + b[2]*x + c[2]*pow2(x)*sqrt(x) + d[2]*log(x) + e[2]*log(x)/pow2(x); } else { fprintf( ioQQQ, " insane levels for Ca XX n=1,2,3 !!!\n" ); cdEXIT(EXIT_FAILURE); } return( Ca20cs123_v ); } /*Ne10cs123 line collision rates for lower levels of hydrogenic neon, n=1,2,3 */ STATIC double Ne10cs123(long int i, long int j) { double Ne10cs123_v, TeUse, t, x; static const double a[3] = {3.346644,151.2435,71.7095}; static const double b[3] = {0.5176036,20.05133,13.1543}; static const double c[3] = {-0.00408072,-0.1311591,-0.1099238}; static const double d[3] = {-3.064742,-129.8303,-71.0617}; static const double e[3] = {-11.87587,-541.8599,-241.2520}; DEBUG_ENTRY( "Ne10cs123()" ); /* These are fits to Table 5 of * >>refer ne10 cs Aggarwal, K.M., & Kingston, A.E. 1991, PhyS, 44, 517 * Ne X collision rates for 1s-3, 2s-3l, and 2p-3l, * principal quantum numbers n and l. * * i is the lower level and runs from 1 to 3 (1s, 2s, 2p) * j is the upper level and runs from 2 to 6 (2s, 2p, 3s, 3p, 3d) * 1s-2s,2p is not done here. * check temperature: fits only good between 3.8 < log Te < 6.2 * */ /* arrays for fits of 3 transitions see the code below for key: */ TeUse = MAX2(phycon.te,6310.); t = MIN2(TeUse,1.6e6); x = log10(t); if( i == 1 && j == 2 ) { /* 1s - 2s (first entry) */ fprintf( ioQQQ, " Neon X 2s-1s not done in Ne10cs123\n" ); cdEXIT(EXIT_FAILURE); } else if( i == 1 && j == 3 ) { /* 1s - 2p (second entry) */ fprintf( ioQQQ, " Neon X 2p-1s not done in Ne10cs123\n" ); cdEXIT(EXIT_FAILURE); } else if( i == 1 && ( j == 4 || j == 5 || j == 6 ) ) { /* 1s - 3 (first entry) */ Ne10cs123_v = a[0] + b[0]*x + c[0]*pow2(x)*sqrt(x) + d[0]*log(x) + e[0]*log(x)/pow2(x); } else if( i == 2 && ( j == 4 || j == 5 || j == 6 ) ) { /* 2s - 3 (second entry) */ Ne10cs123_v = a[1] + b[1]*x + c[1]*pow2(x)*sqrt(x) + d[1]*log(x) + e[1]*log(x)/pow2(x); } else if( i == 3 && ( j == 4 || j == 5 || j == 6 ) ) { /* 2p - 3s (third entry) */ Ne10cs123_v = a[2] + b[2]*x + c[2]*pow2(x)*sqrt(x) + d[2]*log(x) + e[2]*log(x)/pow2(x); } else { fprintf( ioQQQ, " insane levels for Ne X n=1,2,3 !!!\n" ); cdEXIT(EXIT_FAILURE); } return( Ne10cs123_v ); } /*He2cs123 line collision strengths for lower levels of helium ion, n=1,2,3, by K Korista */ STATIC double He2cs123(long int i, long int j) { double He2cs123_v, t; static const double a[11]={0.12176209,0.32916723,0.46546497,0.044501688, 0.040523277,0.5234889,1.4903214,1.4215094,1.0295881,4.769306,9.7226127}; static const double b[11]={0.039936166,2.9711166e-05,-0.020835863,3.0508137e-04, -2.004485e-15,4.41475e-06,1.0622666e-05,2.0538877e-06,0.80638448,2.0967075e-06, 7.6089851e-05}; static const double c[11]={143284.77,0.73158545,-2.159172,0.43254802,2.1338557, 8.9899702e-06,-2.9001451e-12,1.762076e-05,52741.735,-2153.1219,-3.3996921e-11}; DEBUG_ENTRY( "He2cs123()" ); /* These are fits to Table 2 * >>refer he2 cs Aggarwal, K.M., Callaway, J., Kingston, A.E., Unnikrishnan, K. * >>refercon 1992, ApJS, 80, 473 * He II collision rates for 1s-2s, 1s-2p, 1s-3s, 1s-3p, 1s-3d, 2s-3s, 2s-3p, 2s-3d, * 2p-3s, 2p-3p, and 2p-3d. * principal quantum numbers n and l. * * i is the lower level and runs from 1 to 3 (1s, 2s, 2p) * j is the upper level and runs from 2 to 6 (2s, 2p, 3s, 3p, 3d) * check temperature: fits only good between 5,000K and 500,000K * */ /* array for fits of 11 transitions see the code below for key: */ t = phycon.te; if( t < 5000. ) { t = 5000.; } else if( t > 5.0e05 ) { t = 5.0e05; } /**************fits begin here************** * */ if( i == 1 && j == 2 ) { /* 1s - 2s (first entry) */ He2cs123_v = a[0] + b[0]*exp(-t/c[0]); } else if( i == 1 && j == 3 ) { /* 1s - 2p (second entry) */ He2cs123_v = a[1] + b[1]*pow(t,c[1]); } else if( i == 1 && j == 4 ) { /* 1s - 3s (third entry) */ He2cs123_v = a[2] + b[2]*log(t) + c[2]/log(t); } else if( i == 1 && j == 5 ) { /* 1s - 3p (fourth entry) */ He2cs123_v = a[3] + b[3]*pow(t,c[3]); } else if( i == 1 && j == 6 ) { /* 1s - 3d (fifth entry) */ He2cs123_v = a[4] + b[4]*pow(t,c[4]); } else if( i == 2 && j == 4 ) { /* 2s - 3s (sixth entry) */ He2cs123_v = (a[5] + c[5]*t)/(1 + b[5]*t); } else if( i == 2 && j == 5 ) { /* 2s - 3p (seventh entry) */ He2cs123_v = a[6] + b[6]*t + c[6]*t*t; } else if( i == 2 && j == 6 ) { /* 2s - 3d (eighth entry) */ He2cs123_v = (a[7] + c[7]*t)/(1 + b[7]*t); } else if( i == 3 && j == 4 ) { /* 2p - 3s (ninth entry) */ He2cs123_v = a[8] + b[8]*exp(-t/c[8]); } else if( i == 3 && j == 5 ) { /* 2p - 3p (tenth entry) */ He2cs123_v = a[9] + b[9]*t + c[9]/t; } else if( i == 3 && j == 6 ) { /* 2p - 3d (eleventh entry) */ He2cs123_v = a[10] + b[10]*t + c[10]*t*t; } else { /**************fits end here************** */ fprintf( ioQQQ, " insane levels for He II n=1,2,3 !!!\n" ); cdEXIT(EXIT_FAILURE); } return( He2cs123_v ); } /*Fe26cs123 line collision rates for lower levels of hydrogenic iron, n=1,2,3 */ STATIC double Fe26cs123(long int i, long int j) { double Fe26cs123_v, TeUse, t, x; static const double a[3] = {-4.238398,-238.2599,-1211.5237}; static const double b[3] = {-0.4448177,-27.06869,-136.7659}; static const double c[3] = {0.0022861,0.153273,0.7677703}; static const double d[3] = {3.303775,191.7165,972.3731}; static const double e[3] = {15.82689,878.1333,4468.696}; DEBUG_ENTRY( "Fe26cs123()" ); /* These are fits to Table 5 of * >>refer fe26 cs Aggarwal, K.M., & Kingston, A.E. 1993, ApJS, 85, 187 * Fe XXVI collision rates for 1s-3, 2s-3, and 2p-3, * principal quantum numbers n and l. * * i is the lower level and runs from 1 to 3 (1s, 2s, 2p) * j is the upper level and runs from 2 to 6 (2s, 2p, 3s, 3p, 3d) * 1s-2s,2p is not done here. * check temperature: fits only good between 5.2 < log Te < 7.2 * */ /* arrays for fits of 3 transitions see the code below for key: */ TeUse = MAX2(phycon.te,1.585e5); t = MIN2(TeUse,1.585e7); x = log10(t); if( i == 1 && j == 2 ) { /* 1s - 2s (first entry) */ fprintf( ioQQQ, " Fe XXVI 2s-1s not done in Fe26cs123\n" ); cdEXIT(EXIT_FAILURE); } else if( i == 1 && j == 3 ) { /* 1s - 2p (second entry) */ fprintf( ioQQQ, " Fe XXVI 2p-1s not done in Fe26cs123\n" ); cdEXIT(EXIT_FAILURE); } else if( i == 1 && ( j == 4 || j == 5 || j == 6 ) ) { /* 1s - 3 (first entry) */ Fe26cs123_v = a[0] + b[0]*x + c[0]*pow2(x)*sqrt(x) + d[0]*log(x) + e[0]*log(x)/pow2(x); } else if( i == 2 && ( j == 4 || j == 5 || j == 6 ) ) { /* 2s - 3 (second entry) */ Fe26cs123_v = a[1] + b[1]*x + c[1]*pow2(x)*sqrt(x) + d[1]*log(x) + e[1]*log(x)/pow2(x); } else if( i == 3 && ( j == 4 || j == 5 || j == 6 ) ) { /* 2p - 3s (third entry) */ Fe26cs123_v = a[2] + b[2]*x + c[2]*pow2(x)*sqrt(x) + d[2]*log(x) + e[2]*log(x)/pow2(x); } else { fprintf( ioQQQ, " insane levels for Ca XX n=1,2,3 !!!\n" ); cdEXIT(EXIT_FAILURE); } return( Fe26cs123_v ); } STATIC double CS_ThermAve_PR78(long ipISO, long nelem, long nHi, long nLo, double deltaE, double temp ) { double coll_str; DEBUG_ENTRY( "CS_ThermAve_PR78()" ); global_ipISO = ipISO; global_nelem = nelem; global_nHi = nHi; global_nLo = nLo; global_deltaE = deltaE; kTRyd = temp / TE1RYD; if( !iso_ctrl.lgCS_therm_ave[ipISO] ) { /* Must do some thermal averaging for densities greater * than about 10000 and less than about 1e10, * because kT gives significantly different results. * Still, do sparser integration than is done below */ if( (dense.eden > 10000.) && (dense.eden < 1E10 ) ) { coll_str = qg32( 0.0, 6.0, Therm_ave_coll_str_int_PR78); } else { /* Do NOT average over Maxwellian */ coll_str = CS_PercivalRichards78( kTRyd ); } } else { /* DO average over Maxwellian */ coll_str = qg32( 0.0, 1.0, Therm_ave_coll_str_int_PR78); coll_str += qg32( 1.0, 10.0, Therm_ave_coll_str_int_PR78); } return coll_str; } /* The integrand for calculating the thermal average of collision strengths */ STATIC double Therm_ave_coll_str_int_PR78( double EOverKT ) { double integrand; DEBUG_ENTRY( "Therm_ave_coll_str_int_PR78()" ); integrand = exp( -1.*EOverKT ) * CS_PercivalRichards78( EOverKT * kTRyd ); return integrand; } STATIC double CS_PercivalRichards78( double Ebar ) { double cross_section, coll_str; double stat_weight; double A, D, L, F, G, H; double np, n, s, Z_3, xPlus, xMinus, y; long ipISO, nelem; DEBUG_ENTRY( "CS_PercivalRichards78()" ); if( Ebar < global_deltaE ) { DEBUG_ENTRY( "CS_PercivalRichards78()" ); return 0.; } ipISO = global_ipISO; nelem = global_nelem; np = (double)global_nHi; n = (double)global_nLo; s = np - n; ASSERT( s > 0. ); A = (8./3./s) * pow(np/s/n, 3.) * (0.184 - 0.04 * pow( s, -0.66)) * pow( 1. - 0.2*s/n/np, 1.+2.*s); Z_3 = (double)(nelem + 1. - ipISO); D = exp( - Z_3 * Z_3 / n / np / Ebar / Ebar ); L = log( (1. + 0.53 * Ebar * Ebar * n * np / Z_3 / Z_3) / (1. + 0.4*Ebar) ); F = pow( 1. - 0.3 * s * D / n /np, 1. + 2.*s ); G = 0.5* POW3( Ebar * n * n / Z_3 / np ); xPlus = 2. * Z_3 / ( n * n * Ebar * ( sqrt( 2. - n*n/np/np ) + 1. ) ); xMinus = 2. * Z_3 / ( n * n * Ebar * ( sqrt( 2. - n*n/np/np ) - 1. ) ); y = 1. / (1. - D * log ( 18. * s )/ 4. / s); H = POW2( xMinus) * log( 1. + 2.*xMinus/3. ) / ( 2.*y + 1.5*xMinus ); H -= POW2( xPlus ) * log( 1. + 2.* xPlus/3. ) / ( 2.*y + 1.5*xPlus ); /* this is the LHS of equation 1 of PR78 */ cross_section = (A*D*L + F*G*H); /* this is the result after solving equation 1 for the cross section */ cross_section *= PI * POW2( n * n * BOHR_RADIUS_CM / Z_3 ) / Ebar; if( ipISO == ipH_LIKE ) stat_weight = 2. * n * n; else if( ipISO == ipHE_LIKE ) stat_weight = 4. * n * n; else TotalInsanity(); /* convert to collision strength */ coll_str = cross_section * stat_weight * Ebar / ( PI * POW2( BOHR_RADIUS_CM ) ); return coll_str; } #if 0 STATIC void TestPercivalRichards( void ) { double CStemp; /* this reproduces Table 1 of PR78 */ for( long i=0; i<5; i++ ) { double Ebar[5] = {0.1, 0.4, 0.8, 1.0, 10.}; CStemp = CS_PercivalRichards78( 0, 2, 12, 10, Ebar[i] ); } /* this reproduces Table 2 of PR78 */ for( long i=0; i<5; i++ ) { double Ebar[5] = {0.1, 0.4, 0.8, 1.0, 10.}; CStemp = CS_ThermAve_PR78( ipISO, 0, N_(ipHi), N_(ipLo), phycon.te ); } return; } #endif realnum HydroCSInterp(long int nelem, long int ipHi, long int ipLo, long int ipCollider ) { double CStemp; long ipISO = ipH_LIKE; DEBUG_ENTRY( "HydroCSInterp()" ); if( !iso_ctrl.lgColl_excite[ipISO] ) return 0.; /* This set of collision strengths should only be used * if the Storey and Hummer flag is set */ if( opac.lgCaseB_HummerStorey ) { if( N_(ipLo) == N_(ipHi) ) { if( N_(ipHi) <= iso_sp[ipH_LIKE][nelem].n_HighestResolved_max && abs( L_(ipLo) - L_(ipHi) ) != 1 ) { /* if delta L is not +/- 1, set collision strength to zero. */ CStemp = 0.; } else { CStemp = CS_l_mixing_PS64( nelem, iso_sp[ipH_LIKE][nelem].st[ipLo].lifetime(), nelem+1.-ipH_LIKE, iso_sp[ipH_LIKE][nelem].st[ipLo].n(), iso_sp[ipH_LIKE][nelem].st[ipLo].l(), iso_sp[ipH_LIKE][nelem].st[ipHi].g(), ipCollider); } } else { if( N_(ipHi) <= iso_sp[ipH_LIKE][nelem].n_HighestResolved_max && abs( L_(ipLo) - L_(ipHi) ) != 1 ) { /* if delta L is not +/- 1, set collision strength to zero. */ CStemp = 0.; } else if( ipCollider == ipELECTRON ) { CStemp = CS_ThermAve_PR78( ipH_LIKE, nelem, N_(ipHi), N_(ipLo), iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).EnergyErg() / EN1RYD , phycon.te ); } else CStemp = 0.; } } else { /* HCSAR_interp interpolates on a table to return R-matrix collision strengths * for hydrogen only */ if( nelem==ipHYDROGEN && ipCollider==ipELECTRON && N_(ipHi) <= 5 && ( N_(ipHi) != N_(ipLo) ) ) { CStemp = HCSAR_interp(ipLo,ipHi); } else if( nelem==ipHYDROGEN && ipCollider==ipELECTRON && ( N_(ipHi) != N_(ipLo) ) ) { CStemp = hydro_vs_deexcit( ipH_LIKE, nelem, ipHi, ipLo, iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Aul() ); } else if( nelem>ipHYDROGEN && ipCollider==ipELECTRON && N_(ipHi) <= 3 && N_(ipLo) < 3 ) { CStemp = Hydcs123(ipLo + 1,ipHi + 1,nelem,'e'); } else if( nelem>ipHYDROGEN && ipCollider==ipPROTON && ipHi==ipH2p && ipLo==ipH2s ) { CStemp = Hydcs123(ipLo + 1,ipHi + 1,nelem,'p'); } else if( N_(ipLo) == N_(ipHi) ) { if( iso_ctrl.lgCS_Vrinceanu[ipH_LIKE] ) { CStemp = CS_l_mixing_VF01(ipH_LIKE, nelem, iso_sp[ipH_LIKE][nelem].st[ipLo].n(), iso_sp[ipH_LIKE][nelem].st[ipLo].l(), iso_sp[ipH_LIKE][nelem].st[ipHi].l(), iso_sp[ipH_LIKE][nelem].st[ipLo].S(), phycon.te, ipCollider ); } else CStemp = CS_l_mixing_PS64( nelem, iso_sp[ipH_LIKE][nelem].st[ipLo].lifetime(), nelem+1.-ipH_LIKE, iso_sp[ipH_LIKE][nelem].st[ipLo].n(), iso_sp[ipH_LIKE][nelem].st[ipLo].l(), iso_sp[ipH_LIKE][nelem].st[ipHi].g(), ipCollider); } else { ASSERT( N_(ipHi) != N_(ipLo) ); /* highly excited levels */ CStemp = CS_VS80( ipH_LIKE, nelem, ipHi, ipLo, iso_sp[ipH_LIKE][nelem].trans(ipHi,ipLo).Emis().Aul(), phycon.te, ipCollider ); } } return (realnum)CStemp; }