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

Generated for cloudy by doxygen 1.8.1.1