cloudy trunk
Loading...
Searching...
No Matches
heat_sum.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
2 * others. For conditions of distribution and use see copyright notice in license.txt */
3/*HeatSum evaluate heating and secondary ionization for current conditions */
4/*HeatZero is called by ConvBase */
5#include "cddefines.h"
6#include "physconst.h"
7#include "thermal.h"
8#include "heavy.h"
9#include "trace.h"
10#include "secondaries.h"
11#include "conv.h"
12#include "called.h"
13#include "coolheavy.h"
14#include "iso.h"
15#include "mole.h"
16#include "h2.h"
17#include "hmi.h"
18#include "dense.h"
19#include "ionbal.h"
20#include "phycon.h"
21#include "numderiv.h"
22#include "atomfeii.h"
23#include "cooling.h"
24#include "grainvar.h"
25#include "taulines.h"
26#include "deuterium.h"
27
28/* this is the faintest relative heating we will print */
29static const double FAINT_HEAT = 0.02;
30
31static const bool PRT_DERIV = false;
32
33void HeatSum( void )
34{
35 /* use to dim some vectors */
36 const int NDIM = 40;
37
38 /* secondary ionization and excitation due to Compton scattering */
39 double cosmic_ray_ionization_rate ,
40 pair_production_ionization_rate ,
41 fast_neutron_ionization_rate , SecExcitLyaRate;
42
43 /* ionization and excitation rates from hydrogen itself */
44 double SeconIoniz_iso[NISO] ,
45 SeconExcit_iso[NISO];
46
47 long int i,
48 ion,
49 ipnt,
50 ipsave[NDIM],
51 j,
52 jpsave[NDIM],
53 limit;
54
55 double HeatingLo ,
56 HeatingHi ,
57 secmet ,
58 smetla;
59 long ipISO,
60 ns;
61 static long int nhit=0,
62 nzSave=0;
63 double photo_heat_2lev_atoms,
64 photo_heat_ISO_atoms ,
65 photo_heat_UTA ,
66 OtherHeat ,
67 deriv,
68 oldfac,
69 save[NDIM];
70 static double oldheat=0.,
71 oldtemp=0.;
72 double secmetsave[LIMELM];
73
74 realnum SaveOxygen1 , SaveCarbon1;
75
76 /* routine to sum up total heating, and print agents if needed
77 * it also computes the true derivative, dH/dT */
78 double Cosmic_ray_heat_eff ,
79 Cosmic_ray_sec2prim;
80 realnum sec2prim_par_1;
81 realnum sec2prim_par_2;
82
83 DEBUG_ENTRY( "HeatSum()" );
84
85 /*******************************************************************
86 *
87 * reevaluate the secondary ionization and excitation rates
88 *
89 *******************************************************************/
90 /* ElectronFraction is electron fraction
91 * analytic fits to
92 * >>>refer sec ioniz Shull & Van Steenberg (1985; Ap.J. 298, 268).
93 * lgSecOFF turns off secondary ionizations, sets heating efficiency to 100% */
94
95 /* at very low ionization - as per
96 * >>>refer sec ioniz Xu & McCray 1991, Ap.J. 375, 190.
97 * everything goes to asymptote that is not present in Shull and
98 * Van Steenberg - do this by not letting ElecFrac fall below 1e-4 */
99 /* this uses the correct electron density, which will not equal the
100 * current electron density if we are not converged. Using the
101 * correct value aids convergence onto it */
102
103 // ElectronFraction should be the fraction of electrons available
104 // for collisions which are free, as opposed to bound in atoms,
105 // ions or molecules, which seems like the most plausible
106 // generalization of the definition of X in S&vS.
107
108 realnum xValenceDonorDensity, ElectronFraction;
109
110 if (1)
111 {
112 long ipNoble[] = {
114 /*, ipXENON, ipRADON, ipUNUNOCTIUM */
115 };
116
117 long ipFullShell = -1, ipNextFull = 0;
118 xValenceDonorDensity = 0.;
119 for (long nelem=ipHYDROGEN; nelem < LIMELM; ++nelem)
120 {
121 // H -> 1; He -> 2; Li -> 1 (nelem==2,ipFullShell==1), etc.
122 xValenceDonorDensity +=
123 (nelem-ipFullShell)*dense.gas_phase[nelem];
124 if (nelem == ipNoble[ipNextFull])
125 {
126 ipFullShell = nelem;
127 ++ipNextFull;
128 }
129 }
130
131 //fprintf(ioQQQ,"%g %g\n",dense.eden,xValenceDonorDensity);
132
133 ElectronFraction = (realnum)(MAX2(MIN2(1.0,dense.eden/
134 xValenceDonorDensity),1e-4));
135 }
136 else
137 {
138 /* >>chng 03 apr 29, move evaluation of xNeutralParticleDensity to here
139 * from PresTotCurrent since only used for secondary ionization */
140 /* this is total neutral particle density for collisional ionization
141 * for very high ionization conditions it may be zero */
143 realnum xNeutralParticleDensity = 0.;
144 for( long nelem=0; nelem < LIMELM; nelem++ )
145 {
146 xNeutralParticleDensity += dense.xIonDense[nelem][0];
147 }
148
149 xNeutralParticleDensity += deut.xIonDense[0];
150
151 /* now add all the heavy molecules */
152 for( i=0; i < mole_global.num_calc; i++ )
153 {
154 /* add in if real molecule and not counted above */
155 if(mole.species[i].location == NULL && mole_global.list[i]->charge == 0 && mole_global.list[i]->parentLabel.empty())
156 xNeutralParticleDensity += (realnum)mole.species[i].den;
157 }
158
159 {
160 enum {DEBUG_LOC=false};
161 if( DEBUG_LOC )
162 {
163 fprintf(ioQQQ," xIonDense xNeutralParticleDensity tot\t%.3e",xNeutralParticleDensity);
164 for( long nelem=0; nelem < LIMELM; nelem++ )
165 {
166 fprintf(ioQQQ,"\t%.2e",dense.xIonDense[nelem][0]);
167 }
168 fprintf(ioQQQ,"\n");
169 }
170 }
171 xValenceDonorDensity = (xNeutralParticleDensity+dense.EdenTrue);
172 ElectronFraction = (realnum)(MAX2(MIN2(1.0,dense.EdenTrue/
173 xValenceDonorDensity),1e-4));
174 }
175
176 // xBoundValenceDensity must be consistent with ElectronFraction
177 // for expressions for ionizations below to have bounded behaviour
178 // in the high ElectronFraction limit, i.e. have a factor of (1-EF)
179 realnum xBoundValenceDensity = (1.0f-ElectronFraction)*xValenceDonorDensity;
180
181 if( secondaries.lgSecOFF )
182 {
183 secondaries.HeatEfficPrimary = 1.;
184 secondaries.SecIon2PrimaryErg = 0.;
185 secondaries.SecExcitLya2PrimaryErg = 0.;
186 Cosmic_ray_heat_eff = 0.95;
187 Cosmic_ray_sec2prim = 0.05f;
188 }
189 /* >>chng 03 apr 29, only evaluate one time per zone since drove oscillations
190 * in He+ - He0 ionization front in very high Z models, like hizqso */
191 else
192 {
193
194 /* this is heating efficiency for high-energy photo ejections. It is the ratio
195 * of the final heat added to the energy of the photo electron. Fully
196 * ionized, this is unity, and less than unity for neutral media.
197 * It is a scale factor that multiplies the
198 * high energy heating rate */
199 secondaries.HeatEfficPrimary = 0.9971f*(1.f - pow(1.f - pow(ElectronFraction,(realnum)0.2663f),(realnum)1.3163f));
200
201 /* secondary ionizations per primary Rydberg - the number of secondary
202 * ionizations produced for each Rydberg of high energy photo-electron
203 * energy deposited. It multiplies the high-energy heating rate.
204 * It is zero for a fully ionized gas and goes to 0.3908 for very neutral gas */
205 secondaries.SecIon2PrimaryErg = 0.3908f*pow(1.f - pow(ElectronFraction,(realnum)0.4092f),(realnum)1.7592f);
206 /* by dividing by the energy of one Rydberg, converts heating rate
207 * in ergs into secondary ionization rate */
208 secondaries.SecIon2PrimaryErg /= ((realnum)EN1RYD);
209
210 /* This is cosmic ray secondaries per primary (???),
211 * derived to approximate the curve given in Shull and
212 * van Steenberg for cosmic rays at 35 eV */
213 sec2prim_par_1 = -(1.251f + 195.932f * ElectronFraction);
214 sec2prim_par_2 = 1.f + 46.814f * ElectronFraction - 44.969f *
215 ElectronFraction * ElectronFraction;
216
217 Cosmic_ray_sec2prim = (exp(sec2prim_par_1/SDIV( sec2prim_par_2)));
218
219 /* number of secondary excitations of L-alpha per erg of high
220 * energy light - actually all Lyman lines
221 *
222 * Note--This formula is derived for primary energies greater
223 * than 100 eV and is only an approximation. This will
224 * over predict the secondary ionization of L-alpha. We cannot make
225 * fitting functions for this equation at low energies like we did
226 * for the heating efficiency and for the number of secondaries
227 * because the Shull and van Steenberg paper did not publish a similar
228 * curve for L-alpha */
229 secondaries.SecExcitLya2PrimaryErg = 0.4766f*pow(1.f - pow(ElectronFraction,(realnum)0.2735f),(realnum)1.5221f)/((realnum)EN1RYD);
230
231 if( (trace.lgTrace && trace.lgSecIon) )
232 {
233 fprintf( ioQQQ,
234 " csupra H0 %.2e HeatSum eval sec ion effic, ElectronFraction = %.3e HeatEfficPrimary %.3e SecIon2PrimaryErg %.3e\n",
235 secondaries.csupra[ipHYDROGEN][0],
236 ElectronFraction,
237 secondaries.HeatEfficPrimary ,
238 secondaries.SecIon2PrimaryErg );
239 }
240
241 /* cosmic ray heating, counts as non-ionizing heat since already result of cascade,
242 * this was 35 eV per secondary ionization */
243 /* We want to use the heating efficiency that is applicable to a 35 eV
244 * primary electron, the current efficiency used is for 100eV cosmic ray
245 * Here we use the correct value as given in Wolfire et al. 1995 */
246
247 /* In general the equation for the cosmic ray heating rate is:*
248 * *
249 * *
250 * CR_Heat_Rate = (density)*(cosmic ray ionization rate)* *
251 * (energy of primary electron)* *
252 * (Heating efficiency at that energy) *
253 * *
254 * The product of the last two terms gives the amount of heat *
255 * added by the primary electron, and is dependent upon the *
256 * electron fraction. We are using the same average primary *
257 * electron energy as Wolfire et al. (1995). We do not, *
258 * however, use their formula for the heating efficiency. *
259 * This is because the coefficients (f1, f2, and f3) were *
260 * originally intended for primary electron energies >100eV. *
261 * Instead we derived a heating efficiency based on the curves*
262 * given in Shull and van Steenberg (1985). We interpolated *
263 * for an energy of 35 eV the heating efficiency for electron *
264 * fractions between 1e-4 and 1 *
265 **************************************************************/
266
267 /*printf("Here is ElectronFraction:\t%.3e\n", ElectronFraction);*/
268 Cosmic_ray_heat_eff = - 8.189 - 18.394 * ElectronFraction - 6.608 * ElectronFraction * ElectronFraction * log(ElectronFraction)
269 + 8.322 * exp( ElectronFraction ) + 4.961 * sqrt(ElectronFraction);
270 }
271
272 /*******************************************************************
273 *
274 * get total heating from all species
275 *
276 *******************************************************************/
277
278 /* get total heating */
279 photo_heat_2lev_atoms = 0.;
280 photo_heat_ISO_atoms = 0.;
281 photo_heat_UTA = 0.;
282
283 /* add in effects of high-energy opacity of CO using C and O atomic opacities
284 * k-shell and valence photo of C and O in CO is not explicitly counted elsewhere
285 * this trick roughly accounts for it*/
286 SaveOxygen1 = dense.xIonDense[ipOXYGEN][0];
287 SaveCarbon1 = dense.xIonDense[ipCARBON][0];
288
289 /* atomic C and O will include CO during the heating sum loop */
290 fixit(); /* This breaks the invariant for total mass.
291 *
292 * In any case, is CO the only species which contributes?
293 * -- might expect all other molecular species to do so,
294 * i.e. the item added should perhaps be
295 * dense.xMolecules[nelem]) -- code duplicated in
296 * opacity_addtotal
297 */
298 dense.xIonDense[ipOXYGEN][0] += (realnum) findspecieslocal("CO")->den;
299 dense.xIonDense[ipCARBON][0] += (realnum) findspecieslocal("CO")->den;
300
301 /* this will hold cooling due to metal collisional ionization */
302 CoolHeavy.colmet = 0.;
303 /* metals secondary ionization, Lya excitation */
304 secmet = 0.;
305 smetla = 0.;
306 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
307 {
308 SeconIoniz_iso[ipISO] = 0.;
309 SeconExcit_iso[ipISO] = 0.;
310 }
311
312 /* this loop includes hydrogenic, which is special case due to presence
313 * of substantial excited state populations */
314 for( long nelem=0; nelem<LIMELM; ++nelem)
315 {
316 secmetsave[nelem] = 0.;
317 if( dense.lgElmtOn[nelem] )
318 {
319 /* sum heat over all stages of ionization that exist */
320 /* first do the iso sequences,
321 * h-like and he-like are special because full atom always done,
322 * can be substantial, pops in excited states, with little in ground
323 * (true near LTE), these are done in following loop */
324
325 limit = MIN2( dense.IonHigh[nelem] , nelem+1-NISO );
326
327 /* loop over all elements/ions done with as two-level systems */
328 for( ion=dense.IonLow[nelem]; ion < limit; ion++ )
329 {
330 /* this will be heating for a single stage of ionization */
331 HeatingLo = 0.;
332 HeatingHi = 0.;
333 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
334 {
335 /* heating by various sub-shells */
336 HeatingLo += ionbal.PhotoRate_Shell[nelem][ion][ns][1];
337 HeatingHi += ionbal.PhotoRate_Shell[nelem][ion][ns][2];
338 }
339
340 /* total photoelectric heating, all shells, for this stage
341 * of ionization */
342 thermal.heating[nelem][ion] = dense.xIonDense[nelem][ion]*
343 (HeatingLo + HeatingHi*secondaries.HeatEfficPrimary);
344 /* heating due to UTA ionization */
345 photo_heat_UTA += dense.xIonDense[nelem][ion]*
346 ionbal.UTA_heat_rate[nelem][ion];
347
348 /* add to total heating */
349 photo_heat_2lev_atoms += thermal.heating[nelem][ion];
350 /*if( nzone>290 && thermal.heating[nelem][ion]/thermal.htot>0.01 )
351 fprintf(ioQQQ,"buggg\t%li %li %.3f\n", nelem,ion,thermal.heating[nelem][ion]/thermal.htot);*/
352
353 /* Cooling due to collisional ionization of heavy elements by thermal electrons
354 * CollidRate[nelem][ion][1] is cooling, erg/s/atom, eval in ion_collis */
355 CoolHeavy.colmet += dense.xIonDense[nelem][ion]*ionbal.CollIonRate_Ground[nelem][ion][1];
356
357 /* secondary ionization rate, */
358 secmetsave[nelem] += HeatingHi*secondaries.SecIon2PrimaryErg* dense.xIonDense[nelem][ion];
359
360 /* LA excitation rate, =0 if ionized, s-1 cm-3 */
361 smetla += HeatingHi*secondaries.SecExcitLya2PrimaryErg* dense.xIonDense[nelem][ion];
362 }
363 secmet += secmetsave[nelem];
364
365 /* this branch loop over all ions done with full iso sequence */
366 limit = MAX2( limit, dense.IonLow[nelem] );
367 for( ion=MAX2(0,limit); ion < dense.IonHigh[nelem]; ion++ )
368 {
369 /* this is the iso sequence */
370 ipISO = nelem-ion;
371 /* this will be heating for a single stage of ionization */
372 HeatingLo = 0.;
373 HeatingHi = 0.;
374 /* the outer shell contains the Compton recoil part */
375 for( ns=0; ns < Heavy.nsShells[nelem][ion]; ns++ )
376 {
377 /* heating, erg s-1 atom-1, by low energy, and then high energy, photons */
378 HeatingLo += ionbal.PhotoRate_Shell[nelem][ion][ns][1];
379 HeatingHi += ionbal.PhotoRate_Shell[nelem][ion][ns][2];
380 }
381
382 /* net heating, erg cm-3 s-1 */
383 thermal.heating[nelem][ion] = iso_sp[ipISO][nelem].st[0].Pop()*
384 (HeatingLo + HeatingHi*secondaries.HeatEfficPrimary);
385
386 /* heating due to UTA ionization, erg cm-3 s-1 */
387 photo_heat_UTA += dense.xIonDense[nelem][ion]*
388 ionbal.UTA_heat_rate[nelem][ion];
389
390 /* add to total heating */
391 photo_heat_ISO_atoms += thermal.heating[nelem][ion];
392
393 if( secondaries.SecIon2PrimaryErg > SMALLFLOAT && xBoundValenceDensity > SMALLFLOAT )
394 {
395 /* prevent crash in very high ionization conditions
396 * where xBoundValenceDensity is zero */
397 /* secondary ionization rate, */
398 SeconIoniz_iso[ipISO] += HeatingHi*secondaries.SecIon2PrimaryErg*
399 iso_sp[ipISO][nelem].st[0].Pop()/
400 xBoundValenceDensity;
401
402 /* LA excitation rate, =0 if ionized, units excitations s-1 */
403 SeconExcit_iso[ipISO] += HeatingHi*secondaries.SecExcitLya2PrimaryErg*
404 iso_sp[ipISO][nelem].st[0].Pop()/
405 xBoundValenceDensity;
406
407 ASSERT( SeconIoniz_iso[ipISO]>=0. &&
408 SeconExcit_iso[ipISO]>=0.);
409 }
410 }
411
412 /* make sure stages of ionization with no abundances,
413 * also has no heating */
414 for( ion=0; ion<dense.IonLow[nelem]; ion++ )
415 {
416 ASSERT( thermal.heating[nelem][ion] == 0. );
417 }
418 for( ion=dense.IonHigh[nelem]+1; ion<nelem+1; ion++ )
419 {
420 ASSERT( thermal.heating[nelem][ion] == 0. );
421 }
422 }
423 }
424 if( trace.lgTrace && trace.lgSecIon )
425 {
426 double savemax=0.;
427 long int ipsavemax=-1;
428 for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
429 {
430 if( secmetsave[nelem] > savemax )
431 {
432 savemax = secmetsave[nelem];
433 ipsavemax = nelem;
434 }
435 }
436 fprintf( ioQQQ,
437 " HeatSum secmet largest contributor element %li frac of total %.3e, total %.3e\n",
438 ipsavemax,
439 savemax/SDIV(secmet),
440 secmet);
441 }
442
443 /* now reset the abundances */
444 dense.xIonDense[ipOXYGEN][0] = SaveOxygen1;
445 dense.xIonDense[ipCARBON][0] = SaveCarbon1;
446
447 /* convert secmet to proper final units */
448 /*fprintf(ioQQQ,"bugggg\t%li\t%.3e\t%.3e\t%.3e\n", nzone ,
449 secmet , xBoundValenceDensity , secmet / xBoundValenceDensity );*/
450 /* prevent crash when xBoundValenceDensity is zero */
451 if( secondaries.SecIon2PrimaryErg > SMALLFLOAT && xBoundValenceDensity > SMALLFLOAT )
452 {
453 /* convert from s-1 cm-3 to s-1 - a true rate */
454 secmet /= xBoundValenceDensity;
455 smetla /= xBoundValenceDensity;
456 }
457 else
458 {
459 /* zero */
460 secmet = 0.;
461 smetla = 0.;
462 }
463
464 /* >>chng 01 dec 20, do full some over all secondaries */
465 /* bound Compton recoil heating */
466 /* >>chng 02 mar 28, save heating in this var rather than heating[0][18]
467 * since now saved in photo heat
468 * this is only used for a printout and in lines, not as heat since already counted*/
469 ionbal.CompRecoilHeatLocal = 0.;
470 for( long nelem=0; nelem<LIMELM; ++nelem )
471 {
472 for( ion=0; ion<nelem+1; ++ion )
473 {
474 ionbal.CompRecoilHeatLocal +=
475 ionbal.CompRecoilHeatRate[nelem][ion]*secondaries.HeatEfficPrimary*dense.xIonDense[nelem][ion];
476 }
477 }
478 /* >>chng 05 nov 26, include this term - bound ionization of H2
479 * assume total cs is that of two separated H */
480 ionbal.CompRecoilHeatLocal +=
481 2.*ionbal.CompRecoilHeatRate[ipHYDROGEN][0]*secondaries.HeatEfficPrimary*hmi.H2_total;
482
483 /* find heating due to charge transfer
484 * >>chng 05 oct 29, move from here to conv base so that can be treated as cooling if
485 * negative heating */
486 thermal.heating[0][24] = thermal.char_tran_heat;
487
488 /* heating due to pair production */
489 thermal.heating[0][21] =
490 ionbal.PairProducPhotoRate[2]*secondaries.HeatEfficPrimary*(dense.gas_phase[ipHYDROGEN] + 4.F*dense.gas_phase[ipHELIUM]);
491 /* last term above is number of nuclei in helium */
492
493 /* this is heating due to fast neutrons, assumed to secondary ionize */
494 thermal.heating[0][20] =
495 ionbal.xNeutronHeatRate*secondaries.HeatEfficPrimary*dense.gas_phase[ipHYDROGEN];
496
497 /* turbulent heating, assumed to be a non-ionizing heat agent, added here */
498 thermal.heating[0][20] += ionbal.ExtraHeatRate;
499
500 /* UTA heating */
501 thermal.heating[0][7] = photo_heat_UTA;
502 /*fprintf(ioQQQ,"DEBUG UTA heat %.3e\n", photo_heat_UTA/SDIV(thermal.htot ));*/
503
504 // diatomic photoionization heating
505 thermal.heating[0][18] = 0.;
506 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
507 {
508 /*>>KEYWORD H2 photoionization heating */
509 thermal.heating[0][18] += (*diatom)->Abund() *
510 ( (*diatom)->photo_heat_soft + (*diatom)->photo_heat_hard*secondaries.HeatEfficPrimary);
511 }
513
514 /* >>chng 05 nov 27, approximate heating due to H2+, H3+ photoionization
515 * assuming H0 rates */
516 /*>>KEYWORD H2+ photoionization heating; H3+ photoionization heating */
517 thermal.heating[0][26] = (findspecieslocal("H2+")->den+findspecieslocal("H3+")->den) *
518 (ionbal.PhotoRate_Shell[ipHYDROGEN][0][0][1] +
519 ionbal.PhotoRate_Shell[ipHYDROGEN][0][0][2]*secondaries.HeatEfficPrimary);
520
521 /* add on cosmic rays - most important in neutral or molecular gas, use
522 * difference between total H and H+ density as substitute for sum of
523 * H in atoms and all molecules, but only in limit where difference is
524 * significant */
525# if 0
526 double hneut;
527 if( (dense.gas_phase[ipHYDROGEN] - dense.xIonDense[ipHYDROGEN][1])/
528 dense.gas_phase[ipHYDROGEN]<0.001 )
529 {
530 /* limit where most H is ionized - simply use atomic and H2 */
531 hneut = dense.xIonDense[ipHYDROGEN][0] + 2.*(findspecieslocal("H2")->den+findspecieslocal("H2*")->den);
532 }
533 else
534 {
535 /* limit where neutral is significant, use different between total and H+ */
536 hneut = dense.gas_phase[ipHYDROGEN] - dense.xIonDense[ipHYDROGEN][1];
537 }
538# endif
539
540 /* cosmic ray heating */
541 thermal.heating[1][6] =
542 ionbal.CosRayHeatNeutralParticles*
543 xBoundValenceDensity * Cosmic_ray_heat_eff;
544 /* cosmic ray heating of thermal electrons */
545 /* >>refer CR physics Ferland, G.J. & Mushotzky, R.F., 1984, ApJ, 286, 42 */
546 thermal.heating[1][6] += ionbal.CosRayHeatThermalElectrons * dense.eden;
547
548 /* now sum up all heating agents not included in sum over normal bound levels above */
549 OtherHeat = 0.;
550 for( long nelem=0; nelem<LIMELM; ++nelem)
551 {
552 /* we now have ionization solution for this element,sum heat over
553 * all stages of ionization that exist */
554 /* >>>chng 99 may 08, following loop had started at nelem+3, and so missed [1][0],
555 * which is excited states of hydrogenic species. increase this limit */
556 for( i=nelem+1; i < LIMELM; i++ )
557 {
558 OtherHeat += thermal.heating[nelem][i];
559 }
560 }
561
562 thermal.htot = OtherHeat + photo_heat_2lev_atoms + photo_heat_ISO_atoms;
563
564 /* following checks whether total heating is strange, if we are not in search phase */
565 if( called.lgTalk && !conv.lgSearch )
566 {
567 /* print this warning if not constant temperature and neg heat */
568 if( thermal.htot < 0. && !thermal.lgTemperatureConstant)
569 {
570 fprintf( ioQQQ,
571 " Total heating is <0; is this model collisionally ionized? zone is %li\n",
572 nzone );
573 }
574 else if( thermal.htot == 0. )
575 {
576 fprintf( ioQQQ,
577 " Total heating is 0; is the density small? zone is %li\n",
578 nzone);
579 }
580 }
581
582 /* add on line heating to this array, heatl was evaluated in sumcool
583 * not zero, because of roundoff error */
584 if( thermal.heatl/thermal.ctot < -1e-15 )
585 {
586 fprintf( ioQQQ, " HeatSum gets negative line heating,%10.2e%10.2e this is insane.\n",
587 thermal.heatl, thermal.ctot );
588
589 fprintf( ioQQQ, " this is zone%4ld\n", nzone );
590 ShowMe();
592 }
593
594 fixit(); // much or all of this secondary stuff (next ~120 lines) should be updated much earlier
595 // solvers in conv_base use outdated details
596
597 /*******************************************************************
598 *
599 * secondary ionization and excitation rates
600 *
601 *******************************************************************/
602
603 /* the terms cosmic_ray_ionization_rate & SecExcitLyaRate contain energies added in highen.
604 * The only significant one is usually bound Compton heating except when
605 * cosmic rays are present */
606
607 if( secondaries.SecIon2PrimaryErg > SMALLFLOAT && xBoundValenceDensity > SMALLFLOAT )
608 {
609 /* add on ionization rate s-1 cm-3 due to pair production */
610 /* next two terms account for secondary ionization and excitation due to pair productions.
611 * The code integrates over the radiation field where \gamma <-> e- e+ is possible -
612 * that is ionbal.PairProducPhotoRate. For all SEDs of realistic sources (AGN, or x-ray binaries)
613 * the radiation field over this range is small compared with the FUV/x-ray. Pair production will,
614 * as a result, only be important energetically in well-shielded regions where much of the SED
615 * is absorbed. In these regions gas would be neutral or molecular and the pair would produce
616 * the same effects as a cosmic ray - heating, excitation, and ionization. That is why the code is
617 * here. Given the shape of the SED these effects are the largest pair production should produce -
618 * even then it is negligible.
619 */
620
621 /* the photon rates in ionbal.PairProducPhotoRate are the integral over the range that the
622 * do pair production using cross sections from Figure 41 of Experimental Nuclear Physics,
623 * Vol1, E.Segre, ed
624 * factor of four in DensityRatio accounts for additional nucleons in alpha particle */
625 realnum DensityRatio = (dense.gas_phase[ipHYDROGEN] +
626 4.F*dense.gas_phase[ipHELIUM])/xBoundValenceDensity;
627
628 pair_production_ionization_rate =
629 ionbal.PairProducPhotoRate[2]*secondaries.SecIon2PrimaryErg*DensityRatio;
630
631 /* total secondary excitation rate cm-3 s-1 for Lya due to pair production*/
632 SecExcitLyaRate =
633 ionbal.PairProducPhotoRate[2]*secondaries.SecExcitLya2PrimaryErg*1.3333*DensityRatio;
634
635 /* ionization rate of fast neutrons */
636 fast_neutron_ionization_rate =
637 ionbal.xNeutronHeatRate*secondaries.SecIon2PrimaryErg*DensityRatio;
638
639 /* Secondary excitation rate due to neutrons */
640 SecExcitLyaRate +=
641 ionbal.xNeutronHeatRate*secondaries.SecExcitLya2PrimaryErg*1.3333*DensityRatio;
642
643 /* cosmic ray Lya excitation rate */
644 SecExcitLyaRate +=
645 /* multiply by atomic H and He densities */
646 ionbal.CosRayHeatNeutralParticles*secondaries.SecExcitLya2PrimaryErg*1.3333*DensityRatio;
647 }
648 else
649 {
650 /* zero */
651 pair_production_ionization_rate = 0.;
652 SecExcitLyaRate = 0.;
653 fast_neutron_ionization_rate = 0.;
654 }
655
656 /* cosmic ray ionization */
657 /* >>>chng 99 apr 29, term in PhotoRate was not present */
658 cosmic_ray_ionization_rate =
659 /* this term is H^0 cosmic ray ionization, set in highen, did not multiply by
660 * collider density in highen, so do not divide by it here */
661 /* >>chng 99 jun 29, added SecIon2PrimaryErg to cosmic ray rate, also multiply
662 * by number of secondaries per primary*/
663 ionbal.CosRayIonRate*Cosmic_ray_sec2prim +
664 /* this is the cosmic ray heating rate */
665 ionbal.CosRayHeatNeutralParticles*secondaries.SecIon2PrimaryErg;
666
667 /* find total suprathermal collisional ionization rate
668 * CSUPRA is H0 col ionize rate from non-thermal electrons (inverse sec)
669 * x12tot is LA excitation rate, units s-1
670 * SECCMP evaluated in HIGHEN, =ioniz rate for cosmic rays, sec bound */
671
672 /* option to force secondary ionization rate, normally false */
673 if( secondaries.lgCSetOn )
674 {
675 for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
676 {
677 for( ion=0; ion<nelem+1; ++ion )
678 {
679 secondaries.csupra[nelem][ion] = secondaries.SetCsupra*secondaries.csupra_effic[nelem][ion];
680 }
681 }
682 secondaries.x12tot = secondaries.SetCsupra;
683 }
684 else
685 {
686 double csupra = cosmic_ray_ionization_rate +
687 pair_production_ionization_rate +
688 fast_neutron_ionization_rate +
689 SeconIoniz_iso[ipH_LIKE] +
690 SeconIoniz_iso[ipHE_LIKE] +
691 secmet;
692
694 /* now fill in ionization rates for all elements and ion stages */
695 for( long nelem=ipHYDROGEN; nelem<LIMELM; ++nelem )
696 {
697 for( ion=0; ion<nelem+1; ++ion )
698 {
699 secondaries.csupra[nelem][ion] = (realnum)csupra*secondaries.csupra_effic[nelem][ion];
700 }
701 }
702
703 fixit(); // why does x12tot include non Ly-a stuff here, and get used to scale other lines elsewhere?
704 /* secondary suprathermal excitation of Ly-alpha, excitations s-1 */
705 secondaries.x12tot = SecExcitLyaRate +
706 SeconExcit_iso[ipH_LIKE] +
707 SeconExcit_iso[ipHE_LIKE] +
708 smetla;
709
710 if (0)
711 fprintf(ioQQQ,"x12 %ld %15.8g %15.8g %15.8g\n",nzone,
712 secondaries.x12tot,
713 secondaries.SecExcitLya2PrimaryErg,ElectronFraction);
714 }
715
716 if( trace.lgTrace && secondaries.csupra[ipHYDROGEN][0] > 0. )
717 {
718 fprintf( ioQQQ,
719 " HeatSum return CSUPRA %9.2e SECCMP %6.3f SecHI %6.3f SECHE %6.3f SECMET %6.3f efrac= %9.2e \n",
720 secondaries.csupra[ipHYDROGEN][0],
721 (cosmic_ray_ionization_rate+pair_production_ionization_rate+fast_neutron_ionization_rate)/secondaries.csupra[ipHYDROGEN][0],
722 SeconIoniz_iso[ipH_LIKE] / secondaries.csupra[ipHYDROGEN][0],
723 SeconIoniz_iso[ipHE_LIKE]/secondaries.csupra[ipHYDROGEN][0],
724 secmet/secondaries.csupra[ipHYDROGEN][0] ,
725 ElectronFraction );
726 }
727
728 /*******************************************************************
729 *
730 * get derivative of total heating
731 *
732 *******************************************************************/
733
734 /* now get derivative of heating due to photoionization,
735 * >>chng 96 jan 14
736 *>>>>>NB heating decreasing with increasing temp is negative dH/dT */
737 thermal.dHeatdT = -0.7*(photo_heat_2lev_atoms+photo_heat_ISO_atoms+
738 photo_heat_UTA)/phycon.te;
739 /* >>chng 04 feb 28, add this correction factor - when gas totally neutral heating
740 * does not depend on temperature - when ionized depends on recom coefficient - this
741 * smoothly joins two limits */
742 thermal.dHeatdT *= dense.xIonDense[ipHYDROGEN][1]/dense.gas_phase[ipHYDROGEN];
743 if( PRT_DERIV )
744 fprintf(ioQQQ,"DEBUG dhdt 0 %.3e %.3e %.3e \n",
745 thermal.dHeatdT,
746 photo_heat_2lev_atoms,
747 photo_heat_ISO_atoms);
748
749 /* oldfac was factor used in old implementation
750 * any term using it should be rethought */
751 oldfac = -0.7/phycon.te;
752
753 /* carbon monoxide */
754 thermal.dHeatdT += thermal.heating[0][9]*oldfac;
755
756 /* Ly alpha collisional heating */
757 thermal.dHeatdT += thermal.heating[0][10]*oldfac;
758
759 /* line heating */
760 thermal.dHeatdT += thermal.heating[0][22]*oldfac;
761 if( PRT_DERIV )
762 fprintf(ioQQQ,"DEBUG dhdt 1 %.3e \n", thermal.dHeatdT);
763
764 /* free free heating, propto T^-0.5
765 * >>chng 96 oct 30, from heating(1,12) to CoolHeavy.brems_heat_total,
766 * do cooling separately assume CoolHeavy.brems_heat_total goes as T^-3/2
767 * dHTotDT = dHTotDT + heating(1,12) * (-0.5/te) */
768 thermal.dHeatdT += CoolHeavy.brems_heat_total*(-1.5/phycon.te);
769
770 /* >>chng 04 aug 07, use better estimate for heating derivative; needed in PDRs, PvH */
771 /* this includes PE, thermionic, and collisional heating of the gas by the grains */
772 thermal.dHeatdT += gv.dHeatdT;
773
774 /* helium triplets heating */
775 thermal.dHeatdT += thermal.heating[1][2]*oldfac;
776 if( PRT_DERIV )
777 fprintf(ioQQQ,"DEBUG dhdt 2 %.3e \n", thermal.dHeatdT);
778
779 /* hydrogen molecule collisional deexcitation heating */
780 /* >>chng 04 jan 25, add max to prevent cooling from entering here */
781 /*thermal.dHeatdT += MAX2(0.,thermal.heating[0][8])*oldfac;*/
782 if( hmi.HeatH2Dexc_used > 0. )
783 thermal.dHeatdT += hmi.deriv_HeatH2Dexc_used;
784
785 /* >>chng 96 nov 15, had been oldfac below, wrong sign
786 * H- H minus heating, goes up with temp since rad assoc does */
787 thermal.dHeatdT += thermal.heating[0][15]*0.7/phycon.te;
788
789 /* H2+ heating */
790 thermal.dHeatdT += thermal.heating[0][16]*oldfac;
791
792 /* Balmer continuum and all other excited states
793 * - T goes up so does pop and heating
794 * but this all screwed up by change in eden */
795 thermal.dHeatdT += thermal.heating[0][1]*oldfac;
796 if( PRT_DERIV )
797 fprintf(ioQQQ,"DEBUG dhdt 3 %.3e \n", thermal.dHeatdT);
798
799 /* all three body heating, opposite of coll ion cooling */
800 thermal.dHeatdT += thermal.heating[0][3]*oldfac;
801
802 /* bound electron recoil heating
803 thermal.dHeatdT += ionbal.CompRecoilHeatLocal*oldfac; */
804
805 /* Compton heating */
806 thermal.dHeatdT += thermal.heating[0][19]*oldfac;
807
808 /* extra heating source, usually zero */
809 thermal.dHeatdT += thermal.heating[0][20]*oldfac;
810
811 /* pair production */
812 thermal.dHeatdT += thermal.heating[0][21]*oldfac;
813 if( PRT_DERIV )
814 fprintf(ioQQQ,"DEBUG dhdt 4 %.3e \n", thermal.dHeatdT);
815
816 /* derivative of heating due to collisions of H lines, heating(1,24) */
817 for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
818 {
819 for( long nelem=ipISO; nelem<LIMELM; ++nelem)
820 {
821 thermal.dHeatdT += iso_sp[ipISO][nelem].dLTot;
822 }
823 }
824
825 /* heating due to large FeII atom, zero unless atom FeII is entered,
826 * FeII.Fe2_large_heat was entered into thermal.heating[25][27] */
827 if( FeII.Fe2_large_heat > 0. )
828 {
829 thermal.dHeatdT += FeII.ddT_Fe2_large_cool;
830 }
831 if( PRT_DERIV )
832 fprintf(ioQQQ,"DEBUG dhdt 5 %.3e \n", thermal.dHeatdT);
833
834 /* possibility of getting empirical heating derivative
835 * normally false, set true with 'set numerical derivatives' command */
836 if( NumDeriv.lgNumDeriv )
837 {
838 if( ((nzone > 2 && nzone == nzSave) &&
839 ! fp_equal( oldtemp, phycon.te )) && nhit > 4 )
840 {
841 /* hnit is number of tries on this zone - use to stop numerical problems
842 * do not evaluate numerical derivative until well into solution */
843 deriv = (oldheat - thermal.htot)/(oldtemp - phycon.te);
844 thermal.dHeatdT = deriv;
845 }
846 else
847 {
848 deriv = thermal.dHeatdT;
849 }
850
851 /* this happens when new zone starts */
852 if( nzone != nzSave )
853 {
854 nhit = 0;
855 }
856 nzSave = nzone;
857 nhit += 1;
858 oldheat = thermal.htot;
859 oldtemp = phycon.te;
860 }
861
862 if( trace.lgTrace && trace.lgHeatBug )
863 {
864 ipnt = 0;
865 /* this loops through the 2D array that contains all agents counted in htot */
866 for( i=0; i < LIMELM; i++ )
867 {
868 for( j=0; j < LIMELM; j++ )
869 {
870 if( thermal.heating[i][j]/thermal.htot > FAINT_HEAT )
871 {
872 ipsave[ipnt] = i;
873 jpsave[ipnt] = j;
874 save[ipnt] = thermal.heating[i][j]/thermal.htot;
875 /* increment ipnt, but do not let it get too big */
876 ipnt = MIN2((long)NDIM-1,ipnt+1);
877 }
878 }
879 }
880
881 /* now check for possible line heating not counted in 1,23
882 * there should not be any significant heating source here
883 * since they would not be counted in derivative correctly */
884 for( i=0; i < thermal.ncltot; i++ )
885 {
886 if( thermal.heatnt[i]/thermal.htot > FAINT_HEAT )
887 {
888 ipsave[ipnt] = -1;
889 jpsave[ipnt] = i;
890 save[ipnt] = thermal.heatnt[i]/thermal.htot;
891 ipnt = MIN2((long)NDIM-1,ipnt+1);
892 }
893 }
894
895 fprintf( ioQQQ,
896 " HeatSum HTOT %.3e Te:%.3e dH/dT%.2e other %.2e photo 2lev %.2e photo iso %.2e\n",
897 thermal.htot,
898 phycon.te,
899 thermal.dHeatdT ,
900 /* total heating is sum of following three terms
901 * OtherHeat is a sum over all other heating agents */
902 OtherHeat ,
903 photo_heat_2lev_atoms,
904 photo_heat_ISO_atoms);
905
906 fprintf( ioQQQ, " " );
907 for( i=0; i < ipnt; i++ )
908 {
909 /*lint -e644 these three are initialized above */
910 fprintf( ioQQQ, " [%ld][%ld]%6.3f",
911 ipsave[i],
912 jpsave[i],
913 save[i] );
914 /*lint +e644 these three are initialized above */
915 /* throw a cr every n numbers */
916 if( !(i%8) && i>0 && i!=(ipnt-1) )
917 {
918 fprintf( ioQQQ, "\n " );
919 }
920 }
921 fprintf( ioQQQ, "\n" );
922 }
923 return;
924}
925
926/* =============================================================================*/
927/* HeatZero is called by ConvBase */
928void HeatZero( void )
929{
930 long int i , j;
931
932 DEBUG_ENTRY( "HeatZero()" );
933
934 for( i=0; i < LIMELM; i++ )
935 {
936 for( j=0; j < LIMELM; j++ )
937 {
938 thermal.heating[i][j] = 0.;
939 }
940 }
941 return;
942}
t_FeII FeII
Definition atomfeii.cpp:5
t_called called
Definition called.cpp:5
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
const int ipCARBON
Definition cddefines.h:310
#define MIN2
Definition cddefines.h:761
const int ipOXYGEN
Definition cddefines.h:312
const int LIMELM
Definition cddefines.h:258
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int NISO
Definition cddefines.h:261
const int ipHELIUM
Definition cddefines.h:306
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
sys_float SDIV(sys_float x)
Definition cddefines.h:952
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
const int ipNEON
Definition cddefines.h:314
const int ipKRYPTON
Definition cddefines.h:335
const int ipARGON
Definition cddefines.h:322
void ShowMe(void)
Definition service.cpp:181
void fixit(void)
Definition service.cpp:991
double den
Definition mole.h:358
t_conv conv
Definition conv.cpp:5
static const bool PRT_DERIV
Definition cool_eval.cpp:43
t_CoolHeavy CoolHeavy
Definition coolheavy.cpp:5
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
t_deuterium deut
Definition deuterium.cpp:8
GrainVar gv
Definition grainvar.cpp:5
vector< diatomics * > diatoms
Definition h2.cpp:8
vector< diatomics * >::iterator diatom_iter
Definition h2.h:13
void HeatZero(void)
Definition heat_sum.cpp:928
void HeatSum(void)
Definition heat_sum.cpp:33
static const double FAINT_HEAT
Definition heat_sum.cpp:29
t_Heavy Heavy
Definition heavy.cpp:5
t_hmi hmi
Definition hmi.cpp:5
t_ionbal ionbal
Definition ionbal.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipHE_LIKE
Definition iso.h:63
const int ipH_LIKE
Definition iso.h:62
t_mole_global mole_global
Definition mole.cpp:6
t_mole_local mole
Definition mole.cpp:7
molezone * findspecieslocal(const char buf[])
t_NumDeriv NumDeriv
Definition numderiv.cpp:5
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double EN1RYD
Definition physconst.h:179
t_save save
Definition save.cpp:5
t_secondaries secondaries
t_thermal thermal
Definition thermal.cpp:5
t_trace trace
Definition trace.cpp:5