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