cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
iso_cool.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 /*iso_cool compute net cooling due to species in iso-sequences */
4 #include "cddefines.h"
5 #include "physconst.h"
6 #include "taulines.h"
7 #include "hydrogenic.h"
8 #include "elementnames.h"
9 #include "phycon.h"
10 #include "dense.h"
11 #include "thermal.h"
12 #include "cooling.h"
13 #include "iso.h"
14 
15 /* HP cc cannot compile this routine with any optimization */
16 #if defined(__HP_aCC)
17 # pragma OPT_LEVEL 1
18 #endif
19 
20 // set to true to enable debug print of contributors to collisional ionization cooling
21 const bool lgPrintIonizCooling = false;
22 
23 void iso_cool(
24  /* iso sequence, 0 for hydrogenic */
25  long int ipISO ,
26  /* nelem is charge -1, so 0 for H itself */
27  long int nelem)
28 {
29  long int ipHi,
30  ipbig,
31  ipLo,
32  n;
33  double RecCoolExtra,
34  biggest = 0.,
35  dCdT_all,
36  edenHCorr_IonAbund,
37  edenIonAbund,
38  CollisIonizatCoolingTotal,
39  dCollisIonizatCoolingTotalDT,
40  HeatExcited,
41  heat_max,
42  CollisIonizatCooling,
43  CollisIonizatCoolingDT,
44  hlone,
45  thin,
46  ThinCoolingCaseB,
47  ThinCoolingSum,
48  collider;
49 
50  valarray<double> CollisIonizatCoolingArr,
51  CollisIonizatCoolingDTArr,
52  SavePhotoHeat,
53  SaveInducCool,
54  SaveRadRecCool;
55 
56  long int nlo_heat_max , nhi_heat_max;
57 
58  /* place to put string labels for iso lines */
59  char chLabel[NCOLNT_LAB_LEN+1];
60 
61  DEBUG_ENTRY( "iso_cool()" );
62 
63  /* validate the incoming data */
64  ASSERT( nelem >= ipISO );
65  ASSERT( ipISO <NISO );
66  ASSERT( nelem < LIMELM );
67  /* local number of levels may be less than malloced number if continuum
68  * has been lowered due to high density */
69  ASSERT( iso.numLevels_local[ipISO][nelem] <= iso.numLevels_max[ipISO][nelem] );
70 
71  /* for zero abundance or if element has been turned off we need to
72  * set some produced variables to zero */
73  if( dense.xIonDense[nelem][nelem+1-ipISO] <= 0. ||
74  /* >>chng 04 dec 07, add test for recombined species having zero abundance,
75  * this occurs when set ion is in place so very artificial */
76  dense.xIonDense[nelem][nelem-ipISO]<=0. ||
77  !dense.lgElmtOn[nelem] )
78  {
79  /* all global variables must be zeroed here or set below */
80  iso.coll_ion[ipISO][nelem] = 0.;
81  iso.cLya_cool[ipISO][nelem] = 0.;
82  iso.cLyrest_cool[ipISO][nelem] = 0.;
83  iso.cBal_cool[ipISO][nelem] = 0.;
84  iso.cRest_cool[ipISO][nelem] = 0.;
85  iso.xLineTotCool[ipISO][nelem] = 0.;
86  iso.RadRecCool[ipISO][nelem] = 0.;
87  iso.FreeBnd_net_Cool_Rate[ipISO][nelem] = 0.;
88  iso.dLTot[ipISO][nelem] = 0.;
89  iso.RecomInducCool_Rate[ipISO][nelem] = 0.;
90  return;
91  }
92  /* >>chng 05 aug 17, must use real electron density for collisional ionization of H
93  * since in Leiden v4 pdr with its artificial high temperature coll ion can be important
94  * H on H is homonuclear and scaling laws for other elements does not apply
95  * next two were multiplied by dense.EdenHCorr and changed to dense.eden for H,
96  * EdenHCorr for rest */
97  if( nelem==ipHYDROGEN )
98  {
99  /* special version for H0 onto H0 */
100  collider = dense.EdenHontoHCorr;
101  }
102  else
103  {
104  collider = dense.EdenHCorr;
105  }
106 
107  /* create some space, these go to numLevels_local instead of _max
108  * since continuum may have been lowered by density */
109  if( lgPrintIonizCooling )
110  {
111  CollisIonizatCoolingArr.resize( iso.numLevels_local[ipISO][nelem] );
112  CollisIonizatCoolingDTArr.resize( iso.numLevels_local[ipISO][nelem] );
113  }
114  SavePhotoHeat.resize( iso.numLevels_local[ipISO][nelem] );
115  SaveInducCool.resize( iso.numLevels_local[ipISO][nelem] );
116  SaveRadRecCool.resize( iso.numLevels_local[ipISO][nelem] );
117 
118  /***********************************************************************
119  * *
120  * collisional ionization cooling, less three-body recombination heat *
121  * *
122  ***********************************************************************/
123 
124  /* will be net collisional ionization cooling, units erg/cm^3/s */
125  CollisIonizatCoolingTotal = 0.;
126  dCollisIonizatCoolingTotalDT = 0.;
127 
128  /* collisional ionization cooling minus three body heating
129  * depending on how topoff is done, highest level can have large population
130  * and its coupling to continuum can be large, at various times code
131  * had to ignore effects of very highest level, but starting in mid
132  * 20006 all levels have been included
133  * 2008 apr 18, do not include highest - when only 2 collapsed levels
134  * are used several density 13 BLR sims have serious convergence
135  * problems */
136  for( n=0; n < iso.numLevels_local[ipISO][nelem]-1; ++n )
137  {
138  /* total collisional ionization cooling less three body heating */
139  CollisIonizatCooling =
140  iso.xIsoLevNIonRyd[ipISO][nelem][n]*iso.ColIoniz[ipISO][nelem][n]*collider*
141  (StatesElem[ipISO][nelem][n].Pop -iso.PopLTE[ipISO][nelem][n]*dense.eden);
142  CollisIonizatCoolingTotal += CollisIonizatCooling;
143 
144  /* the derivative of the cooling
145  * need extra factor of temp of 1 Ryd since div by square of temp in Ryd */
146  CollisIonizatCoolingDT = CollisIonizatCooling*
147  (iso.xIsoLevNIonRyd[ipISO][nelem][n]/TE1RYD/POW2(phycon.te_ryd));
148 
149  dCollisIonizatCoolingTotalDT += CollisIonizatCoolingDT;
150  // save values for debug printout
151  if( lgPrintIonizCooling )
152  {
153  CollisIonizatCoolingArr[n] = CollisIonizatCooling;
154  CollisIonizatCoolingDTArr[n] = CollisIonizatCoolingDT;
155  }
156  }
157 
158  /* save net collisional ionization cooling less H-3 body heating
159  * for inclusion in printout
160  * convert to physical units, need to convert Ryd to ergs,
161  * and bring to density per vol not per ion */
162  iso.coll_ion[ipISO][nelem] = CollisIonizatCoolingTotal * EN1RYD*
163  dense.xIonDense[nelem][nelem+1-ipISO];
164 
165  /* derivative wrt temp
166  * convert to physical units, need to convert Ryd to ergs,
167  * and bring to density per vol not per ion */
168  double dcool = dCollisIonizatCoolingTotalDT * EN1RYD *
169  dense.xIonDense[nelem][nelem+1-ipISO];
170  fixit();
173  dcool = iso.coll_ion[ipISO][nelem] / phycon.te;
174 
175  /* create a meaningful label */
176  sprintf(chLabel , "IScion%2s%2s" , elementnames.chElementSym[ipISO] ,
177  elementnames.chElementSym[nelem] );
178  /* dump the coolant onto the cooling stack */
179  CoolAdd(chLabel,(realnum)nelem,MAX2(0.,iso.coll_ion[ipISO][nelem]));
180 
181  thermal.dCooldT += dcool;
182 
183  /* heating[0][3] is all three body heating, opposite of collisional
184  * ionization cooling,
185  * would be unusual for this to be non-zero since would require excited
186  * state departure coefficients to be greater than unity */
187  thermal.heating[0][3] += MAX2(0.,-iso.coll_ion[ipISO][nelem]);
188 
189  /* debug block printing individual contributors to collisional ionization cooling */
190  if( lgPrintIonizCooling && nelem==1 && ipISO==1 )
191  {
192  fprintf(ioQQQ,"DEBUG coll ioniz cool contributors:");
193  for( n=0; n < iso.numLevels_local[ipISO][nelem]; n++ )
194  {
195  if( CollisIonizatCoolingArr[n] / SDIV( CollisIonizatCoolingTotal ) > 0.1 )
196  fprintf(ioQQQ," %2li %.1e",
197  n,
198  CollisIonizatCoolingArr[n]/ SDIV( CollisIonizatCoolingTotal ) );
199  }
200  fprintf(ioQQQ,"\n");
201  fprintf(ioQQQ,"DEBUG coll ioniz derivcontributors:");
202  for( n=0; n < iso.numLevels_local[ipISO][nelem]; n++ )
203  {
204  if( CollisIonizatCoolingDTArr[n] / SDIV( dCollisIonizatCoolingTotalDT ) > 0.1 )
205  fprintf(ioQQQ," %2li %.1e",
206  n,
207  CollisIonizatCoolingDTArr[n]/ SDIV( dCollisIonizatCoolingTotalDT ) );
208  }
209  fprintf(ioQQQ,"\n");
210  }
211 
212  /***********************************************************************
213  * *
214  * hydrogen recombination free-bound free bound cooling *
215  * *
216  ***********************************************************************/
217 
218  /* this is the product of the ion abundance times the electron density,
219  * will multiply level pops which are stored relative to ion
220  * EdenHCorr is used in level pops, so should be used here too */
221  edenIonAbund = dense.eden*dense.xIonDense[nelem][nelem+1-ipISO];
222 
223  /* this is the product of the ion abundance times the electron density with
224  * a small correction for the presence of neutrals, which should be used in
225  * reactions that involve collisions between states, but NOT radiative recombination
226  * but should be used in collisional ionization / three body recombination */
227  edenHCorr_IonAbund = collider*dense.xIonDense[nelem][nelem+1-ipISO];
228 
229  /* now do case b sum to compare with exact value below */
230  iso.RadRecCool[ipISO][nelem] = 0.;
231  ThinCoolingSum = 0.;
232 
233  if( ipISO == ipH_LIKE )
234  {
235  /* do ground with special approximate fits to Ferland et al. '92 */
236  thin = HydroRecCool(
237  /* n is the prin quantum number on the physical scale */
238  1 ,
239  /* nelem is the charge on the C scale, 0 is hydrogen */
240  nelem);
241  }
242  else
243  {
244  /* this is the cooling before correction for optical depths */
245  thin = iso.RadRecomb[ipISO][nelem][0][ipRecRad]*
246  /* arg is the scaled temperature, T * n^2 / Z^2,
247  * n is principal quantum number, Z is charge, 1 for H */
248  HCoolRatio(
249  phycon.te * POW2( (double)StatesElem[ipISO][nelem][0].n / (double)(nelem+1-ipISO) ))*
250  /* convert results to energy per unit vol */
251  phycon.te * BOLTZMANN;
252  }
253  /* the cooling, corrected for optical depth */
254  SaveRadRecCool[0] = iso.RadRecomb[ipISO][nelem][0][ipRecNetEsc] * thin;
255  /* this is now total free-bound cooling */
256  iso.RadRecCool[ipISO][nelem] += SaveRadRecCool[0] * edenIonAbund;
257 
258  /* radiative recombination cooling for all excited states */
259  for( n=1; n < iso.numLevels_local[ipISO][nelem]; n++ )
260  {
261  /* this is the cooling before correction for optical depths */
262  thin = iso.RadRecomb[ipISO][nelem][n][ipRecRad]*
263  /* arg is the scaled temperature, T * n^2 / Z^2,
264  * n is principal quantum number, Z is charge, 1 for H */
265  HCoolRatio(
266  phycon.te * POW2( (double)StatesElem[ipISO][nelem][n].n / (double)(nelem+1-ipISO) ))*
267  /* convert results to energy per unit vol */
268  phycon.te * BOLTZMANN;
269 
270  /* the cooling, corrected for optical depth */
271  SaveRadRecCool[n] = iso.RadRecomb[ipISO][nelem][n][ipRecNetEsc] * thin * edenIonAbund;
272  /* this is now total free-bound cooling */
273  iso.RadRecCool[ipISO][nelem] += SaveRadRecCool[n];
274 
275  /* keep track of case b sum for topoff below */
276  ThinCoolingSum += thin;
277  }
278  {
279  /* debug block for state specific recombination cooling */
280  enum {DEBUG_LOC=false};
281  if( DEBUG_LOC )
282  {
283  if( nelem==ipISO )
284  {
285  /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */
286  for( n=0; n < (iso.numLevels_local[ipISO][nelem] - 1); n++ )
287  {
288  fprintf(ioQQQ,"\t%.2f",SaveRadRecCool[n]/ThinCoolingSum);
289  }
290  fprintf(ioQQQ,"\n");
291  }
292  }
293  }
294 
295  /* Case b sum of optically thin radiative recombination cooling.
296  * add any remainder to the sum from above - high precision is needed
297  * to get STE result to converge close to equilibrium - only done for
298  * H-like ions where exact result is known */
299  if( ipISO == ipH_LIKE )
300  {
301  /* these expressions are only valid for hydrogenic sequence */
302  if( nelem == 0 )
303  {
304  /*expansion for hydrogen itself */
305  ThinCoolingCaseB = (-25.859117 +
306  0.16229407*phycon.telogn[0] +
307  0.34912863*phycon.telogn[1] -
308  0.10615964*phycon.telogn[2])/(1. +
309  0.050866793*phycon.telogn[0] -
310  0.014118924*phycon.telogn[1] +
311  0.0044980897*phycon.telogn[2] +
312  6.0969594e-5*phycon.telogn[3]);
313  }
314  else
315  {
316  /* same expansion but for hydrogen ions */
317  ThinCoolingCaseB = (-25.859117 +
318  0.16229407*(phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]) +
319  0.34912863*POW2(phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]) -
320  0.10615964*powi( (phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]),3) )/(1. +
321  0.050866793*(phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]) -
322  0.014118924*POW2(phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]) +
323  0.0044980897*powi( (phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]),3) +
324  6.0969594e-5*powi( (phycon.telogn[0]-phycon.sqlogz[nelem-ipISO]),4) );
325  }
326 
327  /* now convert to linear cooling coefficient */
328  ThinCoolingCaseB = POW3(1.+nelem-ipISO)*pow(10.,ThinCoolingCaseB)/(phycon.te/POW2(1.+nelem-ipISO) );
329 
330  /* this is the error, expect positive since do not include infinite number of levels */
331  RecCoolExtra = ThinCoolingCaseB - ThinCoolingSum;
332  }
333  else
334  {
335  ThinCoolingCaseB = 0.;
336  RecCoolExtra = 0.;
337  }
338 
339  /* don't let the extra be negative - should be positive if H-like, negative for
340  * he-like only due to real difference in recombination coefficients */
341  RecCoolExtra = MAX2(0., RecCoolExtra );
342 
343  /* add error onto total - this is significant for approach to STE */
344  iso.RadRecCool[ipISO][nelem] += RecCoolExtra* edenIonAbund *iso.RadRecomb[ipISO][nelem][iso.numLevels_local[ipISO][nelem]-1][ipRecNetEsc];
345 
346  /***********************************************************************
347  *
348  * heating by photoionization of
349  * excited states of all species
350  *
351  ***********************************************************************/
352 
353  /* photoionization of excited levels */
354  HeatExcited = 0.;
355  ipbig = -1000;
356  for( n=1; n < (iso.numLevels_local[ipISO][nelem] - 1); ++n )
357  {
358  SavePhotoHeat[n] = dense.xIonDense[nelem][nelem+1-ipISO]*StatesElem[ipISO][nelem][n].Pop*
359  iso.PhotoHeat[ipISO][nelem][n];
360  HeatExcited += SavePhotoHeat[n];
361  if( SavePhotoHeat[n] > biggest )
362  {
363  biggest = SavePhotoHeat[n];
364  ipbig = (int)n;
365  }
366  }
367  {
368  /* debug block for heating due to photo of each n */
369  enum {DEBUG_LOC=false};
370  if( DEBUG_LOC && ipISO==0 && nelem==0 && nzone > 700)
371  {
372  /* this was not done above */
373  SavePhotoHeat[ipH1s] = dense.xIonDense[nelem][nelem+1-ipISO]*StatesElem[ipISO][nelem][ipH1s].Pop*
374  iso.PhotoHeat[ipISO][nelem][ipH1s];
375  fprintf(ioQQQ,"ipISO = %li nelem=%li ipbig=%li biggest=%g HeatExcited=%.2e ctot=%.2e\n",
376  ipISO , nelem,
377  ipbig ,
378  biggest,
379  HeatExcited ,
380  thermal.ctot);
381  /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */
382  for(n=ipH1s; n< (iso.numLevels_local[ipISO][nelem] - 1); ++n )
383  {
384  fprintf(ioQQQ,"DEBUG phot heat%2li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
385  n,
386  SavePhotoHeat[n]/HeatExcited,
387  dense.xIonDense[nelem][nelem+1-ipISO],
388  StatesElem[ipISO][nelem][n].Pop,
389  iso.PhotoHeat[ipISO][nelem][n],
390  iso.gamnc[ipISO][nelem][n]);
391  }
392  }
393  }
394 
395  /* FreeBnd_net_Cool_Rate is net cooling due to recombination
396  * RadRecCool is total radiative recombination cooling sum to all levels,
397  * with n>=2 photoionization heating subtracted */
398  iso.FreeBnd_net_Cool_Rate[ipISO][nelem] = iso.RadRecCool[ipISO][nelem] - HeatExcited;
399  /*fprintf(ioQQQ,"DEBUG Hn2\t%.3e\t%.3e\n",
400  -iso.RadRecCool[ipISO][nelem]/SDIV(thermal.htot),
401  HeatExcited/SDIV(thermal.htot));*/
402 
403  /* heating[0][1] is all excited state photoionization heating from ALL
404  * species, this is set to zero in CoolEvaluate before loop where this
405  * is called. */
406  thermal.heating[0][1] += MAX2(0.,-iso.FreeBnd_net_Cool_Rate[ipISO][nelem]);
407 
408  /* net free-bound cooling minus free-free heating */
409  /* create a meaningful label */
410  sprintf(chLabel , "ISrcol%2s%2s" , elementnames.chElementSym[ipISO] ,
411  elementnames.chElementSym[nelem]);
412  CoolAdd(chLabel, (realnum)nelem, MAX2(0.,iso.FreeBnd_net_Cool_Rate[ipISO][nelem]));
413 
414  /* if rec coef goes as T^-.8, but times T, so propto t^.2 */
415  thermal.dCooldT += 0.2*iso.FreeBnd_net_Cool_Rate[ipISO][nelem]*phycon.teinv;
416 
417  /***********************************************************************
418  * *
419  * induced recombination cooling *
420  * *
421  ***********************************************************************/
422 
423  iso.RecomInducCool_Rate[ipISO][nelem] = 0.;
424  /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */
425  for( n=0; n < (iso.numLevels_local[ipISO][nelem] - 1); ++n )
426  {
427  /* >>chng 02 jan 22, removed cinduc, replace with RecomInducCool */
428  SaveInducCool[n] = iso.RecomInducCool_Coef[ipISO][nelem][n]*iso.PopLTE[ipISO][nelem][n]*edenIonAbund;
429  iso.RecomInducCool_Rate[ipISO][nelem] += SaveInducCool[n];
430  thermal.dCooldT += SaveInducCool[n]*
431  (iso.xIsoLevNIonRyd[ipISO][nelem][n]/phycon.te_ryd - 1.5)*phycon.teinv;
432  }
433 
434  {
435  /* print rec cool, induced rec cool, photo heat */
436  enum {DEBUG_LOC=false};
437  if( DEBUG_LOC && ipISO==0 && nelem==5 )
438  {
439  fprintf(ioQQQ," ipISO=%li nelem=%li ctot = %.2e\n",
440  ipISO,
441  nelem,
442  thermal.ctot);
443  fprintf(ioQQQ,"sum\t%.2e\t%.2e\t%.2e\n",
444  HeatExcited,
445  iso.RadRecCool[ipISO][nelem],
446  iso.RecomInducCool_Rate[ipISO][nelem]);
447  fprintf(ioQQQ,"sum\tp ht\tr cl\ti cl\n");
448 
449  /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */
450  for(n=0; n< (iso.numLevels_local[ipISO][nelem] - 1); ++n )
451  {
452  fprintf(ioQQQ,"%li\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e \n",
453  n,
454  SavePhotoHeat[n],
455  SaveRadRecCool[n],
456  SaveInducCool[n] ,
457  iso.RecomInducCool_Coef[ipISO][nelem][n],
458  iso.PopLTE[ipISO][nelem][n],
459  iso.RecomInducRate[ipISO][nelem][n]);
460  }
461  fprintf(ioQQQ," \n");
462  }
463  }
464  /* create a meaningful label - induced rec cooling */
465  sprintf(chLabel , "ISicol%2s%2s" , elementnames.chElementSym[ipISO] ,
466  elementnames.chElementSym[nelem]);
467  /* induced rec cooling */
468  CoolAdd(chLabel,(realnum)nelem,iso.RecomInducCool_Rate[ipISO][nelem]);
469 
470  /* find total collisional energy exchange due to bound-bound */
471  iso.xLineTotCool[ipISO][nelem] = 0.;
472  dCdT_all = 0.;
473  heat_max = 0.;
474  nlo_heat_max = -1;
475  nhi_heat_max = -1;
476 
477  /* loop does not include highest levels - their population may be
478  * affected by topoff */
479  for( ipLo=0; ipLo < iso.numLevels_local[ipISO][nelem]-2; ipLo++ )
480  {
481  for( ipHi=ipLo + 1; ipHi < iso.numLevels_local[ipISO][nelem]-1; ipHi++ )
482  {
483  /* collisional energy exchange between ipHi and ipLo - net cool */
484  hlone =
485  Transitions[ipISO][nelem][ipHi][ipLo].Coll.ColUL*
486  (StatesElem[ipISO][nelem][ipLo].Pop*
487  iso.Boltzmann[ipISO][nelem][ipHi][ipLo]*
488  StatesElem[ipISO][nelem][ipHi].g/StatesElem[ipISO][nelem][ipLo].g -
489  StatesElem[ipISO][nelem][ipHi].Pop)*edenHCorr_IonAbund*
490  Transitions[ipISO][nelem][ipHi][ipLo].EnergyErg;
491 
492  iso.xLineTotCool[ipISO][nelem] += hlone;
493 
494  /* next get derivative */
495  if( hlone > 0. )
496  {
497  /* usual case, this line was a net coolant - for derivative
498  * taking the exponential from ground gave better
499  * representation of effects */
500  dCdT_all += hlone*
501  (Transitions[ipISO][nelem][ipHi][ipH1s].EnergyK*thermal.tsq1 - thermal.halfte);
502  }
503  else
504  {
505  /* this line heats the gas, remember which one it was,
506  * the following three vars never are used, but could be for
507  * debugging */
508  if( hlone < heat_max )
509  {
510  heat_max = hlone;
511  nlo_heat_max = ipLo;
512  nhi_heat_max = ipHi;
513  }
514  dCdT_all -= hlone*thermal.halfte;
515  }
516  }
517  }
518  {
519  /* debug block announcing which line was most important */
520  enum {DEBUG_LOC=false};
521  if( DEBUG_LOC )
522  {
523  if( nelem==ipISO )
524  fprintf(ioQQQ,"%li %li %.2e\n", nlo_heat_max, nhi_heat_max, heat_max );
525  }
526  }
527 
528  /* Lyman line collisional heating/cooling */
529  /* Lya itself */
530  iso.cLya_cool[ipISO][nelem] = 0.;
531  /* lines higher than Lya */
532  iso.cLyrest_cool[ipISO][nelem] = 0.;
533 
534  for( ipHi=1; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ )
535  {
536  hlone = Transitions[ipISO][nelem][ipHi][ipH1s].Coll.ColUL*
537  (StatesElem[ipISO][nelem][0].Pop*iso.Boltzmann[ipISO][nelem][ipHi][0]*
538  StatesElem[ipISO][nelem][ipHi].g/StatesElem[ipISO][nelem][0].g -
539  StatesElem[ipISO][nelem][ipHi].Pop)* edenHCorr_IonAbund*
540  Transitions[ipISO][nelem][ipHi][0].EnergyErg;
541 
542  if( ipHi == iso.nLyaLevel[ipISO] )
543  {
544  iso.cLya_cool[ipISO][nelem] = hlone;
545 
546  }
547  else
548  {
549  /* sum energy in higher lyman lines */
550  iso.cLyrest_cool[ipISO][nelem] += hlone;
551  }
552  }
553 
554  /* Balmer line collisional heating/cooling and derivative
555  * only used for print out, incorrect if not H */
556  iso.cBal_cool[ipISO][nelem] = 0.;
557  for( ipHi=3; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ )
558  {
559  /* single line cooling */
560  ipLo = ipH2s;
561  hlone =
562  Transitions[ipISO][nelem][ipHi][ipLo].Coll.ColUL*(
563  StatesElem[ipISO][nelem][ipLo].Pop*
564  iso.Boltzmann[ipISO][nelem][ipHi][ipLo]*
565  StatesElem[ipISO][nelem][ipHi].g/StatesElem[ipISO][nelem][ipLo].g -
566  StatesElem[ipISO][nelem][ipHi].Pop)*edenHCorr_IonAbund*
567  Transitions[ipISO][nelem][ipHi][ipLo].EnergyErg;
568 
569  ipLo = ipH2p;
570  hlone +=
571  Transitions[ipISO][nelem][ipHi][ipLo].Coll.ColUL*(
572  StatesElem[ipISO][nelem][ipLo].Pop*
573  iso.Boltzmann[ipISO][nelem][ipHi][ipLo]*
574  StatesElem[ipISO][nelem][ipHi].g/StatesElem[ipISO][nelem][ipLo].g -
575  StatesElem[ipISO][nelem][ipHi].Pop)*edenHCorr_IonAbund*
576  Transitions[ipISO][nelem][ipHi][ipLo].EnergyErg;
577 
578  /* this is only used to add to line array */
579  iso.cBal_cool[ipISO][nelem] += hlone;
580  }
581 
582  /* all hydrogen lines from Paschen on up, but not Balmer
583  * only used for printout, incorrect if not H */
584  iso.cRest_cool[ipISO][nelem] = 0.;
585  for( ipLo=3; ipLo < iso.numLevels_local[ipISO][nelem]-1; ipLo++ )
586  {
587  for( ipHi=ipLo + 1; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ )
588  {
589  hlone =
590  Transitions[ipISO][nelem][ipHi][ipLo].Coll.ColUL*(
591  StatesElem[ipISO][nelem][ipLo].Pop*
592  iso.Boltzmann[ipISO][nelem][ipHi][ipLo]*
593  StatesElem[ipISO][nelem][ipHi].g/StatesElem[ipISO][nelem][ipLo].g -
594  StatesElem[ipISO][nelem][ipHi].Pop)*edenHCorr_IonAbund*
595  Transitions[ipISO][nelem][ipHi][ipLo].EnergyErg;
596 
597  iso.cRest_cool[ipISO][nelem] += hlone;
598  }
599  }
600 
601  /* add total line heating or cooling into stacks, derivatives */
602  /* line energy exchange can be either heating or coolant
603  * must add this to total heating or cooling, as appropriate */
604  /* create a meaningful label */
605  sprintf(chLabel , "ISclin%2s%2s" , elementnames.chElementSym[ipISO] ,
606  elementnames.chElementSym[nelem]);
607  if( iso.xLineTotCool[ipISO][nelem] > 0. )
608  {
609  /* species is a net coolant label gives iso sequence, "wavelength" gives element */
610  CoolAdd(chLabel,(realnum)nelem,iso.xLineTotCool[ipISO][nelem]);
611  thermal.dCooldT += dCdT_all;
612  iso.dLTot[ipISO][nelem] = 0.;
613  }
614  else
615  {
616  /* species is a net heat source, thermal.heating[0][23]was set to 0 in CoolEvaluate*/
617  thermal.heating[0][23] -= iso.xLineTotCool[ipISO][nelem];
618  CoolAdd(chLabel,(realnum)nelem,0.);
619  iso.dLTot[ipISO][nelem] = -dCdT_all;
620  }
621 
622  {
623  /* debug print for understanding contributors to HLineTotCool */
624  enum {DEBUG_LOC=false};
625  if( DEBUG_LOC )
626  {
627  if( nelem == 0 )
628  {
629  fprintf(ioQQQ,"%.2e la %.2f restly %.2f barest %.2f hrest %.2f\n",
630  iso.xLineTotCool[ipISO][nelem] ,
631  iso.cLya_cool[ipISO][nelem]/iso.xLineTotCool[ipISO][nelem] ,
632  iso.cLyrest_cool[ipISO][nelem]/iso.xLineTotCool[ipISO][nelem] ,
633  iso.cBal_cool[ipISO][nelem]/iso.xLineTotCool[ipISO][nelem] ,
634  iso.cRest_cool[ipISO][nelem]/iso.xLineTotCool[ipISO][nelem] );
635  }
636  }
637  }
638  {
639  /* this is an option to print out each rate affecting each level in strict ste,
640  * this is a check that the rates are indeed in detailed balance */
641  enum {DEBUG_LOC=false};
642  enum {LTEPOP=true};
643  /* special debug print for gas near STE */
644  if( DEBUG_LOC && (nelem==1 || nelem==0) )
645  {
646  /* LTEPOP is option to assume STE for rates */
647  if( LTEPOP )
648  {
649  /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */
650  for(n=ipH1s; n<iso.numLevels_local[ipISO][nelem]-1; ++n )
651  {
652  fprintf(ioQQQ,"%li\t%li\t%g\t%g\t%g\t%g\tT=\t%g\t%g\t%g\t%g\n", nelem,n,
653  iso.gamnc[ipISO][nelem][n] *iso.PopLTE[ipISO][nelem][n],
654  /* induced recom, intergral times hlte */
655  /*iso.RadRecomb[ipISO][nelem][n][ipRecRad]+iso.rinduc[ipISO][nelem][n] ,*/
656  /* >>chng 02 jan 22, remove rinduc, replace with RecomInducRate */
657  iso.RadRecomb[ipISO][nelem][n][ipRecRad]+
658  iso.RecomInducRate[ipISO][nelem][n]*iso.PopLTE[ipISO][nelem][n] ,
659  /* induced rec */
660  iso.RecomInducRate[ipISO][nelem][n]*iso.PopLTE[ipISO][nelem][n] ,
661  /* spontaneous recombination */
662  iso.RadRecomb[ipISO][nelem][n][ipRecRad] ,
663  /* heating, followed by two processes that must balance it */
664  iso.PhotoHeat[ipISO][nelem][n]*iso.PopLTE[ipISO][nelem][n],
665  iso.RecomInducCool_Coef[ipISO][nelem][n]*iso.PopLTE[ipISO][nelem][n]+SaveRadRecCool[n] ,
666  /* induced rec cooling, integral times hlte */
667  iso.RecomInducCool_Coef[ipISO][nelem][n]*iso.PopLTE[ipISO][nelem][n] ,
668  /* free-bound cooling per unit vol */
669  SaveRadRecCool[n] );
670  }
671  }
672  else
673  {
674  /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */
675  for(n=ipH1s; n<iso.numLevels_local[ipISO][nelem]-1; ++n )
676  {
677  fprintf(ioQQQ,"%li\t%li\t%g\t%g\t%g\t%g\tT=\t%g\t%g\t%g\t%g\n", nelem,n,
678  iso.gamnc[ipISO][nelem][n] *dense.xIonDense[nelem][nelem+1-ipISO]*StatesElem[ipISO][nelem][n].Pop,
679  /* induced recom, intergral times hlte */
680  iso.RadRecomb[ipISO][nelem][n][ipRecRad]*edenIonAbund+
681  iso.RecomInducRate[ipISO][nelem][n]*iso.PopLTE[ipISO][nelem][n] *edenIonAbund ,
682  iso.RadRecomb[ipISO][nelem][n][ipRecRad]*edenIonAbund ,
683  iso.RecomInducRate[ipISO][nelem][n]*iso.PopLTE[ipISO][nelem][n] *edenIonAbund ,
684  /* heating, followed by two processes that must balance it */
685  SavePhotoHeat[n],
686  SaveInducCool[n]+SaveRadRecCool[n]*edenIonAbund ,
687  /* induced rec cooling, integral times hlte */
688  SaveInducCool[n] ,
689  /* free-bound cooling per unit vol */
690  SaveRadRecCool[n]*edenIonAbund );
691  }
692  }
693  }
694  }
695  return;
696 }
697 #if defined(__HP_aCC)
698 #pragma OPTIMIZE OFF
699 #pragma OPTIMIZE ON
700 #endif

Generated for cloudy by doxygen 1.8.1.2