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 /*CoolOxyg evaluate total cooling due to oxygen */ 00004 #include "cddefines.h" 00005 #include "coolheavy.h" 00006 #include "dense.h" 00007 #include "taulines.h" 00008 #include "h2.h" 00009 #include "phycon.h" 00010 #include "embesq.h" 00011 #include "hmi.h" 00012 #include "oxy.h" 00013 #include "colden.h" 00014 #include "mole.h" 00015 #include "ligbar.h" 00016 #include "thermal.h" 00017 #include "lines_service.h" 00018 #include "atoms.h" 00019 #include "cooling.h" 00020 /* oi3pcs make collision strengths with electrons for 3P O I ground term. 00021 *>>refer o1 cs Berrington, K.A. 1988, J.Phys.B, 21, 1083 for Te > 3000K 00022 *>>refer o1 cs Bell, Berrington & Thomas 1998, MNRAS 293, L83 for 50K <= Te <= 3000K 00023 * numbers in cs variables refer to j value: j=0,1,2 (n=3,2,1 in 5 level atom) 00024 * written by Kirk Korista, 29 may 96 00025 * adapted by Peter van Hoof, 30 march 99 (to include Bell et al. data) 00026 * all data above 3000K have been adjusted down by a constant factor to make 00027 * them line up with Bell et al. data. */ 00028 STATIC void oi3Pcs(double *cs21, 00029 double *cs20, 00030 double *cs10); 00031 00032 #if defined (__ICC) && defined(__ia64) && __INTEL_COMPILER < 910 00033 #pragma optimization_level 1 00034 #endif 00035 void CoolOxyg(void) 00036 { 00037 double a21, 00038 a31, 00039 a32, 00040 a41, 00041 a42, 00042 a43, 00043 a51, 00044 a52, 00045 a53, 00046 a54, 00047 a6300, 00048 a6363, 00049 aeff, 00050 cs, 00051 cs2s2p, 00052 cs2s3p , 00053 cs01, 00054 cs02, 00055 cs12, 00056 cs13, 00057 cs21, 00058 cs23, 00059 cs31, 00060 cs32, 00061 cs41, 00062 cs42, 00063 cs43, 00064 cs51, 00065 cs52, 00066 cs53, 00067 cs54, 00068 Te_bounded, 00069 Te_bounded_log , 00070 o3cs23, 00071 p[5], 00072 r12 , 00073 r13, 00074 pump_rate; 00075 /*double cs_ab ,cs_ac;(Commented out- Humeshkar Nemala)*/ 00076 /* these were added by Peter van Hoof to update the collision 00077 * data within the OI ground term */ 00078 double cse01=-1.,cse12=-1.,cse02 =-1., 00079 csh01=-1.,cshe01=-1.,csh201=-1.,csh12=-1.,cshe12=-1.,csh212=-1.,csh02=-1.,cshe02=-1.,csh202 =-1., 00080 csh2o01=-1.,csh2o02=-1.,csh2o12=-1.,csh2p01=-1.,csh2p02=-1.,csh2p12=-1.,csp01=-1.,csp02=-1., 00081 csp12=-1.; 00082 static double go2[5]={4.,6.,4.,4.,2.}; 00083 static double exo2[4]={26808.4,21.0,13637.5,1.5}; 00084 /* these will be used to update change in cs wrt temperature, 00085 * and its affect on cooling derivative */ 00086 static double oi_cs_tsave=-1. , oi_te_tsave=-1. , oi_dcdt_tsave=-1.; 00087 long int i; 00088 static bool lgFirst=true; 00089 static long int *ipO4Pump=NULL, 00090 nO4Pump=0; 00091 double rate_OH_dissoc; 00092 00093 DEBUG_ENTRY( "CoolOxyg()" ); 00094 00095 /* following does the OI Bowen Ly-bet pumping problem */ 00096 atom_oi_calc(&CoolHeavy.coolOi); 00097 CoolAdd("o 1",8446,CoolHeavy.coolOi); 00098 00099 /* [OI] 6300, 6363, 5575, etc 00100 * >>chng 06 oct 02, Humeshkar Nemala incorporate Barklem cs data for OI 00101 * largest difference is 3P - 1D (6300+6363) which is now roughly 3x smaller */ 00102 /* This is the transition from 3P(J=2,1,0;the levels are reversed) to 1D2 00103 * The rate coefficients were converted to CS 00104 *>>refer oi cs Barklem,P.S.,2006,A&A (astroph 0609684) 00105 * Data pts are avilable over 1000,3000,5000,8000,12000,20000 and 50000 00106 * Fits between 1000K and 20000K are reliable;this is not the case 00107 * between 20000K & 50000K*/ 00108 /* this was old calculation of cs12 that was replaced by Humeshkar */ 00109 #if 0 00110 cs12 = 2.68e-5*phycon.te*(1. + 1.67e-6*phycon.te - 2.95e-10*phycon.te* 00111 phycon.te); 00112 cs12 = MAX2(0.1,cs12); 00113 #endif 00114 00115 /* [OI] electron 6300, 6363 collision strength 00116 * >>chng 06 oct 02, Humeshkar Nemala code, based on Barklem paper */ 00117 if(phycon.te < 8E3) 00118 { 00119 cs12 = (4E-08)*(phycon.te*phycon.te70*phycon.te05); 00120 } 00121 else if(phycon.te < 2E4) 00122 { 00123 cs12 = (4.630155E-05)*((phycon.te/phycon.te04)*phycon.te005*phycon.te0001); 00124 } 00125 else 00126 { 00127 /* this branch T > 20,000K */ 00128 cs12 = (1.5394E-03)*(phycon.sqrte*phycon.te10*phycon.te01*phycon.te001*phycon.te0003); 00129 } 00130 00131 /* this block adds on collisional excitation by H0 */ 00132 { 00133 /* >>chng 06 aug 18, add atomic hydrogen collisional processes using rates from 00134 * >>refer O1 coll Krems, R.V., Jameson, M.J. & Dalgarno, A. 2006, ApJ, 647, 1531 00135 * their equation 10 - the deecxiation rate coefficient, cm3 s-1 00136 * >>referold oi cs Federman, S.R., & Shipsey, E.J. 1983, ApJ, 269, 791 */ 00137 double te_scale = phycon.te / 6000.; 00138 double rate_H0 = (1.74*te_scale + 0.6)*1e-12*sexp(0.47*te_scale) / sqrt(te_scale ) * 00139 dense.xIonDense[ipHYDROGEN][0]; 00140 /*fprintf(ioQQQ,"DEBUG OI H0 rate %.2e " , cs12 );*/ 00141 cs12 += ConvRate2CS( 5.f , (realnum)rate_H0 ); 00142 /*fprintf(ioQQQ,"%.2e\n" , cs12 );*/ 00143 } 00144 /* following was old fit to 2-3 5755 transition */ 00145 /*cs23 = 1.05e-3*phycon.sqrte;*/ 00146 /*>>chng 06 oct 02, Humeshkar Nemala's new fit to Barklem paper */ 00147 if(phycon.te < 5E3) 00148 { 00149 cs23 = (7E-08)*(phycon.te*phycon.sqrte*phycon.te10*phycon.te007*phycon.te0001); 00150 } 00151 else if(phycon.te<2E4) 00152 { 00153 cs23 = (1.98479e-04)*((phycon.te70/phycon.te03)*phycon.te003*phycon.te0007); 00154 } 00155 else 00156 { 00157 cs23 = (9.31E-04)*(phycon.sqrte*phycon.te01*phycon.te007*phycon.te0005*phycon.te0001); 00158 } 00159 00160 /* following was old fit to 1,2,3 -> 5 transition */ 00161 /*cs13 = 3.24e-6*phycon.te;*/ 00162 /* following is Humeshkar's new fit to Barklem paper */ 00163 if(phycon.te < 2E4) 00164 { 00165 cs13 = (2E-05)*(phycon.sqrte*phycon.te30*phycon.te05*phycon.te01*(phycon.te004/phycon.te0002)); 00166 } 00167 else 00168 { 00169 cs13 = (1.054E-03)*(phycon.sqrte/phycon.te04)*phycon.te003*phycon.te0005; 00170 } 00171 00172 /* O I 6300, 6363, A from 00173 * >>refer all all Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103, 00174 * >>refercon ed by D.R. Flower, (D. Reidel: Holland), 143 */ 00175 a6300 = TauLines[ipT6300].Emis->Aul*TauLines[ipT6300].Emis->Pesc; 00176 TauLines[ipT6300].Emis->PopOpc = (dense.xIonDense[ipOXYGEN][0]*5./5.); 00177 TauLines[ipT6300].Lo->Pop = (dense.xIonDense[ipOXYGEN][0]*5./5.); 00178 TauLines[ipT6300].Hi->Pop = 0.; 00179 TauLines[ipT6300].Coll.cs = (realnum)(cs12*5./9.); 00180 TauLines[ipT6363].Emis->PopOpc = (dense.xIonDense[ipOXYGEN][0]*5./3.); 00181 TauLines[ipT6363].Lo->Pop = (dense.xIonDense[ipOXYGEN][0]*5./3.); 00182 TauLines[ipT6363].Hi->Pop = 0.; 00183 TauLines[ipT6363].Coll.cs = (realnum)(cs12*3./9.); 00184 a6363 = TauLines[ipT6363].Emis->Aul*TauLines[ipT6363].Emis->Pesc; 00185 00186 a21 = a6300 + a6363 + oxy.d6300; 00187 a32 = TauLines[ipT5577].Emis->Aul*TauLines[ipT5577].Emis->Pesc; 00188 00189 PutCS(cs23,&TauLines[ipT5577]); 00190 00191 /* rate of new populations of O^0 due to dissociation of OH, 00192 * co.rate_OH_dissoc is rate OH -> O + H [cm-3 s-1], 00193 * must make it per unit O atom, so this rate is s-1 excitations per O atom */ 00194 rate_OH_dissoc = CO_findrate("PHOTON,OH=>O,H"); 00195 r12 = rate_OH_dissoc*0.55/SDIV( dense.xIonDense[ipOXYGEN][0] ); 00196 r13 = rate_OH_dissoc*0.05/SDIV( dense.xIonDense[ipOXYGEN][0] ); 00197 00198 /* below is correction for fraction of excitations the produce emission */ 00199 CoolHeavy.c6300_frac_emit = (a6300+a6363)/(a6300+a6363+cs12*dense.cdsqte/5.); 00200 CoolHeavy.c5577_frac_emit = (a32)/(a32+cs23*dense.cdsqte/3.); 00201 00202 /* d6300 is the photoionization rate from the excited level 00203 * was computed when ionization balance done */ 00204 CoolHeavy.c5577 = atom_pop3(9.,5.,1.,cs12,cs13,cs23,a21,7.56e-2,a32, 00205 22590.,25775.7,&oxy.poiexc,dense.xIonDense[ipOXYGEN][0],0.,r12,r13)*a32* 00206 3.57e-12; 00207 00208 TauLines[ipT5577].Emis->PopOpc = oxy.poiexc; 00209 TauLines[ipT5577].Lo->Pop = oxy.poiexc; 00210 TauLines[ipT5577].Hi->Pop = 0.; 00211 00212 /* convert poiexc to relative abundances */ 00213 /* >>chng 04 apr 20, include correction for fraction of 6300 due to OH pump */ 00214 CoolHeavy.c5577 *= (1.-r13/(SDIV(atoms.c13))); 00215 00216 CoolHeavy.c6300 = oxy.poiexc*a6300*TauLines[ipT6300].EnergyErg * 00217 (1.-r12/(SDIV(atoms.c12))); 00218 CoolHeavy.c6363 = oxy.poiexc*a6363*TauLines[ipT6363].EnergyErg * 00219 (1.-r12/(SDIV(atoms.c12))); 00220 /* must introduce correction for fraction of 6300 that is photo produced */ 00221 thermal.dCooldT += CoolHeavy.c6300*(2.28e4*thermal.tsq1 + thermal.halfte) * 00222 /* note that atoms.c12 has all ra tes 1->2 */ 00223 (1.-r12/(SDIV(atoms.c12))); 00224 00225 oxy.poiexc /= (realnum)MAX2(1e-20,dense.xIonDense[ipOXYGEN][0]); 00226 CoolAdd("o 1",5577,CoolHeavy.c5577); 00227 CoolAdd("o 1",6300,CoolHeavy.c6300); 00228 CoolAdd("o 1",6363,CoolHeavy.c6363); 00229 00230 /* O I fine structure lines rad data from 00231 * >>refer all cs Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103, 00232 * >>refercon ed by D.R. Flower, (D. Reidel: Holland), 143 00233 * >>refer o1 as Berrington, K.A. 1988, J.Phys. B, 21, 1083 00234 * hydrogen collisions from 00235 * >>refer oi cs Tielens, A.G.G., & Hollenbach, D. 1985, ApJ, 291, 722 00236 * factor in () is their rate coef 00237 * assume H2 and H are same 00238 * CDSQTE = 8.629E-6*EDEN/SQRTE 00239 * cs01 = 9.8e-6*te + (4.2e-12*te70/te03) / cdsqte * 3. * hdcor 00240 * cs12 = 2.6e-6*te + (1.5e-10*sqrte/te03/te03)/cdsqte*hdcor 00241 * cs02 = 2.9e-6*te + (1.1e-12*te70*te10)/cdsqte*hdcor 00242 * evaluate fits to OI electron rates, indices on var do not agree 00243 * with Kirk's in sub, but are OK */ 00244 /*==============================================================*/ 00245 /* >>>chng 99 jun 01, 00246 * following changed to be parallel to Peter van Hoof's changes 00247 * in the Fortran C90.05 version 00248 * following is collisions with electrons */ 00249 oi3Pcs(&cse01,&cse02,&cse12); 00250 /* remember the electron part of the cs */ 00251 cs = cse01; 00252 00253 /* rate coefficients for collisional de-excitation of O^o(3P) with H^o(2S1/2) 00254 * NOTE: due to lack of data these relations are extrapolated to higher Te ! 00255 * >>refer o1 cs Launay & Roueff 1977, AA 56, 289 00256 * the first fit is for Te <= 300K, the second for higher temps 00257 * these data are valid for 50K <= Te <= 1000K*/ 00258 csh12 = MAX2(2.00e-11*phycon.te30*phycon.te05*phycon.te02, 00259 7.67e-12*phycon.sqrte*phycon.te03); 00260 00261 /* these data are valid for 100K <= Te <= 1000K */ 00262 csh01 = MIN2(3.12e-12*phycon.te70*phycon.te02*phycon.te02, 00263 7.51e-12*phycon.sqrte*phycon.te05*phycon.te03); 00264 00265 /* these data are valid for 150K <= Te <= 1000K*/ 00266 csh02 = MIN2(6.96e-13*phycon.te/phycon.te10/phycon.te02, 00267 1.49e-12*phycon.te70*phycon.te05); 00268 00269 /* rate coefficients for collisional de-excitation of O^o(3P) with He^o(1S) 00270 * NOTE: due to lack of data these relations are extrapolated to higher Te ! 00271 * >>refer oi cs Monteiro & Flower 1987, MNRAS 228, 101 00272 * the first fit is for Te <= 300K, the second for higher temps 00273 * these data are valid for 100K <= Te <= 1000K */ 00274 cshe12 = MIN2(8.09e-16*phycon.te32*phycon.te10*phycon.te03, 00275 9.72e-15*phycon.te*phycon.te20); 00276 00277 cshe01 = 1.57e-12*phycon.te70/phycon.te01; 00278 00279 cshe02 = MIN2(1.80e-12*phycon.te70*phycon.te05, 00280 4.45e-12*phycon.te70/phycon.te10); 00281 00282 if( phycon.te<=1.5e3 ) 00283 { 00284 /* >>chng 04 mar 15, use explicit ortho-para densities */ 00285 double ortho_frac = h2.ortho_density/SDIV(hmi.H2_total); 00286 /* rate coefficients for collisional de-excitation of O^o(3P) with H2(J=1,0) 00287 * >>refer oi cs Jaquet et al. 1992, J.Phys.B 25, 285 00288 * these data are valid for 40K <= Te <= 1500K 00289 * the first entry is contribution from ortho H2, the second para H2.*/ 00290 csh2o12 = 2.21e-14*phycon.te*phycon.te10/phycon.te01; 00291 csh2p12 = 9.45e-15*phycon.te*phycon.te20/phycon.te01; 00292 csh212 = ortho_frac*csh2o12 + (1.-ortho_frac)*csh2p12; 00293 00294 csh2o01 = 2.37e-11*phycon.te30*phycon.te10/phycon.te02; 00295 csh2p01 = 3.40e-11*phycon.te30*phycon.te02; 00296 csh201 = ortho_frac*csh2o01 + (1.-ortho_frac)*csh2p01; 00297 00298 csh2o02 = 4.86e-11*phycon.te30*phycon.te02*phycon.te02; 00299 csh2p02 = 6.82e-11*phycon.te30/phycon.te02; 00300 csh202 = ortho_frac*csh2o02 + (1.-ortho_frac)*csh2p02; 00301 } 00302 else 00303 { 00304 csh212 = 0.; 00305 csh201 = 0.; 00306 csh202 = 0.; 00307 } 00308 00309 /* effective collision strength of O^o(3P) with p 00310 * >>refer oi cs Pequignot, D. 1990, A&A 231, 499 00311 * original data: 00312 * >>refer oi cs Chambaud et al., 1980, J.Phys.B, 13, 4205 (upto 5000K) 00313 * >>refer oi cs Roueff, private communication (10,000K and 20,000K)*/ 00314 if( phycon.te<=1000. ) 00315 { 00316 csp01 = MAX2(2.22e-5*phycon.te/phycon.te10, 00317 /* >>chng 05 jul 05, eden to dense.EdenHCorr */ 00318 /*2.69e-6*phycon.te*phycon.te30)*dense.xIonDense[ipHYDROGEN][1]/dense.eden;*/ 00319 2.69e-6*phycon.te*phycon.te30)*dense.xIonDense[ipHYDROGEN][1]/dense.EdenHCorr; 00320 00321 csp02 = MIN2(7.07e-8*phycon.te32*phycon.te10,2.46e-7* 00322 /* >>chng 05 jul 05, eden to dense.EdenHCorr */ 00323 /*phycon.te32/phycon.te10)*dense.xIonDense[ipHYDROGEN][1]/dense.eden;*/ 00324 phycon.te32/phycon.te10)*dense.xIonDense[ipHYDROGEN][1]/dense.EdenHCorr; 00325 } 00326 else 00327 { 00328 csp01 = MIN2(2.69e-6*phycon.te*phycon.te30,9.21e-5*phycon.te/phycon.te10/ 00329 /* >>chng 05 jul 05, eden to dense.EdenHCorr */ 00330 /*phycon.te03)*dense.xIonDense[ipHYDROGEN][1]/dense.eden;*/ 00331 phycon.te03)*dense.xIonDense[ipHYDROGEN][1]/dense.EdenHCorr; 00332 00333 csp02 = MIN2(2.46e-7*phycon.te32/phycon.te10,5.21e-5*phycon.te/ 00334 /* >>chng 05 jul 05, eden to dense.EdenHCorr */ 00335 /*phycon.te20)*dense.xIonDense[ipHYDROGEN][1]/dense.eden;*/ 00336 phycon.te20)*dense.xIonDense[ipHYDROGEN][1]/dense.EdenHCorr; 00337 } 00338 00339 csp12 = MIN2(2.35e-6*phycon.te*phycon.te05*phycon.te01,3.98e-5* 00340 /* >>chng 05 jul 05, eden to dense.EdenHCorr */ 00341 /*phycon.te20)*dense.xIonDense[ipHYDROGEN][1]/dense.eden;*/ 00342 /*phycon.te70/phycon.te01)*dense.xIonDense[ipHYDROGEN][1]/dense.eden;*/ 00343 phycon.te70/phycon.te01)*dense.xIonDense[ipHYDROGEN][1]/dense.EdenHCorr; 00344 00345 /*cs01 = cse01+csp01+3.*(csh01*dense.xIonDense[ipHYDROGEN][0] + cshe01*dense.xIonDense[ipHELIUM][0] + csh201*hmi.Hmolec[ipMH2g])/ 00346 dense.cdsqte; 00347 cs12 = cse12+csp12+(csh12*dense.xIonDense[ipHYDROGEN][0] + cshe12*dense.xIonDense[ipHELIUM][0] + csh212*hmi.Hmolec[ipMH2g])/ 00348 dense.cdsqte; 00349 cs02 = cse02+csp02+(csh02*dense.xIonDense[ipHYDROGEN][0] + cshe02*dense.xIonDense[ipHELIUM][0] + csh202*hmi.Hmolec[ipMH2g])/ 00350 dense.cdsqte;*/ 00351 00352 cs01 = cse01+csp01+3.*(csh01*dense.xIonDense[ipHYDROGEN][0] + cshe01*dense.xIonDense[ipHELIUM][0] + csh201*hmi.H2_total)/ 00353 dense.cdsqte; 00354 cs12 = cse12+csp12+(csh12*dense.xIonDense[ipHYDROGEN][0] + cshe12*dense.xIonDense[ipHELIUM][0] + csh212*hmi.H2_total)/ 00355 dense.cdsqte; 00356 cs02 = cse02+csp02+(csh02*dense.xIonDense[ipHYDROGEN][0] + cshe02*dense.xIonDense[ipHELIUM][0] + csh202*hmi.H2_total)/ 00357 dense.cdsqte; 00358 00359 /* end changes */ 00360 /*==============================================================*/ 00361 PutCS(cs01,&TauLines[ipT63]); 00362 PutCS(cs12,&TauLines[ipT146]); 00363 PutCS(cs02,&TauDummy); 00364 atom_level3(&TauLines[ipT63],&TauLines[ipT146],&TauDummy); 00365 00366 /* now save pops to add col den in radinc */ 00367 for( i=0; i<3; ++i) 00368 { 00369 /* >>chng 02 oct 23, bug - had been O1Colden rather than O1Pops */ 00370 colden.O1Pops[i] = (realnum)atoms.PopLevels[i]; 00371 } 00372 /*fprintf(ioQQQ,"colllll\t%li\t%.3e\t%.3e\n", 00373 nzone, 00374 colden.O1Pops[0], 00375 colden.O1Pops[1] );*/ 00376 00377 /* >>>chng 99 apr 29, for orion pdr.in calc, temp oscillated due to 00378 * bad dC/dT estimate, due to rapid change in cs with T. added following*/ 00379 /* when neutral collisions onto OI dominate the cooling the derivative 00380 * produced by atom_level3 is vastly too small because of the strong temperature 00381 * dependence of the neutral collision strength. Increase cooling derivative by 00382 * ratio of cs due to electron (nearly constant) and total */ 00383 /* following was correction for this between 00 apr 29 and 02 jul 25 */ 00384 # if 0 00385 if( cs01 > SMALLFLOAT ) 00386 { 00387 thermal.dCooldT += 00388 (cs01-cs)/cs01*TauLines[ipT63].cool* 00389 TauLines[ipT63].EnergyK*thermal.tsq1; 00390 } 00391 # endif 00392 00393 /* >>chng 02 jul 25, keep track of change in cs for 63 micron line */ 00394 if( fabs(phycon.te - oi_te_tsave)/phycon.te > 1e-4 ) 00395 { 00396 /* very first time we come through, previous values are -1 */ 00397 if(oi_te_tsave>0. ) 00398 { 00399 oi_dcdt_tsave = (cs01-oi_cs_tsave) / (phycon.te-oi_te_tsave); 00400 } 00401 else 00402 { 00403 oi_dcdt_tsave = 0.; 00404 } 00405 oi_cs_tsave = cs01; 00406 oi_te_tsave = phycon.te; 00407 00408 /* >>chng 03 jan 21, this factor could become very large - it should 00409 * always be positive since neutral cs's are, and not much bigger than 00410 * the usual derivative */ 00411 /* can't be negative */ 00412 oi_dcdt_tsave = MAX2( 0. , oi_dcdt_tsave); 00413 /* can't be bigger than several times normal dC/dT */ 00414 oi_dcdt_tsave = MIN2( TauLines[ipT63].EnergyK*thermal.tsq1*4. , oi_dcdt_tsave); 00415 } 00416 /* this is only derivative due to change in collision strength, which is capped to be 00417 * less than 4x the thermal derivative just above */ 00418 thermal.dCooldT += TauLines[ipT63].Coll.cool*oi_dcdt_tsave; 00419 /* this is a bebug print statement for the numerical cs derivative */ 00420 { 00421 /*@-redef@*/ 00422 enum{DEBUG_LOC=false}; 00423 /*@+redef@*/ 00424 if( DEBUG_LOC ) 00425 { 00426 fprintf(ioQQQ,"DEBUG OI\t%.2f\tte\t%.5e\tcool\t%.5e\tcs\t%.4e\told\t%.4e\tnew\t%.4e\n", 00427 fnzone, 00428 phycon.te, 00429 TauLines[ipT63].Coll.cool , 00430 TauLines[ipT63].Coll.cs , 00431 TauLines[ipT63].Coll.cool*TauLines[ipT63].EnergyK*thermal.tsq1, 00432 TauLines[ipT63].Coll.cool*oi_dcdt_tsave ); 00433 } 00434 } 00435 /* [O II] 00436 * two sets of A's, Zeippen 1982 being the default 00437 * >>chng 97 feb 19, change in last sig fig as per Bob Rubin e-mail feb 11 97 00438 * A from NIST compilation 00439 * >>refer o2 as Wiese, W.L., Fuhr, J.R., Deters, T.M. 1996, J Phys Chem Ref Data, 00440 * >>refercon Monograph 7 00441 * these are selected with the set atomic data ion oxygen 2 command */ 00442 if( dense.lgAsChoose[ipOXYGEN][1] ) 00443 { 00444 /* not used unless requested with set atomic data command */ 00445 a21 = 3.058e-5; 00446 a31 = 1.776e-4; 00447 a41 = 5.22e-2; 00448 a51 = 2.12e-2; 00449 a32 = 1.30e-7; 00450 a42 = 0.09907; 00451 a52 = 0.0519; 00452 a43 = 0.0534; 00453 a53 = 0.08672; 00454 a54 = 1.41e-10; 00455 } 00456 else 00457 { 00458 /* default - the Zeippen (1982) values 00459 *>>refer O2 As Zeippen, C. J. 1982, MNRAS, 198, 111 00460 * these are recommended by 00461 *>>refer O2 As Wang et al. 2004 A&A, 427, 873 00462 >>refer O2 As Chen, Qing & Li Phys Rev A, 2007, 76, 042507 00463 */ 00464 a21 = 3.82e-5; 00465 a31 = 1.65e-4; 00466 a41 = 5.64e-2; 00467 a51 = 2.3202e-2; 00468 a32 = 1.20e-7; 00469 a42 = 0.11678; 00470 a52 = 0.0615; 00471 a43 = 0.0614; 00472 a53 = 0.10175; 00473 a54 = 2.08e-11; 00474 } 00475 /* >>chng 01 may 18, update cs to 00476 * >>referold o2 cs McLaughlin, B.M., & Bell, K.L. 1998, J Phys B 31, 4317 00477 * >>chng 02 mar 13, go back to older values as per Seaton/Osterbrock correspondence 00478 * this is a mix of the McLauthlin Bell 93 and 98 results, with the expected 00479 * branching to levels within the term 00480 * >>chng 04 nov 01, statistical weights were reversed, caught by Kevin Blagrave 00481 * >>refer o2 cs A.K.Pradhan, M.Montenegro,S.N.Nahar & W.Eissner,2006,MNRAS,366,L6-L9 00482 */ 00483 cs21 = (realnum)(0.8229*(phycon.te007*phycon.te0005*phycon.te0001)); 00484 cs31 = (realnum)(0.5976/phycon.te002); 00485 cs32 = (realnum)(2.6084/(phycon.te07/(phycon.te003*phycon.te0001))); 00486 cs41 = (realnum)(0.2471*phycon.te02*(phycon.te007/phycon.te0004)); 00487 cs42 = (realnum)(0.7075*phycon.te03*phycon.te004*phycon.te0002); 00488 cs43 = (realnum)(0.415*phycon.te04*(phycon.te004/phycon.te0004)); 00489 cs51 = (realnum)(0.1294*(phycon.te02/(phycon.te001*phycon.te0004))); 00490 cs52 = (realnum)(0.2841*phycon.te04*phycon.te0004); 00491 cs53 = (realnum)(0.2737*phycon.te04*phycon.te003*phycon.te0001); 00492 cs54 = (realnum)(0.2037*phycon.te04*(phycon.te002/phycon.te0004)); 00493 atom_pop5(go2,exo2,cs21,cs31,cs41,cs51,cs32,cs42,cs52,cs43,cs53,cs54, 00494 a21,a31,a41,a51,a32,a42,a52,a43,a53,a54,p,dense.xIonDense[ipOXYGEN][1]); 00495 00496 CoolHeavy.O3730 = (realnum)(p[1]*a21*5.34e-12); 00497 CoolHeavy.O3726 = (realnum)(p[2]*a31*5.34e-12); 00498 CoolHeavy.O2471 = (realnum)((p[3]*a41 + p[4]*a51)*8.05e-12); 00499 CoolHeavy.O7323 = (realnum)((p[3]*a42 + p[4]*a52)*2.72e-12); 00500 CoolHeavy.O7332 = (realnum)((p[3]*a43 + p[4]*a53)*2.71e-12); 00501 CoolHeavy.c3727 = CoolHeavy.O3730 + CoolHeavy.O3726; 00502 CoolHeavy.c7325 = CoolHeavy.O7323 + CoolHeavy.O7332; 00503 thermal.dCooldT += CoolHeavy.c3727*(3.86e4*thermal.tsq1 - thermal.halfte); 00504 CoolAdd("O 2",3727,CoolHeavy.c3727); 00505 CoolAdd("O 2",7325,CoolHeavy.c7325); 00506 CoolAdd("O 2",2470,CoolHeavy.c7325*0.78); 00507 00508 /* remember ratio of radiative to total decays, to use for estimating 00509 * recombination contribution in lines_lv1_li_ne */ 00510 if( (p[3] + p[4]) > SMALLFLOAT ) 00511 { 00512 CoolHeavy.O2_A3_tot = (p[3]*(a41+a42+a43) + p[4]*(a51+a52+a53) ) / 00513 ( (p[3]*(a41+a42+a43) + p[4]*(a51+a52+a53) ) + 00514 ( p[3]*(cs41+cs42+cs43)/go2[3] + p[4]*(cs51+cs52+cs53)/go2[4]) * 00515 dense.cdsqte ); 00516 } 00517 else 00518 { 00519 CoolHeavy.O2_A3_tot = 0.; 00520 } 00521 00522 if( (p[1] + p[2]) > SMALLFLOAT ) 00523 { 00524 CoolHeavy.O2_A2_tot = (p[1]*a21 + p[2]*a31 ) / 00525 ( (p[1]*a21 + p[2]*a31 ) + 00526 ( p[1]*cs21/go2[1] + p[2]*cs31/go2[2]) * 00527 dense.cdsqte ); 00528 } 00529 else 00530 { 00531 CoolHeavy.O2_A2_tot = 0.; 00532 } 00533 00534 /* O II 833.9A, CS 00535 * >>refer o2 cs McLaughlin, B.M., & Bell, K.L. 1993, ApJ, 408, 753 */ 00536 /* >>chng 01 aug 10, turn this line back on - no reason given for turning it off. */ 00537 cs = 0.355*phycon.te10*phycon.te10*phycon.te003*phycon.te001* 00538 phycon.te001; 00539 PutCS(cs,&TauLines[ipT834]); 00540 atom_level2(&TauLines[ipT834]); 00541 00542 /* O III 1666, crit den=2.6+10, A from 00543 * >>refer all cs Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103, 00544 * >>refercon ed by D.R. Flower, (D. Reidel: Holland), 143 00545 * c.s. 00546 * >>referold o3 cs Lennon, D.J. Burke, V.M. 1994, A&AS, 103, 273 */ 00547 /* >>chng 06 jun 30- Humeshkar Nemala */ 00548 /*Refer o3 cs Aggarwal,K.M. & Keenan,F. P.1999,ApJS,123,311*/ 00549 /*Data available in the temperature range 2500 K to 2E6 K*/ 00550 /*fits at temperatures below 30000K and at temperatures above 30000K*/ 00551 if(phycon.te < 30000) 00552 { 00553 cs = (realnum)(0.2519*(phycon.te07*phycon.te02*phycon.te007*phycon.te001*phycon.te0002)); 00554 } 00555 else 00556 { 00557 cs = (realnum)(6.166388*(1/(phycon.te20*phycon.te01*phycon.te002))); 00558 } 00559 /*PutCS(0.756,&TauLines[ipT1666]);*/ 00560 PutCS(cs ,&TauLines[ipT1666]); 00561 /*>>chng 06 jun 30- Humeshkar Nemala*/ 00562 /*>>Refer o3 cs Aggarwal,K.M. & Keenan,F. P.1999,ApJS,123,311*/ 00563 /*Data available in the temperature range 2500 K to 2E6 K*/ 00564 /*fits at temperatures below 30000K and at temperatures above 30000K*/ 00565 if(phycon.te < 30000) 00566 { 00567 cs=(realnum)((0.1511)*(phycon.te07*phycon.te02*phycon.te007*phycon.te001*phycon.te0002)); 00568 } 00569 else 00570 { 00571 cs=(realnum)((3.668474)*(1/(phycon.te20*phycon.te01*phycon.te001*phycon.te0002))); 00572 } 00573 /*PutCS(0.454,&TauLines[ipT1661]);*/ 00574 PutCS(cs,&TauLines[ipT1661]); 00575 /*this line correspond to the transition 3P0-5S^o2*/ 00576 /*>>chng 06 jun 30- Humeshkar Nemala*/ 00577 /*>>Refer o3 cs Aggarwal,K.M. & Keenan,F. P.1999,ApJS,123,311*/ 00578 /*Data available in the temperature range 2500 K to 2E6 K*/ 00579 /*fits at temperatures below 30000K and at temperatures above 30000K*/ 00580 if(phycon.te < 30000) 00581 { 00582 cs = (realnum)(0.0504*((phycon.te10/phycon.te01)*phycon.te005*phycon.te003*phycon.te0002)); 00583 } 00584 else 00585 { 00586 cs = (realnum)(1.223633/(phycon.te20*phycon.te01*phycon.te001*phycon.te0002)); 00587 } 00588 /*PutCS(1.3,&TauDummy);*/ 00589 PutCS(cs,&TauDummy); 00590 atom_level3(&TauDummy,&TauLines[ipT1666],&TauLines[ipT1661]); 00591 /* if( nLev3Fail.gt.0 )write(qq,'('' nlev3fail='',i3)') nLev3Fail 00592 * 00593 * OIII 304 not in cooling function */ 00595 TauLines[ipT304].Emis->PopOpc = dense.xIonDense[ipOXYGEN][2]; 00596 TauLines[ipT304].Lo->Pop = dense.xIonDense[ipOXYGEN][2]; 00597 TauLines[ipT304].Hi->Pop = 0.; 00598 00599 /* o iii 5007+4959, As 96 NIST */ 00600 /* >>chng 01 may 04, As updated to 00601 * >>refer o3 as Storey, P.J., & Zeippen, C.J., 2000, MNRAS, 312, 813-816, 00602 * changed by 10 percent! */ 00603 aeff = 0.027205 + oxy.d5007r; 00604 /* following data used in routine that deduces OIII temp from spectrum 00605 * 00606 * >>refer o3 cs Lennon, D.J. Burke, V.M. 1994, A&AS, 103, 273 00607 * IP paper Y(ki) differ significantly from those 00608 * calculated by 00609 * >>refer o3 cs Burke, V.M., Lennon, D.J., & Seaton, M.J. 1989, MNRAS, 236, 353 00610 * especially for 1D2-1S0. 00611 * BLS89 is adopted for 1D2-1S0 and LB94 for 3P2,1-1D2. 00612 * NB!! these cs's must be kept parallel with those in Peimbert analysis */ 00613 00614 /* >>refer o3 cs Aggarwal,K.M. & Keenan,F. P.1999,ApJS,123,311*/ 00615 /* >>chng 06 jul 24- Humeshkar Nemala */ 00616 /*cs are available over two temperature ranges: below 30000K and above 30000K*/ 00617 /*The old values seemed to saturate at around 100,000K.The new values are around 00618 10-15% different from the old values*/ 00619 /*The cs of the transitions 3P0,1,2 to 1D2 are added together to give oxy.o3cs12 */ 00620 /*the cs of the transition 1D2-1S0 is mentioned as o3cs23*/ 00621 /*The cs of the transitions 3P0,1,2 to 1S0 are added together to give oxy.o3cs13*/ 00622 if(phycon.te < 3E4) 00623 { 00624 oxy.o3cs12 = (realnum)(0.9062*(phycon.te10/phycon.te002)); 00625 o3cs23 = (realnum)(0.0995*phycon.te10*phycon.te07*(phycon.te007/phycon.te0002)); 00626 oxy.o3cs13 = (realnum)(0.1237*phycon.te07*phycon.te02*phycon.te005*phycon.te0005); 00627 00628 } 00629 else if(phycon.te < 6E4) 00630 { 00631 oxy.o3cs12 = (realnum)(1.710073*(phycon.te04/phycon.te004)*phycon.te0004); 00632 oxy.o3cs13 = (realnum)(0.1963109*phycon.te05*phycon.te0007); 00633 o3cs23 = (realnum)(0.781266/(phycon.te02*phycon.te003*phycon.te0001)); 00634 } 00635 else 00636 { 00637 oxy.o3cs12 = (realnum)(6.452638/((phycon.te10/phycon.te02)*phycon.te004*phycon.te0003)); 00638 oxy.o3cs13 = (realnum)(0.841578/(phycon.te07*phycon.te01*(phycon.te002/phycon.te0004))); 00639 o3cs23 = (realnum)(0.781266/(phycon.te02*phycon.te003*phycon.te0001)); 00640 } 00641 # if 0 00642 /*Old code commented -Humeshkar Nemala */ 00643 if( phycon.te > 3981 && phycon.te <= 1e5 ) 00644 { 00645 oxy.o3cs12 = (realnum)(3.0211144 - 101.57536/phycon.sqrte + 355.06905* 00646 /*really is base e log of temperature */ 00647 log(phycon.te)/phycon.te); 00648 o3cs23 = 0.32412181 + 79.051672/phycon.sqrte - 4374.7816/ 00649 phycon.te; 00650 } 00651 else if( phycon.te < 3981. ) 00652 { 00653 oxy.o3cs12 = 2.15f; 00654 o3cs23 = 0.4781; 00655 } 00656 else 00657 { 00658 oxy.o3cs12 = 2.6594f; 00659 o3cs23 = 0.5304; 00660 } 00661 oxy.o3cs13 = 00662 (realnum)(MIN2(0.36,0.0784*phycon.te10*phycon.te03*phycon.te01* 00663 phycon.te003)); 00664 /*End :Old Code Commented Out-Humeshkar Nemala*/ 00665 /* >>chng 01 may 04, derive a21 from aeff, had been 0.02431 00666 a21 = 0.02431;*/ 00667 # endif 00668 00669 a21 = aeff - oxy.d5007r; 00670 a31 = .215634; 00671 a32 = 1.71; 00672 oxy.o3ex23 = 32947.; 00673 oxy.o3br32 = (realnum)(a32/(a31 + a32)); 00674 oxy.o3enro = (realnum)(4.56e-12/3.98e-12); 00675 /* POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2) */ 00676 oxy.poiii3 = 00677 (realnum)(atom_pop3(9.,5.,1.,oxy.o3cs12,oxy.o3cs13,o3cs23, 00678 a21,a31,a32,28990.,oxy.o3ex23,&oxy.poiii2,dense.xIonDense[ipOXYGEN][2], 00679 oxy.d5007r,0.,0.)); 00680 CoolHeavy.c4363 = oxy.poiii3*a32*4.56e-12; 00681 CoolHeavy.c5007 = oxy.poiii2*a21*3.98e-12; 00682 oxy.d5007t = (realnum)(CoolHeavy.c5007*oxy.d5007r/aeff); 00683 thermal.dCooldT += CoolHeavy.c5007*(2.88e4*thermal.tsq1 - thermal.halfte); 00684 thermal.dCooldT += CoolHeavy.c4363*6.20e4*thermal.tsq1; 00685 CoolAdd("O 3",5007,CoolHeavy.c5007); 00686 CoolAdd("O 3",4363,CoolHeavy.c4363); 00687 CoolAdd("O 3",2331,CoolHeavy.c4363*0.236); 00688 /* >>chng 02 jun 27, prevent crash on Sun with gcc 3.1, PvH */ 00689 /* if( xIonDense[ipOXYGEN][2] > 0. ) */ 00690 if( MAX2(oxy.poiii2,oxy.poiii3) > 0.f && dense.xIonDense[ipOXYGEN][2]>SMALLFLOAT) 00691 { 00692 oxy.poiii3 /= dense.xIonDense[ipOXYGEN][2]; 00693 oxy.poiii2 /= dense.xIonDense[ipOXYGEN][2]; 00694 } 00695 00696 /* OIII IR lines, col str from iron project, 00697 * >>referold o3 cs Lennon, D.J. Burke, V.M. 1994, A&AS, 103, 273 */ 00698 /*>>refer o3 cs Aggarwal,K.M. & Keenan,F. P.1999,ApJS,123,311*/ 00699 /*88 microns refers to the 3P0-3P1 transition*/ 00700 if(phycon.te < 3E4) 00701 { 00702 cs = (realnum)(0.3993*(phycon.te03/phycon.te001)*phycon.te0001); 00703 } 00704 else if(phycon.te < 1E5) 00705 { 00706 cs = (realnum)(0.245712*phycon.te07*phycon.te005*phycon.te001*phycon.te0002); 00707 } 00708 else 00709 { 00710 cs = (realnum)(1.310467/((phycon.te07/phycon.te001)*phycon.te0002)); 00711 } 00712 /*PutCS(0.5590,&TauLines[ipTO88]);*/ 00713 PutCS(cs,&TauLines[ipTO88]); 00714 /*52 microns refers to the 3P1-3P2 transition*/ 00715 if(phycon.te < 3E4) 00716 { 00717 cs = (realnum)(0.7812*(phycon.te05/phycon.te0001)); 00718 } 00719 else if(phycon.te < 1.2E5) 00720 { 00721 cs = (realnum)(0.516157*phycon.te07*phycon.te02*phycon.te0001); 00722 } 00723 else 00724 { 00725 cs = (realnum)(4.038402/(phycon.te05*phycon.te03*phycon.te005*phycon.te0007*phycon.te0001)); 00726 } 00727 /*PutCS(1.335,&TauLines[ipT52]);*/ 00728 PutCS(cs,&TauLines[ipT52]); 00729 /*TauDummy refers to the 3P0-3P2 transition*/ 00730 if(phycon.te < 3E4) 00731 { 00732 cs = (realnum)(0.1337*phycon.te07*phycon.te002*phycon.te0002); 00733 } 00734 else if(phycon.te < 1.2E5) 00735 { 00736 cs = (realnum)(0.086446*phycon.te10*phycon.te01*phycon.te004*phycon.te0005); 00737 } 00738 else 00739 { 00740 cs = (realnum)(0.82031/(phycon.te07*phycon.te007*phycon.te0007*phycon.te0002)); 00741 } 00742 /*PutCS(0.2832,&TauDummy);*/ 00743 PutCS(cs,&TauDummy); 00744 atom_level3(&TauLines[ipTO88],&TauLines[ipT52],&TauDummy); 00745 00746 /* O III 834A, CS */ 00747 /* >>referold o3 cs Aggarwal, K.M., 1985 A&A 146, 149. */ 00748 /* >>chng 06 jul 22- Humeshkar Nemala */ 00749 /*>>refer o3 cs Aggarwal,K.M. & Keenan,F. P.1999,ApJS,123,311*/ 00750 /*the cs of OIII 834A was obtained by summing the cs of the six transitions: 00751 3P0-3D^o1,3P1-3D^o1,3P2-3D^o1,3P1-3D^o2,3P2-3D^o2,3P2-3D^o3*/ 00752 if(phycon.te < 3E4) 00753 { 00754 cs = (realnum)(0.9595*phycon.te04*phycon.te0002); 00755 } 00756 else 00757 { 00758 cs = (realnum)(0.10798*phycon.te20*phycon.te05*phycon.te002*phycon.te0001); 00759 } 00760 /*PutCS(5.,&TauLines[ipT835]);*/ 00761 PutCS(cs,&TauLines[ipT835]); 00762 atom_level2(&TauLines[ipT835]); 00763 /* call DumpLine( t835 ); */ 00764 00765 /* these cs data only extend over a modest Temp range */ 00766 Te_bounded = MAX2(phycon.te,4500.); 00767 Te_bounded = MIN2(Te_bounded,450000.); 00768 Te_bounded_log = log(Te_bounded); 00769 00770 /* O IV 789A, CS from 00771 * >>refer o4 cs Zhang, H.L., Graziani, M., Pradhan, A.K. 1994, A&A, 283, 319 */ 00772 cs = -3.0102462 + 109.22973/Te_bounded_log - 18666.357/Te_bounded; 00773 PutCS(cs,&TauLines[ipT789]); 00774 atom_level2(&TauLines[ipT789]); 00775 00776 /* O IV 26 micron, CS 00777 * >>referold o4 cs Blum, R.D., & Pradhan, A.K. 1992, ApJS 80, 425 00778 * A= 00779 * >>refer o4 as Brage, T., Judge, P.G., & Brekke, P. 1996, ApJ. 464, 1030 */ 00780 /* O IV 26 micron */ 00781 /*cs = 9.2141614 - 5.1727759e-3*sqrt(Te_bounded) - 58.116447/Te_bounded_log;*/ 00782 00783 /* >>chng 06 jun 30- Humeshkar Nemala */ 00784 /*>>refer o4 cs Zhang,H.L.,Graziani,M. & Pradhan,A.K. 1994,A&A,283,319*/ 00785 /*cs given over the temp range :900K to 450000K 00786 00787 if(phycon.te < 1.8E4) 00788 { 00789 cs = (realnum)(0.5363*phycon.te10*phycon.te05*(phycon.te01/phycon.te0004)); 00790 } 00791 else if(phycon.te < 5.4E4) 00792 { 00793 cs = (realnum)(1.554285*phycon.te05*phycon.te001); 00794 } 00795 else 00796 { 00797 cs = (realnum)(95.373669/(phycon.te30*phycon.te02*(phycon.te007/phycon.te0002))); 00798 } 00799 00800 * >>chng 06 nov 08 - NPA. Update collision strength to new data from: 00801 * >>refer o4 cs Tayal, S. 2006, ApJS 166, 634 00802 * Equation derived by using TableCurve, and goes to zero as 00803 * T => 0 and T => infinity */ 00804 00805 cs = (realnum)(exp(3.265723 - 0.00014686984*phycon.alnte*phycon.sqrte 00806 -22.035066/phycon.alnte)); 00807 /* cs goes to zero at very high T, which will cause an assert in the 00808 * n-level solver - don't let this happen */ 00809 cs = MAX2( cs , 1e-4 ); 00810 PutCS(cs,&TauLines[ipT26]); 00811 00812 /* one time initialization if first call, and level 2 lines are on */ 00813 if( lgFirst && nWindLine ) 00814 { 00815 lgFirst = false; 00816 nO4Pump = 0; 00817 for( i=0; i<nWindLine; ++i ) 00818 { 00819 /* don't test on nelem==ipIRON since lines on physics, not C, scale */ 00820 if( TauLine2[i].Hi->nelem ==8 && TauLine2[i].Hi->IonStg==4 ) 00821 { 00822 ++nO4Pump; 00823 } 00824 } 00825 if( nO4Pump<0 ) 00826 TotalInsanity(); 00827 else if( nO4Pump > 0 ) 00828 /* create the space - can't malloc 0 bytes */ 00829 ipO4Pump = (long *)MALLOC((unsigned)(nO4Pump)*sizeof(long) ); 00830 nO4Pump = 0; 00831 for( i=0; i<nWindLine; ++i ) 00832 { 00833 /* don't test on nelem==ipIRON since lines on physics, not C, scale */ 00834 if( TauLine2[i].Hi->nelem ==8 && TauLine2[i].Hi->IonStg==4 ) 00835 { 00836 # if 0 00837 DumpLine( &TauLine2[i] ); 00838 # endif 00839 ipO4Pump[nO4Pump] = i; 00840 ++nO4Pump; 00841 } 00842 } 00843 } 00844 else 00845 /* level 2 lines are not enabled */ 00846 nO4Pump = 0; 00847 00848 /* now sum pump rates */ 00849 pump_rate = 0.; 00850 for( i=0; i<nO4Pump; ++i ) 00851 { 00852 pump_rate += TauLine2[ipO4Pump[i]].Emis->pump; 00853 # if 0 00854 fprintf(ioQQQ,"DEBUG C %li %.3e %.3e\n", 00855 i, 00856 TauLine2[ipO4Pump[i]].WLAng , TauLine2[ipO4Pump[i]].pump ); 00857 # endif 00858 } 00859 00860 /*atom_level2(&TauLines[ipT26]);*/ 00861 /*AtomSeqBoron compute cooling from 5-level boron sequence model atom */ 00862 /* >>refer o4 cs Blum, R.D., & Pradhan, A.K., 1992, ApJS 80, 425 00863 * >>refer o4 cs Zhang, H.L., Graziani, M., Pradhan, A.K. 1994, A&A, 283, 319 */ 00864 AtomSeqBoron(&TauLines[ipT26], 00865 &TauLines[ipO4_1400], 00866 &TauLines[ipO4_1397], 00867 &TauLines[ipO4_1407], 00868 &TauLines[ipO4_1405], 00869 &TauLines[ipO4_1401], 00870 0.1367 , 0.1560 , 1.1564 , 0.0457 , pump_rate,"O 4"); 00871 00872 /* O V 630, CS from 00873 * >>refer o5 cs Berrington, K.A., Burke, P.G., Dufton, P.L., Kingston, A.E. 1985, 00874 * >>refercon At. Data Nucl. Data Tables, 33, 195 */ 00875 cs = MIN2(4.0,1.317*phycon.te10/phycon.te03); 00876 PutCS(cs,&TauLines[ipT630]); 00877 atom_level2(&TauLines[ipT630]); 00878 00879 /* O V 1218; coll data from 00880 * >>refer o4 cs Berrington, K.A., Burke, P.G., Dufton, P.L., Kingston, A.E. 1985, 00881 * >>refercon At. Data Nucl. Data Tables, 33, 195 */ 00882 if( phycon.te <= 3.16e4 ) 00883 { 00884 cs = 3.224/(phycon.te10*phycon.te03*phycon.te03*phycon.te003); 00885 } 00886 else 00887 { 00888 cs = 10.549/(phycon.te10*phycon.te10*phycon.te10/phycon.te03); 00889 } 00890 /* >>chng 01 sep 09, AtomSeqBeryllium will reset this to 1/3 so critical density correct */ 00891 PutCS(cs,&TauLines[ipT1214]); 00892 /* c1214 = AtomSeqBeryllium( .87,1.05,3.32, t1214, .0216) * 1.64E-11 00893 * AtomSeqBeryllium(CS23,CS24,CS34,tarray,A41) 00894 * A's 00895 * >>refer o5 cs Fleming, J., Bell, K.L, Hibbert, A., Vaeck, N., Godefroid, M.R. 00896 * >>refercon 1996, MNRAS, 279, 1289 */ 00897 AtomSeqBeryllium(.87,0.602,2.86,&TauLines[ipT1214],.02198); 00898 embesq.em1218 = (realnum)(atoms.PopLevels[3]*0.0216*1.64e-11); 00899 00900 /* O VI 1035 li seq 00901 * generate collision strengths, then stuff them in 00902 * >>refer o6 vs Cochrane, D.M., & McWhirter, R.W.P. 1983, PhyS, 28, 25 */ 00903 ligbar(8,&TauLines[ipT1032],&TauLines[ipT150],&cs2s2p,&cs2s3p); 00904 PutCS(cs2s2p,&TauLines[ipT1032]); 00905 PutCS(cs2s2p*0.5,&TauLines[ipT1037]); 00906 /* no data for the 2-3 transition */ 00907 PutCS(1.0,&TauDummy); 00908 /* solve the 3 level atom */ 00909 atom_level3(&TauLines[ipT1037],&TauDummy,&TauLines[ipT1032]); 00910 00911 PutCS(cs2s3p,&TauLines[ipT150]); 00912 atom_level2(&TauLines[ipT150]); 00913 return; 00914 } 00915 00916 /*"oi3pcs.h" make collision strengths with electrons for 3P O I ground term. 00917 *>>refer o1 cs Berrington, K.A. 1988, J.Phys.B, 21, 1083 for Te > 3000K 00918 *>>refer o1 cs Bell, Berrington & Thomas 1998, MNRAS 293, L83 for 50K <= Te <= 3000K 00919 * numbers in cs variables refer to j value: j=0,1,2 (n=3,2,1 in 5 level atom) 00920 * written by Kirk Korista, 29 may 96 00921 * adapted by Peter van Hoof, 30 march 99 (to include Bell et al. data) 00922 * all data above 3000K have been adjusted down by a constant factor to make 00923 * them line up with Bell et al. data. */ 00924 STATIC void oi3Pcs(double *cs21, 00925 double *cs20, 00926 double *cs10) 00927 { 00928 double a, 00929 b, 00930 c, 00931 d; 00932 00933 DEBUG_ENTRY( "oi3Pcs()" ); 00934 /* local variables */ 00935 00936 /* 3P2 - 3P1 */ 00937 if( phycon.te <= 3.0e3 ) 00938 { 00939 *cs21 = 1.49e-4*phycon.sqrte/phycon.te02/phycon.te02; 00940 } 00941 else if( phycon.te <= 1.0e4 ) 00942 { 00943 a = -5.5634127e-04; 00944 b = 8.3458102e-08; 00945 c = 2.3068232e-04; 00946 *cs21 = 0.228f*(a + b*phycon.te32 + c*phycon.sqrte); 00947 } 00948 else 00949 { 00950 *cs21 = 0.228*MIN2(0.222,5.563e-06*phycon.te*phycon.te05* 00951 phycon.te02); 00952 } 00953 00954 /* 3P2 - 3P0 */ 00955 if( phycon.te <= 3.0e3 ) 00956 { 00957 *cs20 = 4.98e-5*phycon.sqrte; 00958 } 00959 else if( phycon.te <= 1.0e4 ) 00960 { 00961 a = -3.7178028e-04; 00962 b = 2.0569267e-08; 00963 c = 1.1898539e-04; 00964 *cs20 = 0.288*(a + b*phycon.te32 + c*phycon.sqrte); 00965 } 00966 else 00967 { 00968 *cs20 = 0.288*MIN2(0.0589,1.015e-05*phycon.te/phycon.te10/ 00969 phycon.te02/phycon.te005); 00970 } 00971 00972 /* 3P1 - 3P0 */ 00973 if( phycon.te <= 3.0e3 ) 00974 { 00975 *cs10 = 1.83e-9*phycon.te*phycon.te30*phycon.te05; 00976 } 00977 else if( phycon.te <= 1.0e4 ) 00978 { 00979 a = -5.9364373e-04; 00980 b = 0.02946867; 00981 c = 10768.675; 00982 d = 3871.9826; 00983 *cs10 = 0.0269*(a + b*exp(-0.5*POW2((phycon.te-c)/d))); 00984 } 00985 else 00986 { 00987 *cs10 = 0.0269*MIN2(0.074,7.794e-08*phycon.te32/phycon.te10/ 00988 phycon.te01); 00989 } 00990 00991 return; 00992 }