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 /*ChargTranEval fill in the HCharExcIonOf and Rec arrays with Kingdon's fitted CT with H */ 00004 /*ChargTranSumHeat sum net heating due to charge transfer, called in HeatSum */ 00005 /*HCTIon H charge transfer ionization*/ 00006 /*HCTRecom H charge transfer recombination */ 00007 /*ChargTranPun, punch charge transfer coefficient */ 00008 /*MakeHCTData holds data for charge transfer fits */ 00009 #include "cddefines.h" 00010 #include "phycon.h" 00011 #include "physconst.h" 00012 #include "abund.h" 00013 #include "dense.h" 00014 #include "iso.h" 00015 #include "thermal.h" 00016 #include "mole.h" 00017 #include "elementnames.h" 00018 #include "heavy.h" 00019 #include "trace.h" 00020 #include "conv.h" 00021 #include "atmdat.h" 00022 00023 /*HCTion H charge transfer ionization, H+ + A => H + A+ */ 00024 STATIC double HCTIon(long int ion, 00025 long int nelem); 00026 00027 /*HCTRecom H charge transfer recombination, H + A+ => H+ + A */ 00028 STATIC double HCTRecom(long int ion, 00029 long int nelem); 00030 00031 /* the translated block data */ 00032 STATIC void MakeHCTData(void); 00033 00034 /* the structures for storing the charge transfer data, these are filled in 00035 * at the end of this file, in what used to be a block data */ 00036 static double CTIonData[LIMELM][4][8]; 00037 static double CTRecombData[LIMELM][4][7]; 00038 /* this will be flag for whether or not charge transfer data 00039 * have been initialized */ 00040 static bool lgCTDataDefined = false; 00041 00042 /*ChargTranEval update charge trans rate coefficients if temperature has changed */ 00043 void ChargTranEval( 00044 /* return value is H ionization rate (s-1) due to O+ charge transfer */ 00045 double *O_HIonRate ) 00046 { 00047 long int ion, 00048 nelem; 00049 double a, b, c, a_op, b_op, c_op, d_op, e_op, f_op, a_o, 00050 b_o, c_o, d_o, e_o, f_o, g_o; 00051 00052 static double TeUsed = -1.; 00053 00054 DEBUG_ENTRY( "ChargTranEval()" ); 00055 00056 /* first is to force reevaluation on very first call */ 00057 if( !conv.nTotalIoniz || fabs(phycon.te-TeUsed)/phycon.te > 0.01 ) 00058 { 00059 /* refresh hydrogen charge transfer arrays */ 00060 /* >>chng 01 apr 25, lower limit had been 2, lithium, changed to 1, helium */ 00061 for( nelem=ipHELIUM; nelem < LIMELM; nelem++ ) 00062 { 00063 for( ion=0; ion <= nelem; ion++ ) 00064 { 00065 /* >>chng 01 apr 28, add factors to turn off ct, 00066 * had previously been handled downstream */ 00067 00068 /* HCharExcIonOf[nelem][ion]*hii is ion => ion+1 for nelem */ 00069 /* charge transfer ionization O^0 + H+ -> O+ + H0 00070 * is HCharExcIonOf[ipOXYGEN][0]*dense.xIonDense[ipHYDROGEN][1] 00071 * charge transfer recombination of atomic hydrogen is 00072 * HCharExcIonOf[ipOXYGEN][0]*dense.xIonDense[ipOXYGEN][0] */ 00073 atmdat.HCharExcIonOf[nelem][ion] = HCTIon(ion,nelem); 00074 00075 /* HCharExcRecTo[nelem][ion]*hi is ion+1 => ion of nelem */ 00076 /* charge transfer recombination O+ + H0 -> O^0 + H+ is 00077 * HCharExcRecTo[ipOXYGEN][0]*dense.xIonDense[ipHYDROGEN][0] 00078 * charge transfer ionization of atomic hydrogen is 00079 * HCharExcRecTo[ipOXYGEN][0]*dense.xIonDense[ipOXYGEN][1] */ 00080 atmdat.HCharExcRecTo[nelem][ion] = HCTRecom(ion,nelem); 00081 } 00082 00083 /* zero out helium charge transfer arrays */ 00084 for( ion=0; ion < LIMELM; ion++ ) 00085 { 00086 atmdat.HeCharExcIonOf[nelem][ion] = 0.; 00087 atmdat.HeCharExcRecTo[nelem][ion] = 0.; 00088 } 00089 } 00090 00091 /* The above included only the radiative charge transfer from 00092 * Stancil et al 1998. must explicitly add on the ct fitted by Kingdon & Ferland, 00093 * The process H0 + He+ -> He0 + H+ */ 00094 if( phycon.te > 6000. ) 00095 atmdat.HCharExcRecTo[ipHELIUM][0] += 7.47e-15*pow(phycon.te/1e4,2.06)* 00096 (1.+9.93*sexp(3.89e-4*phycon.te) ); 00097 00098 /* >>chng 04 jun 21 -- NPA. Put in the charge transfer rate for: 00099 He+ + H => He + H+ as defined in the UMIST database. This is only 00100 used if the "set UMIST rates" command is used */ 00101 if(!co.lgUMISTrates) 00102 { 00103 atmdat.HCharExcRecTo[ipHELIUM][0] = 4.85e-15*pow(phycon.te/300, 0.18); 00104 } 00105 /* >>chng 04 jun 21 -- NPA. update the charge transfer rates between hydrogen 00106 and oxygen to: 00107 >>refer O CT Stancil et al. 1999, A&AS, 140, 225-234 00108 Instead of using the UMIST rate, the program TableCurve was used 00109 to generate a fit to the data listed in Tables 2, 3, and 4 of the 00110 aforementioned reference. The following fitting equations agree 00111 very well with the published data. */ 00112 00113 /* At or below 10K, just use the value the formula's below give 00114 at 10K.*/ 00115 /* do both O+ -> O and O -> O+ for low T limit */ 00116 if(phycon.te <= 10. ) 00117 { 00118 atmdat.HCharExcRecTo[ipOXYGEN][0] = 3.744e-10; 00119 atmdat.HCharExcIonOf[ipOXYGEN][0] = 4.749e-20; 00120 } 00121 00122 /* this does O+ -> O for all higher temperatures */ 00123 if( phycon.te > 10.) 00124 { 00125 a_op = 2.3344302e-10; 00126 b_op = 2.3651505e-10; 00127 c_op = -1.3146803e-10; 00128 d_op = 2.9979994e-11; 00129 e_op = -2.8577012e-12; 00130 f_op = 1.1963502e-13; 00131 00132 /* equation rank 53 of TableCurve */ 00133 atmdat.HCharExcRecTo[ipOXYGEN][0] = a_op + b_op*phycon.alnte + c_op*pow(phycon.alnte,2) + d_op*pow(phycon.alnte,3) 00134 + e_op*pow(phycon.alnte,4) + f_op*pow(phycon.alnte, 5); 00135 } 00136 00137 /* now do O -> O+ 00138 * The next two equations were chosen to match up at 200K, so that there 00139 *are no discontinuities */ 00140 if((phycon.te > 10.) && (phycon.te <= 190.)) 00141 { 00142 a = -21.134531; 00143 b = -242.06831; 00144 c = 84.761441; 00145 00146 /* equation rank 2 of TableCurve */ 00147 atmdat.HCharExcIonOf[ipOXYGEN][0] = exp((a + b/SDIV(phycon.te) + c/SDIV(phycon.tesqrd))); 00148 } 00149 00150 else if((phycon.te > 190.) && (phycon.te <= 200.)) 00151 { 00152 00153 /* We are using two fitting formula's for this rate, in order to assure no 00154 sudden "jumps" in the rate, the rate between 190-200K is made to 00155 increase linearly. The formula below gets the same answer as the equation 00156 above at 190K, and gets the same answer as the the formula below this one 00157 at 200K*/ 00158 atmdat.HCharExcIonOf[ipOXYGEN][0] = 2.18733e-12*(phycon.te-190) + 1.85823e-10; 00159 } 00160 00161 else if(phycon.te > 200.) 00162 { 00163 00164 a_o = -7.6767404e-14; 00165 b_o = -3.7282001e-13; 00166 c_o = -1.488594e-12; 00167 d_o = -3.6606214e-12; 00168 e_o = 2.0699463e-12; 00169 f_o = -2.6139493e-13; 00170 g_o = 1.1580844e-14; 00171 00172 /* equation rank 120 of TableCurve */ 00173 atmdat.HCharExcIonOf[ipOXYGEN][0] = a_o + b_o*phycon.alnte + c_o*pow(phycon.alnte,2) + d_o*pow(phycon.alnte,3) 00174 + e_o*pow(phycon.alnte,4) + f_o*pow(phycon.alnte, 5) + g_o*pow(phycon.alnte,6); 00175 } 00176 00177 /* use UMIST rates for charge transfer if UMIST command is used - disagree 00178 * by 20% at 5000K and diverge at high temperature */ 00179 if(!co.lgUMISTrates) 00180 { 00181 atmdat.HCharExcIonOf[ipOXYGEN][0] = HMRATE(7.31e-10,0.23,225.9); 00182 atmdat.HCharExcRecTo[ipOXYGEN][0] = HMRATE(5.66e-10,0.36,-8.60); 00183 } 00184 00185 /* >>chng 01 may 07, following had all been +=, ok if above was zero. 00186 * changed to = and added HCTOn */ 00187 /* >>chng 01 jan 30, add following block of CT reactions */ 00188 /* ionization, added as per Phillip Stancil OK, number comes from 00189 * >>refer Fe CT Tielens, A.G.G.M., & Hollenbach, D., 1985a, ApJ, 294, 722-746 */ 00190 /* >>refer Fe CT Prasad, S.S., & Huntress, W.T., 1980, ApJS, 43, 1-35 */ 00191 /* the actual rate comes from the following paper: */ 00192 /* >>refer Fe CT Pequignot, D., & Aldrovandi, S.M.V., 1986, A&A, 161, 169-176 */ 00193 /* Fe0 + H+ => Fe+ + H0 */ 00194 /*>>chng 05 sep 15, GS, old rate had problem in predicting observed Fe I column density along HD185418. 00195 *>> refer Private communication with Stancil, data taken from ORNL web site, 00196 * "There is a well known problem with the Fe charge transfer rate coefficients: i.e., there are no accurate calculations nor or there 00197 any experiments. For Fe + H+ -> Fe+ + H, I estimated for Gary a few years ago the value of 5.4e-9. So mid way between the two values 00198 you are using. I have some notes on it in my office, but not with me. See: http://cfadc.phy.ornl.gov/astro/ps/data/home.html 00199 value changed from 3e-9 to 5.4e-9 */ 00200 00201 atmdat.HCharExcIonOf[ipIRON][0] = 5.4e-9; 00202 /*>>chng 06 sep 20 - following sets removes Fe ionization ct to be similar to Mg */ 00203 /*atmdat.HCharExcIonOf[ipIRON][0] = 0.;broken();rm this line */ 00204 00205 /* all remaining entries are from Pequignot & Aldrovandi*/ 00206 /* >>refer Al CT Pequignot, D., & Aldrovandi, S.M.V., 1986, A&A, 161, 169-176 */ 00207 /* Al0 + H+ => Al+ + H0 */ 00208 atmdat.HCharExcIonOf[ipALUMINIUM][0] = 3e-9; 00209 00210 /* >>refer P CT Pequignot, D., & Aldrovandi, S.M.V., 1986, A&A, 161, 169-176 */ 00211 /* P0 + H+ => P+ + H0 */ 00212 atmdat.HCharExcIonOf[ipPHOSPHORUS][0] = 1e-9; 00213 00214 /* >>refer Cl CT Pequignot, D., & Aldrovandi, S.M.V., 1986, A&A, 161, 169-176 */ 00215 /* Cl0 + H+ => Cl+ + H0 */ 00216 atmdat.HCharExcIonOf[ipCHLORINE][0] = 1e-9; 00217 00218 /* >>refer Ti CT Pequignot, D., & Aldrovandi, S.M.V., 1986, A&A, 161, 169-176 */ 00219 /* Ti0 + H+ => Cl+ + H0 */ 00220 atmdat.HCharExcIonOf[ipTITANIUM][0] = 3e-9; 00221 00222 /* >>refer Mn CT Pequignot, D., & Aldrovandi, S.M.V., 1986, A&A, 161, 169-176 */ 00223 /* Mn0 + H+ => Mn+ + H0 */ 00224 atmdat.HCharExcIonOf[ipMANGANESE][0] = 3e-9; 00225 00226 /* >>refer Ni CT Pequignot, D., & Aldrovandi, S.M.V., 1986, A&A, 161, 169-176 */ 00227 /* Ni0 + H+ => Ni+ + H0 */ 00228 atmdat.HCharExcIonOf[ipNICKEL][0] = 3e-9; 00229 00230 /* >>chng 01 feb 02, add following CT reaction from */ 00231 /* >>refer Na0 CT Dutta, C.M., Nordlander, P., Kimura, M., & Dalgarno, A., 2001, Phys REv A, 63, 022709 */ 00232 /* this is roughly their number around 500K - they do not give explicit values, rather 00233 * a small figure. Previous calculations were 5 orders of mag smaller at this temp. 00234 * ND this deposits electron into n=2 */ 00235 /* Na0 + H+ => Na+ + H0(n=2) */ 00236 atmdat.HCharExcIonOf[ipSODIUM][0] = 7e-12; 00237 00238 /* >>chng 05 sep 15,GS, add following CT reaction from */ 00239 /* >>refer Na0 CT Watanabe, Dutta, C.M., Nordlander, P., Kimura, M., & Dalgarno, A., 2002, Phys REv A, 66, 044701 */ 00240 /* this is roughly their number around 50K - they do not give explicit values, rather 00241 * a small figure. this deposits electron into n=1 00242 * Na0 + H+ => Na+ + H0(n=1) 00243 * add to previous rate which was for population of n=2 */ 00244 atmdat.HCharExcIonOf[ipSODIUM][0] += 0.7e-12; 00245 00246 /* >>chng 05 sep 15, GS, add following CT reaction from 00247 * >>refer K0 CT Watanabe, Dutta, C.M., Nordlander, P., Kimura, M., & Dalgarno, A., 2002, Phys REv A, 66, 044701 00248 * this is roughly their number around 50K - they do not give explicit values, rather 00249 * a small figure. 00250 * K0 + H+ => K+ + H0(n=1) */ 00251 atmdat.HCharExcIonOf[ipPOTASSIUM][0] = 1.25e-12; 00252 00253 /* >>chng 05 sep 15, GS, add following CT reaction from 00254 * >>refer S0 CT ORNL data base for charge transfer 00255 * This rate is valid for 1e3 to 1e4. Due to the small value, I did not put any limit on temp. 00256 * Earlier, other reactions also assume constant value 00257 * S0 + H+ => H + S+ */ 00258 atmdat.HCharExcIonOf[ipSULPHUR][0] = 1.e-14; 00259 00260 if( phycon.te < 1e5 ) 00261 { 00262 00263 /* >>chng 05 sep 15, GS, add following CT reaction from 00264 * >>refer Mg0 CT ORNL data base for charge transfer, 00265 * this rate is valid for temp 5e3 to 3e4, The rate goes down very fast in low temp. So I did not put a lower cut of for temp 00266 * Mg0 + H+ => H + Mg+ */ 00267 atmdat.HCharExcIonOf[ipMAGNESIUM][0] = 9.76e-12*pow((phycon.te/1e4),3.14)*(1. + 55.54*sexp(1.12*phycon.te/1e4)); 00268 /*>>chng 06 jul 20, do not allow this to fall below UMIST rate - above fit not intended for 00269 * very low temperatures */ 00270 /*>>chng 06 aug 01, UMIST is bogus - email exchange with Phillip Stancil, late July 2006 */ 00271 /*atmdat.HCharExcIonOf[ipMAGNESIUM][0] = MAX2( 1.1e-9 , atmdat.HCharExcIonOf[ipMAGNESIUM][0]);*/ 00272 /*>>chng 06 sep 20 - following sets Mg ionization ct to Fe */ 00273 /*atmdat.HCharExcIonOf[ipMAGNESIUM][0] = 5.4e-9;broken(); rm this line */ 00274 00275 /* >>chng 05 sep 15, GS, add following CT reaction from 00276 * >>refer Si0 CT ORNL data base for charge transfer 00277 * this rate is valid for temp 1e3 to 2e5, The rate goes down very fast in low temp. So I did not put a lower cut of for temp 00278 * Si0 + H+ => H + Si+ */ 00279 atmdat.HCharExcIonOf[ipSILICON][0] = 0.92e-12*pow((phycon.te/1e4),1.15)*(1. + 0.80*sexp(0.24*phycon.te/1e4)); 00280 /*>>chng 06 jul 20, do not allow this to fall below UMIST rate - above fit not intended for 00281 * very low temperatures */ 00283 /*>>chng 06 aug 01, UMIST rate is bogus as per Phillip Stancil emails of late July 2006 */ 00284 /*atmdat.HCharExcIonOf[ipSILICON][0] = MAX2( 9.9e-10 , atmdat.HCharExcIonOf[ipSILICON][0]);*/ 00285 00286 /* >>chng 05 sep 15, GS, add following CT reaction from 00287 * >>refer Li0 CT ORNL data base for charge transfer 00288 * this rate is valid for temp 1e2 to 1e4, The rate goes down very fast in low temp. So I did not put a lower cut of for temp 00289 * Li0 + H+ => H + Li+ */ 00290 atmdat.HCharExcIonOf[ipLITHIUM][0] = 2.84e-12*pow((phycon.te/1e4),1.99)*(1. + 375.54*sexp(54.07*phycon.te/1e4)); 00293 } 00294 else 00295 { 00297 atmdat.HCharExcIonOf[ipMAGNESIUM][0] = 0.; 00298 atmdat.HCharExcIonOf[ipSILICON][0] = 0.; 00299 atmdat.HCharExcIonOf[ipLITHIUM][0] = 0.; 00300 } 00301 00302 { 00303 /*>>chng 06 jul 07, Terry Yun add these charge transfer reactions */ 00304 /*>>refer N0 ct Lin, C.Y., Stancil, P.C., Gu, J.P., Buenker, R.J. & Kimura, M., 2005, Phys Rev A, 71, 062708 00305 * and combined with data from 00306 *>>refer N0 ct Butler, S.E. & Dalgarno, A. 1979, ApJ, 234, 765 */ 00307 00308 /* natural log of te */ 00309 double tefac = phycon.te * phycon.alnte; 00310 00311 /* N(4S) + H+ -> N+(3P) + H */ 00312 /* >>chng 06 jul 10, add exp for endoergic reaction */ 00313 double ct_from_n0grhp_to_npgrh0 = (1.64e-16*phycon.te - 8.76e-17*tefac + 2.41e-20*phycon.tesqrd + 9.83e-13*phycon.alnte )* 00314 sexp( 10985./phycon.te ); 00315 00316 /* N(2D) + H+ -> N+(3P) + H */ 00318 /*double ct_from_n0exhp_to_npgrh0 = 1.51e-15*phycon.te -1.61e-16*tefac + 7.74e-21*phycon.tesqrd + 1.34e-16*phycon.alnte;*/ 00319 00320 /* N+(3P) + H -> N(4S) + H+ endoergic */ 00321 double ct_from_npgrh0_to_n0grhp = (1.56e-15*phycon.te - 1.79e-16*tefac + 1.15e-20*phycon.tesqrd + 1.08e-13*phycon.alnte); 00322 00323 /* N+(3P) + H0 -> N(2D) + H+ */ 00324 /* >>chng 06 jul 10, add exp for endoergic reaction */ 00325 atmdat.HCharExcRecTo_N0_2D = (6.83e-16*phycon.te - 7.40e-17*tefac + 3.73e-21*phycon.tesqrd + 1.75e-15*phycon.alnte)* 00326 sexp( 16680./phycon.te ); 00327 00328 /* these rates are from the ground state into all possible states of the 00329 * species that is produced */ 00330 atmdat.HCharExcIonOf[ipNITROGEN][0] = ct_from_n0grhp_to_npgrh0; 00331 atmdat.HCharExcRecTo[ipNITROGEN][0] = ct_from_npgrh0_to_n0grhp + atmdat.HCharExcRecTo_N0_2D; 00332 } 00333 00334 /*>>chng 06 aug 01, update O++ and N++ -- H0 CT recombination 00335 *>>refer O3 CT Barragan, P. Errea, L.F., Mendez, L., Rabadan, I. & Riera, A. 00336 *>>refercon 2006, ApJ, 636, 544 */ 00337 /* O+2 + H -> O+ + H+ */ 00338 if( phycon.te <= 1500. ) 00339 { 00340 atmdat.HCharExcRecTo[ipOXYGEN][1] = 0.5337e-9*pow( (phycon.te/100.) ,-0.076); 00341 } 00342 else 00343 { 00344 atmdat.HCharExcRecTo[ipOXYGEN][1] = 0.4344e-9 + 00345 0.6340e-9*pow( log10(phycon.te/1500.) ,2.06 ); 00346 } 00347 00348 /* N+2 + H -> N+ + H+ */ 00349 if( phycon.te <= 1500. ) 00350 { 00351 atmdat.HCharExcRecTo[ipNITROGEN][1] = 0.8692e-9*pow( (phycon.te/1500.) ,0.17); 00352 } 00353 else if( phycon.te <= 20000. ) 00354 { 00355 atmdat.HCharExcRecTo[ipNITROGEN][1] = 0.9703e-9*pow( (phycon.te/10000.) ,0.058); 00356 } 00357 else 00358 { 00359 atmdat.HCharExcRecTo[ipNITROGEN][1] = 1.0101e-9 + 00360 1.4589e-9*pow( log10(phycon.te/20000.) ,2.06 ); 00361 } 00362 00363 /* ===================== helium charge transfer ====================*/ 00364 00365 /* atmdat.HeCharExcIonOf is ionization, */ 00366 /* [0] is Atom^0 + He+ => Atom+1 + He0 00367 * [n] is Atom^+n + He+ => Atom^+n-1 + He0 */ 00368 00369 /* atmdat.HeCharExcRecTo is recombination */ 00370 /* [0] is Atom^+1 + He0 => Atom^0 + He^+ 00371 * [n] is Atom^+n+1 + He0 => Atom^+n + He^+ */ 00372 00373 /* Carbon */ 00374 /* recombination */ 00375 /* C+3 + He => C+2 + He+ */ 00376 /* >>refer C+3 CT Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */ 00377 atmdat.HeCharExcRecTo[ipCARBON][2] = 4.6e-19*phycon.tesqrd; 00378 00379 /* C+4 + He => C+3 + He+ */ 00380 /* >>refer C+4 CT Butler, S.E., & Dalgarno, 1980b */ 00381 atmdat.HeCharExcRecTo[ipCARBON][3] = 1e-14; 00382 00383 /* ionization */ 00384 /* C0 + He+ => C+ + He0 */ 00385 /* unknown reference, from older version of the code */ 00386 /*atmdat.HeCharExcIonOf[ipCARBON][0] = 4.17e-17*(phycon.te/phycon.te03);*/ 00387 00388 /* >>chng 04 jun 21 -- update this rate to that given in the UMIST database - NPA */ 00389 atmdat.HeCharExcIonOf[ipCARBON][0] = 6.3e-15*pow((phycon.te/300),0.75); 00390 00391 /* C+1 + He+ => C+2 + He */ 00392 /* >>refer C+1 CT Butler, S.E., Heil,T.G., & Dalgarno, A. 1980, ApJ, 241, 442*/ 00393 atmdat.HeCharExcIonOf[ipCARBON][1] = 00394 5e-20*phycon.tesqrd*sexp(0.07e-4*phycon.te)*sexp(6.29/phycon.te_eV); 00395 00396 /* nitrogen */ 00397 /* recombination */ 00398 /* N+2 => N+ Butler and Dalgarno 1980B 00399 * ct with update 00400 * >>refer N+2 CT Sun, Sadeghpour, Kirby, Dalgarno, and Lafyatis, cfa preprint 4208 00401 * this agrees with exp 00402 * >>refer N+2 CT Fand&Kwong, ApJ 474, 529 */ 00403 atmdat.HeCharExcRecTo[ipNITROGEN][1] = 0.8e-10; 00404 00405 /* N+3 => N+2 */ 00406 /* >>refer N+3 CT Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */ 00407 atmdat.HeCharExcRecTo[ipNITROGEN][2] = 1.5e-10; 00408 00409 /* ce rate from quantal calc of ce, 00410 * >>refer N+4 CT Feickert, Blint, Surratt, and Watson, (preprint Sep 84). Ap.J. in press, 00411 * >>refer N+4 CT Rittby et al J Phys B 17, L677, 1984. 00412 * CR = 1.E-9 + 8E-12 * TE10 * SQRTE */ 00413 atmdat.HeCharExcRecTo[ipNITROGEN][3] = 2e-9; 00414 00415 /* ionization */ 00416 /* N+1 + He+ => N+2 + He */ 00417 /* >>refer N+1 CT Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */ 00418 atmdat.HeCharExcIonOf[ipNITROGEN][1] = 00419 3.7e-20*phycon.tesqrd*sexp(0.063e-4*phycon.te)*sexp(1.44/phycon.te_eV); 00420 00421 /* oxygen */ 00422 /* recombination */ 00423 /* O+2 + He => O+1 + He+ */ 00424 /* >>refer O+2 CT Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */ 00425 atmdat.HeCharExcRecTo[ipOXYGEN][1] = 3.2e-14*phycon.te/phycon.te05; 00426 /* O+3 + He => O+2 + He+ */ 00427 /* >>refer O+3 CT Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */ 00428 atmdat.HeCharExcRecTo[ipOXYGEN][2] = 1e-9; 00429 /* O+4 + He => O+3 + He+ */ 00430 /* >>refer O+4 CT Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */ 00431 atmdat.HeCharExcRecTo[ipOXYGEN][3] = 6e-10; 00432 00433 /* ionization */ 00434 /* O0 + He+ => O+ + He0 */ 00435 /* >>refer O0 CT Zhao et al., ApJ, 615, 1063 */ 00436 atmdat.HeCharExcIonOf[ipOXYGEN][0] = 00437 4.991E-15 * pow( phycon.te / 1e4, 0.3794 )* sexp( phycon.te/1.121E6 ) + 00438 2.780E-15 * pow( phycon.te / 1e4, -0.2163 )* exp( -1. * MIN2(1e7, phycon.te)/(-8.158E5) ); 00439 00440 /* neon */ 00441 /* recombination */ 00442 /* Ne+2 + He => Ne+1 + He+ */ 00443 /* >>refer Ne+2 CT Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */ 00444 atmdat.HeCharExcRecTo[ipNEON][1] = 1e-14; 00445 /* Ne+3 + He => Ne+2 + He+ */ 00446 /* >>refer Ne+3 CT Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */ 00447 atmdat.HeCharExcRecTo[ipNEON][2] = 1e-16*phycon.sqrte; 00448 /* Ne+4 + He => Ne+3 + He+ */ 00449 /* >>refer Ne+4 CT Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */ 00450 atmdat.HeCharExcRecTo[ipNEON][3] = 1.7e-11*phycon.sqrte; 00451 00452 /* magnesium */ 00453 /* recombination */ 00454 /* Mg+3 + Heo => Mg+2 */ 00455 /* >>refer Mg+3 CT Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */ 00456 atmdat.HeCharExcRecTo[ipMAGNESIUM][2] = 7.5e-10; 00457 /* Mg+4 + Heo => Mg+3 */ 00458 /* >>refer Mg+4 CT Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */ 00459 atmdat.HeCharExcRecTo[ipMAGNESIUM][3] = 1.4e-10*phycon.te30; 00460 00461 00462 /* silicon */ 00463 /* recombination */ 00464 /* Si+3 +He => Si+2 + He+ */ 00465 /* >>refer Si+3 CT Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 00466 * scale by 1.3 to bring into agreement with 00467 * >>refer Si+3 CT Fang, Z., & Kwong, H.S. 1997 ApJ 483, 527 */ 00468 atmdat.HeCharExcRecTo[ipSILICON][2] += phycon.sqrte*phycon.te10*phycon.te10* 00469 1.3*1.5e-12; 00470 00471 /* Si+4 + Heo => Si+3 00472 * >>refer Si+4 CT Opradolce et al., 1985, A&A, 148, 229 */ 00473 atmdat.HeCharExcRecTo[ipSILICON][3] = 2.54e-11*phycon.sqrte/phycon.te03/ 00474 phycon.te01/phycon.te01; 00475 00476 /* ionization */ 00477 /* Si0 + He+ => Si+ + He0 */ 00478 /* >>refer Si0 CT Prasad, S.S., & Huntress, W.T., 1980, ApJS, 43, 1-35 */ 00479 atmdat.HeCharExcIonOf[ipSILICON][0] = 3.3e-9; 00480 00481 /* Si+1 + He+ => Si+2 + He */ 00482 /* >>refer Si+1 CT Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */ 00483 atmdat.HeCharExcIonOf[ipSILICON][1] = 00484 1.5e-11*phycon.te20*phycon.te05*sexp(6.91/phycon.te_eV); 00485 00486 /* Si+2 + He+ => Si+3 + He */ 00487 /* >>refer Si+2 CT Gargaud, M., McCarroll, R., & Valiron, P. 1982, A&ASup, 45, 603 */ 00488 atmdat.HeCharExcIonOf[ipSILICON][2] = 00489 1.15e-11*phycon.sqrte*sexp(8.88/phycon.te_eV); 00490 00491 /* sulphur */ 00492 /* recombination */ 00493 /* S+3 + Heo => S+2 */ 00494 /* >>refer S+3 CT Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */ 00495 atmdat.HeCharExcRecTo[ipSULPHUR][2] = phycon.sqrte*1.1e-11; 00496 00497 /* S+4 + Heo => S+3 */ 00498 /* >>refer S+4 CT Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */ 00499 /* >>chng 04 jul 01, from [ipSULPHUR][2] to [3] - bug */ 00500 atmdat.HeCharExcRecTo[ipSULPHUR][3] = 4.8e-14*phycon.te30; 00501 00502 /* ionization */ 00503 /* S+1 + He+ => S+2 + He */ 00504 /* >>refer S+1 CT Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */ 00505 atmdat.HeCharExcIonOf[ipSULPHUR][1] = 00506 4.4e-16*phycon.te*phycon.te20*sexp(0.036e-4*phycon.te)*sexp(9.2/phycon.te_eV); 00507 00508 /* S+2 + He+ => S+3 + He */ 00509 /* >>refer S+2 CT Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */ 00510 atmdat.HeCharExcIonOf[ipSULPHUR][2] = 00511 5.5e-18*phycon.te*phycon.sqrte*phycon.te10*sexp(0.046e-4*phycon.te)*sexp(10.5/phycon.te_eV); 00512 00513 /* Argon */ 00514 /* recombination */ 00515 /* >>refer Ar+2 CT Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */ 00516 atmdat.HeCharExcRecTo[ipARGON][1] = 1.3e-10; 00517 00518 /* >>refer Ar+3 CT Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */ 00519 atmdat.HeCharExcRecTo[ipARGON][2] = 1.e-14; 00520 00521 /* >>refer Ar+4 CT Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */ 00522 atmdat.HeCharExcRecTo[ipARGON][3] = 1.6e-8/phycon.te30; 00523 00524 /* ionization */ 00525 /* Ar+1 + He+ => Ar+2 + He0 */ 00526 /* >>refer Ar+1 CT Butler, S.E., & Dalgarno, A., 1980b, ApJ, 241, 838 */ 00527 atmdat.HeCharExcIonOf[ipARGON][1] = 1.1e-10; 00528 00529 TeUsed = phycon.te; 00530 00531 if(!co.lgUMISTrates) 00532 { 00533 /* Set all charge transfer rates equal to zero that do not appear 00534 in the UMIST database. This if statement is only performed 00535 if the "set UMIST rates" command is used */ 00536 00537 atmdat.HCharExcIonOf[ipHELIUM][0] = 0; 00538 atmdat.HCharExcIonOf[ipCARBON][0] = 0; 00539 atmdat.HCharExcRecTo[ipCARBON][0] = 0; 00540 00541 atmdat.HeCharExcRecTo[ipCARBON][0] = 0; 00542 atmdat.HeCharExcIonOf[ipOXYGEN][0] = 0; 00543 atmdat.HeCharExcRecTo[ipOXYGEN][0] = 0; 00544 } 00545 00546 00547 /* this is set false with the no charge transfer command */ 00548 if( !atmdat.lgCTOn ) 00549 { 00550 for( nelem=0; nelem< LIMELM; ++nelem ) 00551 { 00552 for( ion=0; ion<LIMELM; ++ion ) 00553 { 00554 atmdat.HCharExcIonOf[nelem][ion] = 0.; 00555 atmdat.HCharExcRecTo[nelem][ion] = 0.; 00556 atmdat.HeCharExcIonOf[nelem][ion] = 0.; 00557 atmdat.HeCharExcRecTo[nelem][ion] = 0.; 00558 } 00559 } 00560 } 00561 } 00562 00563 /* return value is H ionization rate (s-1) due to O+ charge transfer */ 00564 *O_HIonRate = atmdat.HCharExcRecTo[ipOXYGEN][0]*dense.xIonDense[ipOXYGEN][1]; 00565 return; 00566 } 00567 00568 /*================================================================================* 00569 *================================================================================*/ 00570 double ChargTranSumHeat(void) 00571 { 00572 long int ion, 00573 nelem; 00574 double SumCTHeat_v; 00575 00576 DEBUG_ENTRY( "ChargTranSumHeat()" ); 00577 00578 /* second dimension is ionization stage, 00579 * 1=+1 for parent, etc 00580 * third dimension is atomic weight of atom */ 00581 00582 /* make sure data are initialized */ 00583 ASSERT( lgCTDataDefined ); 00584 00585 SumCTHeat_v = 0.; 00586 /* >>chng 01 apr 25, lower limit had been 0 should have been 1 (helium) */ 00587 for( nelem=ipHELIUM; nelem < LIMELM; nelem++ ) 00588 { 00589 /* >>chng >>01 apr 25, loops had been to LIMELM, which may have done no harm 00590 * since extra array elements were set to zero, but is incorrect since the physical 00591 * limit is the number of stages of ionization */ 00592 int limit = MIN2(4, nelem); 00593 /* this first group of lower stages of ionization have exact rate coefficients */ 00594 for( ion=0; ion < limit; ion++ ) 00595 { 00596 /* CTIonData[nelem][ion][7] and CTRecombData[nelem][ion][6] are the energy deficits in eV, 00597 * atmdat.HCharExcIonOf[nelem][ion] and atmdat.HCharExcIonOf[nelem][ion] 00598 * save the rate coefficients 00599 * this is sum of heat exchange in eV s^-1 cm^-3 */ 00600 SumCTHeat_v += 00601 00602 /* heating due to ionization of heavy element, recombination of hydrogen */ 00603 atmdat.HCharExcIonOf[nelem][ion]*CTIonData[nelem][ion][7]* 00604 (double)dense.xIonDense[ipHYDROGEN][1]*(double)dense.xIonDense[nelem][ion] + 00605 00606 /* heating due to recombination of heavy element, ionization of hydrogen */ 00607 atmdat.HCharExcRecTo[nelem][ion]*CTRecombData[nelem][ion][6]* 00608 StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop*dense.xIonDense[ipHYDROGEN][1]* 00609 (double)dense.xIonDense[nelem][ion+1]; 00610 } 00611 00612 /* >>chng >>01 apr 25, following loop had been to LIMELM, change to nelem */ 00613 /* following do not have exact energies, so use 2.86*(Z-1) */ 00614 for( ion=4; ion < nelem; ion++ ) 00615 { 00616 SumCTHeat_v += 00617 atmdat.HCharExcRecTo[nelem][ion]* 2.86*(double)ion * 00618 (double)dense.xIonDense[ipHYDROGEN][0]*(double)dense.xIonDense[nelem][ion+1]; 00619 } 00620 } 00621 00622 /* convert from eV to ergs, HCharHeatOn usually 1, set to 0 with no CTHeat, 00623 * EN1EV is ergs in 1 eV, 1.602176e-012*/ 00624 SumCTHeat_v *= EN1EV * atmdat.HCharHeatOn; 00625 00626 if( thermal.htot > 1e-35 ) 00627 { 00628 /* remember largest fractions of heating and cooling for comment */ 00629 atmdat.HCharHeatMax = MAX2(atmdat.HCharHeatMax, 00630 SumCTHeat_v/thermal.htot ); 00631 00632 atmdat.HCharCoolMax = MAX2(atmdat.HCharCoolMax, 00633 -SumCTHeat_v/thermal.htot); 00634 } 00635 00636 /* debug code to print out the contributors to total CT heating */ 00637 { 00638 /*@-redef@*/ 00639 enum {DEBUG_LOC=false}; 00640 /*@+redef@*/ 00641 if( DEBUG_LOC) 00642 { 00643 # define FRAC 0.1 00644 for( nelem=ipHELIUM; nelem < LIMELM; nelem++ ) 00645 { 00646 /* >>chng >>01 apr 25, loops had been to LIMELM, which may have done no harm 00647 * since extra array elements were set to zero, but is incorrect since the physical 00648 * limit is the number of stages of ionization */ 00649 int limit = MIN2(4, nelem); 00650 /* this first group of lower stages of ionization have exact rate coefficients */ 00651 for( ion=0; ion < limit; ion++ ) 00652 { 00653 if( 00654 /* heating due to ionization of heavy element, recombination of hydrogen */ 00655 (atmdat.HCharExcIonOf[nelem][ion]*CTIonData[nelem][ion][7]* 00656 (double)dense.xIonDense[ipHYDROGEN][1]*(double)dense.xIonDense[nelem][ion] + 00657 00658 /* heating due to recombination of heavy element, ionization of hydrogen */ 00659 atmdat.HCharExcRecTo[nelem][ion]*CTRecombData[nelem][ion][6]* 00660 StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop*dense.xIonDense[ipHYDROGEN][1]* 00661 (double)dense.xIonDense[nelem][ion+1])/SumCTHeat_v> FRAC ) 00662 00663 fprintf(ioQQQ,"DEBUG ct %li %li %.3f\n", nelem, ion, 00664 (atmdat.HCharExcIonOf[nelem][ion]*CTIonData[nelem][ion][7]* 00665 (double)dense.xIonDense[ipHYDROGEN][1]*(double)dense.xIonDense[nelem][ion] + 00666 00667 /* heating due to recombination of heavy element, ionization of hydrogen */ 00668 atmdat.HCharExcRecTo[nelem][ion]*CTRecombData[nelem][ion][6]* 00669 StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop*dense.xIonDense[ipHYDROGEN][1]* 00670 (double)dense.xIonDense[nelem][ion+1]) ); 00671 } 00672 00673 for( ion=4; ion < nelem; ion++ ) 00674 { 00675 if( 00676 (atmdat.HCharExcRecTo[nelem][ion]* 2.86*(double)ion * 00677 (double)dense.xIonDense[ipHYDROGEN][0]*(double)dense.xIonDense[nelem][ion+1])/SumCTHeat_v> FRAC ) 00678 fprintf(ioQQQ,"DEBUG ct %li %li %.3f\n", nelem, ion, 00679 (atmdat.HCharExcRecTo[nelem][ion]* 2.86*(double)ion * 00680 (double)dense.xIonDense[ipHYDROGEN][0]*(double)dense.xIonDense[nelem][ion+1]) ); 00681 } 00682 } 00683 # undef FRAC 00684 fprintf(ioQQQ,"DEBUT ct tot %.e3\n", SumCTHeat_v / thermal.htot ); 00685 } 00686 } 00687 return( SumCTHeat_v ); 00688 } 00689 00690 /*================================================================================* 00691 *================================================================================*/ 00692 STATIC double HCTIon( 00693 /* ion is stage of ionization on C scale, 0 for atom */ 00694 long int ion, 00695 /* nelem is atomic number of element on C scale, 1 to 29 */ 00696 /* HCTIon(1,5) is C+ + H+ => C++ + H */ 00697 long int nelem) 00698 { 00699 double HCTIon_v, 00700 tused; 00701 00702 DEBUG_ENTRY( "HCTIon()" ); 00703 /* H charge transfer ionization, using Jim Kingdon's ctdata.for */ 00704 00705 /* set up the rate coefficients if this is first call */ 00706 if( !lgCTDataDefined ) 00707 { 00708 if( trace.lgTrace ) 00709 { 00710 fprintf( ioQQQ," HCTIon doing 1-time init of charge transfer data\n"); 00711 } 00712 lgCTDataDefined = true; 00713 MakeHCTData(); 00714 } 00715 00716 /* check that data have been linked in, 00717 * error here would mean that data below had not been loaded */ 00718 ASSERT( CTIonData[2][0][0] > 0. ); 00719 00720 /* return zero for highly ionized species */ 00721 if( ion >= 3 ) 00722 { 00723 HCTIon_v = 0.; 00724 return( HCTIon_v ); 00725 } 00726 00727 /*begin sanity checks */ 00728 /* check that ionization stage is ok for this element*/ 00729 ASSERT( ion >= 0); 00730 ASSERT( ion <= nelem ); 00731 00732 /* now check the element is valid, this is routine HCTIon */ 00733 ASSERT( nelem > 0 ); 00734 ASSERT( nelem < LIMELM ); 00735 00736 /* may be no entry, first coefficient is zero in this case */ 00737 if( CTIonData[nelem][ion][0] <= 0. ) 00738 { 00739 HCTIon_v = 0.; 00740 00741 } 00742 00743 else 00744 { 00745 /* Make sure te is between temp. boundaries; set constant outside of range */ 00746 tused = MAX2((double)phycon.te,CTIonData[nelem][ion][4]); 00747 tused = MIN2(tused,CTIonData[nelem][ion][5]); 00748 tused *= 1e-4; 00749 00750 /* the interpolation equation */ 00751 HCTIon_v = CTIonData[nelem][ion][0]*1e-9*(pow(tused,CTIonData[nelem][ion][1]))* 00752 (1. + CTIonData[nelem][ion][2]*exp(CTIonData[nelem][ion][3]*tused))* 00753 exp(-CTIonData[nelem][ion][6]/tused); 00754 } 00755 return( HCTIon_v ); 00756 } 00757 00758 /*================================================================================* 00759 *================================================================================*/ 00760 STATIC double HCTRecom( 00761 /* ion is stage of ionization on C scale, 0 for rec to atom */ 00762 long int ion, 00763 /* nelem is atomic number of element on C scale, 1 = he up to LIMELM */ 00764 /* HCTRecom(1,5) would be C+2 + H => C+ + H+ */ 00765 long int nelem) 00766 { 00767 double HCTRecom_v, 00768 tused; 00769 00770 DEBUG_ENTRY( "HCTRecom()" ); 00771 /* 00772 * H charge transfer recombination using Jim Kingdon's block ctdata.for 00773 */ 00774 00775 /* set up the rate coefficients if this is first call */ 00776 if( !lgCTDataDefined ) 00777 { 00778 if( trace.lgTrace ) 00779 { 00780 fprintf( ioQQQ," HCTIon doing 1-time init of charge transfer data\n"); 00781 } 00782 lgCTDataDefined = true; 00783 MakeHCTData(); 00784 } 00785 00786 /* this is check that data have been set up properly, will 00787 * fail if arrays are not initialized properly */ 00788 ASSERT( CTRecombData[1][0][0] > 0. ); 00789 00790 /* use Dalgarno estimate for highly ionized species, number reset with 00791 * set charge transfer command */ 00792 if( ion > 3 ) 00793 { 00794 /* >>chng 96 nov 25, added this option, default is 1.92e-9 00795 * Dalgarno's charge transfer */ 00796 HCTRecom_v = atmdat.HCTAlex*((double)ion+1.); 00797 return( HCTRecom_v ); 00798 } 00799 00800 /* check that ion stage within bound for this atom */ 00801 ASSERT( ion >= 0 && ion <= nelem ); 00802 00803 /* now check the element is valid, this is routine HCTIon */ 00804 ASSERT( nelem > 0 && nelem < LIMELM ); 00805 00806 tused = MAX2((double)phycon.te,CTRecombData[nelem][ion][4]); 00807 tused = MIN2(tused,CTRecombData[nelem][ion][5]); 00808 tused *= 1e-4; 00809 00810 if( tused == 0. ) 00811 { 00812 HCTRecom_v = 0.; 00813 return( HCTRecom_v ); 00814 } 00815 00816 /* the interpolation equation */ 00817 HCTRecom_v = CTRecombData[nelem][ion][0]*1e-9*(pow(tused,CTRecombData[nelem][ion][1]))* 00818 (1. + CTRecombData[nelem][ion][2]*sexp(-CTRecombData[nelem][ion][3]*tused)); 00819 00820 /* in sexp negative sign not typo - there are negative signs already 00821 * in coefficient, and sexp has implicit negative sign */ 00822 return( HCTRecom_v ); 00823 } 00824 00825 /*================================================================================* 00826 *================================================================================*/ 00827 /*block data with Jim Kingdon's charge transfer data */ 00828 /* >>refer H CT Kingdon, J, B., & Ferland, G.J. 1996, ApJS, 106, 205 */ 00829 /* 00830 * first dimension is atomic number of atom, 0 for H 00831 * second dimension is ionization stage, 00832 * 1=+0 for parent, etc 00833 * third dimension is atomic number of atom 00834 * second dimension is ionization stage, 00835 * 1=+1 for parent, etc 00836 */ 00837 00838 /* digital form of the fits to the charge transfer 00839 * ionization rate coefficients 00840 * 00841 * Note: First parameter is in units of 1e-9! 00842 * Note: Seventh parameter is in units of 1e4 K */ 00843 00844 /* digital form of the fits to the charge transfer 00845 * recombination rate coefficients (total) 00846 * 00847 * Note: First parameter is in units of 1e-9! 00848 * recombination 00849 */ 00850 00851 /* holds data for charge transfer fits */ 00852 STATIC void MakeHCTData(void) 00853 { 00854 long int i, 00855 j, 00856 nelem, 00857 _r; 00858 00859 DEBUG_ENTRY( "MakeHCTData()" ); 00860 00861 /* >>chng 01 apr 24, zero out this block, as per PvH comments that 00862 * translated block data's do not fully initialize arrays */ 00863 /* first zero out entire arrays, since some may not have charge transfer data */ 00864 for( nelem=0; nelem<LIMELM; ++nelem ) 00865 { 00866 for( i=0; i<4; ++i ) 00867 { 00868 for( j=0; j<7; ++j ) 00869 { 00870 CTIonData[nelem][i][j] = 0.; 00871 CTRecombData[nelem][i][j] = 0.; 00872 } 00873 CTIonData[nelem][i][7] = 0.; 00874 } 00875 } 00876 00877 /* 00878 * following are coefficients for charge transfer ionization, 00879 * H+ + A => H + A+ 00880 */ 00881 /* Lithium +0 */ 00882 { static double _itmp0[] = {2.84e-3 , 1.99 , 375.54 , -54.07 , 1e2 , 1e4 , 0., 00883 -10.}; 00884 00885 for( i=1, _r = 0; i <= 8; i++ ) 00886 { 00887 CTIonData[2][0][i-1] = _itmp0[_r++]; 00888 } 00889 } 00890 00891 /* C+0 ionization */ 00892 { static double _itmp1[] = {1.07e-6 , 3.15 , 176.43 , -4.29 , 1e3 , 1e5 , 0. ,2.34}; 00893 for( i=1, _r = 0; i <= 8; i++ ) 00894 { 00895 CTIonData[5][0][i-1] = _itmp1[_r++]; 00896 } 00897 } 00898 { static double _itmp2[] = {4.55e-3,-0.29,-0.92,-8.38,1e2,5e4, 00899 1.086,-0.94}; 00900 for( i=1, _r = 0; i <= 8; i++ ) 00901 { 00902 CTIonData[6][0][i-1] = _itmp2[_r++]; 00903 } 00904 } 00905 /* oxygen */ 00906 { static double _itmp3[] = {7.40e-2,0.47,24.37,-0.74,1e1,1e4, 00907 0.023,-0.02}; 00908 for( i=1, _r = 0; i <= 8; i++ ) 00909 { 00910 CTIonData[7][0][i-1] = _itmp3[_r++]; 00911 } 00912 } 00913 { static double _itmp4[] = {3.34e-6,9.31,2632.31,-3.04,1e3, 00914 2e4,0.0,-1.74}; 00915 for( i=1, _r = 0; i <= 8; i++ ) 00916 { 00917 CTIonData[10][0][i-1] = _itmp4[_r++]; 00918 } 00919 } 00920 { static double _itmp5[] = {9.76e-3,3.14,55.54,-1.12,5e3,3e4, 00921 0.0,1.52}; 00922 for( i=1, _r = 0; i <= 8; i++ ) 00923 { 00924 CTIonData[11][0][i-1] = _itmp5[_r++]; 00925 } 00926 } 00927 { static double _itmp6[] = {7.60e-5,0.00,-1.97,-4.32,1e4,3e5, 00928 1.670,-1.44}; 00929 for( i=1, _r = 0; i <= 8; i++ ) 00930 { 00931 CTIonData[11][1][i-1] = _itmp6[_r++]; 00932 } 00933 } 00934 { static double _itmp7[] = {0.92,1.15,0.80,-0.24,1e3,2e5,0.0, 00935 0.12}; 00936 for( i=1, _r = 0; i <= 8; i++ ) 00937 { 00938 CTIonData[13][0][i-1] = _itmp7[_r++]; 00939 } 00940 } 00941 /* Si+1 ionization */ 00942 { static double _itmp8[] = {2.26 , 7.36e-2 , -0.43 , -0.11 , 2e3 , 1e5 , 3.031 00943 ,-2.72}; 00944 for( i=1, _r = 0; i <= 8; i++ ) 00945 { 00946 CTIonData[13][1][i-1] = _itmp8[_r++]; 00947 } 00948 } 00949 { static double _itmp9[] = {1.00e-5,0.00,0.00,0.00,1e3,1e4, 00950 0.0,-3.24}; 00951 for( i=1, _r = 0; i <= 8; i++ ) 00952 { 00953 CTIonData[15][0][i-1] = _itmp9[_r++]; 00954 } 00955 } 00956 { static double _itmp10[] = {4.39,0.61,-0.89,-3.56,1e3,3e4, 00957 3.349,-2.89}; 00958 for( i=1, _r = 0; i <= 8; i++ ) 00959 { 00960 CTIonData[23][1][i-1] = _itmp10[_r++]; 00961 } 00962 } 00963 { static double _itmp11[] = {2.83e-1,6.80e-3,6.44e-2,-9.70, 00964 1e3,3e4,2.368,-2.04}; 00965 for( i=1, _r = 0; i <= 8; i++ ) 00966 { 00967 CTIonData[24][1][i-1] = _itmp11[_r++]; 00968 } 00969 } 00970 { static double _itmp12[] = {2.10,7.72e-2,-0.41,-7.31,1e4,1e5, 00971 3.005,-2.56}; 00972 for( i=1, _r = 0; i <= 8; i++ ) 00973 { 00974 CTIonData[25][1][i-1] = _itmp12[_r++]; 00975 } 00976 } 00977 { static double _itmp13[] = {1.20e-2,3.49,24.41,-1.26,1e3,3e4, 00978 4.044,-3.49}; 00979 for( i=1, _r = 0; i <= 8; i++ ) 00980 { 00981 CTIonData[26][1][i-1] = _itmp13[_r++]; 00982 } 00983 } 00984 /* CT recombination, A+n + H => A+n-1 + H+ */ 00985 /* >>chng 01 may 03, first coefficient multiplied by 0.25, as per comment in 00986 * >>refer Li CT Stancil, P.C., & Zygelman, B., 1996, ApJ, 472, 102 00987 * which corrected the error in 00988 * >>refer He CT Zygelman, B., Dalgarno, A., Kimura, M., & Lane, N.F., 00989 * >>refercon 1989, Phys. Rev. A, 40, 2340 00990 * this was used in the original Kingdon & Ferland paper so no correction required 00991 * >>chng 04 apr 27, He was in error above as well, factor of 4, noted in 00992 * >>refer He CT Stancil, P.C., Lepp, S., & Dalgarno, A. 1998, ApJ, 509, 1 00993 */ 00994 { static double _itmp14[] = {/*7.47e-6*/1.87e-6,2.06,9.93,-3.89,6e3,1e5, 00995 10.99}; 00996 for( i=1, _r = 0; i <= 7; i++ ) 00997 { 00998 CTRecombData[1][0][i-1] = _itmp14[_r++]; 00999 } 01000 } 01001 { static double _itmp15[] = {1.00e-5,0.,0.,0.,1e3,1e7,-40.81}; 01002 for( i=1, _r = 0; i <= 7; i++ ) 01003 { 01004 CTRecombData[1][1][i-1] = _itmp15[_r++]; 01005 } 01006 } 01007 for( i=1; i <= 7; i++ ) 01008 { 01009 CTRecombData[2][0][i-1] = 0.f; 01010 } 01011 { static double _itmp16[] = {1.26,0.96,3.02,-0.65,1e3,3e4,3.02}; 01012 for( i=1, _r = 0; i <= 7; i++ ) 01013 { 01014 CTRecombData[2][1][i-1] = _itmp16[_r++]; 01015 } 01016 } 01017 { static double _itmp17[] = {1.00e-5,0.,0.,0.,2e3,5e4,-108.83}; 01018 for( i=1, _r = 0; i <= 7; i++ ) 01019 { 01020 CTRecombData[2][2][i-1] = _itmp17[_r++]; 01021 } 01022 } 01023 for( i=1; i <= 7; i++ ) 01024 { 01025 CTRecombData[3][0][i-1] = 0.f; 01026 } 01027 { static double _itmp18[] = {1.00e-5,0.,0.,0.,2e3,5e4,-4.61}; 01028 for( i=1, _r = 0; i <= 7; i++ ) 01029 { 01030 CTRecombData[3][1][i-1] = _itmp18[_r++]; 01031 } 01032 } 01033 { static double _itmp19[] = {1.00e-5,0.,0.,0.,2e3,5e4,-140.26}; 01034 for( i=1, _r = 0; i <= 7; i++ ) 01035 { 01036 CTRecombData[3][2][i-1] = _itmp19[_r++]; 01037 } 01038 } 01039 { static double _itmp20[] = {5.17,0.82,-0.69,-1.12,2e3,5e4, 01040 10.59}; 01041 for( i=1, _r = 0; i <= 7; i++ ) 01042 { 01043 CTRecombData[3][3][i-1] = _itmp20[_r++]; 01044 } 01045 } 01046 for( i=1; i <= 7; i++ ) 01047 { 01048 CTRecombData[4][0][i-1] = 0.f; 01049 } 01050 { static double _itmp21[] = {2.00e-2,0.,0.,0.,1e3,1e9,2.46}; 01051 for( i=1, _r = 0; i <= 7; i++ ) 01052 { 01053 CTRecombData[4][1][i-1] = _itmp21[_r++]; 01054 } 01055 } 01056 { static double _itmp22[] = {1.00e-5,0.,0.,0.,2e3,5e4,-24.33}; 01057 for( i=1, _r = 0; i <= 7; i++ ) 01058 { 01059 CTRecombData[4][2][i-1] = _itmp22[_r++]; 01060 } 01061 } 01062 /* B+4 recombinatino */ 01063 { static double _itmp23[] = {2.74 , 0.93 , -0.61 , -1.13 , 2e3 , 5e4 , 01064 11.}; 01065 for( i=1, _r = 0; i <= 7; i++ ) 01066 { 01067 CTRecombData[4][3][i-1] = _itmp23[_r++]; 01068 } 01069 } 01070 /* C+1 recombinatino */ 01071 { static double _itmp24[] = {4.88e-7 , 3.25 , -1.12 , -0.21 , 5.5e3 , 1e5 , 01072 -2.34}; 01073 for( i=1, _r = 0; i <= 7; i++ ) 01074 { 01075 CTRecombData[5][0][i-1] = _itmp24[_r++]; 01076 } 01077 } 01078 { static double _itmp25[] = {1.67e-4,2.79,304.72,-4.07,5e3, 01079 5e4,4.01}; 01080 for( i=1, _r = 0; i <= 7; i++ ) 01081 { 01082 CTRecombData[5][1][i-1] = _itmp25[_r++]; 01083 } 01084 } 01085 { static double _itmp26[] = {3.25,0.21,0.19,-3.29,1e3,1e5,5.73}; 01086 for( i=1, _r = 0; i <= 7; i++ ) 01087 { 01088 CTRecombData[5][2][i-1] = _itmp26[_r++]; 01089 } 01090 } 01091 { static double _itmp27[] = {332.46,-0.11,-9.95e-1,-1.58e-3, 01092 1e1,1e5,11.30}; 01093 for( i=1, _r = 0; i <= 7; i++ ) 01094 { 01095 CTRecombData[5][3][i-1] = _itmp27[_r++]; 01096 } 01097 } 01098 { static double _itmp28[] = {1.01e-3,-0.29,-0.92,-8.38,1e2, 01099 5e4,0.94}; 01100 for( i=1, _r = 0; i <= 7; i++ ) 01101 { 01102 CTRecombData[6][0][i-1] = _itmp28[_r++]; 01103 } 01104 } 01105 { static double _itmp29[] = {3.05e-1,0.60,2.65,-0.93,1e3,1e5, 01106 4.56}; 01107 for( i=1, _r = 0; i <= 7; i++ ) 01108 { 01109 CTRecombData[6][1][i-1] = _itmp29[_r++]; 01110 } 01111 } 01112 { static double _itmp30[] = {4.54,0.57,-0.65,-0.89,1e1,1e5, 01113 6.40}; 01114 for( i=1, _r = 0; i <= 7; i++ ) 01115 { 01116 CTRecombData[6][2][i-1] = _itmp30[_r++]; 01117 } 01118 } 01119 /* N+4 recombination */ 01120 { static double _itmp31[] = { 2.95 , 0.55 , -0.39 , -1.07 , 1e3 , 1e6 , 01121 11.}; 01122 for( i=1, _r = 0; i <= 7; i++ ) 01123 { 01124 CTRecombData[6][3][i-1] = _itmp31[_r++]; 01125 } 01126 } 01127 { static double _itmp32[] = {1.04,3.15e-2,-0.61,-9.73,1e1,1e4, 01128 0.02}; 01129 for( i=1, _r = 0; i <= 7; i++ ) 01130 { 01131 CTRecombData[7][0][i-1] = _itmp32[_r++]; 01132 } 01133 } 01134 { static double _itmp33[] = {1.04,0.27,2.02,-5.92,1e2,1e5,6.65}; 01135 for( i=1, _r = 0; i <= 7; i++ ) 01136 { 01137 CTRecombData[7][1][i-1] = _itmp33[_r++]; 01138 } 01139 } 01140 { static double _itmp34[] = {3.98,0.26,0.56,-2.62,1e3,5e4,5.}; 01141 for( i=1, _r = 0; i <= 7; i++ ) 01142 { 01143 CTRecombData[7][2][i-1] = _itmp34[_r++]; 01144 } 01145 } 01146 { static double _itmp35[] = {2.52e-1,0.63,2.08,-4.16,1e3,3e4, 01147 8.47}; 01148 for( i=1, _r = 0; i <= 7; i++ ) 01149 { 01150 CTRecombData[7][3][i-1] = _itmp35[_r++]; 01151 } 01152 } 01153 for( i=1; i <= 7; i++ ) 01154 { 01155 CTRecombData[8][0][i-1] = 0.f; 01156 } 01157 { static double _itmp36[] = {1.00e-5,0.,0.,0.,2e3,5e4,-21.37}; 01158 for( i=1, _r = 0; i <= 7; i++ ) 01159 { 01160 CTRecombData[8][1][i-1] = _itmp36[_r++]; 01161 } 01162 } 01163 { static double _itmp37[] = {9.86,0.29,-0.21,-1.15,2e3,5e4, 01164 5.6}; 01165 for( i=1, _r = 0; i <= 7; i++ ) 01166 { 01167 CTRecombData[8][2][i-1] = _itmp37[_r++]; 01168 } 01169 } 01170 { static double _itmp38[] = {7.15e-1,1.21,-0.70,-0.85,2e3,5e4, 01171 11.8}; 01172 for( i=1, _r = 0; i <= 7; i++ ) 01173 { 01174 CTRecombData[8][3][i-1] = _itmp38[_r++]; 01175 } 01176 } 01177 for( i=1; i <= 7; i++ ) 01178 { 01179 CTRecombData[9][0][i-1] = 0.f; 01180 } 01181 { static double _itmp39[] = {1.00e-5,0.,0.,0.,5e3,5e4,-27.36}; 01182 for( i=1, _r = 0; i <= 7; i++ ) 01183 { 01184 CTRecombData[9][1][i-1] = _itmp39[_r++]; 01185 } 01186 } 01187 { static double _itmp40[] = {14.73,4.52e-2,-0.84,-0.31,5e3, 01188 5e4,5.82}; 01189 for( i=1, _r = 0; i <= 7; i++ ) 01190 { 01191 CTRecombData[9][2][i-1] = _itmp40[_r++]; 01192 } 01193 } 01194 { static double _itmp41[] = {6.47,0.54,3.59,-5.22,1e3,3e4,8.60}; 01195 for( i=1, _r = 0; i <= 7; i++ ) 01196 { 01197 CTRecombData[9][3][i-1] = _itmp41[_r++]; 01198 } 01199 } 01200 for( i=1; i <= 7; i++ ) 01201 { 01202 CTRecombData[10][0][i-1] = 0.f; 01203 } 01204 { static double _itmp42[] = {1.00e-5,0.,0.,0.,2e3,5e4,-33.68}; 01205 for( i=1, _r = 0; i <= 7; i++ ) 01206 { 01207 CTRecombData[10][1][i-1] = _itmp42[_r++]; 01208 } 01209 } 01210 { static double _itmp43[] = {1.33,1.15,1.20,-0.32,2e3,5e4,6.25}; 01211 for( i=1, _r = 0; i <= 7; i++ ) 01212 { 01213 CTRecombData[10][2][i-1] = _itmp43[_r++]; 01214 } 01215 } 01216 { static double _itmp44[] = {1.01e-1,1.34,10.05,-6.41,2e3,5e4, 01217 11.}; 01218 for( i=1, _r = 0; i <= 7; i++ ) 01219 { 01220 CTRecombData[10][3][i-1] = _itmp44[_r++]; 01221 } 01222 } 01223 for( i=1; i <= 7; i++ ) 01224 { 01225 CTRecombData[11][0][i-1] = 0.f; 01226 } 01227 { static double _itmp45[] = {8.58e-5,2.49e-3,2.93e-2,-4.33, 01228 1e3,3e4,1.44}; 01229 for( i=1, _r = 0; i <= 7; i++ ) 01230 { 01231 CTRecombData[11][1][i-1] = _itmp45[_r++]; 01232 } 01233 } 01234 { static double _itmp46[] = {6.49,0.53,2.82,-7.63,1e3,3e4,5.73}; 01235 for( i=1, _r = 0; i <= 7; i++ ) 01236 { 01237 CTRecombData[11][2][i-1] = _itmp46[_r++]; 01238 } 01239 } 01240 { static double _itmp47[] = {6.36,0.55,3.86,-5.19,1e3,3e4,8.60}; 01241 for( i=1, _r = 0; i <= 7; i++ ) 01242 { 01243 CTRecombData[11][3][i-1] = _itmp47[_r++]; 01244 } 01245 } 01246 for( i=1; i <= 7; i++ ) 01247 { 01248 CTRecombData[12][0][i-1] = 0.f; 01249 } 01250 { static double _itmp48[] = {1.00e-5,0.,0.,0.,1e3,3e4,-5.23}; 01251 for( i=1, _r = 0; i <= 7; i++ ) 01252 { 01253 CTRecombData[12][1][i-1] = _itmp48[_r++]; 01254 } 01255 } 01256 { static double _itmp49[] = {7.11e-5,4.12,1.72e4,-22.24,1e3, 01257 3e4,8.17}; 01258 for( i=1, _r = 0; i <= 7; i++ ) 01259 { 01260 CTRecombData[12][2][i-1] = _itmp49[_r++]; 01261 } 01262 } 01263 { static double _itmp50[] = {7.52e-1,0.77,6.24,-5.67,1e3,3e4, 01264 8.}; 01265 for( i=1, _r = 0; i <= 7; i++ ) 01266 { 01267 CTRecombData[12][3][i-1] = _itmp50[_r++]; 01268 } 01269 } 01270 for( i=1; i <= 7; i++ ) 01271 { 01272 CTRecombData[13][0][i-1] = 0.f; 01273 } 01274 /* Si+2 recombination */ 01275 { static double _itmp51[] = {6.77 , 7.36e-2 , -0.43 , -0.11 , 5e2 , 1e5 , 01276 2.72}; 01277 for( i=1, _r = 0; i <= 7; i++ ) 01278 { 01279 CTRecombData[13][1][i-1] = _itmp51[_r++]; 01280 } 01281 } 01282 { static double _itmp52[] = {4.90e-1,-8.74e-2,-0.36,-0.79,1e3, 01283 3e4,4.23}; 01284 for( i=1, _r = 0; i <= 7; i++ ) 01285 { 01286 CTRecombData[13][2][i-1] = _itmp52[_r++]; 01287 } 01288 } 01289 { static double _itmp53[] = {7.58,0.37,1.06,-4.09,1e3,5e4,7.49}; 01290 for( i=1, _r = 0; i <= 7; i++ ) 01291 { 01292 CTRecombData[13][3][i-1] = _itmp53[_r++]; 01293 } 01294 } 01295 for( i=1; i <= 7; i++ ) 01296 { 01297 CTRecombData[14][0][i-1] = 0.f; 01298 } 01299 { static double _itmp54[] = {1.74e-4,3.84,36.06,-0.97,1e3,3e4, 01300 3.45}; 01301 for( i=1, _r = 0; i <= 7; i++ ) 01302 { 01303 CTRecombData[14][1][i-1] = _itmp54[_r++]; 01304 } 01305 } 01306 { static double _itmp55[] = {9.46e-2,-5.58e-2,0.77,-6.43,1e3, 01307 3e4,7.29}; 01308 for( i=1, _r = 0; i <= 7; i++ ) 01309 { 01310 CTRecombData[14][2][i-1] = _itmp55[_r++]; 01311 } 01312 } 01313 { static double _itmp56[] = {5.37,0.47,2.21,-8.52,1e3,3e4,9.71}; 01314 for( i=1, _r = 0; i <= 7; i++ ) 01315 { 01316 CTRecombData[14][3][i-1] = _itmp56[_r++]; 01317 } 01318 } 01319 { static double _itmp57[] = {3.82e-7,11.10,2.57e4,-8.22,1e3, 01320 1e4,-3.24}; 01321 for( i=1, _r = 0; i <= 7; i++ ) 01322 { 01323 CTRecombData[15][0][i-1] = _itmp57[_r++]; 01324 } 01325 } 01326 { static double _itmp58[] = {1.00e-5,0.,0.,0.,1e3,3e4,-9.73}; 01327 for( i=1, _r = 0; i <= 7; i++ ) 01328 { 01329 CTRecombData[15][1][i-1] = _itmp58[_r++]; 01330 } 01331 } 01332 { static double _itmp59[] = {2.29,4.02e-2,1.59,-6.06,1e3,3e4, 01333 5.73}; 01334 for( i=1, _r = 0; i <= 7; i++ ) 01335 { 01336 CTRecombData[15][2][i-1] = _itmp59[_r++]; 01337 } 01338 } 01339 { static double _itmp60[] = {6.44,0.13,2.69,-5.69,1e3,3e4,8.60}; 01340 for( i=1, _r = 0; i <= 7; i++ ) 01341 { 01342 CTRecombData[15][3][i-1] = _itmp60[_r++]; 01343 } 01344 } 01345 for( i=1; i <= 7; i++ ) 01346 { 01347 CTRecombData[16][0][i-1] = 0.f; 01348 } 01349 { static double _itmp61[] = {1.00e-5,0.,0.,0.,1e3,3e4,-10.21}; 01350 for( i=1, _r = 0; i <= 7; i++ ) 01351 { 01352 CTRecombData[16][1][i-1] = _itmp61[_r++]; 01353 } 01354 } 01355 { static double _itmp62[] = {1.88,0.32,1.77,-5.70,1e3,3e4,8.}; 01356 for( i=1, _r = 0; i <= 7; i++ ) 01357 { 01358 CTRecombData[16][2][i-1] = _itmp62[_r++]; 01359 } 01360 } 01361 { static double _itmp63[] = {7.27,0.29,1.04,-10.14,1e3,3e4, 01362 9.}; 01363 for( i=1, _r = 0; i <= 7; i++ ) 01364 { 01365 CTRecombData[16][3][i-1] = _itmp63[_r++]; 01366 } 01367 } 01368 for( i=1; i <= 7; i++ ) 01369 { 01370 CTRecombData[17][0][i-1] = 0.f; 01371 } 01372 { static double _itmp64[] = {1.00e-5,0.,0.,0.,1e3,3e4,-14.03}; 01373 for( i=1, _r = 0; i <= 7; i++ ) 01374 { 01375 CTRecombData[17][1][i-1] = _itmp64[_r++]; 01376 } 01377 } 01378 { static double _itmp65[] = {4.57,0.27,-0.18,-1.57,1e3,3e4, 01379 5.73}; 01380 for( i=1, _r = 0; i <= 7; i++ ) 01381 { 01382 CTRecombData[17][2][i-1] = _itmp65[_r++]; 01383 } 01384 } 01385 { static double _itmp66[] = {6.37,0.85,10.21,-6.22,1e3,3e4, 01386 8.60}; 01387 for( i=1, _r = 0; i <= 7; i++ ) 01388 { 01389 CTRecombData[17][3][i-1] = _itmp66[_r++]; 01390 } 01391 } 01392 for( i=1; i <= 7; i++ ) 01393 { 01394 CTRecombData[18][0][i-1] = 0.f; 01395 } 01396 { static double _itmp67[] = {1.00e-5,0.,0.,0.,1e3,3e4,-18.02}; 01397 for( i=1, _r = 0; i <= 7; i++ ) 01398 { 01399 CTRecombData[18][1][i-1] = _itmp67[_r++]; 01400 } 01401 } 01402 { static double _itmp68[] = {4.76,0.44,-0.56,-0.88,1e3,3e4, 01403 6.}; 01404 for( i=1, _r = 0; i <= 7; i++ ) 01405 { 01406 CTRecombData[18][2][i-1] = _itmp68[_r++]; 01407 } 01408 } 01409 { static double _itmp69[] = {1.00e-5,0.,0.,0.,1e3,3e4,-47.3}; 01410 for( i=1, _r = 0; i <= 7; i++ ) 01411 { 01412 CTRecombData[18][3][i-1] = _itmp69[_r++]; 01413 } 01414 } 01415 for( i=1; i <= 7; i++ ) 01416 { 01417 CTRecombData[19][0][i-1] = 0.f; 01418 } 01419 { static double _itmp70[] = {0.,0.,0.,0.,1e1,1e9,0.}; 01420 for( i=1, _r = 0; i <= 7; i++ ) 01421 { 01422 CTRecombData[19][1][i-1] = _itmp70[_r++]; 01423 } 01424 } 01425 { static double _itmp71[] = {3.17e-2,2.12,12.06,-0.40,1e3,3e4, 01426 6.6}; 01427 for( i=1, _r = 0; i <= 7; i++ ) 01428 { 01429 CTRecombData[19][2][i-1] = _itmp71[_r++]; 01430 } 01431 } 01432 { static double _itmp72[] = {2.68,0.69,-0.68,-4.47,1e3,3e4, 01433 9.9}; 01434 for( i=1, _r = 0; i <= 7; i++ ) 01435 { 01436 CTRecombData[19][3][i-1] = _itmp72[_r++]; 01437 } 01438 } 01439 for( i=1; i <= 7; i++ ) 01440 { 01441 CTRecombData[20][0][i-1] = 0.f; 01442 } 01443 { static double _itmp73[] = {0.,0.,0.,0.,1e1,1e9,0.}; 01444 for( i=1, _r = 0; i <= 7; i++ ) 01445 { 01446 CTRecombData[20][1][i-1] = _itmp73[_r++]; 01447 } 01448 } 01449 { static double _itmp74[] = {7.22e-3,2.34,411.50,-13.24,1e3, 01450 3e4,3.5}; 01451 for( i=1, _r = 0; i <= 7; i++ ) 01452 { 01453 CTRecombData[20][2][i-1] = _itmp74[_r++]; 01454 } 01455 } 01456 { static double _itmp75[] = {1.20e-1,1.48,4.00,-9.33,1e3,3e4, 01457 10.61}; 01458 for( i=1, _r = 0; i <= 7; i++ ) 01459 { 01460 CTRecombData[20][3][i-1] = _itmp75[_r++]; 01461 } 01462 } 01463 for( i=1; i <= 7; i++ ) 01464 { 01465 CTRecombData[21][0][i-1] = 0.f; 01466 } 01467 { static double _itmp76[] = {0.,0.,0.,0.,1e1,1e9,0.}; 01468 for( i=1, _r = 0; i <= 7; i++ ) 01469 { 01470 CTRecombData[21][1][i-1] = _itmp76[_r++]; 01471 } 01472 } 01473 { static double _itmp77[] = {6.34e-1,6.87e-3,0.18,-8.04,1e3, 01474 3e4,4.3}; 01475 for( i=1, _r = 0; i <= 7; i++ ) 01476 { 01477 CTRecombData[21][2][i-1] = _itmp77[_r++]; 01478 } 01479 } 01480 { static double _itmp78[] = {4.37e-3,1.25,40.02,-8.05,1e3,3e4, 01481 5.3}; 01482 for( i=1, _r = 0; i <= 7; i++ ) 01483 { 01484 CTRecombData[21][3][i-1] = _itmp78[_r++]; 01485 } 01486 } 01487 for( i=1; i <= 7; i++ ) 01488 { 01489 CTRecombData[22][0][i-1] = 0.f; 01490 } 01491 { static double _itmp79[] = {1.00e-5,0.,0.,0.,1e3,3e4,-1.05}; 01492 for( i=1, _r = 0; i <= 7; i++ ) 01493 { 01494 CTRecombData[22][1][i-1] = _itmp79[_r++]; 01495 } 01496 } 01497 { static double _itmp80[] = {5.12,-2.18e-2,-0.24,-0.83,1e3, 01498 3e4,4.7}; 01499 for( i=1, _r = 0; i <= 7; i++ ) 01500 { 01501 CTRecombData[22][2][i-1] = _itmp80[_r++]; 01502 } 01503 } 01504 { static double _itmp81[] = {1.96e-1,-8.53e-3,0.28,-6.46,1e3, 01505 3e4,6.2}; 01506 for( i=1, _r = 0; i <= 7; i++ ) 01507 { 01508 CTRecombData[22][3][i-1] = _itmp81[_r++]; 01509 } 01510 } 01511 for( i=1; i <= 7; i++ ) 01512 { 01513 CTRecombData[23][0][i-1] = 0.f; 01514 } 01515 { static double _itmp82[] = {5.27e-1,0.61,-0.89,-3.56,1e3,3e4, 01516 2.89}; 01517 for( i=1, _r = 0; i <= 7; i++ ) 01518 { 01519 CTRecombData[23][1][i-1] = _itmp82[_r++]; 01520 } 01521 } 01522 { static double _itmp83[] = {10.90,0.24,0.26,-11.94,1e3,3e4, 01523 5.4}; 01524 for( i=1, _r = 0; i <= 7; i++ ) 01525 { 01526 CTRecombData[23][2][i-1] = _itmp83[_r++]; 01527 } 01528 } 01529 { static double _itmp84[] = {1.18,0.20,0.77,-7.09,1e3,3e4,6.6}; 01530 for( i=1, _r = 0; i <= 7; i++ ) 01531 { 01532 CTRecombData[23][3][i-1] = _itmp84[_r++]; 01533 } 01534 } 01535 for( i=1; i <= 7; i++ ) 01536 { 01537 CTRecombData[24][0][i-1] = 0.f; 01538 } 01539 { static double _itmp85[] = {1.65e-1,6.80e-3,6.44e-2,-9.70, 01540 1e3,3e4,2.04}; 01541 for( i=1, _r = 0; i <= 7; i++ ) 01542 { 01543 CTRecombData[24][1][i-1] = _itmp85[_r++]; 01544 } 01545 } 01546 { static double _itmp86[] = {14.20,0.34,-0.41,-1.19,1e3,3e4, 01547 6.}; 01548 for( i=1, _r = 0; i <= 7; i++ ) 01549 { 01550 CTRecombData[24][2][i-1] = _itmp86[_r++]; 01551 } 01552 } 01553 { static double _itmp87[] = {4.43e-1,0.91,10.76,-7.49,1e3,3e4, 01554 7.}; 01555 for( i=1, _r = 0; i <= 7; i++ ) 01556 { 01557 CTRecombData[24][3][i-1] = _itmp87[_r++]; 01558 } 01559 } 01560 for( i=1; i <= 7; i++ ) 01561 { 01562 CTRecombData[25][0][i-1] = 0.f; 01563 } 01564 { static double _itmp88[] = {1.26,7.72e-2,-0.41,-7.31,1e3,1e5, 01565 2.56}; 01566 for( i=1, _r = 0; i <= 7; i++ ) 01567 { 01568 CTRecombData[25][1][i-1] = _itmp88[_r++]; 01569 } 01570 } 01571 { static double _itmp89[] = {3.42,0.51,-2.06,-8.99,1e3,1e5, 01572 6.3}; 01573 for( i=1, _r = 0; i <= 7; i++ ) 01574 { 01575 CTRecombData[25][2][i-1] = _itmp89[_r++]; 01576 } 01577 } 01578 { static double _itmp90[] = {14.60,3.57e-2,-0.92,-0.37,1e3, 01579 3e4,10.}; 01580 for( i=1, _r = 0; i <= 7; i++ ) 01581 { 01582 CTRecombData[25][3][i-1] = _itmp90[_r++]; 01583 } 01584 } 01585 for( i=1; i <= 7; i++ ) 01586 { 01587 CTRecombData[26][0][i-1] = 0.f; 01588 } 01589 { static double _itmp91[] = {5.30,0.24,-0.91,-0.47,1e3,3e4, 01590 2.9}; 01591 for( i=1, _r = 0; i <= 7; i++ ) 01592 { 01593 CTRecombData[26][1][i-1] = _itmp91[_r++]; 01594 } 01595 } 01596 { static double _itmp92[] = {3.26,0.87,2.85,-9.23,1e3,3e4,6.}; 01597 for( i=1, _r = 0; i <= 7; i++ ) 01598 { 01599 CTRecombData[26][2][i-1] = _itmp92[_r++]; 01600 } 01601 } 01602 { static double _itmp93[] = {1.03,0.58,-0.89,-0.66,1e3,3e4, 01603 10.51}; 01604 for( i=1, _r = 0; i <= 7; i++ ) 01605 { 01606 CTRecombData[26][3][i-1] = _itmp93[_r++]; 01607 } 01608 } 01609 for( i=1; i <= 7; i++ ) 01610 { 01611 CTRecombData[27][0][i-1] = 0.f; 01612 } 01613 { static double _itmp94[] = {1.05,1.28,6.54,-1.81,1e3,1e5,3.0}; 01614 for( i=1, _r = 0; i <= 7; i++ ) 01615 { 01616 CTRecombData[27][1][i-1] = _itmp94[_r++]; 01617 } 01618 } 01619 { static double _itmp95[] = {9.73,0.35,0.90,-5.33,1e3,3e4,5.2}; 01620 for( i=1, _r = 0; i <= 7; i++ ) 01621 { 01622 CTRecombData[27][2][i-1] = _itmp95[_r++]; 01623 } 01624 } 01625 { static double _itmp96[] = {6.14,0.25,-0.91,-0.42,1e3,3e4, 01626 10.}; 01627 for( i=1, _r = 0; i <= 7; i++ ) 01628 { 01629 CTRecombData[27][3][i-1] = _itmp96[_r++]; 01630 } 01631 } 01632 for( i=1; i <= 7; i++ ) 01633 { 01634 CTRecombData[28][0][i-1] = 0.f; 01635 } 01636 { static double _itmp97[] = {1.47e-3,3.51,23.91,-0.93,1e3,3e4, 01637 3.44}; 01638 for( i=1, _r = 0; i <= 7; i++ ) 01639 { 01640 CTRecombData[28][1][i-1] = _itmp97[_r++]; 01641 } 01642 } 01643 { static double _itmp98[] = {9.26,0.37,0.40,-10.73,1e3,3e4, 01644 5.6}; 01645 for( i=1, _r = 0; i <= 7; i++ ) 01646 { 01647 CTRecombData[28][2][i-1] = _itmp98[_r++]; 01648 } 01649 } 01650 { static double _itmp99[] = {11.59,0.20,0.80,-6.62,1e3,3e4, 01651 9.}; 01652 for( i=1, _r = 0; i <= 7; i++ ) 01653 { 01654 CTRecombData[28][3][i-1] = _itmp99[_r++]; 01655 } 01656 } 01657 for( i=1; i <= 7; i++ ) 01658 { 01659 CTRecombData[29][0][i-1] = 0.f; 01660 } 01661 { static double _itmp100[] = {1.00e-5,0.,0.,0.,1e3,3e4,-4.37}; 01662 for( i=1, _r = 0; i <= 7; i++ ) 01663 { 01664 CTRecombData[29][1][i-1] = _itmp100[_r++]; 01665 } 01666 } 01667 { static double _itmp101[] = {6.96e-4,4.24,26.06,-1.24,1e3, 01668 3e4,7.8}; 01669 for( i=1, _r = 0; i <= 7; i++ ) 01670 { 01671 CTRecombData[29][2][i-1] = _itmp101[_r++]; 01672 } 01673 } 01674 { static double _itmp102[] = {1.33e-2,1.56,-0.92,-1.20,1e3, 01675 3e4,11.73}; 01676 for( i=1, _r = 0; i <= 7; i++ ) 01677 { 01678 CTRecombData[29][3][i-1] = _itmp102[_r++]; 01679 } 01680 } 01681 } 01682 01683 /* ChargTranPun, punch charge transfer coefficient */ 01684 void ChargTranPun( FILE* ipPnunit , char* chPunch ) 01685 { 01686 long int j, jj; 01687 /* restore temp when done with this routine */ 01688 double TempSave = phycon.te; 01689 01690 DEBUG_ENTRY( "ChargTranPun()" ); 01691 01692 /* this is the usual charge transfer option */ 01693 if( strcmp( chPunch,"CHAR") == 0 ) 01694 { 01695 /* charge exchange rate coefficients, entered with 01696 * punch charge transfer command. Queries Jim Kingdon's 01697 * CT fits and routines to get H - He and higher CT rates 01698 * rates are evaluated at the current temperature, which can 01699 * be specified with the constant temperature command */ 01700 /* first group is charge transfer recombination, 01701 * the process A+x + H => A+x-1 + H+ */ 01702 fprintf( ipPnunit, "#element\tion\n"); 01703 for( j=1; j < LIMELM; j++ ) 01704 { 01705 /* first number is atomic number, so add 1 to get onto physical scale */ 01706 /* >>chng 00 may 09, caught by Jon Slavin */ 01707 fprintf( ipPnunit, "%s\t", elementnames.chElementSym[j] ); 01708 /*fprintf( ipPnunit, "%3ld\t", j+1 );*/ 01709 01710 for( jj=0; jj < j; jj++ ) 01711 { 01712 fprintf( ipPnunit, "%.2e\t", 01713 HCTRecom(jj,j) ); 01714 } 01715 fprintf( ipPnunit, "\n" ); 01716 } 01717 01718 /* second group is charge transfer ionization, 01719 * the process A+x + H+ => A+x+1 + H */ 01720 fprintf( ipPnunit, "\n#ionization rates, atomic number\n"); 01721 for( j=1; j < LIMELM; j++ ) 01722 { 01723 fprintf( ipPnunit, "%s\t", elementnames.chElementSym[j] ); 01724 for( jj=0; jj < j; jj++ ) 01725 { 01726 fprintf( ipPnunit, "%.2e\t", 01727 HCTIon(jj,j) ); 01728 } 01729 fprintf( ipPnunit, "\n" ); 01730 } 01731 } 01732 01733 /* this is the charge transfer to produce output for AGN3, 01734 * invoked with the punch charge transfer AGN command */ 01735 else if( strcmp( chPunch,"CHAG") == 0 ) 01736 { 01737 /* this is boundary between two tables */ 01738 double BreakEnergy = 100./13.0; 01739 double OHRate; 01740 realnum teinit = 5000.; 01741 realnum tefinal = 20000.; 01742 realnum te1; 01743 long int nelem, ion; 01744 01745 te1 = teinit; 01746 fprintf(ipPnunit,"H ioniz\n X+i\\Te"); 01747 while( te1 <= tefinal ) 01748 { 01749 fprintf(ipPnunit,"\t%.0f K",te1); 01750 te1 *= 2.; 01751 } 01752 fprintf(ipPnunit,"\n"); 01753 01754 /* make sure rates already evaluated at least one time */ 01755 ChargTranEval(&OHRate); 01756 01757 /* loop over all elements, H charge transfer ionization */ 01758 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem ) 01759 { 01760 /* this list of elements included in the AGN tables is defined in zeroabun.c */ 01761 if( abund.lgAGN[nelem] ) 01762 { 01763 for( ion=0; ion<=nelem; ++ion ) 01764 { 01765 /* skip high ionization CT */ 01766 if( Heavy.Valence_IP_Ryd[nelem][ion] > BreakEnergy ) 01767 break; 01768 /* most of these are actually zero */ 01769 if( atmdat.HCharExcIonOf[nelem][ion] == 0 ) 01770 continue; 01771 01772 /* print chemical symbol */ 01773 fprintf(ipPnunit,"%s", 01774 elementnames.chElementSym[nelem]); 01775 /* now ionization stage */ 01776 if( ion==0 ) 01777 { 01778 fprintf(ipPnunit,"0 "); 01779 } 01780 else if( ion==1 ) 01781 { 01782 fprintf(ipPnunit,"+ "); 01783 } 01784 else 01785 { 01786 fprintf(ipPnunit,"+%li",ion); 01787 } 01788 01789 /* fully define the new temperature */ 01790 TempChange(teinit , false); 01791 01792 while( phycon.te <= tefinal ) 01793 { 01794 dense.IonLow[nelem] = 0; 01795 dense.IonHigh[nelem] = nelem+1; 01796 ChargTranEval(&OHRate); 01797 01798 fprintf(ipPnunit,"\t%.2e",atmdat.HCharExcIonOf[nelem][ion]); 01799 TempChange(phycon.te *2.f , false); 01800 } 01801 fprintf(ipPnunit,"\n"); 01802 } 01803 fprintf(ipPnunit,"\n"); 01804 } 01805 } 01806 01807 te1 = teinit; 01808 fprintf(ipPnunit,"H recom\n X+i\\Te"); 01809 while( te1 <= tefinal ) 01810 { 01811 fprintf(ipPnunit,"\t%.0f K",te1); 01812 te1 *= 2.; 01813 } 01814 fprintf(ipPnunit,"\n"); 01815 01816 /* loop over all elements, H charge transfer recombination */ 01817 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem ) 01818 { 01819 /* this list of elements included in the AGN tables is defined in zeroabun.c */ 01820 if( abund.lgAGN[nelem] ) 01821 { 01822 for( ion=0; ion<=nelem; ++ion ) 01823 { 01824 /* skip high ionization CT */ 01825 if( Heavy.Valence_IP_Ryd[nelem][ion] > BreakEnergy ) 01826 break; 01827 /* most of these are actually zero */ 01828 if( atmdat.HCharExcRecTo[nelem][ion] == 0 ) 01829 continue; 01830 01831 /* print chemical symbol */ 01832 fprintf(ipPnunit,"%s", 01833 elementnames.chElementSym[nelem]); 01834 /* now ionization stage */ 01835 if( ion==0 ) 01836 { 01837 fprintf(ipPnunit,"0 "); 01838 } 01839 else if( ion==1 ) 01840 { 01841 fprintf(ipPnunit,"+ "); 01842 } 01843 else 01844 { 01845 fprintf(ipPnunit,"+%li",ion); 01846 } 01847 01848 /* fully define the new temperature */ 01849 TempChange(teinit , false); 01850 while( phycon.te <= tefinal ) 01851 { 01852 dense.IonLow[nelem] = 0; 01853 dense.IonHigh[nelem] = nelem+1; 01854 ChargTranEval(&OHRate); 01855 01856 fprintf(ipPnunit,"\t%.2e",atmdat.HCharExcRecTo[nelem][ion]); 01857 TempChange(phycon.te *2.f , false); 01858 } 01859 fprintf(ipPnunit,"\n"); 01860 } 01861 fprintf(ipPnunit,"\n"); 01862 } 01863 } 01864 01865 # if 0 01866 te1 = teinit; 01867 fprintf(ipPnunit,"He recom\n Elem\\Te"); 01868 while( te1 <= tefinal ) 01869 { 01870 fprintf(ipPnunit,"\t%.0f",te1); 01871 te1 *= 2.; 01872 } 01873 fprintf(ipPnunit,"\n"); 01874 01875 /* loop over all elements, H charge transfer recombination */ 01876 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem ) 01877 { 01878 /* this list of elements included in the AGN tables is defined in zeroabun.c */ 01879 if( abund.lgAGN[nelem] ) 01880 { 01881 for( ion=0; ion<=nelem; ++ion ) 01882 { 01883 /* most of these are actually zero */ 01884 if( atmdat.HeCharExcRecTo[nelem][ion] == 0 ) 01885 continue; 01886 fprintf(ipPnunit,"%.2s%.2s", 01887 elementnames.chElementSym[nelem], 01888 elementnames.chIonStage[ion]); 01889 01890 /* fully define the new temperature */ 01891 TempChange(teinit , false); 01892 while( phycon.te <= tefinal ) 01893 { 01894 dense.IonLow[nelem] = 0; 01895 dense.IonHigh[nelem] = nelem+1; 01896 ChargTranEval(); 01897 01898 fprintf(ipPnunit,"\t%.2e",atmdat.HeCharExcRecTo[nelem][ion]); 01899 TempChange(phycon.te *2.fprintf , false); 01900 } 01901 fprintf(ipPnunit,"\n"); 01902 } 01903 fprintf(ipPnunit,"\n"); 01904 } 01905 } 01906 01907 01908 te1 = teinit; 01909 fprintf(ipPnunit,"He ioniz\n Elem\\Te"); 01910 while( te1 <= tefinal ) 01911 { 01912 fprintf(ipPnunit,"\t%.0f",te1); 01913 te1 *= 2.; 01914 } 01915 fprintf(ipPnunit,"\n"); 01916 01917 /* loop over all elements, H charge transfer recombination */ 01918 for( nelem=ipHELIUM; nelem<LIMELM; ++nelem ) 01919 { 01920 /* this list of elements included in the AGN tables is defined in zeroabun.c */ 01921 if( abund.lgAGN[nelem] ) 01922 { 01923 for( ion=0; ion<=nelem; ++ion ) 01924 { 01925 /* most of these are actually zero */ 01926 if( atmdat.HeCharExcIonOf[nelem][ion] == 0 ) 01927 continue; 01928 fprintf(ipPnunit,"%.2s%.2s", 01929 elementnames.chElementSym[nelem], 01930 elementnames.chIonStage[ion]); 01931 01932 /* fully define the new temperature */ 01933 TempChange(teinit , false); 01934 while( phycon.te <= tefinal ) 01935 { 01936 dense.IonLow[nelem] = 0; 01937 dense.IonHigh[nelem] = nelem+1; 01938 ChargTranEval(); 01939 01940 fprintf(ipPnunit,"\t%.2e",atmdat.HeCharExcIonOf[nelem][ion]); 01941 TempChange(phycon.te*2.f , true); 01942 } 01943 fprintf(ipPnunit,"\n"); 01944 } 01945 fprintf(ipPnunit,"\n"); 01946 } 01947 } 01948 # endif 01949 } 01950 else 01951 { 01952 fprintf( ioQQQ, " punch charge keyword insane\n" ); 01953 cdEXIT(EXIT_FAILURE); 01954 } 01955 01956 TempChange(TempSave , false); 01957 return; 01958 }