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 #include "cddefines.h" 00004 #include "atmdat.h" 00005 #include "conv.h" 00006 #include "dense.h" 00007 #include "heavy.h" 00008 #include "helike_cs.h" 00009 #include "hydrogenic.h" 00010 #include "hydro_vs_rates.h" 00011 #include "ionbal.h" 00012 #include "iso.h" 00013 #include "opacity.h" 00014 #include "phycon.h" 00015 #include "physconst.h" 00016 #include "rfield.h" 00017 #include "secondaries.h" 00018 #include "trace.h" 00019 #include "taulines.h" 00020 00021 /* These are masses relative to the proton mass of the electron, proton, he+, and alpha particle. */ 00022 static double ColliderMass[4] = {ELECTRON_MASS/PROTON_MASS, 1.0, 4.0, 4.0}; 00023 00024 void iso_collisional_ionization( long ipISO, long nelem ) 00025 { 00026 ASSERT( ipISO <= NISO ); 00027 00028 DEBUG_ENTRY( "iso_collisional_ionization()" ); 00029 00030 /* collisional ionization from ground */ 00031 iso.ColIoniz[ipISO][nelem][0] = iso.lgColl_ionize[ipISO] * 00032 t_ADfA::Inst().coll_ion( nelem+1, 1+ipISO, phycon.te ); 00033 00034 for( long ipHi=1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00035 { 00036 if( nelem == ipISO ) 00037 { 00038 /* use routine from Vriens and Smeets (1981). */ 00039 /* >>refer iso neutral col.ion. Vriens, L., & Smeets, A.H.M. 1980, Phys Rev A 22, 940 */ 00040 iso.ColIoniz[ipISO][nelem][ipHi] = hydro_vs_ioniz( ipISO, nelem, ipHi ); 00041 } 00042 else 00043 { 00044 /* ions */ 00045 /* use hydrogenic ionization rates for ions 00046 * >>refer iso ions col.ion. Allen 1973, Astro. Quan. for low Te. 00047 * >>refer iso ions col.ion. Sampson and Zhang 1988, ApJ, 335, 516 for High Te. 00048 * */ 00049 iso.ColIoniz[ipISO][nelem][ipHi] = Hion_coll_ioniz_ratecoef( ipISO, nelem, ipHi ); 00050 } 00051 00052 /* iso.lgColl_ionize is option to turn off collisions, "atom XX-like collis off" comnd */ 00053 /* always leave highest level coupled to continuum */ 00054 if( ipHi < iso.numLevels_max[ipISO][nelem] - 1 ) 00055 iso.ColIoniz[ipISO][nelem][ipHi] *= iso.lgColl_ionize[ipISO]; 00056 } 00057 00058 /* Here we arbitrarily scale the highest level ionization to account for the fact 00059 * that, if the atom is not full size, this level should be interacting with higher 00060 * levels and not just the continuum. We did add on collisional excitation terms instead 00061 * but this caused a problem at low temperatures because the collisional ionization was 00062 * a sum of terms with different Boltzmann factors, while PopLTE had just one Boltzmann 00063 * factor. The result was a collisional recombination that had residual exponentials of 00064 * the form exp(x/kT), which blew up at small T. */ 00065 if( !iso.lgLevelsLowered[ipISO][nelem] ) 00066 { 00067 iso.ColIoniz[ipISO][nelem][iso.numLevels_max[ipISO][nelem]-1] *= 10.;/*1.E5;*/ 00068 iso.ColIoniz[ipISO][nelem][iso.numLevels_max[ipISO][nelem]-1] *= 10.;/*1.E5;*/ 00069 } 00070 00071 return; 00072 } 00073 00074 void iso_suprathermal( long ipISO, long nelem ) 00075 { 00076 DEBUG_ENTRY( "iso_suprathermal()" ); 00077 00078 /* check that we were called with valid parameters */ 00079 ASSERT( ipISO < NISO ); 00080 ASSERT( nelem >= ipISO ); 00081 ASSERT( nelem < LIMELM ); 00082 00083 /*********************************************************************** 00084 * * 00085 * get secondary excitation by suprathermal electrons * 00086 * * 00087 ***********************************************************************/ 00088 00089 for( long i=1; i < iso.numLevels_max[ipISO][nelem]; i++ ) 00090 { 00091 if( Transitions[ipISO][nelem][i][0].ipCont > 0 ) 00092 { 00093 /* get secondaries for all iso lines by scaling LyA 00094 * excitation by ratio of cross section (oscillator strength/energy) 00095 * Born approximation or plane-wave approximation */ 00096 secondaries.Hx12[ipISO][nelem][i] = secondaries.x12tot * 00097 (Transitions[ipISO][nelem][i][0].Emis->gf/ 00098 Transitions[ipISO][nelem][i][0].EnergyWN) / 00099 (Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][0].Emis->gf/ 00100 Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][0].EnergyWN); 00101 } 00102 else 00103 secondaries.Hx12[ipISO][nelem][i] = 0.; 00104 } 00105 00106 return; 00107 } 00108 00109 /*============================*/ 00110 /* evaluate collisional rates */ 00111 void iso_collide( long ipISO, long nelem ) 00112 { 00113 long ipHi, ipLo; 00114 double factor, ConvLTEPOP, edenCorr; 00115 00116 /* this array stores the last temperature at which collision data were evaluated for 00117 * each species of the isoelectronic sequence. */ 00118 static double TeUsed[NISO][LIMELM]={ 00119 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0}, 00120 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0} }; 00121 00122 /* collision strengths are assumed to be roughly constant for small changes in temperature 00123 * and are not recalculated as often as other data. This array stores the last temperature 00124 * at which collision strengths were evaluated for each species of the isoelectronic sequence. */ 00125 static double TeUsedForCS[NISO][LIMELM]={ 00126 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0}, 00127 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0} }; 00128 00129 DEBUG_ENTRY( "iso_collide()" ); 00130 00131 /* check that we were called with valid parameters */ 00132 ASSERT( ipISO < NISO ); 00133 ASSERT( nelem >= ipISO ); 00134 ASSERT( nelem < LIMELM ); 00135 00136 /* skip most of this routine if temperature has not changed, 00137 * the check on conv.nTotalIoniz is to make sure that we redo this 00138 * on the very first call in a grid calc - it is 0 on the first call */ 00139 if( fp_equal( TeUsed[ipISO][nelem], phycon.te ) && conv.nTotalIoniz ) 00140 { 00141 ASSERT( Transitions[ipISO][nelem][ iso.nLyaLevel[ipISO] ][0].Coll.ColUL >= 0. ); 00142 00143 if( trace.lgTrace ) 00144 { 00145 fprintf( ioQQQ, 00146 " iso_collide called %s nelem %li - no reeval Boltz fac, LTE dens\n", 00147 iso.chISO[ipISO], nelem ); 00148 } 00149 } 00150 else 00151 { 00152 TeUsed[ipISO][nelem] = phycon.te; 00153 00154 if( trace.lgTrace ) 00155 { 00156 fprintf( ioQQQ, 00157 " iso_collide called %s nelem %li - will reeval Boltz fac, LTE dens\n", 00158 iso.chISO[ipISO], nelem ); 00159 } 00160 00161 /********************************************************** 00162 * * 00163 * Boltzmann factors for all levels, and * 00164 * collisional ionization and excitation * 00165 * * 00166 **********************************************************/ 00167 00168 for( ipHi=1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00169 { 00170 for( ipLo=0; ipLo<ipHi; ipLo++ ) 00171 { 00172 iso.Boltzmann[ipISO][nelem][ipHi][ipLo] = 00173 sexp( Transitions[ipISO][nelem][ipHi][ipLo].EnergyK / phycon.te ); 00174 } 00175 } 00176 00177 /* HION_LTE_POP is planck^2 / (2 pi m_e k ), raised to 3/2 next */ 00178 factor = HION_LTE_POP*dense.AtomicWeight[nelem]/ 00179 (dense.AtomicWeight[nelem]+ELECTRON_MASS/ATOMIC_MASS_UNIT); 00180 00181 /* term in () is stat weight of electron * ion */ 00182 ConvLTEPOP = pow(factor,1.5)/(2.*iso.stat_ion[ipISO])/phycon.te32; 00183 00184 iso.lgPopLTE_OK[ipISO][nelem] = true; 00185 00186 /* fully define Boltzmann factors to continuum for model levels */ 00187 for( ipLo=0; ipLo<iso.numLevels_max[ipISO][nelem]; ipLo++ ) 00188 { 00189 /* this Boltzmann factor is exp( +ioniz energy / Te ) */ 00190 StatesElem[ipISO][nelem][ipLo].ConBoltz = 00191 dsexp(iso.xIsoLevNIonRyd[ipISO][nelem][ipLo]/phycon.te_ryd); 00192 00193 /*************************************** 00194 * * 00195 * LTE abundances for all levels * 00196 * * 00197 ***************************************/ 00198 00199 if( StatesElem[ipISO][nelem][ipLo].ConBoltz > SMALLDOUBLE ) 00200 { 00201 /* LTE population of given level. */ 00202 iso.PopLTE[ipISO][nelem][ipLo] = 00203 StatesElem[ipISO][nelem][ipLo].g / StatesElem[ipISO][nelem][ipLo].ConBoltz * ConvLTEPOP; 00204 ASSERT( iso.PopLTE[ipISO][nelem][ipLo] < BIGDOUBLE ); 00205 } 00206 else 00207 { 00208 iso.PopLTE[ipISO][nelem][ipLo] = 0.; 00209 } 00210 00211 /* now check for any zeros - if present then matrix cannot be used */ 00212 if( iso.PopLTE[ipISO][nelem][ipLo] <= 0. ) 00213 { 00214 iso.lgPopLTE_OK[ipISO][nelem] = false; 00215 } 00216 } 00217 00218 iso_collisional_ionization( ipISO, nelem ); 00219 00220 /*********************************************************** 00221 * * 00222 * collision strengths for all lines in iso sequence * 00223 * * 00224 ***********************************************************/ 00225 00226 /* Calculate initial value of collision strengths, OR 00227 * Reevaluate collision strengths if the temperature has changed by more than 15%. */ 00228 if( (TeUsedForCS[ipISO][nelem] == 0.) || 00229 ( TeUsedForCS[ipISO][nelem]/phycon.te > 1.15 ) || 00230 ( TeUsedForCS[ipISO][nelem]/phycon.te < 0.85 ) ) 00231 { 00232 /* Update temperature at which collision strengths are evaluated. */ 00233 TeUsedForCS[ipISO][nelem] = phycon.te; 00234 00235 for( ipHi=1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00236 { 00237 for( ipLo=0; ipLo < ipHi; ipLo++ ) 00238 { 00239 if( ipISO == ipH_LIKE ) 00240 { 00241 /* collision strengths due to electron impact */ 00242 Transitions[ipISO][nelem][ipHi][ipLo].Coll.cs = 00243 HydroCSInterp( nelem , ipHi , ipLo, ipELECTRON ); 00244 /* proton impact */ 00245 Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipPROTON] = 00246 HydroCSInterp( nelem , ipHi , ipLo, ipPROTON ); 00247 /* He+ impact */ 00248 Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipHE_PLUS] = 00249 HydroCSInterp( nelem , ipHi , ipLo, ipHE_PLUS ); 00250 /* and He++ impact */ 00251 Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipALPHA] = 00252 HydroCSInterp( nelem , ipHi , ipLo, ipALPHA ); 00253 } 00254 else if( ipISO == ipHE_LIKE ) 00255 { 00256 /* collision strengths due to electron impact */ 00257 Transitions[ipISO][nelem][ipHi][ipLo].Coll.cs = 00258 HeCSInterp( nelem , ipHi , ipLo, ipELECTRON ); 00259 /* proton impact */ 00260 Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipPROTON] = 00261 HeCSInterp( nelem , ipHi , ipLo, ipPROTON ); 00262 /* and He+ impact */ 00263 Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipHE_PLUS] = 00264 HeCSInterp( nelem , ipHi , ipLo, ipHE_PLUS ); 00265 00266 Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipALPHA] = 0.; 00267 } 00268 else 00269 TotalInsanity(); 00270 00271 if( N_(ipHi) != N_(ipLo) ) 00272 { 00273 /* kill all collisional excitation if this set with 00274 * atom XX-like collisional excitation */ 00275 Transitions[ipISO][nelem][ipHi][ipLo].Coll.cs *= iso.lgColl_excite[ipISO]; 00276 Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipPROTON] *= iso.lgColl_excite[ipISO]; 00277 Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipHE_PLUS] *= iso.lgColl_excite[ipISO]; 00278 Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipALPHA] *= iso.lgColl_excite[ipISO]; 00279 00280 } 00281 else 00282 { 00283 Transitions[ipISO][nelem][ipHi][ipLo].Coll.cs *= iso.lgColl_l_mixing[ipISO]; 00284 Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipPROTON] *= iso.lgColl_l_mixing[ipISO]; 00285 Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipHE_PLUS] *= iso.lgColl_l_mixing[ipISO]; 00286 Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipALPHA] *= iso.lgColl_excite[ipISO]; 00287 } 00288 00289 /* check for sanity */ 00290 ASSERT( Transitions[ipISO][nelem][ipHi][ipLo].Coll.cs >= 0. ); 00291 ASSERT( Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipPROTON] >= 0. ); 00292 ASSERT( Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipHE_PLUS] >= 0. ); 00293 ASSERT( Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipALPHA] >= 0. ); 00294 } 00295 } 00296 } 00297 00298 /* the case b hummer and storey option, 00299 * this kills collisional excitation and ionization from n=1 and n=2 */ 00300 if( opac.lgCaseB_HummerStorey && ipISO==ipH_LIKE ) 00301 { 00302 for( ipLo=0; ipLo<iso.numLevels_max[ipISO][nelem]-1; ipLo++ ) 00303 { 00304 if( N_(ipLo)>=3 ) 00305 break; 00306 00307 iso.ColIoniz[ipISO][nelem][ipLo] = 0.; 00308 00309 for( ipHi=ipLo+1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00310 { 00311 Transitions[ipISO][nelem][ipHi][ipLo].Coll.cs = 0.; 00312 Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipPROTON] = 0.; 00313 Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipHE_PLUS] = 0.; 00314 Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipALPHA] = 0.; 00315 } 00316 } 00317 } 00318 00319 /*********************************************************************** 00320 * * 00321 * collisional deexcitation for all lines * 00322 * * 00323 ***********************************************************************/ 00324 for( ipHi=1; ipHi<iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00325 { 00326 for( ipLo=0; ipLo<ipHi; ipLo++ ) 00327 { 00328 double reduced_mass_proton = dense.AtomicWeight[nelem]*ColliderMass[ipPROTON]/ 00329 (dense.AtomicWeight[nelem]+ColliderMass[ipPROTON])*ATOMIC_MASS_UNIT; 00330 00331 double reduced_mass_heplus = dense.AtomicWeight[nelem]*ColliderMass[ipHE_PLUS]/ 00332 (dense.AtomicWeight[nelem]+ColliderMass[ipHE_PLUS])*ATOMIC_MASS_UNIT; 00333 00334 /******************************************************** 00335 ******************************************************** 00336 * NB - the collision strengths for proton and helium * 00337 * ion impact are multiplied by the ratio of collider * 00338 * to electron densities to precorrect for getting * 00339 * multiplied by the electron density when being put * 00340 * into rate matrix later. They are also multiplied * 00341 * by the pow(m_e/reduced mass)^1.5 to correct for the * 00342 * fact that COLL_CONST contains the electron mass. * 00343 * Here, reduced mass means the reduced mass of the * 00344 * collider-target system. * 00345 ******************************************************** 00346 ********************************************************/ 00347 Transitions[ipISO][nelem][ipHi][ipLo].Coll.ColUL = (realnum)( 00348 ( 00349 /* due to electron impact */ 00350 Transitions[ipISO][nelem][ipHi][ipLo].Coll.cs + 00351 /* due to proton impact */ 00352 Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipPROTON]* 00353 (realnum)(dense.xIonDense[ipHYDROGEN][1]/dense.EdenHCorr)* 00354 pow(ELECTRON_MASS/reduced_mass_proton, 1.5) + 00355 /* due to he+ impact */ 00356 Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipHE_PLUS]* 00357 (realnum)(dense.xIonDense[ipHELIUM][1]/dense.EdenHCorr)* 00358 pow(ELECTRON_MASS/reduced_mass_heplus, 1.5) + 00359 /* due to he++ impact */ 00360 Transitions[ipISO][nelem][ipHi][ipLo].Coll.csi[ipALPHA]* 00361 (realnum)(dense.xIonDense[ipHELIUM][2]/dense.EdenHCorr)* 00362 pow(ELECTRON_MASS/reduced_mass_heplus, 1.5) 00363 ) / phycon.sqrte*COLL_CONST/(double)StatesElem[ipISO][nelem][ipHi].g ); 00364 00365 if( ipISO == ipH_LIKE ) 00366 { 00367 if( N_(ipHi) > iso.n_HighestResolved_max[ipISO][nelem] && 00368 N_(ipLo) <= iso.n_HighestResolved_max[ipISO][nelem] ) 00369 { 00370 Transitions[ipISO][nelem][ipHi][ipLo].Coll.ColUL *= 00371 (8.f/3.f)*(log((realnum)N_(ipHi))+2.f); 00372 } 00373 } 00374 else if( ipISO == ipHE_LIKE ) 00375 { 00376 fixit(); 00377 /* This is intended to be a trick to get the correct collisional excitation from 00378 * collapsed levels to resolved levels. Is it needed or does the stat weight used 00379 * above handle it automatically? If it is needed, is this correct? */ 00380 if( N_(ipHi) > iso.n_HighestResolved_max[ipISO][nelem] && 00381 N_(ipLo) <= iso.n_HighestResolved_max[ipISO][nelem] ) 00382 { 00383 Transitions[ipISO][nelem][ipHi][ipLo].Coll.ColUL *= 00384 (2.f/3.f)*(log((realnum)N_(ipHi))+2.f); 00385 } 00386 } 00387 } 00388 } 00389 00390 if( (trace.lgTrace && trace.lgIsoTraceFull[ipISO]) && (nelem == trace.ipIsoTrace[ipISO]) ) 00391 { 00392 fprintf( ioQQQ, " iso_collide: %s Z=%li de-excitation rates coefficients\n", iso.chISO[ipISO], nelem + 1 ); 00393 for( ipHi=1; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00394 { 00395 fprintf( ioQQQ, " %li\t", ipHi ); 00396 for( ipLo=0; ipLo < ipHi; ipLo++ ) 00397 { 00398 fprintf( ioQQQ,PrintEfmt("%10.3e", Transitions[ipISO][nelem][ipHi][ipLo].Coll.ColUL )); 00399 } 00400 fprintf( ioQQQ, "\n" ); 00401 } 00402 00403 fprintf( ioQQQ, " iso_collide: %s Z=%li collisional ionization coefficients\n", iso.chISO[ipISO], nelem + 1 ); 00404 for( ipHi=0; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00405 { 00406 fprintf( ioQQQ,PrintEfmt("%10.3e", iso.ColIoniz[ipISO][nelem][ipHi] )); 00407 } 00408 fprintf( ioQQQ, "\n" ); 00409 00410 fprintf( ioQQQ, " iso_collide: %s Z=%li continuum boltzmann factor\n", iso.chISO[ipISO], nelem + 1 ); 00411 for( ipHi=0; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00412 { 00413 fprintf( ioQQQ,PrintEfmt("%10.3e", StatesElem[ipISO][nelem][ipHi].ConBoltz )); 00414 } 00415 fprintf( ioQQQ, "\n" ); 00416 00417 fprintf( ioQQQ, " iso_collide: %s Z=%li continuum boltzmann factor\n", iso.chISO[ipISO], nelem + 1 ); 00418 for( ipHi=0; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ ) 00419 { 00420 fprintf( ioQQQ,PrintEfmt("%10.3e", iso.PopLTE[ipISO][nelem][ipHi] )); 00421 } 00422 fprintf( ioQQQ, "\n" ); 00423 } 00424 } 00425 00426 iso_suprathermal( ipISO, nelem ); 00427 00428 /* >>chng 05 aug 17, must use real electron density for collisional ionization of H 00429 * since in Leiden v4 pdr with its artificial high temperature coll ion can be important 00430 * H on H is homonuclear and scaling laws for other elements does not apply 00431 * next two were multiplied by dense.EdenHCorr and changed to dense.eden for H, 00432 * EdenHCorr for rest */ 00433 if( nelem==ipHYDROGEN ) 00434 { 00435 /* special version for H0 onto H0 */ 00436 edenCorr = dense.EdenHontoHCorr; 00437 } 00438 else 00439 { 00440 edenCorr = dense.EdenHCorr; 00441 } 00442 00443 /* this must be reevaluated every time since eden can change when Te does not */ 00444 /* save into main array - collisional ionization by thermal electrons */ 00445 ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0] = 00446 iso.ColIoniz[ipISO][nelem][0]*edenCorr; 00447 00448 /* cooling due to collisional ionization, which only includes thermal electrons */ 00449 ionbal.CollIonRate_Ground[nelem][nelem-ipISO][1] = 00450 ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0]* 00451 rfield.anu[Heavy.ipHeavy[nelem][nelem-ipISO]-1]*EN1RYD; 00452 00453 return; 00454 }