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