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 /*CoolCalc compute calcium cooling */ 00004 #include "cddefines.h" 00005 #include "taulines.h" 00006 #include "doppvel.h" 00007 #include "phycon.h" 00008 #include "ca.h" 00009 #include "dense.h" 00010 #include "thermal.h" 00011 #include "opacity.h" 00012 #include "rfield.h" 00013 #include "ligbar.h" 00014 #include "lines_service.h" 00015 #include "atoms.h" 00016 #include "cooling.h" 00017 00018 void CoolCalc(void) 00019 { 00020 realnum p2; 00021 double a21, 00022 a31, 00023 a41, 00024 a42, 00025 a51, 00026 a52, 00027 a53, 00028 c21, 00029 Ca2pop[5] , 00030 cs, 00031 cs2s2p, 00032 cs2s3p , 00033 cs01, 00034 cs02, 00035 cs12, 00036 cs14, 00037 cs15, 00038 d41, 00039 d42, 00040 d51, 00041 d52, 00042 d53, 00043 hlgam, 00044 op41, 00045 op51, 00046 opckh, 00047 opcxyz, 00048 PhotoRate2, 00049 p3, 00050 PhotoRate3, 00051 PhotoRate4, 00052 PhotoRate5, 00053 r21, 00054 r31, 00055 r41, 00056 r42, 00057 r51, 00058 r52, 00059 r53; 00060 static double gCa2[5]={2.,4.,6.,2.,4.}; 00061 static double exCa2[4]={13650.2,60.7,11480.6,222.9}; 00062 static realnum opCax = 0.f; 00063 static realnum opCay = 0.f; 00064 static realnum opCaz = 0.f; 00065 00066 DEBUG_ENTRY( "CoolCalc()" ); 00067 00068 /* Ca I 4228 */ 00069 MakeCS(&TauLines[ipCaI4228]); 00070 atom_level2(&TauLines[ipCaI4228]); 00071 00072 /* photoionization of evcited levels by Ly-alpha */ 00073 hlgam = rfield.otslin[ Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont -1]; 00074 PhotoRate5 = 1.7e-18*hlgam; 00075 PhotoRate4 = 8.4e-19*hlgam; 00076 PhotoRate3 = 7.0e-18*hlgam; 00077 PhotoRate2 = 4.8e-18*hlgam; 00078 00079 /* spontaneous decays 00080 * frist two trans prob from 00081 * >>refer ca2 as Zeippen, C.J. 1990, A&A, 229, 248 */ 00082 a21 = 1.02*TauLines[ipT7324].Emis->Pesc; 00083 a31 = 1.05*TauLines[ipT7291].Emis->Pesc; 00084 a41 = 1.4e8*TauLines[ipT3969].Emis->Pesc; 00085 a51 = 1.4e8*TauLines[ipT3934].Emis->Pesc; 00086 a42 = 7.9e6*TauLines[ipT8662].Emis->Pesc; 00087 a52 = 8.2e5*TauLines[ipT8498].Emis->Pesc; 00088 a53 = 7.48e6*TauLines[ipT8542].Emis->Pesc; 00089 00090 /* destruction of IR triplet by continuous opacities */ 00091 opcxyz = opac.opacity_abs[ TauLines[ipT7324].ipCont -1]; 00092 00093 /* opcxyz = opac(icaxyz) */ 00094 if( opcxyz > 0. ) 00095 { 00096 d52 = 5.6*opcxyz/(opcxyz + opCax)*(1. - TauLines[ipT8498].Emis->Pesc); 00097 d53 = 5.6*opcxyz/(opcxyz + opCay)*(1. - TauLines[ipT8542].Emis->Pesc); 00098 d42 = 5.6*opcxyz/(opcxyz + opCaz)*(1. - TauLines[ipT8662].Emis->Pesc); 00099 } 00100 else 00101 { 00102 d52 = 0.; 00103 d53 = 0.; 00104 d42 = 0.; 00105 } 00106 00107 /* near UV dest of KH by background continuum */ 00108 opckh = opac.opacity_abs[ TauLines[ipT3969].ipCont -1]; 00109 00110 /* opckh = opac(icakh) */ 00111 if( opckh > 0. ) 00112 { 00113 op51 = dense.xIonDense[ipCALCIUM][1]*3.89e-7/DoppVel.doppler[19]; 00114 d51 = 5.6*opckh/(opckh + op51); 00115 op41 = dense.xIonDense[ipCALCIUM][1]*1.96e-7/DoppVel.doppler[19]; 00116 d41 = 5.6*opckh/(opckh + op41); 00117 } 00118 else 00119 { 00120 op51 = 0.; 00121 d51 = 0.; 00122 op41 = 0.; 00123 d41 = 0.; 00124 } 00125 /* net rates */ 00126 r21 = PhotoRate2 + a21; 00127 r31 = PhotoRate3 + a31; 00128 r41 = a41 + PhotoRate4 + d41; 00129 r51 = a51 + PhotoRate5 + d51; 00130 r42 = a42 + d42; 00131 r52 = a52 + d52; 00132 r53 = a53 + d53; 00133 cs14 = 0.923*phycon.te10*phycon.te10; 00134 cs15 = cs14*2.; 00135 TauLines[ipT3969].Coll.cs = (realnum)cs14; 00136 TauLines[ipT3934].Coll.cs = (realnum)cs15; 00137 00138 /* following used to correct rec contribution 00139 * fcakh = a51 / ( a51 + eden*1.5e-5 / sqrte ) 00140 * cs 1-2 from 00141 * >>refer ca2 cs Saraph, H.E. 1970, J.Phys. B, 3, 952 00142 * other 00143 * >>refer ca2 cs Chidichimo, M.C. 1981, J.Phys. B, 14, 4149 */ 00144 atom_pop5(gCa2,exCa2,5.8,8.6,cs14,cs15,20.6,22.9,9.8,3.4,44.4,1.0, 00145 r21,r31,r41,r51,0.,r42,r52,0.,r53,0.,Ca2pop,dense.xIonDense[ipCALCIUM][1]); 00146 00147 /* CDSQTE = 8.629E-6*EDEN/SQRTE */ 00148 c21 = 5.8/4.*dense.cdsqte; 00149 00150 /* remember largest ratio of Ly-al removal to total */ 00151 if( dense.xIonDense[ipCALCIUM][1] > 0. ) 00152 { 00153 ca.Ca2RmLya = MAX2(ca.Ca2RmLya,(realnum)(PhotoRate2/(PhotoRate2+a21+c21))); 00154 } 00155 ca.Cak = (realnum)(Ca2pop[4]*a51*5.06e-12); 00156 ca.Cah = (realnum)(Ca2pop[3]*a41*5.01e-12); 00157 ca.Cax = (realnum)(Ca2pop[4]*a52*2.34e-12); 00158 ca.Cay = (realnum)(Ca2pop[4]*a53*2.33e-12); 00159 ca.Caz = (realnum)(Ca2pop[3]*a42*2.30e-12); 00160 ca.Caf1 = (realnum)(Ca2pop[2]*a31*2.73e-12); 00161 ca.Caf2 = (realnum)(Ca2pop[1]*a21*2.72e-12); 00162 ca.popca2ex = (realnum)(Ca2pop[1] + Ca2pop[2] + Ca2pop[3] + Ca2pop[4]); 00163 00164 /* this is the total cooling due to the model atom */ 00165 ca.Cair = ca.Cax + ca.Cay + ca.Caz; 00166 ca.c7306 = ca.Caf1 + ca.Caf2; 00167 ca.Cakh = ca.Cak + ca.Cah; 00168 00169 /* add the CaII lines to the cooling stack */ 00170 /* >>chng 04 feb 28, next three, deriv wrt temp had not been present 00171 * extra factor of three in second term of deriv is due to extra T-1 00172 * dependence of cooling on recombination coefficient. Almost always 00173 * CaII cooling is most important when Ca+ is a trace stage of ionization 00174 * so changes in temp change cooling due to change in Ca+ abundance */ 00175 CoolAdd("Ca 2",7306,ca.c7306); 00176 thermal.dCooldT -= ca.c7306*4.*thermal.halfte; 00177 00178 CoolAdd("Ca 2",8400,ca.Cair); 00179 /* >>chng 04 feb 28, this is emperical fit to cooling deriv - does down with 00180 * increasing temp */ 00181 thermal.dCooldT -= ca.Cair*4.*thermal.halfte; 00182 00183 /* >>chng 04 feb 28 wl had been 3800, chng to 3950, closer to mean of K and R */ 00184 CoolAdd("Ca 2",3950,ca.Cakh); 00185 /* >>chng 04 feb 28, this is emperical fit to cooling deriv - does down with 00186 * increasing temp */ 00187 thermal.dCooldT -= ca.Cakh*4.*thermal.halfte; 00188 00189 /*fprintf(ioQQQ,"DEBUG ca2\t%.2f\t%.5e\t%.4e\t%.4e\n", 00190 fnzone, phycon.te,ca.Cakh,dense.xIonDense[ipCALCIUM][1]);*/ 00191 00192 /* level populations that will be used for excited state photoionization */ 00193 ca.dstCala = (realnum)(Ca2pop[4]*PhotoRate5 + Ca2pop[3]*PhotoRate4); 00194 ca.dCakh = (realnum)(ca.dstCala*5.03e-12); 00195 ca.dCaf12 = (realnum)((Ca2pop[2]*PhotoRate3 + Ca2pop[1]*PhotoRate2)*2.31e-12); 00196 opCax = (realnum)(Ca2pop[1]*1.13e-8/DoppVel.doppler[19]); 00197 opCay = (realnum)(Ca2pop[2]*6.87e-8/DoppVel.doppler[19]); 00198 opCaz = (realnum)(Ca2pop[1]*5.74e-8/DoppVel.doppler[19]); 00199 00200 /* total rate Lalpha destroys CaII, 00201 * this is only used in ioncali to increase ionization rate by 00202 * adding it to the ct vector */ 00203 if( dense.xIonDense[ipCALCIUM][1] > 0. ) 00204 { 00205 ca.dstCala = (realnum)( 00206 (ca.dstCala + ca.dCaf12/2.31e-12)/dense.xIonDense[ipCALCIUM][1]); 00207 /* take average of old and new */ 00208 ca.dstCala = (realnum)((ca.dstCala + ca.oldcla)/2.); 00209 ca.oldcla = ca.dstCala; 00210 { 00211 /*@-redef@*/ 00212 enum {DEBUG_LOC=false}; 00213 /*@+redef@*/ 00214 if( DEBUG_LOC ) 00215 { 00216 fprintf(ioQQQ," hlgam is %e\n", hlgam); 00217 } 00218 } 00219 } 00220 else 00221 { 00222 ca.dstCala = 0.; 00223 } 00224 ca.Ca3d = (realnum)(Ca2pop[1] + Ca2pop[2]); 00225 ca.Ca4p = (realnum)(Ca2pop[3] + Ca2pop[4]); 00226 00227 /* incl stimulated emission for Calcium II 5-level atom */ 00228 TauLines[ipT3934].Emis->PopOpc = (Ca2pop[0] - Ca2pop[4]/2.); 00229 TauLines[ipT3934].Hi->Pop = Ca2pop[4]; 00230 TauLines[ipT3934].Lo->Pop = Ca2pop[0]; 00231 TauLines[ipT3969].Emis->PopOpc = (Ca2pop[0] - Ca2pop[3]); 00232 TauLines[ipT3969].Hi->Pop = Ca2pop[3]; 00233 TauLines[ipT3969].Lo->Pop = Ca2pop[0]; 00234 00235 opac.tpcah[0] = TauLines[ipT3969].Emis->TauIn; 00236 TauLines[ipT8498].Emis->PopOpc = (Ca2pop[1] - Ca2pop[4]); 00237 TauLines[ipT8498].Hi->Pop = Ca2pop[4]; 00238 TauLines[ipT8498].Lo->Pop = Ca2pop[1]; 00239 TauLines[ipT8542].Emis->PopOpc = (Ca2pop[2] - Ca2pop[4]*1.5); 00240 TauLines[ipT8542].Hi->Pop = Ca2pop[4]; 00241 TauLines[ipT8542].Lo->Pop = Ca2pop[2]; 00242 TauLines[ipT8662].Emis->PopOpc = (Ca2pop[1] - Ca2pop[3]*2.); 00243 TauLines[ipT8662].Hi->Pop = Ca2pop[3]; 00244 TauLines[ipT8662].Lo->Pop = Ca2pop[1]; 00245 TauLines[ipT7291].Emis->PopOpc = dense.xIonDense[ipCALCIUM][1]; 00246 TauLines[ipT7291].Hi->Pop = 0.; 00247 TauLines[ipT7291].Lo->Pop = dense.xIonDense[ipCALCIUM][1]; 00248 TauLines[ipT7324].Emis->PopOpc = dense.xIonDense[ipCALCIUM][1]; 00249 TauLines[ipT7324].Hi->Pop = 0.; 00250 TauLines[ipT7324].Lo->Pop = dense.xIonDense[ipCALCIUM][1]; 00251 00252 /* Ca IV 3.2 micron; data from 00253 * >>refer ca4 as Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103, 00254 * >>refercon ed by D.R. Flower, (D. Reidel: Holland), 143 00255 * Y(ik) from 00256 * >>refer ca4 cs Pelan, J., & Berrington, K.A. 1995, A&A Suppl, 110, 209 */ 00257 if( phycon.te <= 1e5 ) 00258 { 00259 cs = MAX2(1.0,8.854e-3*phycon.sqrte); 00260 } 00261 else if( phycon.te < 2.512e5 ) 00262 { 00263 cs = 2.8; 00264 } 00265 else 00266 { 00267 cs = 641.1/(phycon.te30*phycon.te10*phycon.te02*phycon.te02/ 00268 phycon.te003); 00269 } 00270 PutCS(cs,&TauLines[ipTCa3]); 00271 atom_level2(&TauLines[ipTCa3]); 00272 00273 /* [Ca V] IR 4.16, 11.47 micron; A from 00274 * >>refer ca5 as Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103, 00275 * >>refercon ed by D.R. Flower, (D. Reidel: Holland), 143 00276 * cs from 00277 * >>refer ca5 cs Galavis, M.E., Mendoza, C., & Zeippen, C.J. 1995, A&AS, 111, 347 00278 * >>chng 96 jul 16, big changes in cs */ 00279 cs = MIN2(3.3,0.392*phycon.te20/phycon.te005/phycon.te003); 00280 cs = MAX2(2.2,cs); 00281 PutCS(cs,&TauLines[ipTCa4]); 00282 00283 /* >>chng 96 aug 02, following had error in te dep in ver 90.01 */ 00284 cs = MIN2(0.93,0.162*phycon.te10*phycon.te05*phycon.te003* 00285 phycon.te001); 00286 cs = MAX2(0.67,cs); 00287 PutCS(cs,&TauLines[ipTCa12]); 00288 00289 cs = MIN2(0.97,0.0894*phycon.te20*phycon.te01*phycon.te005); 00290 cs = MAX2(0.60,cs); 00291 PutCS(cs,&TauDummy); 00292 00293 atom_level3(&TauLines[ipTCa4],&TauLines[ipTCa12],&TauDummy); 00294 00295 /* Ca V lines from 1d, 1s; A from 00296 * >>refer ca5 as Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103, 00297 * >>refercon ed by D.R. Flower, (D. Reidel: Holland), 143 00298 * cs from 00299 * >>refer ca5 cs Galavis, M.E., Mendoza, C., & Zeippen, C.J. 1995, A&AS, 111, 347 00300 * POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2) */ 00301 cs01 = MIN2(4.1,0.533*phycon.te20/phycon.te01); 00302 cs01 = MAX2(2.8,cs01); 00303 cs02 = MIN2(0.87,5.22e-03*phycon.sqrte); 00304 p3 = atom_pop3(9.,5.,1.,cs01,cs02,1.35,2.326,23.2,3.73,2.57e4,3.60e4, 00305 &p2,dense.xIonDense[ipCALCIUM][4],0.,0.,0.); 00306 00307 ca.c3997 = p3*3.73*4.98e-12; 00308 ca.c2414 = p3*23.1*8.245e-12; 00309 ca.Ca6087 = p2*0.426*3.268e-12; 00310 ca.c5311 = p2*1.90*3.747e-12; 00311 00312 CoolAdd("Ca 5",3997,ca.c3997); 00313 CoolAdd("Ca 5",2414,ca.c2414); 00314 CoolAdd("Ca 5",6087,ca.Ca6087); 00315 CoolAdd("Ca 5",5311,ca.c5311); 00316 00317 /* Ca VII lines from 1d, 1s 00318 * all cs from 00319 * >>refer ca72 cs Galavis, M.E., Mendoza, C., & Zeippen, C.J. 1995, A&AS, 111, 347 00320 * POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2) */ 00321 cs01 = MIN2(4.4,22.25/(phycon.te20/phycon.te02/phycon.te02)); 00322 cs01 = MAX2(3.5,cs01); 00323 00324 cs12 = MIN2(1.20,0.303*phycon.te30*phycon.te03); 00325 cs12 = MAX2(0.62,cs12); 00326 00327 cs02 = MIN2(0.959,7.889/(phycon.te20*phycon.te05/phycon.te01)); 00328 cs02 = MAX2(0.50,cs02); 00329 00330 p3 = atom_pop3(9.,5.,1.,cs01,cs02,cs12,3.124,30.4,6.81,2.91e4,3.90e4, 00331 &p2,dense.xIonDense[ipCALCIUM][6],0.,0.,0.); 00332 00333 ca.Ca3688 = p3*6.81*5.40e-12; 00334 ca.Ca2112 = p3*30.4*9.42e-12; 00335 ca.Ca5620 = p2*2.15*3.548e-12; 00336 ca.Ca4941 = p2*0.974*4.037e-12; 00337 CoolAdd("Ca 7",3688,ca.Ca3688); 00338 CoolAdd("Ca 7",2112,ca.Ca2112); 00339 CoolAdd("Ca 7",5620,ca.Ca5620); 00340 CoolAdd("Ca 7",4941,ca.Ca4941); 00341 00342 /* all cs from 00343 * >>refer ca7 cs Galavis, M.E., Mendoza, C., & Zeippen, C.J. 1995, A&AS, 111, 347 00344 * [Ca VII] 4.09, 6.15 mic 3P lines */ 00345 cs = MIN2(5.354,0.406*phycon.te20*phycon.te03*phycon.te01); 00346 cs = MAX2(3.702,cs); 00347 PutCS(cs,&TauLines[ipCa0741]); 00348 00349 cs = MIN2(1.59,0.183*phycon.te20); 00350 cs = MAX2(1.153,cs); 00351 PutCS(cs,&TauLines[ipCa0761]); 00352 00353 cs = MIN2(1.497,0.0917*phycon.te20*phycon.te05* phycon.te01); 00354 cs = MAX2(1.005,cs); 00355 PutCS(cs,&TauDummy); 00356 00357 /* atom_level3( t10,t21,t20) */ 00358 atom_level3(&TauLines[ipCa0761],&TauLines[ipCa0741],&TauDummy); 00359 00360 /* [Ca VIII] 2.32 microns, cs 00361 * >>refer ca8 cs Saraph, H.E., & Storey, P.J. A&AS, 115, 151 */ 00362 cs = MIN2(6.75,22.04/(phycon.te10*phycon.te02*phycon.te005)); 00363 PutCS(cs,&TauLines[ipCa08232]); 00364 atom_level2(&TauLines[ipCa08232]); 00365 00366 /* [Ca 12] 3328.78A, cs from 00367 * >>refer ca12 cs Saraph, H.E. & Tully, J.A. 1994, A&AS, 107, 29 */ 00368 cs = MIN2(0.172,0.0118*phycon.te20*phycon.te01); 00369 cs = MAX2(0.10,cs); 00370 PutCS(cs,&TauLines[ipCa12333]); 00371 atom_level2(&TauLines[ipCa12333]); 00372 00373 /* Li seq Ca 18 2s2p 2s3p, 2s2p as two separate lines 00374 * >>refer ca18 cs Cochrane, D.M., & McWhirter, R.W.P. 1983, PhyS, 28, 25 */ 00375 ligbar(20,&TauLines[ipTCa302],&TauLines[ipTCa19],&cs2s2p,&cs2s3p); 00376 00377 PutCS(cs2s2p,&TauLines[ipTCa302]); 00378 atom_level2(&TauLines[ipTCa302]); 00379 00380 /* funny factor (should have been 0.5) due to energy change */ 00381 PutCS(cs2s2p*0.439,&TauLines[ipTCa345]); 00382 atom_level2(&TauLines[ipTCa345]); 00383 00384 PutCS(cs2s3p,&TauLines[ipTCa19]); 00385 atom_level2(&TauLines[ipTCa19]); 00386 return; 00387 }