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 /*HeatSum evaluate heating and secondary ionization for current conditions */ 00004 /*HeatZero is called by ConvBase */ 00005 #include "cddefines.h" 00006 #include "physconst.h" 00007 #include "thermal.h" 00008 #include "heavy.h" 00009 #include "trace.h" 00010 #include "secondaries.h" 00011 #include "conv.h" 00012 #include "called.h" 00013 #include "coolheavy.h" 00014 #include "iso.h" 00015 #include "mole.h" 00016 #include "hmi.h" 00017 #include "dense.h" 00018 #include "ionbal.h" 00019 #include "phycon.h" 00020 #include "numderiv.h" 00021 #include "atomfeii.h" 00022 #include "cooling.h" 00023 #include "grainvar.h" 00024 /* this is the faintest relative heating we will print */ 00025 static const double FAINT_HEAT = 0.02; 00026 00027 static const bool PRT_DERIV = false; 00028 00029 void HeatSum( void ) 00030 { 00031 /* use to dim some vectors */ 00032 const int NDIM = 40; 00033 00034 /* secondary ionization and excitation due to Compton scattering */ 00035 double cosmic_ray_ionization_rate , 00036 pair_production_ionization_rate , 00037 fast_neutron_ionization_rate , SecExcitLyaRate; 00038 00039 /* ionization and excitation rates from hydrogen itself */ 00040 double SeconIoniz_iso[NISO] , 00041 SeconExcit_iso[NISO]; 00042 00043 long int i, 00044 ion, 00045 ipnt, 00046 ipsave[NDIM], 00047 j, 00048 jpsave[NDIM], 00049 limit , 00050 nelem; 00051 double HeatingLo , 00052 HeatingHi , 00053 secmet , 00054 smetla; 00055 long ipISO, 00056 ns; 00057 static long int nhit=0, 00058 nzSave=0; 00059 double photo_heat_2lev_atoms, 00060 photo_heat_ISO_atoms , 00061 photo_heat_UTA , 00062 OtherHeat , 00063 deriv, 00064 oldfac, 00065 save[NDIM]; 00066 static double oldheat=0., 00067 oldtemp=0.; 00068 double secmetsave[LIMELM]; 00069 00070 realnum SaveOxygen1 , SaveCarbon1; 00071 00073 realnum xNeutralParticleDensity; 00074 00075 /* routine to sum up total heating, and print agents if needed 00076 * it also computes the true derivative, dH/dT */ 00077 realnum ElectronFraction; 00078 double Cosmic_ray_heat_eff , 00079 Cosmic_ray_sec2prim; 00080 realnum sec2prim_par_1; 00081 realnum sec2prim_par_2; 00082 00083 DEBUG_ENTRY( "HeatSum()" ); 00084 00085 /******************************************************************* 00086 * 00087 * reevaluate the secondary ionization and excitation rates 00088 * 00089 *******************************************************************/ 00090 /* >>chng 03 apr 29, move evaluation of xNeutralParticleDensity to here 00091 * from PresTotCurrent since only used for secondary ionization */ 00092 /* this is total neutral particle density for collisional ionization 00093 * for very high ionization conditions it may be zero */ 00094 xNeutralParticleDensity = 0.; 00095 for( nelem=0; nelem < LIMELM; nelem++ ) 00096 { 00097 xNeutralParticleDensity += dense.xIonDense[nelem][0]; 00098 } 00099 00100 /* now add all the heavy molecules */ 00101 for( i=0; i < mole.num_comole_calc; i++ ) 00102 { 00103 /* add in if real molecule and not counted above */ 00104 if(COmole[i]->n_nuclei > 1) 00105 xNeutralParticleDensity += COmole[i]->hevmol; 00106 } 00107 /* hydrogen molecules */ 00108 xNeutralParticleDensity += hmi.Hmolec[ipMH2p] + hmi.Hmolec[ipMHm] + hmi.Hmolec[ipMH3p] + 00109 (realnum)2.*hmi.H2_total; 00110 00111 { 00112 enum {DEBUG_LOC=false}; 00113 if( DEBUG_LOC ) 00114 { 00115 fprintf(ioQQQ," xIonDense xNeutralParticleDensity tot\t%.3e",xNeutralParticleDensity); 00116 for( nelem=0; nelem < LIMELM; nelem++ ) 00117 { 00118 fprintf(ioQQQ,"\t%.2e",dense.xIonDense[nelem][0]); 00119 } 00120 fprintf(ioQQQ,"\n"); 00121 } 00122 } 00123 00124 00125 /* ElectronFraction is electron fraction 00126 * analytic fits to 00127 * >>>refer sec ioniz Shull & Van Steenberg (1985; Ap.J. 298, 268). 00128 * lgSecOFF turns off secondary ionizations, sets heating efficiency to 100% */ 00129 00130 /* at very low ionization - as per 00131 * >>>refer sec ioniz Xu & McCray 1991, Ap.J. 375, 190. 00132 * everything goes to asymptote that is not present in Shull and 00133 * Van Steenberg - do this by not letting ElecFrac fall below 1e-4 */ 00134 /*ElectronFraction = (realnum)(MAX2(dense.eden/dense.gas_phase[ipHYDROGEN],1e-4));*/ 00135 /* this uses the correct electron density, which will not equal the 00136 * current electron density if we are not converged. Using the 00137 * correct value aids convergence onto it */ 00138 ElectronFraction = (realnum)(MAX2(dense.EdenTrue/dense.gas_phase[ipHYDROGEN],1e-4)); 00139 00140 /* damp out case where electron fraction is oscillating, this 00141 * happens in neutral CR ionized clouds due to feedback between 00142 * ionization efficiency and elecron fraction */ 00143 static realnum OldElectronFraction = 0, 00144 OlderElectronFraction = 0; 00145 if( !conv.nTotalIoniz ) 00146 { 00147 OldElectronFraction = 0; 00148 OlderElectronFraction = 0; 00149 } 00150 if( (ElectronFraction-OldElectronFraction)* 00151 (OldElectronFraction-OlderElectronFraction) < 0. ) 00152 ElectronFraction = ( ElectronFraction+OldElectronFraction)/2.f; 00153 00154 OlderElectronFraction = OldElectronFraction; 00155 OldElectronFraction = ElectronFraction; 00156 00157 if( secondaries.lgSecOFF || ElectronFraction > 0.95 ) 00158 { 00159 secondaries.HeatEfficPrimary = 1.; 00160 secondaries.SecIon2PrimaryErg = 0.; 00161 secondaries.SecExcitLya2PrimaryErg = 0.; 00162 Cosmic_ray_heat_eff = 0.95; 00163 Cosmic_ray_sec2prim = 0.05f; 00164 } 00165 /* >>chng 03 apr 29, only evaluate one time per zone since drove oscillations 00166 * in He+ - He0 ionization front in very high Z models, like hizqso */ 00167 else 00168 { 00169 00170 /* this is heating efficiency for high-energy photo ejections. It is the ratio 00171 * of the final heat added to the energy of the photo electron. Fully 00172 * ionized, this is unity, and less than unity for neutral media. 00173 * It is a scale factor that multiplies the 00174 * high energy heating rate */ 00175 secondaries.HeatEfficPrimary = 0.9971f*(1.f - pow(1.f - pow(ElectronFraction,(realnum)0.2663f),(realnum)1.3163f)); 00176 00177 /* secondary ionizations per primary Rydberg - the number of secondary 00178 * ionizations produced for each Rydberg of high energy photo-electron 00179 * energy deposited. It multiplies the high-energy heating rate. 00180 * It is zero for a fully ionized gas and goes to 0.3908 for very neutral gas */ 00181 secondaries.SecIon2PrimaryErg = 0.3908f*pow(1.f - pow(ElectronFraction,(realnum)0.4092f),(realnum)1.7592f); 00182 /* by dividing by the energy of one Rydberg, converts heating rate 00183 * in ergs into secondary ionization rate */ 00184 secondaries.SecIon2PrimaryErg = secondaries.SecIon2PrimaryErg/((realnum)EN1RYD); 00185 00186 /* This is cosmic ray secondaries per primary (???), 00187 * derived to approximate the curve given in Shull and 00188 * van Steenberg for cosmic rays at 35 eV */ 00189 sec2prim_par_1 = -(1.251f + 195.932f * ElectronFraction); 00190 sec2prim_par_2 = 1.f + 46.814f * ElectronFraction - 44.969f * 00191 ElectronFraction * ElectronFraction; 00192 00193 Cosmic_ray_sec2prim = (exp(sec2prim_par_1/SDIV( sec2prim_par_2))); 00194 00195 /* number of secondary excitations of L-alpha per erg of high 00196 * energy light - actually all Lyman lines 00197 * 00198 * Note--This formula is derived for primary energies greater 00199 * than 100 eV and is only an approximation. This will 00200 * over predict the secondary ionization of L-alpha. We cannot make 00201 * fitting functions for this equation at low energies like we did 00202 * for the heating efficiency and for the number of secondaries 00203 * because the Shull and van Steenberg paper did not publish a similar 00204 * curve for L-alpha */ 00205 secondaries.SecExcitLya2PrimaryErg = 0.4766f*pow(1.f - pow(ElectronFraction,(realnum)0.2735f),(realnum)1.5221f)/((realnum)EN1RYD); 00206 00207 if( (trace.lgTrace && trace.lgSecIon) ) 00208 { 00209 fprintf( ioQQQ, 00210 " csupra H0 %.2e HeatSum eval sec ion effic, ElectronFraction = %.3e HeatEfficPrimary %.3e SecIon2PrimaryErg %.3e\n", 00211 secondaries.csupra[ipHYDROGEN][0], 00212 ElectronFraction, 00213 secondaries.HeatEfficPrimary , 00214 secondaries.SecIon2PrimaryErg ); 00215 } 00216 00217 /* cosmic ray heating, counts as non-ionizing heat since already result of cascade, 00218 * this was 35 eV per secondary ionization */ 00219 /* We want to use the heating efficiency that is applicable to a 35 eV 00220 * primary electron, the current efficiency used is for 100eV cosmic ray 00221 * Here we use the correct value as given in Wolfire et al. 1995 */ 00222 00223 /* In general the equation for the cosmic ray heating rate is:* 00224 * * 00225 * * 00226 * CR_Heat_Rate = (density)*(cosmic ray ionization rate)* * 00227 * (energy of primary electron)* * 00228 * (Heating efficiency at that energy) * 00229 * * 00230 * The product of the last two terms gives the amount of heat * 00231 * added by the primary electron, and is dependent upon the * 00232 * electron fraction. We are using the same average primary * 00233 * electron energy as Wolfire et al. (1995). We do not, * 00234 * however, use their formula for the heating efficiency. * 00235 * This is because the coefficients (f1, f2, and f3) were * 00236 * originally intended for primary electron energies >100eV. * 00237 * Instead we derived a heating efficiency based on the curves* 00238 * given in Shull and van Steenberg (1985). We interpolated * 00239 * for an energy of 35 eV the heating efficiency for electron * 00240 * fractions between 1e-4 and 1 * 00241 **************************************************************/ 00242 00243 /*printf("Here is ElectronFraction:\t%.3e\n", ElectronFraction);*/ 00244 Cosmic_ray_heat_eff = - 8.189 - 18.394 * ElectronFraction - 6.608 * ElectronFraction * ElectronFraction * log(ElectronFraction) 00245 + 8.322 * exp( ElectronFraction ) + 4.961 * sqrt(ElectronFraction); 00246 } 00247 00248 /******************************************************************* 00249 * 00250 * get total heating from all species 00251 * 00252 *******************************************************************/ 00253 00254 /* get total heating */ 00255 photo_heat_2lev_atoms = 0.; 00256 photo_heat_ISO_atoms = 0.; 00257 photo_heat_UTA = 0.; 00258 00259 /* add in effects of high-energy opacity of CO using C and O atomic opacities 00260 * k-shell and valence photo of C and O in CO is not explicitly counted elsewhere 00261 * this trick roughly accounts for it*/ 00262 SaveOxygen1 = dense.xIonDense[ipOXYGEN][0]; 00263 SaveCarbon1 = dense.xIonDense[ipCARBON][0]; 00264 00265 /* atomic C and O will include CO during the heating sum loop */ 00266 dense.xIonDense[ipOXYGEN][0] += findspecies("CO")->hevmol; 00267 dense.xIonDense[ipCARBON][0] += findspecies("CO")->hevmol; 00268 00269 /* this will hold cooling due to metal collisional ionization */ 00270 CoolHeavy.colmet = 0.; 00271 /* metals secondary ionization, Lya excitation */ 00272 secmet = 0.; 00273 smetla = 0.; 00274 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00275 { 00276 SeconIoniz_iso[ipISO] = 0.; 00277 SeconExcit_iso[ipISO] = 0.; 00278 } 00279 00280 /* this loop includes hydrogenic, which is special case due to presence 00281 * of substantial excited state populations */ 00282 for( nelem=0; nelem<LIMELM; ++nelem) 00283 { 00284 secmetsave[nelem] = 0.; 00285 if( dense.lgElmtOn[nelem] ) 00286 { 00287 /* sum heat over all stages of ionization that exist */ 00288 /* first do the iso sequences, 00289 * h-like and he-like are special because full atom always done, 00290 * can be substantial, pops in excited states, with little in ground 00291 * (true near LTE), these are done in following loop */ 00292 00293 limit = MIN2( dense.IonHigh[nelem] , nelem+1-NISO ); 00294 00295 /* loop over all elements/ions done with as two-level systems */ 00296 for( ion=dense.IonLow[nelem]; ion < limit; ion++ ) 00297 { 00298 /* this will be heating for a single stage of ionization */ 00299 HeatingLo = 0.; 00300 HeatingHi = 0.; 00301 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ ) 00302 { 00303 /* heating by various sub-shells */ 00304 HeatingLo += ionbal.PhotoRate_Shell[nelem][ion][ns][1]; 00305 HeatingHi += ionbal.PhotoRate_Shell[nelem][ion][ns][2]; 00306 } 00307 00308 /* total photoelectric heating, all shells, for this stage 00309 * of ionization */ 00310 thermal.heating[nelem][ion] = dense.xIonDense[nelem][ion]* 00311 (HeatingLo + HeatingHi*secondaries.HeatEfficPrimary); 00312 /* heating due to UTA ionization */ 00313 photo_heat_UTA += dense.xIonDense[nelem][ion]* 00314 ionbal.UTA_heat_rate[nelem][ion]; 00315 00316 /* add to total heating */ 00317 photo_heat_2lev_atoms += thermal.heating[nelem][ion]; 00318 /*if( nzone>290 && thermal.heating[nelem][ion]/thermal.htot>0.01 ) 00319 fprintf(ioQQQ,"buggg\t%li %li %.3f\n", nelem,ion,thermal.heating[nelem][ion]/thermal.htot);*/ 00320 00321 /* Cooling due to collisional ionization of heavy elements by thermal electrons 00322 * CollidRate[nelem][ion][1] is cooling, erg/s/atom, eval in ion_collis */ 00323 CoolHeavy.colmet += dense.xIonDense[nelem][ion]*ionbal.CollIonRate_Ground[nelem][ion][1]; 00324 00325 /* secondary ionization rate, */ 00326 secmetsave[nelem] += HeatingHi*secondaries.SecIon2PrimaryErg* dense.xIonDense[nelem][ion]; 00327 00328 /* LA excitation rate, =0 if ionized, s-1 cm-3 */ 00329 smetla += HeatingHi*secondaries.SecExcitLya2PrimaryErg* dense.xIonDense[nelem][ion]; 00330 } 00331 secmet += secmetsave[nelem]; 00332 00333 /* this branch loop over all ions done with full iso sequence */ 00334 limit = MAX2( limit, dense.IonLow[nelem] ); 00335 for( ion=MAX2(0,limit); ion < dense.IonHigh[nelem]; ion++ ) 00336 { 00337 /* this is the iso sequence */ 00338 ipISO = nelem-ion; 00339 /* this will be heating for a single stage of ionization */ 00340 HeatingLo = 0.; 00341 HeatingHi = 0.; 00342 /* the outer shell contains the Compton recoil part */ 00343 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ ) 00344 { 00345 /* heating, erg s-1 atom-1, by low energy, and then high energy, photons */ 00346 HeatingLo += ionbal.PhotoRate_Shell[nelem][ion][ns][1]; 00347 HeatingHi += ionbal.PhotoRate_Shell[nelem][ion][ns][2]; 00348 } 00349 00350 /* net heating, erg cm-3 s-1 */ 00351 thermal.heating[nelem][ion] = dense.xIonDense[nelem][ion+1]* 00352 StatesElem[ipISO][nelem][0].Pop*(HeatingLo + HeatingHi*secondaries.HeatEfficPrimary); 00353 00354 /* heating due to UTA ionization, erg cm-3 s-1 */ 00355 photo_heat_UTA += dense.xIonDense[nelem][ion]* 00356 ionbal.UTA_heat_rate[nelem][ion]; 00357 00358 /* add to total heating */ 00359 photo_heat_ISO_atoms += thermal.heating[nelem][ion]; 00360 00361 if( secondaries.SecIon2PrimaryErg > SMALLFLOAT && xNeutralParticleDensity > SMALLFLOAT ) 00362 { 00363 /* prevent crash in very high ionization conditions 00364 * where xNeutralParticleDensity is zero */ 00365 /* secondary ionization rate, */ 00366 SeconIoniz_iso[ipISO] += HeatingHi*secondaries.SecIon2PrimaryErg* 00367 StatesElem[ipISO][nelem][0].Pop*dense.xIonDense[nelem][ion+1]/ 00368 xNeutralParticleDensity; 00369 00370 /* LA excitation rate, =0 if ionized, units excitations s-1 */ 00371 SeconExcit_iso[ipISO] += HeatingHi*secondaries.SecExcitLya2PrimaryErg* 00372 StatesElem[ipISO][nelem][0].Pop*dense.xIonDense[nelem][ion+1]/ 00373 xNeutralParticleDensity; 00374 00375 ASSERT( SeconIoniz_iso[ipISO]>=0. && 00376 SeconExcit_iso[ipISO]>=0.); 00377 } 00378 } 00379 00380 /* make sure stages of ionization with no abundances, 00381 * also has no heating */ 00382 for( ion=0; ion<dense.IonLow[nelem]; ion++ ) 00383 { 00384 ASSERT( thermal.heating[nelem][ion] == 0. ); 00385 } 00386 for( ion=dense.IonHigh[nelem]+1; ion<nelem+1; ion++ ) 00387 { 00388 ASSERT( thermal.heating[nelem][ion] == 0. ); 00389 } 00390 } 00391 } 00392 if( trace.lgTrace && trace.lgSecIon ) 00393 { 00394 double savemax=0.; 00395 long int ipsavemax=-1; 00396 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem ) 00397 { 00398 if( secmetsave[nelem] > savemax ) 00399 { 00400 savemax = secmetsave[nelem]; 00401 ipsavemax = nelem; 00402 } 00403 } 00404 fprintf( ioQQQ, 00405 " HeatSum secmet largest contributor element %li frac of total %.3e, total %.3e\n", 00406 ipsavemax, 00407 savemax/SDIV(secmet), 00408 secmet); 00409 } 00410 00411 /* now reset the abundances */ 00412 dense.xIonDense[ipOXYGEN][0] = SaveOxygen1; 00413 dense.xIonDense[ipCARBON][0] = SaveCarbon1; 00414 00415 /* convert secmet to proper final units */ 00416 /*fprintf(ioQQQ,"bugggg\t%li\t%.3e\t%.3e\t%.3e\n", nzone , 00417 secmet , xNeutralParticleDensity , secmet / xNeutralParticleDensity );*/ 00418 /* prevent crash when xNeutralParticleDensity is zero */ 00419 if( secondaries.SecIon2PrimaryErg > SMALLFLOAT && xNeutralParticleDensity > SMALLFLOAT ) 00420 { 00421 /* convert from s-1 cm-3 to s-1 - a true rate */ 00422 secmet /= xNeutralParticleDensity; 00423 smetla /= xNeutralParticleDensity; 00424 } 00425 else 00426 { 00427 /* zero */ 00428 secmet = 0.; 00429 smetla = 0.; 00430 } 00431 00432 /* >>chng 01 dec 20, do full some over all secondaries */ 00433 /* bound Compton recoil heating */ 00434 /* >>chng 02 mar 28, save heating in this var rather than heating[0][18] 00435 * since now saved in photo heat 00436 * this is only used for a printout and in lines, not as heat since already counted*/ 00437 ionbal.CompRecoilHeatLocal = 0.; 00438 for( nelem=0; nelem<LIMELM; ++nelem ) 00439 { 00440 for( ion=0; ion<nelem+1; ++ion ) 00441 { 00442 ionbal.CompRecoilHeatLocal += 00443 ionbal.CompRecoilHeatRate[nelem][ion]*secondaries.HeatEfficPrimary*dense.xIonDense[nelem][ion]; 00444 } 00445 } 00446 /* >>chng 05 nov 26, include this term - bound ionization of H2 00447 * assume total cs is that of two separated H */ 00448 ionbal.CompRecoilHeatLocal += 00449 2.*ionbal.CompRecoilHeatRate[ipHYDROGEN][0]*secondaries.HeatEfficPrimary*hmi.H2_total; 00450 00451 /* find heating due to charge transfer 00452 * >>chng 05 oct 29, move from here to conv base so that can be treated as cooling if 00453 * negative heating */ 00454 thermal.heating[0][24] = thermal.char_tran_heat; 00455 00456 /* heating due to pair production */ 00457 thermal.heating[0][21] = 00458 ionbal.PairProducPhotoRate[2]*secondaries.HeatEfficPrimary*(dense.gas_phase[ipHYDROGEN] + 4.F*dense.gas_phase[ipHELIUM]); 00459 /* last term above is number of nuclei in helium */ 00460 00461 /* this is heating due to fast neutrons, assumed to secondary ionize */ 00462 thermal.heating[0][20] = 00463 ionbal.xNeutronHeatRate*secondaries.HeatEfficPrimary*dense.gas_phase[ipHYDROGEN]; 00464 00465 /* turbulent heating, assumed to be a non-ionizing heat agent, added here */ 00466 thermal.heating[0][20] += ionbal.ExtraHeatRate; 00467 00468 /* UTA heating */ 00469 thermal.heating[0][7] = photo_heat_UTA; 00470 /*fprintf(ioQQQ,"DEBUG UTA heat %.3e\n", photo_heat_UTA/SDIV(thermal.htot ));*/ 00471 00472 /* >>chng 05 nov 26, include heating due to H2 photoionization */ 00473 /*>>KEYWORD H2 photoionization heating */ 00474 thermal.heating[0][18] = hmi.H2_total * 00475 (hmi.H2_photo_heat_soft + hmi.H2_photo_heat_hard*secondaries.HeatEfficPrimary); 00478 /* >>chng 05 nov 27, approximate heating due to H2+, H3+ photoionization 00479 * assuming H0 rates */ 00480 /*>>KEYWORD H2+ photoionization heating; H3+ photoionization heating */ 00481 thermal.heating[0][26] = (hmi.Hmolec[ipMH2p]+hmi.Hmolec[ipMH3p]) * 00482 (ionbal.PhotoRate_Shell[ipHYDROGEN][0][0][1] + 00483 ionbal.PhotoRate_Shell[ipHYDROGEN][0][0][2]*secondaries.HeatEfficPrimary); 00484 00485 /* add on cosmic rays - most important in neutral or molecular gas, use 00486 * difference between total H and H+ density as substitute for sum of 00487 * H in atoms and all molecules, but only in limit where difference is 00488 * significant */ 00489 # if 0 00490 double hneut; 00491 if( (dense.gas_phase[ipHYDROGEN] - dense.xIonDense[ipHYDROGEN][1])/ 00492 dense.gas_phase[ipHYDROGEN]<0.001 ) 00493 { 00494 /* limit where most H is ionized - simply use atomic and H2 */ 00495 hneut = dense.xIonDense[ipHYDROGEN][0] + 2.*(hmi.Hmolec[ipMH2g]+hmi.Hmolec[ipMH2s]); 00496 } 00497 else 00498 { 00499 /* limit where neutral is significant, use different between total and H+ */ 00500 hneut = dense.gas_phase[ipHYDROGEN] - dense.xIonDense[ipHYDROGEN][1]; 00501 } 00502 # endif 00503 00504 /* cosmic ray heating */ 00505 thermal.heating[1][6] = 00506 ionbal.CosRayHeatNeutralParticles* 00507 xNeutralParticleDensity * Cosmic_ray_heat_eff; 00508 /* cosmic ray heating of thermal electrons */ 00509 /* >>refer CR physics Ferland, G.J. & Mushotzky, R.F., 1984, ApJ, 286, 42 */ 00510 thermal.heating[1][6] += ionbal.CosRayHeatThermalElectrons * dense.eden; 00511 00512 /* now sum up all heating agents not included in sum over normal bound levels above */ 00513 OtherHeat = 0.; 00514 for( nelem=0; nelem<LIMELM; ++nelem) 00515 { 00516 /* we now have ionization solution for this element,sum heat over 00517 * all stages of ionization that exist */ 00518 /* >>>chng 99 may 08, following loop had started at nelem+3, and so missed [1][0], 00519 * which is excited states of hydrogenic species. increase this limit */ 00520 for( i=nelem+1; i < LIMELM; i++ ) 00521 { 00522 OtherHeat += thermal.heating[nelem][i]; 00523 } 00524 } 00525 00526 thermal.htot = OtherHeat + photo_heat_2lev_atoms + photo_heat_ISO_atoms; 00527 00528 /* following checks whether total heating is strange, if we are not in search phase */ 00529 if( called.lgTalk && !conv.lgSearch ) 00530 { 00531 /* print this warning if not constant temperature and neg heat */ 00532 if( thermal.htot < 0. && !thermal.lgTemperatureConstant) 00533 { 00534 fprintf( ioQQQ, 00535 " Total heating is <0; is this model collisionally ionized? zone is %li\n", 00536 nzone ); 00537 } 00538 else if( thermal.htot == 0. ) 00539 { 00540 fprintf( ioQQQ, 00541 " Total heating is 0; is the density small? zone is %li\n", 00542 nzone); 00543 } 00544 } 00545 00546 /* add on line heating to this array, heatl was evaluated in sumcool 00547 * not zero, because of roundoff error */ 00548 if( thermal.heatl/thermal.ctot < -1e-15 ) 00549 { 00550 fprintf( ioQQQ, " HeatSum gets negative line heating,%10.2e%10.2e this is insane.\n", 00551 thermal.heatl, thermal.ctot ); 00552 00553 fprintf( ioQQQ, " this is zone%4ld\n", nzone ); 00554 ShowMe(); 00555 cdEXIT(EXIT_FAILURE); 00556 } 00557 00558 /******************************************************************* 00559 * 00560 * secondary ionization and excitation rates 00561 * 00562 *******************************************************************/ 00563 00564 /* the terms cosmic_ray_ionization_rate & SecExcitLyaRate contain energies added in highen. 00565 * The only significant one is usually bound Compton heating except when 00566 * cosmic rays are present */ 00567 00568 if( secondaries.SecIon2PrimaryErg > SMALLFLOAT && xNeutralParticleDensity > SMALLFLOAT ) 00569 { 00570 realnum DensityRatio = (dense.gas_phase[ipHYDROGEN] + 00571 4.F*dense.gas_phase[ipHELIUM])/xNeutralParticleDensity; 00572 00573 /* add on ionization rate s-1 cm-3 due to pair production */ 00574 pair_production_ionization_rate = 00575 ionbal.PairProducPhotoRate[2]*secondaries.SecIon2PrimaryErg*DensityRatio; 00576 00577 /* total secondary excitation rate cm-3 s-1 for Lya due to pair production*/ 00578 SecExcitLyaRate = 00579 ionbal.PairProducPhotoRate[2]*secondaries.SecExcitLya2PrimaryErg*1.3333*DensityRatio; 00580 00581 /* ionization rate of fast neutrons */ 00582 fast_neutron_ionization_rate = 00583 ionbal.xNeutronHeatRate*secondaries.SecIon2PrimaryErg*DensityRatio; 00584 00585 /* secondary excitation rate due to neutrons */ 00586 SecExcitLyaRate += 00587 ionbal.xNeutronHeatRate*secondaries.SecExcitLya2PrimaryErg*1.3333*DensityRatio; 00588 00589 /* cosmic ray Lya excitation rate */ 00590 SecExcitLyaRate += 00591 /* multiply by atomic H and He densities */ 00592 ionbal.CosRayHeatNeutralParticles*secondaries.SecExcitLya2PrimaryErg*1.3333*DensityRatio; 00593 } 00594 else 00595 { 00596 /* zero */ 00597 pair_production_ionization_rate = 0.; 00598 SecExcitLyaRate = 0.; 00599 fast_neutron_ionization_rate = 0.; 00600 } 00601 00602 /* cosmic ray ionization */ 00603 /* >>>chng 99 apr 29, term in PhotoRate was not present */ 00604 cosmic_ray_ionization_rate = 00605 /* this term is cosmic ray ionization, set in highen, did not multiply by 00606 * collider density in highen, so do not divide by it here */ 00607 /* >>chng 99 jun 29, added SecIon2PrimaryErg to cosmic ray rate, also multiply 00608 * by number of secondaries per primary*/ 00609 ionbal.CosRayIonRate*Cosmic_ray_sec2prim + 00610 /* this is the cosmic ray heating rate */ 00611 ionbal.CosRayHeatNeutralParticles*secondaries.SecIon2PrimaryErg; 00612 00613 /* find total suprathermal collisional ionization rate 00614 * CSUPRA is H0 col ionize rate from non-thermal electrons (inverse sec) 00615 * x12tot is LA excitation rate, units s-1 00616 * SECCMP evaluated in HIGHEN, =ioniz rate for cosmic rays, sec bound */ 00617 00618 /* option to force secondary ionization rate, normally false */ 00619 if( secondaries.lgCSetOn ) 00620 { 00621 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem ) 00622 { 00623 for( ion=0; ion<nelem+1; ++ion ) 00624 { 00625 secondaries.csupra[nelem][ion] = secondaries.SetCsupra*secondaries.csupra_effic[nelem][ion]; 00626 } 00627 } 00628 secondaries.x12tot = secondaries.SetCsupra; 00629 } 00630 else 00631 { 00632 double csupra; 00633 double facold , facnew; 00634 /* xNeutralParticleDensity is total neutral particle density */ 00635 if( secondaries.csupra[ipHYDROGEN][0] / SDIV( iso.RateLevel2Cont[ipH_LIKE][ipHYDROGEN][ipH1s] ) > 0.1 ) 00636 { 00637 /* supra are dominant H destruction - make small changes */ 00638 facold = 0.9; 00639 } 00640 else 00641 { 00642 /* secondaries are not important - only use new */ 00643 facold = 0.; 00644 } 00645 facnew = 1. - facold; 00646 csupra = (secondaries.csupra[ipHYDROGEN][0]* facold + facnew* 00647 (cosmic_ray_ionization_rate + 00648 pair_production_ionization_rate + 00649 fast_neutron_ionization_rate + 00650 SeconIoniz_iso[ipH_LIKE] + 00651 SeconIoniz_iso[ipHE_LIKE] + 00652 secmet )); 00653 00655 /* now fill in ionization rates for all elements and ion stages */ 00656 for( nelem=ipHYDROGEN; nelem<LIMELM; ++nelem ) 00657 { 00658 for( ion=0; ion<nelem+1; ++ion ) 00659 { 00660 secondaries.csupra[nelem][ion] = (realnum)csupra*secondaries.csupra_effic[nelem][ion]; 00661 } 00662 } 00663 00664 /* secondary suprathermal excitation of Ly-alpha, excitations s-1 */ 00665 secondaries.x12tot = (realnum)( secondaries.x12tot*facold + facnew* 00666 /* these terms have units excitations s-1 */ 00667 (SecExcitLyaRate + 00668 SeconExcit_iso[ipH_LIKE] + 00669 SeconExcit_iso[ipHE_LIKE] + 00670 smetla)); 00671 } 00672 00673 if( trace.lgTrace && secondaries.csupra[ipHYDROGEN][0] > 0. ) 00674 { 00675 fprintf( ioQQQ, 00676 " HeatSum return CSUPRA %9.2e SECCMP %6.3f SecHI %6.3f SECHE %6.3f SECMET %6.3f efrac= %9.2e \n", 00677 secondaries.csupra[ipHYDROGEN][0], 00678 (cosmic_ray_ionization_rate+pair_production_ionization_rate+fast_neutron_ionization_rate)/secondaries.csupra[ipHYDROGEN][0], 00679 SeconIoniz_iso[ipH_LIKE] / secondaries.csupra[ipHYDROGEN][0], 00680 SeconIoniz_iso[ipHE_LIKE]/secondaries.csupra[ipHYDROGEN][0], 00681 secmet/secondaries.csupra[ipHYDROGEN][0] , 00682 ElectronFraction ); 00683 } 00684 00685 /******************************************************************* 00686 * 00687 * get derivative of total heating 00688 * 00689 *******************************************************************/ 00690 00691 /* now get derivative of heating due to photoionization, 00692 * >>chng 96 jan 14 00693 *>>>>>NB heating decreasing with increasing temp is negative dH/dT */ 00694 thermal.dHeatdT = -0.7*(photo_heat_2lev_atoms+photo_heat_ISO_atoms+ 00695 photo_heat_UTA)/phycon.te; 00696 /* >>chng 04 feb 28, add this correction factor - when gas totally neutral heating 00697 * does not depend on temperature - when ionized depends on recom coefficient - this 00698 * smoothly joins two limits */ 00699 thermal.dHeatdT *= dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN]; 00700 if( PRT_DERIV ) 00701 fprintf(ioQQQ,"DEBUG dhdt 0 %.3e %.3e %.3e \n", 00702 thermal.dHeatdT, 00703 photo_heat_2lev_atoms, 00704 photo_heat_ISO_atoms); 00705 00706 /* oldfac was factor used in old implementation 00707 * any term using it should be rethought */ 00708 oldfac = -0.7/phycon.te; 00709 00710 /* carbon monoxide */ 00711 thermal.dHeatdT += thermal.heating[0][9]*oldfac; 00712 00713 /* Ly alpha collisional heating */ 00714 thermal.dHeatdT += thermal.heating[0][10]*oldfac; 00715 00716 /* line heating */ 00717 thermal.dHeatdT += thermal.heating[0][22]*oldfac; 00718 if( PRT_DERIV ) 00719 fprintf(ioQQQ,"DEBUG dhdt 1 %.3e \n", thermal.dHeatdT); 00720 00721 /* free free heating, propto T^-0.5 00722 * >>chng 96 oct 30, from heating(1,12) to CoolHeavy.brems_heat_total, 00723 * do cooling separately assume CoolHeavy.brems_heat_total goes as T^-3/2 00724 * dHTotDT = dHTotDT + heating(1,12) * (-0.5/te) */ 00725 thermal.dHeatdT += CoolHeavy.brems_heat_total*(-1.5/phycon.te); 00726 00727 /* >>chng 04 aug 07, use better estimate for heating derivative; needed in PDRs, PvH */ 00728 /* this includes PE, thermionic, and collisional heating of the gas by the grains */ 00729 thermal.dHeatdT += gv.dHeatdT; 00730 00731 /* helium triplets heating */ 00732 thermal.dHeatdT += thermal.heating[1][2]*oldfac; 00733 if( PRT_DERIV ) 00734 fprintf(ioQQQ,"DEBUG dhdt 2 %.3e \n", thermal.dHeatdT); 00735 00736 /* hydrogen molecule collisional deexcitation heating */ 00737 /* >>chng 04 jan 25, add max to prevent cooling from entering here */ 00738 /*thermal.dHeatdT += MAX2(0.,thermal.heating[0][8])*oldfac;*/ 00739 if( hmi.HeatH2Dexc_used > 0. ) 00740 thermal.dHeatdT += hmi.deriv_HeatH2Dexc_used; 00741 00742 /* >>chng 96 nov 15, had been oldfac below, wrong sign 00743 * H- H minus heating, goes up with temp since rad assoc does */ 00744 thermal.dHeatdT += thermal.heating[0][15]*0.7/phycon.te; 00745 00746 /* H2+ heating */ 00747 thermal.dHeatdT += thermal.heating[0][16]*oldfac; 00748 00749 /* Balmer continuum and all other excited states 00750 * - T goes up so does pop and heating 00751 * but this all screwed up by change in eden */ 00752 thermal.dHeatdT += thermal.heating[0][1]*oldfac; 00753 if( PRT_DERIV ) 00754 fprintf(ioQQQ,"DEBUG dhdt 3 %.3e \n", thermal.dHeatdT); 00755 00756 /* all three body heating, opposite of coll ion cooling */ 00757 thermal.dHeatdT += thermal.heating[0][3]*oldfac; 00758 00759 /* bound electron recoil heating 00760 thermal.dHeatdT += ionbal.CompRecoilHeatLocal*oldfac; */ 00761 00762 /* Compton heating */ 00763 thermal.dHeatdT += thermal.heating[0][19]*oldfac; 00764 00765 /* extra heating source, usually zero */ 00766 thermal.dHeatdT += thermal.heating[0][20]*oldfac; 00767 00768 /* pair production */ 00769 thermal.dHeatdT += thermal.heating[0][21]*oldfac; 00770 if( PRT_DERIV ) 00771 fprintf(ioQQQ,"DEBUG dhdt 4 %.3e \n", thermal.dHeatdT); 00772 00773 /* derivative of heating due to collisions of H lines, heating(1,24) */ 00774 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO ) 00775 { 00776 for( nelem=ipISO; nelem<LIMELM; ++nelem) 00777 { 00778 thermal.dHeatdT += iso.dLTot[ipISO][nelem]; 00779 } 00780 } 00781 00782 /* heating due to large FeII atom, zero unless atom FeII is entered, 00783 * FeII.Fe2_large_heat was entered into thermal.heating[25][27] */ 00784 if( FeII.Fe2_large_heat > 0. ) 00785 { 00786 thermal.dHeatdT += FeII.ddT_Fe2_large_cool; 00787 } 00788 if( PRT_DERIV ) 00789 fprintf(ioQQQ,"DEBUG dhdt 5 %.3e \n", thermal.dHeatdT); 00790 00791 /* possibility of getting empirical heating derivative 00792 * normally false, set true with 'set numerical derivatives' command */ 00793 if( NumDeriv.lgNumDeriv ) 00794 { 00795 if( ((nzone > 2 && nzone == nzSave) && 00796 ! fp_equal( oldtemp, phycon.te )) && nhit > 4 ) 00797 { 00798 /* hnit is number of tries on this zone - use to stop numerical problems 00799 * do not evaluate numerical derivative until well into solution */ 00800 deriv = (oldheat - thermal.htot)/(oldtemp - phycon.te); 00801 thermal.dHeatdT = deriv; 00802 } 00803 else 00804 { 00805 deriv = thermal.dHeatdT; 00806 } 00807 00808 /* this happens when new zone starts */ 00809 if( nzone != nzSave ) 00810 { 00811 nhit = 0; 00812 } 00813 nzSave = nzone; 00814 nhit += 1; 00815 oldheat = thermal.htot; 00816 oldtemp = phycon.te; 00817 } 00818 00819 if( trace.lgTrace && trace.lgHeatBug ) 00820 { 00821 ipnt = 0; 00822 /* this loops through the 2D array that contains all agents counted in htot */ 00823 for( i=0; i < LIMELM; i++ ) 00824 { 00825 for( j=0; j < LIMELM; j++ ) 00826 { 00827 if( thermal.heating[i][j]/thermal.htot > FAINT_HEAT ) 00828 { 00829 ipsave[ipnt] = i; 00830 jpsave[ipnt] = j; 00831 save[ipnt] = thermal.heating[i][j]/thermal.htot; 00832 /* increment ipnt, but do not let it get too big */ 00833 ipnt = MIN2((long)NDIM-1,ipnt+1); 00834 } 00835 } 00836 } 00837 00838 /* now check for possible line heating not counted in 1,23 00839 * there should not be any significant heating source here 00840 * since they would not be counted in derivative correctly */ 00841 for( i=0; i < thermal.ncltot; i++ ) 00842 { 00843 if( thermal.heatnt[i]/thermal.htot > FAINT_HEAT ) 00844 { 00845 ipsave[ipnt] = -1; 00846 jpsave[ipnt] = i; 00847 save[ipnt] = thermal.heatnt[i]/thermal.htot; 00848 ipnt = MIN2((long)NDIM-1,ipnt+1); 00849 } 00850 } 00851 00852 fprintf( ioQQQ, 00853 " HeatSum HTOT %.3e Te:%.3e dH/dT%.2e other %.2e photo 2lev %.2e photo iso %.2e\n", 00854 thermal.htot, 00855 phycon.te, 00856 thermal.dHeatdT , 00857 /* total heating is sum of following three terms 00858 * OtherHeat is a sum over all other heating agents */ 00859 OtherHeat , 00860 photo_heat_2lev_atoms, 00861 photo_heat_ISO_atoms); 00862 00863 fprintf( ioQQQ, " " ); 00864 for( i=0; i < ipnt; i++ ) 00865 { 00866 /*lint -e644 these three are initialized above */ 00867 fprintf( ioQQQ, " [%ld][%ld]%6.3f", 00868 ipsave[i], 00869 jpsave[i], 00870 save[i] ); 00871 /*lint +e644 these three are initialized above */ 00872 /* throw a cr every n numbers */ 00873 if( !(i%8) && i>0 && i!=(ipnt-1) ) 00874 { 00875 fprintf( ioQQQ, "\n " ); 00876 } 00877 } 00878 fprintf( ioQQQ, "\n" ); 00879 } 00880 return; 00881 } 00882 00883 /* =============================================================================*/ 00884 /* HeatZero is called by ConvBase */ 00885 void HeatZero( void ) 00886 { 00887 long int i , j; 00888 00889 DEBUG_ENTRY( "HeatZero()" ); 00890 00891 for( i=0; i < LIMELM; i++ ) 00892 { 00893 for( j=0; j < LIMELM; j++ ) 00894 { 00895 thermal.heating[i][j] = 0.; 00896 } 00897 } 00898 return; 00899 }