cloudy  trunk
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cool_eval.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 /*CoolEvaluate main routine to call others, to evaluate total cooling */
4 #include "cddefines.h"
5 #include "cool_eval.h"
6 #include "taulines.h"
7 #include "wind.h"
8 #include "coolheavy.h"
9 #include "radius.h"
10 #include "conv.h"
11 #include "h2.h"
12 #include "rt.h"
13 #include "opacity.h"
14 #include "ionbal.h"
15 #include "trace.h"
16 #include "dynamics.h"
17 #include "grainvar.h"
18 #include "atmdat.h"
19 #include "atoms.h"
20 #include "called.h"
21 #include "hmi.h"
22 #include "numderiv.h"
23 #include "magnetic.h"
24 #include "phycon.h"
25 #include "hyperfine.h"
26 #include "iso.h"
27 #include "thermal.h"
28 #include "cooling.h"
29 #include "pressure.h"
30 #include "mole.h"
31 #include "rfield.h"
32 #include "doppvel.h"
33 #include "freebound.h"
34 #include "dense.h"
35 #include "atmdat_gaunt.h"
36 #include "iterations.h"
37 #include "vectorize.h"
38 
39 
40 /* H2 cooling according to
41  * >>refer h2 cool Glover & Abel 2008, MNRAS, 388, 1627; Sections 2.3.1-2.3.6 */
42 STATIC double CoolH2_GA08 ( double Temp );
43 
44 /*fndneg search cooling array to find negative values */
45 STATIC void fndneg(void);
46 
47 /*fndstr search cooling stack to find strongest values */
48 STATIC void fndstr(double tot, double dc);
49 
50 STATIC void CoolHyperfine (void);
51 
52 STATIC double eeBremsCooling(double Te, // electron temperature in K
53  double eden); // electron density in cm^-3
54 
55 /* set true to debug derivative of heating and cooling */
56 static const bool PRT_DERIV = false;
57 
58 void CoolEvaluate(double *tot)
59 {
60  static long int nhit = 0,
61  nzSave=0;
62 
63  static double oltcool=0.,
64  oldtemp=0.;
65 
66  DEBUG_ENTRY( "CoolEvaluate()" );
67 
68  /* returns tot, the total cooling,
69  * and dc, the derivative of the cooling */
70 
71  if( trace.lgTrace )
72  fprintf( ioQQQ, " COOLR TE:%.4e zone %li %li Cool:%.4e Heat:%.4e eden:%.4e edenTrue:%.4e\n",
73  phycon.te,
76 
77  /* must call TempChange since ionization has changed, there are some
78  * terms that affect collision rates (H0 term in electron collision) */
79  TempChange(phycon.te , false);
80 
81  /* now zero out the cooling stack */
82  CoolZero();
83  if( PRT_DERIV )
84  fprintf(ioQQQ,"DEBUG dCdT 0 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
85  if( gv.lgGrainPhysicsOn )
86  {
87  /* grain heating and cooling */
88  /* grain recombination cooling, evaluated elsewhere
89  * can either heat or cool the gas, do cooling here */
90  CoolAdd("dust",0,MAX2(0.,gv.GasCoolColl));
91 
92  /* grain cooling proportional to temperature ^3/2 */
93  thermal.dCooldT += MAX2(0.,gv.GasCoolColl)*3./(2.*phycon.te);
94 
95  /* these are the various heat agents from grains */
96  /* options to force gas heating or cooling by grains to zero - for tests only ! */
97  if( gv.lgDustOn() && gv.lgDHetOn )
98  {
99  /* rate dust heats gas by photoelectric effect */
101 
102  /* if grains hotter than gas then collisions with gas act
103  * to heat the gas, add this in here
104  * a symmetric statement appears in COOLR, where cooling is added on */
105  thermal.setHeating(0,14, MAX2(0.,-gv.GasCoolColl) );
106 
107  /* this is gas heating due to thermionic emissions */
109  }
110  else
111  {
112  thermal.setHeating(0,13,0.);
113  thermal.setHeating(0,14,0.);
114  thermal.setHeating(0,25,0.);
115  }
116  }
117  else if( gv.lgBakesPAH_heat )
118  {
119  /* >>chng 06 jul 21, option to include Bakes PAH hack with grain physics off,
120  * needed to test dynamics models */
122  }
123 
124  if( PRT_DERIV )
125  fprintf(ioQQQ,"DEBUG dCdT 1 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
126 
127  /* molecular molecules molecule cooling */
128  if( mole_global.lgNoMole )
129  {
130  /* this branch - do not include molecules */
131  hmi.hmicol = 0.;
133  /* line cooling within simple H2 molecule - zero when big used */
134  CoolHeavy.h2line = 0.;
135  /* H + H+ => H2+ cooling */
136  CoolHeavy.H2PlsCool = 0.;
137  CoolHeavy.HD = 0.;
138 
139  /* thermal.heating(0,8) is heating due to collisions within X of H2 */
140  thermal.setHeating(0,8,0.);
141  /* thermal.heating(0,15) is H minus heating*/
142  thermal.setHeating(0,15,0.);
143  /* thermal.heating(0,16) is H2+ heating */
144  thermal.setHeating(0,16,0.);
145  hmi.HeatH2Dish_used = 0.;
146  hmi.HeatH2Dexc_used = 0.;
148  }
149 
150  else
151  {
152  /* save various molecular heating/cooling agent */
153  thermal.setHeating(0,15, hmi.hmihet);
155  /* now get heating from H2 molecule, either simple or from big one */
157  {
158  if( h2.lgEvaluated )
159  {
160  /* these are explicitly from big H2 molecule,
161  * first is heating due to radiative pump of excited states, followed by
162  * radiative decay into continuum of X, followed by dissociation of molecule
163  * with kinetic energy, typically 0.25 - 0.5 eV per event */
166  if (0)
167  fprintf(ioQQQ,"DEBUG big %.2f\t%.5e\t%.2e\t%.2e\t%.2e\n",
170  /* negative sign because right term is really deriv of heating,
171  * but will be used below as deriv of cooling */
173  }
174  else
175  {
176  hmi.HeatH2Dish_used = 0;
177  hmi.HeatH2Dexc_used = 0;
179  }
180  }
181 
182  else if( hmi.chH2_small_model_type == 'T' )
183  {
184  /* TH85 dissociation heating */
185  /* these come from approximations in TH85, see comments above */
189  }
190  else if( hmi.chH2_small_model_type == 'H' )
191  {
192  /* Burton et al. 1990 */
196  }
197  else if( hmi.chH2_small_model_type == 'B' )
198  {
199  /* Bertoldi & Draine */
203  }
204  else if( hmi.chH2_small_model_type == 'E' )
205  {
206  /* this is the default when small H2 used */
210  }
211  else
212  TotalInsanity();
213 
214  /* heating due to photodissociation heating */
216 
217  /* heating due to continuum photodissociation */
218  thermal.setHeating(0,28,0.);
219  {
220  double heat = 0.;
221  for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
222  {
223  if( (*diatom)->lgEnabled && mole_global.lgStancil )
224  {
225  heat += (*diatom)->Cont_Diss_Heat_Rate();
226  }
227  }
228  thermal.setHeating(0,28,MAX2( 0., heat ));
229  CoolAdd("H2cD",0,MAX2(0.,-heat));
230  }
231 
232  /* heating (usually cooling in big H2) due to collisions within X */
233  /* add to heating is net heating is positive */
235  // HD heating, cooling is CoolHeavy.HD
236  MAX2(0.,hd.HeatDexc) + MAX2(0.,hd.HeatDiss));
237 
238  /* add to cooling if net heating is negative */
239  CoolAdd("H2cX",0,MAX2(0.,-hmi.HeatH2Dexc_used));
240  /*fprintf(ioQQQ,"DEBUG coolh2\t%.2f\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
241  fnzone, phycon.te, dense.eden, hmi.H2_total, thermal.ctot, -hmi.HeatH2Dexc_used );*/
242  /* add to net derivative */
243  /*thermal.dCooldT += MAX2(0.,-hmi.HeatH2Dexc_used)* ( 30172. * thermal.tsq1 - thermal.halfte );*/
244  /* >>chng 04 jan 25, check sign to prevent cooling from entering here,
245  * also enter neg sign since going into cooling stack (bug), in heatsum
246  * same term adds to deriv of heating */
247  if( hmi.HeatH2Dexc_used < 0. )
249 
250  /* H + H+ => H2+ cooling */
251  CoolHeavy.H2PlsCool = (realnum)(MAX2((2.325*phycon.te-1875.)*1e-20,0.)*
253 
254  if( h2.lgEnabled )
255  {
256  /* this is simplified approximation to H2 rotation cooling,
257  * big molecule does this far better */
258  CoolHeavy.h2line = 0.;
259  }
260  else
261  {
263 
264  {
265  enum {DEBUG_LOC=false};
266  if( DEBUG_LOC && nzone>187&& iteration > 1)
267  {
268  fprintf(ioQQQ,"h2coolbug\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
269  phycon.te,
270  CoolHeavy.h2line,
271  hmi.H2_total,
272  findspecieslocal("H-")->den,
274  }
275  }
276 
277  }
278 
279  if( hd.lgEnabled )
280  {
281  // heating was thermal.setHeating(0,8 above, with H2
282  CoolHeavy.HD = -MIN2(0.,hd.HeatDexc) - MIN2(0.,hd.HeatDiss);
283  }
284  else
285  {
286  /* >>chng 22 aug 13, use Flower et al. (2000, MNRAS, 314, 753)
287  * HD cooling function
288  * http://ccp7.dur.ac.uk/cooling_by_HD/node3.html */
289 #if 0
290  factor = 0.25*pow2( log10(double(dense.gas_phase[ipHYDROGEN])))+
291  (0.283978*log10(double(dense.gas_phase[ipHYDROGEN]))-1.27333)*
292  (sin(2.03657*log10(phycon.te)+4.63258))-2.08189*
293  log10(double(dense.gas_phase[ipHYDROGEN]))+4.66288;
294 
295  CoolHeavy.HD = hmi.HD_total*exp10((0.5*log10(double(dense.gas_phase[ipHYDROGEN]))) +
296  (-26.2982*pow(log10(phycon.te), -0.215807)) - sqrt(factor));
297 #endif
298 
299  double aa=-26.2982, bb=-0.215807, omeg=2.03657, phi=4.63258,
300  c1=0.283978, c2=-1.27333, d1=-2.08189, d2=4.66288;
301 
302  // the sum of hydrogen nuclei in all forms
303  double y = log10(dense.gas_phase[ipHYDROGEN]);
304  double x = phycon.alogte;
305 
306  double w = 0.5 * y + aa * pow(x,bb)
307  - sqrt( 0.25 * y*y
308  + (c1*y+c2) * sin(omeg*x+phi) + (d1*y+d2) );
309 
310  CoolHeavy.HD = hmi.HD_total * exp10(w);
311 
312  }
313  }
314 
315  fixit("test and enable chemical heating by default");
316 #if 0
317  double chemical_heating = mole.chem_heat();
318  thermal.setHeating(0,29, MAX2(0.,chemical_heating) );
319  /* add to cooling if net heating is negative */
320  CoolAdd("Chem",0,MAX2(0.,-chemical_heating));
321 #endif
322 
323  /* cooling due to charge transfer ionization / recombination */
324  CoolAdd("CT C" , 0. , thermal.char_tran_cool );
325 
326  /* H- FB; H + e -> H- + hnu */
327  /* H- FF is in with H ff */
328  CoolAdd("H-fb",0,hmi.hmicol);
329 
330  /* >>chng 96 nov 15, fac of 2 in deriv to help convergence in very dense
331  * models where H- is important, this takes change in eden into
332  * partial account */
334  if( PRT_DERIV )
335  fprintf(ioQQQ,"DEBUG dCdT 2 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
336 
337  CoolAdd("H2ln",0,CoolHeavy.h2line);
338  /* >>chng 00 oct 21, added coef of 3.5, sign had been wrong */
339  /*thermal.dCooldT += CoolHeavy.h2line*phycon.teinv;*/
340  /* >>chng 03 mar 17, change 3.5 to 15 as per behavior in primal.in */
341  /*thermal.dCooldT += 3.5*CoolHeavy.h2line*phycon.teinv;*/
342  /* >>chng 03 may 18, from 15 to 30 as per behavior in primal.in - overshoots happen */
343  /*thermal.dCooldT += 15.*CoolHeavy.h2line*phycon.teinv;*/
344  /*>>chng 03 oct 03, from 30 to 3, no more overshoots in primalin */
345  /*thermal.dCooldT += 30.*CoolHeavy.h2line*phycon.teinv;*/
347 
348  {
349  /* problems with H2 cooling */
350  enum {DEBUG_LOC=false};
351  if( DEBUG_LOC /*&& nzone>300 && iteration > 1*/)
352  {
353  fprintf(ioQQQ,"CoolEvaluate debuggg\t%.2f\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\n",
354  fnzone,
355  phycon.te,
356  hmi.H2_total ,
358  findspecieslocal("H-")->den ,
359  dense.eden);
360  }
361  }
362 
363  CoolAdd("HDro",0,CoolHeavy.HD);
365 
366  CoolAdd("H2+ ",0,CoolHeavy.H2PlsCool);
368 
369  /* heating due to three-body, will be incremented in iso_cool*/
370  thermal.setHeating(0,3,0.);
371  /* heating due to hydrogen lines */
372  thermal.setHeating(0,23,0.);
373  /* heating due to photoionization of all excited states of hydrogen species */
374  thermal.setHeating(0,1,0.);
375 
376  /* isoelectronic species cooling, mainly lines, and level ionization */
377  for( long int ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
378  {
379  for( long int nelem=ipISO; nelem < LIMELM; nelem++ )
380  {
381  /* must always call iso_cool since we must zero those variables
382  * that would have been set had the species been present */
383  iso_cool( ipISO , nelem );
384  }
385  }
386 
387  /* free-free free free brems emission for all ions */
388  /* highest frequency where we have non-zero Boltzmann factors */
389  long limit = min( rfield.ipMaxBolt, rfield.nflux );
390 
391  CoolHeavy.brems_cool_h = 0.;
396 
397  if( CoolHeavy.lgFreeOn )
398  {
400  ASSERT(limit < rfield.nflux_with_check);
401 
402  t_brems_den sum;
404  double bfac = dense.eden * FREE_FREE_EMIS / phycon.sqrte * EN1RYD;
405 
406  /* do hydrogen first, before main loop since want to break out as separate
407  * coolant, and what to add on H- brems */
410 
411  /* now do helium, both He+ and He++ */
413  sum.den_Hepp*bfac*t_gaunt::Inst().brems_cool( 2, phycon.te );
414 
415  /* heavy elements */
417  for( long ion=1; ion < LIMELM+1; ++ion )
418  if( sum.den_ion[ion] > 0. )
420  sum.den_ion[ion]*bfac*t_gaunt::Inst().brems_cool( ion, phycon.te );
421 
422  /* ipEnergyBremsThin is index to energy where gas becomes optically thin to brems,
423  * so this loop is over optically thin frequencies
424  * do not include optically thick part as net emission since self absorbed */
426  for( long int i=rfield.ipEnergyBremsThin; i < limit; i++ )
427  {
428  /* the total heating due to bremsstrahlung */
430  opac.FreeFreeOpacity[i]*rfield.flux[0][i]*rfield.anu(i)*EN1RYD;
431  }
432 
433  enum {DEBUG_LOC=false};
434  if( DEBUG_LOC && nzone>60 /*&& iteration > 1*/)
435  {
436  double sumfield = 0., sumtot=0., sum1=0., sum2=0.;
437  for( long int i=rfield.ipEnergyBremsThin; i<limit; i++ )
438  {
439  sumtot += opac.FreeFreeOpacity[i]*rfield.flux[0][i]*rfield.anu(i);
440  sumfield += rfield.flux[0][i]*rfield.anu(i);
441  sum1 += opac.FreeFreeOpacity[i]*rfield.flux[0][i]*rfield.anu(i);
442  sum2 += opac.FreeFreeOpacity[i]*rfield.flux[0][i];
443  }
444  fprintf(ioQQQ,"DEBUG brems heat\t%.2f\t%.3e\t%.3e\t%.3e\t%e\t%.3e\t%.3e\n",
445  fnzone,
447  sumtot/SDIV(sumfield) ,
448  sum1/SDIV(sum2),
449  phycon.te ,
451  opac.FreeFreeOpacity[1218]);
452  }
453  }
454 
455  /* these two terms are both large, nearly canceling, near lte */
462  /*fprintf(ioQQQ,"DEBUG brems\t%.2f\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
463  fnzone,
464  phycon.te,
465  CoolHeavy.brems_cool_net,
466  CoolHeavy.brems_cool_h ,
467  CoolHeavy.brems_cool_he ,
468  CoolHeavy.brems_cool_hminus,
469  CoolHeavy.brems_cool_metals ,
470  CoolHeavy.brems_heat_total);*/
471 
472  /* net free free brems cooling, count as cooling if positive */
473  CoolAdd( "FF c" , 0, MAX2(0.,CoolHeavy.brems_cool_net) );
474 
475  /* now stuff into heating array if negative */
477 
478  /* >>chng 96 oct 30, from HFFNet to just FreeFreeCool,
479  * since HeatSum picks up CoolHeavy.brems_heat_total */
484  if( PRT_DERIV )
485  fprintf(ioQQQ,"DEBUG dCdT 3 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
486 
487  /* >>chng 02 jun 21, net cooling already includes this */
488  /* end of brems cooling */
489 
490  /* heavy element recombination cooling, do not count hydrogenic since
491  * already done above, also helium singlets have been done */
492  /* >>chng 02 jul 21, put in charge dependent rec term */
493  CoolHeavy.heavfb = 0.;
494  for( long int nelem=ipLITHIUM; nelem < LIMELM; nelem++ )
495  {
496  if( dense.lgElmtOn[nelem] )
497  {
498  /* note that detailed iso seq atoms are done in iso_cool */
499  long limit_lo = MAX2( 1 , dense.IonLow[nelem] );
500  long limit_hi = MIN2( nelem-NISO+1, dense.IonHigh[nelem] );
501  for( long int ion=limit_lo; ion<=limit_hi; ++ion )
502  {
503  /* factor of 0.9 is roughly correct for nebular conditions, see
504  * >>refer H rec cooling LaMothe, J., & Ferland, G.J., 2001, PASP, 113, 165 */
505  /* note that ionbal.RR_rate_coef_used is the rate coef, cm3 s-1, needs eden */
506  /* >>chng 02 nov 07, move rec arrays around, this now has ONLY rad rec,
507  * previously had included grain rec and three body */
508  /* recombination cooling for iso-seq atoms are done in iso_cool */
509  double one = dense.xIonDense[nelem][ion] * ionbal.RR_rate_coef_used[nelem][ion-1]*
510  dense.eden * phycon.te * BOLTZMANN;
511  /*fprintf(ioQQQ,"debugggfb\t%li\t%li\t%.3e\t%.3e\t%.3e\n", nelem, ion, one,
512  dense.xIonDense[nelem][ion] , ionbal.RR_rate_coef_used[nelem][ion]);*/
513  CoolHeavy.heavfb += one;
514  thermal.elementcool[nelem] += one;
515  }
516  }
517  }
518 
519  /*fprintf(ioQQQ,"debuggg hvFB\t%i\t%.2f\t%.2e\t%.2e\n",iteration, fnzone,CoolHeavy.heavfb, dense.eden);*/
520 
521  CoolAdd("hvFB",0,CoolHeavy.heavfb);
523  if( PRT_DERIV )
524  fprintf(ioQQQ,"DEBUG dCdT 4 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
525 
527 
528  /* >>chng 97 mar 12, added deriv */
530  CoolAdd("eeff",0,CoolHeavy.eebrm);
531 
532  /* add advective heating and cooling */
533  /* this is cooling due to loss of matter from this region */
534  CoolAdd("adve",0,dynamics.Cool() );
535  /* >>chng02 dec 04, rm factor of 8 in front of dCooldT */
537  /* local heating due to matter moving into this location */
538  thermal.setHeating(1,5, dynamics.Heat() );
540 
541  /* total Compton cooling */
543  CoolAdd("Comp",0,CoolHeavy.tccool);
545 
546  /* option for "extra" cooling, expressed as power-law in temperature, these
547  * are set with the CEXTRA command */
548  if( thermal.lgCExtraOn )
549  {
550  CoolHeavy.cextxx =
551  (realnum)(thermal.CoolExtra*pow(phycon.te/1e4,(double)thermal.cextpw));
552  }
553  else
554  {
555  CoolHeavy.cextxx = 0.;
556  }
557  CoolAdd("Extr",0,CoolHeavy.cextxx);
558 
559  /* cooling due to wind expansion, only for winds expansion cooling */
560  if( wind.lgBallistic() )
561  {
562  realnum dDensityDT = -(realnum)(wind.AccelTotalOutward/wind.windv + 2.*wind.windv/
563  radius.Radius);
564  CoolHeavy.expans = -2.5*pressure.PresGasCurr*dDensityDT;
565  }
566  else if( dynamics.lgTimeDependentStatic &&
568  {
569  realnum dens = scalingDensity();
570  realnum dDensityDT =
573  // pdV work term
574  CoolHeavy.expans = -pressure.PresGasCurr*dDensityDT;
575  }
576  else
577  {
578  CoolHeavy.expans = 0.;
579  }
580  CoolAdd("Expn",0,CoolHeavy.expans);
582 
583  /* cyclotron cooling */
584  /* coef is 4/3 /8pi /c * vtherm(elec) */
585  CoolHeavy.cyntrn = 4.5433e-25f*magnetic.pressure*PI8*dense.eden*phycon.te;
586  CoolAdd("Cycl",0,CoolHeavy.cyntrn);
588 
589  /* heavy element collisional ionization
590  * derivative should be zero since increased coll ion rate
591  * decreases neutral fraction by proportional amount */
592  CoolAdd("Hvin",0,CoolHeavy.colmet);
593 
594 
595  double xIonDenseSave[LIMELM][LIMELM+1];
597  {
598  for( int nelem=0; nelem < LIMELM; nelem++ )
599  {
600  for( int ion=0; ion<=nelem+1; ++ion )
601  {
602  xIonDenseSave[nelem][ion] = dense.xIonDense[nelem][ion];
603  // zero abundance of species if we are using Chianti for this ion
604  if( dense.lgIonChiantiOn[nelem][ion] || dense.lgIonStoutOn[nelem][ion] )
605  dense.xIonDense[nelem][ion] = 0.;
606  }
607  }
608  }
609 
610  (*TauDummy).Zero();
611  (*(*TauDummy).Hi()).g() = 0.;
612  (*(*TauDummy).Lo()).g() = 0.;
613  (*(*TauDummy).Hi()).IonStg() = 0;
614  (*(*TauDummy).Lo()).IonStg() = 0;
615  (*(*TauDummy).Hi()).nelem() = 0;
616  (*(*TauDummy).Lo()).nelem() = 0;
617  (*TauDummy).Emis().Aul() = 0.;
618  (*TauDummy).EnergyWN() = 0.;
619  (*TauDummy).WLAng() = 0.;
620 
621  // reset abundances to original values, may have been set zero to protect against old cloudy lines
623  {
624  // this clause, first reset abundances set to zero when Chianti included
625  for( int nelem=0; nelem < LIMELM; nelem++ )
626  {
627  for( int ion=0; ion<=nelem+1; ++ion )
628  {
629  dense.xIonDense[nelem][ion] = xIonDenseSave[nelem][ion];
630  }
631  }
632  }
633 
634  /* opacity project lines Dima Verner added with g-bar approximation */
635  CoolDima();
636 
637  /* do external database lines */
638  dBaseTrim();
640  dBase_solve();
641 
642  /* Hyperfine line cooling */
643  CoolHyperfine();
644 
645  if( PRT_DERIV )
646  fprintf(ioQQQ,"DEBUG dCdT 6 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
647 
648  /* Print number of levels for each species */
649  {
650  enum {DEBUG_LOC=false};
651  if( DEBUG_LOC )
652  {
653  static bool lgMustPrintHeader=true;
654  if( lgMustPrintHeader )
655  {
656  lgMustPrintHeader = false;
657  printf("DEBUG Levels\t%s",dBaseSpecies[0].chLabel );
658  for( long ipSpecies=1; ipSpecies<nSpecies; ipSpecies++ )
659  {
660  printf("\t%s",dBaseSpecies[ipSpecies].chLabel );
661  }
662  printf("\n" );
663  printf("DEBUG Max\t%li" ,dBaseSpecies[0].numLevels_max );
664  for( long ipSpecies=1; ipSpecies<nSpecies; ipSpecies++ )
665  {
666  printf( "\t%li" ,dBaseSpecies[ipSpecies].numLevels_max );
667  }
668  printf("\n");
669  }
670  printf("DEBUG Local\t%li" ,dBaseSpecies[0].numLevels_local );
671  for( long ipSpecies=1; ipSpecies<nSpecies; ipSpecies++ )
672  {
673  printf("\t%li" ,dBaseSpecies[ipSpecies].numLevels_local );
674  }
675  printf("\n");
676  }
677  }
678 
679  /* now add up all the coolants */
680  CoolSum(tot);
681  if( PRT_DERIV )
682  fprintf(ioQQQ,"DEBUG dCdT 14 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
683 
684  /* negative cooling */
685  if( *tot <= 0. )
686  {
687  fprintf( ioQQQ, " COOLR; cooling is <=0, this is impossible.\n" );
688  ShowMe();
690  }
691 
692  /* bad derivative */
693  if( thermal.dCooldT == 0. )
694  {
695  fprintf( ioQQQ, " COOLR; cooling slope <=0, this is impossible.\n" );
696  if( *tot > 0. && dense.gas_phase[ipHYDROGEN] < 1e-4 )
697  {
698  fprintf( ioQQQ, " Probably due to very low density.\n" );
699  }
700  ShowMe();
702  }
703 
704  if( trace.lgTrace )
705  {
706  fndstr(*tot,thermal.dCooldT);
707  }
708 
709  /* lgTSetOn true for constant temperature model */
710  if( (((((!thermal.lgTemperatureConstant) && *tot < 0.) && called.lgTalk) &&
711  !conv.lgSearch) && thermal.lgCNegChk) && nzone > 0 )
712  {
713  fprintf( ioQQQ,
714  " NOTE Negative cooling, zone %4ld, =%10.2e coola=%10.2e CHION=%10.2e Te=%10.2e\n",
715  nzone,
716  *tot,
717  iso_sp[ipH_LIKE][ipHYDROGEN].cLya_cool,
718  iso_sp[ipH_LIKE][ipHYDROGEN].coll_ion,
719  phycon.te );
720  fndneg();
721  }
722 
723  /* possibility of getting empirical cooling derivative
724  * normally false, set true with 'set numerical derivatives' command */
725  if( NumDeriv.lgNumDeriv )
726  {
727  if( ((nzone > 2 && nzone == nzSave) && ! fp_equal( oldtemp, phycon.te )) && nhit > 4 )
728  {
729  /* hnit is number of tries on this zone - use to stop numerical problems
730  * do not evaluate numerical deriv until well into solution */
731  thermal.dCooldT = (oltcool - *tot)/(oldtemp - phycon.te);
732  }
733  if( nzone != nzSave )
734  nhit = 0;
735 
736  nzSave = nzone;
737  nhit += 1;
738  oltcool = *tot;
739  oldtemp = phycon.te;
740  }
741  return;
742 }
743 
744 
745 /*******************************************************************************
746  *
747  * GLOVER & ABEL (2008) H2 LOW-DENSITY COOLING
748  *
749  ******************************************************************************/
750 
751 /* Implementation of log10 of H2 cooling rate, given by
752  * >>refer eqn(27) Glover & Abel 2008, MNRAS, 388, 1627 */
753 STATIC double ga08_sum ( double Temp, double *coeff, int ncoeff )
754 {
755  DEBUG_ENTRY( "ga08_sum()" );
756  double logT = log10( Temp * 0.001 );
757  double logTn = 1.0;
758  double sum = 0.;
759  for( int i = 0; i < ncoeff; i++ )
760  {
761  sum += coeff[i] * logTn;
762  logTn *= logT;
763  }
764 
765  return sum;
766 }
767 
768 /* Cooling rate of ortho-H2 due to collisions with atomic H at T <= 100 K
769  * >>refer eqn(28) Glover & Abel 2008, MNRAS, 388, 1627 */
770 STATIC double ga08_oH2_H_b100 ( double Temp )
771 {
772  DEBUG_ENTRY( "ga08_oH2_H_b100()" );
773  const double DE = 852.5; // energy in K
774  return 5.09e-27 * pow( Temp * 0.001, 0.5) * dsexp( DE / Temp );
775 }
776 
777 /* Cooling rate of ortho-H2 due to collisions with atomic H at 100 K < T < 1000K
778  * >>refer eqn(27) & Table 1 second column Glover & Abel 2008, MNRAS, 388, 1627 */
779 STATIC double ga08_oH2_H_b1000 ( double Temp )
780 {
781  DEBUG_ENTRY( "ga08_oH2_H_b1000()" );
782  double coeff_oH2_H_cool[] = { -24.330855, 4.4404496, -4.0460989, -1.1390725, 9.8094223, 8.6273872 };
783  int ncoeff = sizeof( coeff_oH2_H_cool ) / sizeof( coeff_oH2_H_cool[0] );
784  return ga08_sum( Temp, coeff_oH2_H_cool, ncoeff );
785 }
786 
787 /* Cooling rate of ortho-H2 due to collisions with atomic H at 1000 K <= T < 6000K
788  * >>refer eqn(27) & Table 1 third column Glover & Abel 2008, MNRAS, 388, 1627 */
789 STATIC double ga08_oH2_H_b6000 ( double Temp )
790 {
791  DEBUG_ENTRY( "ga08_oH2_H_b6000()" );
792  /* >>refer h2 cool H Glover & Abel 2008, MNRAS, 388, 1627; Table 1, third column */
793  double coeff_oH2_H_warm[] = { -24.329086, 4.6105087, -3.9505350, 12.363818, -32.403165, 48.853562, -38.542008, 12.066770 };
794  int ncoeff = sizeof( coeff_oH2_H_warm ) / sizeof( coeff_oH2_H_warm[0] );
795  return ga08_sum( Temp, coeff_oH2_H_warm, ncoeff );
796 }
797 
798 /* Cooling rate of ortho-H2 due to collisions with atomic H at T > 6000K
799  * Artificial tail to give the cooling rate a finite value across temperature range */
800 /* Modify the ortho-H2 cooling rate due to atomic H collisions for a small
801  * temperature range above 100 K to obtain a continuous function across 100 K.
802  * Returns log10 of correction. */
803 STATIC double ga08_oH2_H_stitch_100 ( double Temp )
804 {
805  DEBUG_ENTRY( "ga08_oH2_H_stitch_100()" );
806  const double Tbreak = 100.;
807  static double Temp_c = -1.;
808  if( Temp_c < 0. )
809  {
810  // _bb -> below break; _ab -> above break
811  double crate_bb = log10( ga08_oH2_H_b100( Tbreak ) );
812  double crate_ab = ga08_oH2_H_b1000( Tbreak ); // this is log10( rate )
813  Temp_c = Tbreak / ( 1. + (crate_ab - crate_bb) / LOG10_E );
814  if( 0 )
815  printf("oH2 H, at 100K:\t %f (%.10e vs %.10e)\n", Temp_c, crate_bb, crate_ab);
816  }
817 
818  double corr = 0.;
819  if( Temp > Tbreak && Temp <= Temp_c )
820  corr = LOG10_E * ( 1. - Temp / Temp_c );
821  return corr;
822 }
823 
824 /* Modify the ortho-H2 cooling rate due to atomic H collisions for a small
825  * temperature range below 1000 K to obtain a continuous function across 1000 K.
826  * Returns log10 of correction. */
827 STATIC double ga08_oH2_H_stitch_1000 ( double Temp )
828 {
829  DEBUG_ENTRY( "ga08_oH2_H_stitch_1000()" );
830  const double Tbreak = 1000.;
831  static double Temp_c = -1.;
832  if( Temp_c < 0. )
833  {
834  // _bb -> below break; _ab -> above break
835  double crate_bb = ga08_oH2_H_b1000( Tbreak );
836  double crate_ab = ga08_oH2_H_b6000( Tbreak ); // this is log10( rate )
837  Temp_c = Tbreak / ( 1. + (crate_ab - crate_bb) / LOG10_E );
838  if( 0 )
839  printf("oH2 H, at 1000K:\t %f (%.10e vs %.10e)\n", Temp_c, crate_bb, crate_ab);
840  }
841 
842  double corr = 0.;
843  if( Temp >= Temp_c && Temp_c < Tbreak )
844  corr = LOG10_E * ( Temp / Temp_c - 1. );
845  return corr;
846 }
847 
848 /* Cooling rate of ortho-H2 due to collisions with atomic H at any T.
849  * Extrapolate below 100 K. */
850 STATIC double ga08_oH2_H ( double Temp )
851 {
852  DEBUG_ENTRY( "ga08_oH2_H()" );
853  static double cool_above_6000 = -1.;
854  if( cool_above_6000 < 0. )
855  {
856  cool_above_6000 = pow( 10., ga08_oH2_H_b6000( 6000.0 ) );
857  }
858 
859  double cool = 0.;
860  if( Temp <= 100. )
861  {
862  cool = ga08_oH2_H_b100( Temp );
863  }
864  else if( Temp < 1000. )
865  {
866  cool = ga08_oH2_H_b1000( Temp );
867  cool += ga08_oH2_H_stitch_100( Temp );
868  cool += ga08_oH2_H_stitch_1000( Temp );
869  cool = pow(10., cool);
870  }
871  else if( Temp < 6000. )
872  {
873  cool = pow( 10., ga08_oH2_H_b6000( Temp ) );
874  }
875  else
876  {
877  cool = cool_above_6000;
878  }
879  return cool;
880 }
881 
882 
883 
884 /* Cooling rate of para-H2 due to collisions with atomic H at T <= 100 K
885  * >>refer eqn(29) Glover & Abel 2008, MNRAS, 388, 1627 */
886 STATIC double ga08_pH2_H_b100 ( double Temp )
887 {
888  DEBUG_ENTRY( "ga08_pH2_H_b100()" );
889  const double E20 = 509.85; // energy in K
890  return 8.16e-26 * pow(Temp * 0.001, 0.5) * dsexp( E20 / Temp );
891 }
892 
893 /* Cooling rate of para-H2 due to collisions with atomic H at 100 K < T < 1000K
894  * >>refer eqn(27) & Table 2 second column Glover & Abel 2008, MNRAS, 388, 1627 */
895 STATIC double ga08_pH2_H_b1000 ( double Temp )
896 {
897  DEBUG_ENTRY( "ga08_pH2_H_b1000()" );
898  double coeff_pH2_H_cool[] = { -24.216387, 3.3237480, -11.642384, -35.553366, -35.105689, -10.922078 };
899  int ncoeff = sizeof( coeff_pH2_H_cool ) / sizeof( coeff_pH2_H_cool[0] );
900  return ga08_sum( Temp, coeff_pH2_H_cool, ncoeff );
901 }
902 
903 /* Cooling rate of para-H2 due to collisions with atomic H at 1000 K <= T < 6000K
904  * >>refer eqn(27) & Table 2 third column Glover & Abel 2008, MNRAS, 388, 1627 */
905 STATIC double ga08_pH2_H_b6000 ( double Temp )
906 {
907  DEBUG_ENTRY( "ga08_pH2_H_b6000()" );
908  double coeff_pH2_H_warm[] = { -24.216387, 4.2046488, -1.3155285, -1.6552763, 4.1780102, -0.56949697, -3.3824407, 1.0904027 };
909  int ncoeff = sizeof( coeff_pH2_H_warm ) / sizeof( coeff_pH2_H_warm[0] );
910  return ga08_sum( Temp, coeff_pH2_H_warm, ncoeff );
911 }
912 
913 /* Modify the para-H2 cooling rate due to atomic H collisions for a small
914  * temperature range above 100 K to obtain a continuous function across 100 K.
915  * Returns log10 of correction. */
916 STATIC double ga08_pH2_H_stitch_100 ( double Temp )
917 {
918  DEBUG_ENTRY( "ga08_pH2_H_stitch_100()" );
919  const double Tbreak = 100.;
920  static double Temp_c = -1.;
921  if( Temp_c < 0. )
922  {
923  // _bb -> below break; _ab -> above break
924  double crate_bb = log10( ga08_pH2_H_b100( Tbreak ) );
925  double crate_ab = ga08_pH2_H_b1000( Tbreak ); // this is log10( rate )
926  Temp_c = Tbreak / ( 1. + (crate_ab - crate_bb) / LOG10_E );
927  if( 0 )
928  printf("pH2 H, at 100K:\t %f\n", Temp_c);
929  }
930 
931  double corr = 0.;
932  if( Temp > Tbreak && Temp <= Temp_c )
933  corr = LOG10_E * ( 1. - Temp / Temp_c );
934  return corr;
935 }
936 
937 /* Cooling rate of para-H2 due to collisions with atomic H at any T.
938  * Extrapolate below 100 K. */
939 STATIC double ga08_pH2_H ( double Temp )
940 {
941  DEBUG_ENTRY( "ga08_pH2_H()" );
942  static double cool_above_6000 = -1.;
943  if( cool_above_6000 < 0. )
944  {
945  cool_above_6000 = pow( 10., ga08_pH2_H_b6000( 6000.0 ) );
946  }
947 
948  double cool = 0.;
949  if( Temp <= 100. )
950  {
951  cool = ga08_pH2_H_b100( Temp );
952  }
953  else
954  if( Temp < 1000. )
955  {
956  cool = ga08_pH2_H_b1000( Temp );
957  cool += ga08_pH2_H_stitch_100( Temp );
958  cool = pow(10., cool);
959  }
960  else if( Temp < 6000. )
961  {
962  cool = pow( 10., ga08_pH2_H_b6000( Temp ) );
963  }
964  else
965  {
966  cool = cool_above_6000;
967  }
968  return cool;
969 }
970 
971 
972 
973 /* Cooling rate of para-H2 due to collisions with para-H2 at 100 K <= T <= 6000 K
974  * >>refer h2 cool para-para Glover & Abel 2008, MNRAS, 388, 1627; Table 3, second column */
975 STATIC double ga08_pH2_pH2_a100_b6000 ( double Temp )
976 {
977  DEBUG_ENTRY( "ga08_pH2_pH2_a100_b6000()" );
978  double coeff_pH2_pH2[] = { -23.889798, 1.8550774, -0.55593388, 0.28429361, -0.20581113, 0.13112378 };
979  int ncoeff = sizeof( coeff_pH2_pH2 ) / sizeof( coeff_pH2_pH2[0] );
980  return ga08_sum( Temp, coeff_pH2_pH2, ncoeff );
981 }
982 
983 /* Cooling rate of para-H2 due to collisions with para-H2 at any T.
984  * Extrapolate below 100 K. */
985 STATIC double ga08_pH2_pH2 ( double Temp )
986 {
987  DEBUG_ENTRY( "ga08_pH2_pH2()" );
988  static double cool_above_6000 = -1.;
989  if( cool_above_6000 < 0. )
990  {
991  cool_above_6000 = pow( 10., ga08_pH2_pH2_a100_b6000( 6000. ) );
992  }
993 
994  double cool = 0.;
995  if( Temp <= 6000. )
996  cool = pow( 10., ga08_pH2_pH2_a100_b6000( Temp ) );
997  else
998  cool = cool_above_6000;
999  return cool;
1000 }
1001 
1002 /* Cooling rate of para-H2 due to collisions with para-H2 at 100 K <= T <= 6000 K
1003  * >>refer h2 cool para-ortho Glover & Abel 2008, MNRAS, 388, 1627; Table 3, third column */
1004 STATIC double ga08_pH2_oH2_a100_b6000 ( double Temp )
1005 {
1006  DEBUG_ENTRY( "ga08_pH2_oH2_a100_b6000()" );
1007  double coeff_pH2_oH2[] = { -23.748534, 1.76676480, -0.58634325, 0.31074159, -0.17455629, 0.18530758 };
1008  int ncoeff = sizeof( coeff_pH2_oH2 ) / sizeof( coeff_pH2_oH2[0] );
1009  return ga08_sum( Temp, coeff_pH2_oH2, ncoeff );
1010 }
1011 
1012 /* Cooling rate of para-H2 due to collisions with ortho-H2 at any T.
1013  * Extrapolate below 100 K. */
1014 STATIC double ga08_pH2_oH2 ( double Temp )
1015 {
1016  DEBUG_ENTRY( "ga08_pH2_oH2()" );
1017  static double cool_above_6000 = -1.;
1018  if( cool_above_6000 < 0. )
1019  {
1020  cool_above_6000 = pow( 10., ga08_pH2_oH2_a100_b6000( 6000.0 ) );
1021  }
1022 
1023  double cool = 0.;
1024  if( Temp <= 6000. )
1025  cool = pow( 10., ga08_pH2_oH2_a100_b6000( Temp ) );
1026  else
1027  cool = cool_above_6000;
1028  return cool;
1029 }
1030 
1031 
1032 
1033 /* Cooling rate of ortho-H2 due to collisions with para-H2 at 100 K <= T <= 6000 K
1034  * >>refer h2 cool ortho-para Glover & Abel 2008, MNRAS, 388, 1627; Table 4, second column */
1035 STATIC double ga08_oH2_pH2_a100_b6000 ( double Temp )
1036 {
1037  DEBUG_ENTRY( "ga08_oH2_pH2_a100_b6000()" );
1038  double coeff_oH2_pH2[] = { -24.126177, 2.3258217, -1.0082491, 0.54823768, -0.33679759, 0.20771406 };
1039  int ncoeff = sizeof( coeff_oH2_pH2 ) / sizeof( coeff_oH2_pH2[0] );
1040  return ga08_sum( Temp, coeff_oH2_pH2, ncoeff );
1041 }
1042 
1043 /* Cooling rate of ortho-H2 due to collisions with para-H2 at any T.
1044  * Extrapolate below 100 K. */
1045 STATIC double ga08_oH2_pH2 ( double Temp )
1046 {
1047  DEBUG_ENTRY( "ga08_oH2_pH2()" );
1048  static double cool_above_6000 = -1.;
1049  if( cool_above_6000 < 0. )
1050  {
1051  cool_above_6000 = pow( 10., ga08_oH2_pH2_a100_b6000( 6000. ) );
1052  }
1053  double cool = 0.;
1054  if( Temp <= 6000. )
1055  cool = pow( 10., ga08_oH2_pH2_a100_b6000( Temp ) );
1056  else
1057  cool = cool_above_6000;
1058  return cool;
1059 }
1060 
1061 
1062 
1063 /* Cooling rate of ortho-H2 due to collisions with ortho-H2 at 100 K <= T <= 6000 K
1064  * >>refer h2 cool ortho-ortho Glover & Abel 2008, MNRAS, 388, 1627; Table 4, third column */
1065 STATIC double ga08_oH2_oH2_a100_b6000 ( double Temp )
1066 {
1067  DEBUG_ENTRY( "ga08_oH2_oH2_a100_b6000()" );
1068  double coeff_oH2_oH2[] = { -24.020047, 2.2687566, -1.0200304, 0.83561432, -0.40772247, 0.096025713 };
1069  int ncoeff = sizeof( coeff_oH2_oH2 ) / sizeof( coeff_oH2_oH2[0] );
1070  return ga08_sum( Temp, coeff_oH2_oH2, ncoeff );
1071 }
1072 
1073 /* Cooling rate of ortho-H2 due to collisions with ortho-H2 at any T.
1074  * Extrapolate below 100 K. */
1075 STATIC double ga08_oH2_oH2 ( double Temp )
1076 {
1077  DEBUG_ENTRY( "ga08_oH2_oH2()" );
1078  static double cool_above_6000 = -1.;
1079  if( cool_above_6000 < 0. )
1080  {
1081  cool_above_6000 = pow( 10., ga08_oH2_oH2_a100_b6000( 6000.0 ) );
1082  }
1083 
1084  double cool = 0.;
1085  if( Temp <= 6000. )
1086  cool = pow( 10., ga08_oH2_oH2_a100_b6000( Temp ) );
1087  else
1088  cool = cool_above_6000;
1089  return cool;
1090 }
1091 
1092 
1093 
1094 /* Cooling rate of para-H2 due to collisions with He at T <= 6000 K
1095  * >>refer h2 cool He-para Glover & Abel 2008, MNRAS, 388, 1627; Table 5, second column */
1096 STATIC double ga08_pH2_He_b6000 ( double Temp )
1097 {
1098  DEBUG_ENTRY( "ga08_pH2_He_b6000()" );
1099  double coeff_pH2_He[] = { -23.489029 , 1.8210825 , -0.59110559 , 0.42280623 , -0.30171138 , 0.12872839 };
1100  int ncoeff = sizeof( coeff_pH2_He ) / sizeof( coeff_pH2_He[0] );
1101  return ga08_sum( Temp, coeff_pH2_He, ncoeff );
1102 }
1103 
1104 /* Cooling rate of para-H2 due to collisions with He at any T */
1105 STATIC double ga08_pH2_He ( double Temp )
1106 {
1107  DEBUG_ENTRY( "ga08_pH2_He()" );
1108  static double cool_above_6000 = -1.;
1109  if( cool_above_6000 < 0. )
1110  {
1111  cool_above_6000 = pow( 10., ga08_pH2_He_b6000( 6000.0 ) );
1112  }
1113 
1114  double cool = 0.;
1115  if( Temp <= 6000. )
1116  cool = pow( 10., ga08_pH2_He_b6000( Temp ) );
1117  else
1118  cool = cool_above_6000;
1119  return cool;
1120 }
1121 
1122 
1123 /* Cooling rate of ortho-H2 due to collisions with He at T <= 6000 K
1124  * >>refer h2 cool He-ortho Glover & Abel 2008, MNRAS, 388, 1627; Table 5, third column */
1125 STATIC double ga08_oH2_He_b6000 ( double Temp )
1126 {
1127  DEBUG_ENTRY( "ga08_oH2_He_b6000()" );
1128  double coeff_oH2_He[] = { -23.7749 , 2.40654 , -1.23449 , 0.739874 , -0.258940 , 0.120573 };
1129  int ncoeff = sizeof( coeff_oH2_He ) / sizeof( coeff_oH2_He[0] );
1130  return ga08_sum( Temp, coeff_oH2_He, ncoeff );
1131 }
1132 
1133 
1134 /* Cooling rate of ortho-H2 due to collisions with He at any T */
1135 STATIC double ga08_oH2_He ( double Temp )
1136 {
1137  DEBUG_ENTRY( "ga08_oH2_He()" );
1138  static double cool_above_6000 = -1.;
1139  if( cool_above_6000 < 0. )
1140  {
1141  cool_above_6000 = pow( 10., ga08_oH2_He_b6000( 6000.0 ) );
1142  }
1143 
1144  double cool = 0.;
1145  if( Temp <= 6000. )
1146  cool = pow( 10., ga08_oH2_He_b6000( Temp ) );
1147  else
1148  cool = cool_above_6000;
1149  return cool;
1150 }
1151 
1152 
1153 
1154 /* Cooling rate of para-H2 due to collisions with He at 10 K <= T <= 10000 K
1155  * >>refer h2 cool p-para Glover & Abel 2008, MNRAS, 388, 1627; Table 6, second column */
1156 STATIC double ga08_pH2_p_a10_b1e4 ( double Temp )
1157 {
1158  DEBUG_ENTRY( "ga08_pH2_p_a10_b1e4()" );
1159  double coeff_pH2_p[] = { -21.757160 , 1.3998367 , -0.37209530 , 0.061554519 , -0.37238286 , 0.23314157 };
1160  int ncoeff = sizeof( coeff_pH2_p ) / sizeof( coeff_pH2_p[0] );
1161  return ga08_sum( Temp, coeff_pH2_p, ncoeff );
1162 }
1163 
1164 /* Cooling rate of para-H2 due to collisions with p at any T */
1165 STATIC double ga08_pH2_p ( double Temp )
1166 {
1167  DEBUG_ENTRY( "ga08_pH2_p()" );
1168  static double cool_above_10000 = -1.;
1169  if( cool_above_10000 < 0. )
1170  {
1171  cool_above_10000 = pow( 10., ga08_pH2_p_a10_b1e4( 1e4 ) );
1172  }
1173 
1174  double cool = 0.;
1175  if( Temp <= 1e4 )
1176  cool = pow( 10., ga08_pH2_p_a10_b1e4( Temp ) );
1177  else
1178  cool = cool_above_10000;
1179  return cool;
1180 }
1181 
1182 
1183 /* Cooling rate of ortho-H2 due to collisions with p at 10 K <= T <= 10000 K
1184  * >>refer h2 cool p-ortho Glover & Abel 2008, MNRAS, 388, 1627; Table 6, third column */
1185 STATIC double ga08_oH2_p_a10_b1e4 ( double Temp )
1186 {
1187  DEBUG_ENTRY( "ga08_oH2_p_a10_b1e4()" );
1188  double coeff_oH2_p[] = { -21.706641 , 1.3901283 , -0.34993699 , 0.075402398 , -0.23170723 , 0.068938876 };
1189  int ncoeff = sizeof( coeff_oH2_p ) / sizeof( coeff_oH2_p[0] );
1190  return ga08_sum( Temp, coeff_oH2_p, ncoeff );
1191 }
1192 
1193 /* Cooling rate of ortho-H2 due to collisions with p at any T */
1194 STATIC double ga08_oH2_p ( double Temp )
1195 {
1196  DEBUG_ENTRY( "ga08_oH2_p()" );
1197  static double cool_above_10000 = -1.;
1198  if( cool_above_10000 < 0. )
1199  {
1200  cool_above_10000 = pow( 10., ga08_oH2_p_a10_b1e4( 1e4 ) );
1201  }
1202 
1203  double cool = 0.;
1204  if( Temp <= 1e4 )
1205  cool = pow( 10., ga08_oH2_p_a10_b1e4( Temp ) );
1206  else
1207  cool = cool_above_10000;
1208  return cool;
1209 }
1210 
1211 
1212 
1213 /* Cooling rate of para-H2 due to collisions with e at 10 K <= T <= 1000 K
1214  * >>refer h2 cool e-para Glover & Abel 2008, MNRAS, 388, 1627; Table 7, second column */
1215 STATIC double ga08_pH2_e_a10_b1000 ( double Temp )
1216 {
1217  DEBUG_ENTRY( "ga08_pH2_e_a10_b1000()" );
1218  double coeff_pH2_e_cool[] = { -22.817869 , 0.95653474 , 0.79283462 , 0.56811779 , 0.27895033 , 0.056049813 };
1219  int ncoeff = sizeof( coeff_pH2_e_cool ) / sizeof( coeff_pH2_e_cool[0] );
1220  return ga08_sum( Temp, coeff_pH2_e_cool, ncoeff );
1221 }
1222 
1223 /* Cooling rate of para-H2 due to collisions with e at 1000 K < T <= 10000 K
1224  * >>refer h2 cool e-para Glover & Abel 2008, MNRAS, 388, 1627; Table 7, third column */
1225 STATIC double ga08_pH2_e_a1000_b1e4 ( double Temp )
1226 {
1227  DEBUG_ENTRY( "ga08_pH2_e_a1000_b1e4()" );
1228  double coeff_pH2_e_warm[] = { -22.817869 , 0.66916141 , 7.1191428 , -11.176835 , 7.0467275 , -1.6471816 };
1229  int ncoeff = sizeof( coeff_pH2_e_warm ) / sizeof( coeff_pH2_e_warm[0] );
1230  return ga08_sum( Temp, coeff_pH2_e_warm, ncoeff );
1231 }
1232 
1233 /* Cooling rate of para-H2 due to collisions with e at any T */
1234 STATIC double ga08_pH2_e ( double Temp )
1235 {
1236  DEBUG_ENTRY( "ga08_pH2_e()" );
1237  static double cool_above_10000 = -1.;
1238  if( cool_above_10000 < 0. )
1239  {
1240  cool_above_10000 = pow( 10., ga08_pH2_e_a1000_b1e4( 1e4 ) );
1241  }
1242 
1243  const double E20 = 509.85; // energy in K
1244  double cool = 0.;
1245 
1246  if( Temp <= 1000. )
1247  cool = pow( 10., ga08_pH2_e_a10_b1000( Temp ) );
1248  else if( Temp <= 1e4 )
1249  cool = pow( 10., ga08_pH2_e_a1000_b1e4( Temp ) );
1250  else
1251  cool = cool_above_10000;
1252 
1253  return dsexp( E20 / Temp ) * cool;
1254 }
1255 
1256 
1257 
1258 /* Cooling rate of ortho-H2 due to collisions with e at 10 K <= T <= 10000 K
1259  * >>refer h2 cool e-para Glover & Abel 2008, MNRAS, 388, 1627; Table 7, second column */
1260 STATIC double ga08_oH2_e_a10_b1e4 ( double Temp )
1261 {
1262  DEBUG_ENTRY( "ga08_oH2_e_a10_b1e4()" );
1263  double coeff_oH2_e[] = { -21.703215 , 0.76059565 , 0.50644890 , 0.050371349 , -0.10372467 , -0.035709409 };
1264  int ncoeff = sizeof( coeff_oH2_e ) / sizeof( coeff_oH2_e[0] );
1265  return ga08_sum( Temp, coeff_oH2_e, ncoeff );
1266 }
1267 
1268 /* Cooling rate of ortho-H2 due to collisions with e at any T */
1269 STATIC double ga08_oH2_e ( double Temp )
1270 {
1271  DEBUG_ENTRY( "ga08_oH2_e()" );
1272  static double cool_above_10000 = -1.;
1273  if( cool_above_10000 < 0. )
1274  {
1275  cool_above_10000 = pow( 10., ga08_oH2_e_a10_b1e4( 1e4 ) );
1276  }
1277 
1278  double cool = 0.;
1279  if( Temp <= 1e4 )
1280  cool = pow( 10., ga08_oH2_e_a10_b1e4( Temp ) );
1281  else
1282  cool = cool_above_10000;
1283 
1284  const double E31 = 845.;
1285  return dsexp( E31 / Temp ) * cool;
1286 }
1287 
1288 
1289 
1290 /* H2 cooling per H2 molecule according to
1291  * >>refer h2 cool Glover & Abel 2008, MNRAS, 388, 1627; Sections 2.3.1-2.3.6 */
1292 STATIC double CoolH2_GA08 ( double Temp )
1293 {
1294  DEBUG_ENTRY( "CoolH2_GA08()" );
1295 
1296  double x_para = h2.para_density / dense.gas_phase[ ipHYDROGEN ];
1297  double x_ortho = h2.ortho_density / dense.gas_phase[ ipHYDROGEN ];
1298  double f_para = h2.para_density / (h2.para_density + h2.ortho_density);
1299  double f_ortho = 1. - f_para;
1300 
1301  /* Collisions with atomic H */
1302  double cool_H = dense.xIonDense[ipHYDROGEN][0] *
1303  ( f_ortho * ga08_oH2_H( Temp )
1304  + f_para * ga08_pH2_H( Temp ) );
1305 
1306  /* Collisions with H2 */
1307  double cool_H2 = hmi.H2_total *
1308  ( x_para * x_para * ga08_pH2_pH2( Temp )
1309  + x_para * x_ortho * ga08_pH2_oH2( Temp )
1310  + x_ortho * x_para * ga08_oH2_pH2( Temp )
1311  + x_ortho * x_ortho * ga08_oH2_oH2( Temp ) );
1312 
1313  /* Collisions with He */
1314  double cool_He = dense.xIonDense[ipHELIUM][0] *
1315  ( f_para * ga08_pH2_He( Temp )
1316  + f_ortho * ga08_oH2_He( Temp ) );
1317 
1318  /* Collisions with protons */
1319  double cool_p = dense.xIonDense[ipHYDROGEN][1] *
1320  ( f_para * ga08_pH2_p( Temp )
1321  + f_ortho * ga08_oH2_p( Temp ) );
1322 
1323  /* non-equilibrium cooling */
1324  // /* >>refer h2 cool eq(35) J=1<->0 Glover & Abel 2008, MNRAS, 388, 1627; Table 6, third column */
1325  // const double E10 = 170.5;
1326  // cool_rate_per_H2 += 4.76e-24 * ( 9.*dsexp( E10 / Temp ) * x_para - x_ortho ) * dense.xIonDense[ipHYDROGEN][1];
1327 
1328  /* Collisions with electrons */
1329  double cool_e = dense.eden *
1330  ( f_para * ga08_pH2_e( Temp )
1331  + f_ortho * ga08_oH2_e( Temp ) );
1332 
1333  double cooling_low = cool_H + cool_H2 + cool_He + cool_p + cool_e;
1334  double cooling_high = h2.interpolate_LTE_Cooling( Temp );
1335 
1336  {
1337  enum { DEBUG_LOC = false };
1338  if( DEBUG_LOC && iterations.lgLastIt )
1339  {
1340  printf( "%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\t%e\n",
1341  Temp,
1343  hmi.H2_total,
1344  cool_H,
1345  cool_H2,
1346  cool_He,
1347  cool_p ,
1348  cool_e ,
1349  cooling_low,
1350  cooling_high,
1351  cooling_high / (1. + cooling_high / cooling_low),
1352  cooling_high / (1. + cooling_high / cooling_low) * hmi.H2_total
1353  );
1354  }
1355  }
1356 
1357  /* NB NB NB
1358  *
1359  * The interpolation below will fail if the high density cooling is 0,
1360  * that is, the total cooling will be 0, even if the low-density cooling
1361  * is not. This may happen if permitted by the function that computes
1362  * LTE cooling, and/or the tabulated data (see definition of
1363  * interpolate_LTE_Cooling() and associated data file). */
1364  return cooling_high / (1. + cooling_high / cooling_low);
1365 }
1366 
1368 {
1369  DEBUG_ENTRY( "CoolHyperfine()" );
1370  static double TeEvalCS = 0., TeEvalCS_21cm=0.;
1371  static double electron_rate_21cm,
1372  atomic_rate_21cm,
1373  proton_rate_21cm;
1374 
1375 
1376  for( int ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
1377  {
1378  for( int nelem=ipISO; nelem < LIMELM; nelem++ )
1379  {
1380  if ( ! dense.lgElmtOn[nelem] )
1381  continue;
1382  ionbal.ExcitationGround[nelem][nelem-ipISO] +=
1383  iso_sp[ipISO][nelem].fb[0].RateLevel2Cont;
1384  }
1385  }
1386 
1387  {
1388  enum {DEBUG_LOC=false};
1389  if ( DEBUG_LOC )
1390  {
1391  for (int nelem = 0; nelem < LIMELM; nelem++)
1392  {
1393  fprintf(ioQQQ, "ionbal.ExcitationGround: nelem = %d", nelem);
1394  for (int ion = 0; ion < nelem+1; ion++)
1395  {
1396  fprintf(ioQQQ,"\t%.6e", ionbal.ExcitationGround[nelem][ion]);
1397  }
1398  fprintf(ioQQQ, "\n");
1399  }
1400  }
1401  }
1402 
1403  /* evaluate H 21 cm spin changing collisions */
1404  if( !fp_equal(phycon.te,TeEvalCS_21cm) )
1405  {
1406  {
1407  /* this prints table of rates at points given in original data paper */
1408  enum {DEBUG_LOC=false};
1409  if( DEBUG_LOC )
1410  {
1411 # define N21CM_TE 16
1412  int n;
1413  double teval[N21CM_TE]={2.,5.,10.,20.,50.,100.,200.,500.,1000.,
1414  2000.,3000.,5000.,7000.,10000.,15000.,20000.};
1415  for( n = 0; n<N21CM_TE; ++n )
1416  {
1417  fprintf(
1418  ioQQQ,"DEBUG 21 cm deex Te=\t%.2e\tH0=\t%.2e\tp=\t%.2e\te=\t%.2e\n",
1419  teval[n],
1420  H21cm_H_atom( teval[n] ),
1421  H21cm_proton( teval[n] ),
1422  H21cm_electron( teval[n] ) );
1423  }
1425 # undef N21CM_TE
1426  }
1427  }
1428  /*only evaluate T dependent part when Te changes, but update
1429  * density part below since densities may constantly change */
1430  atomic_rate_21cm = H21cm_H_atom( phycon.te );
1431  proton_rate_21cm = H21cm_proton( phycon.te );
1432  electron_rate_21cm = H21cm_electron( phycon.te );
1433  TeEvalCS_21cm = phycon.te;
1434  }
1435  /* H 21 cm emission/population,
1436  * cs will be sum of e cs and H cs converted from rate */
1437  double cs = (electron_rate_21cm * dense.eden +
1438  atomic_rate_21cm * dense.xIonDense[ipHYDROGEN][0] +
1439  proton_rate_21cm * dense.xIonDense[ipHYDROGEN][1] ) *
1440  3./ dense.cdsqte;
1441  PutCS( cs , HFLines[0] );
1442 
1443  /* do level pops for H 21 cm which is a special case since Lya pumping in included */
1444  RT_line_one_escape( HFLines[0], true,0.f, GetDopplerWidth(dense.AtomicWeight[(*HFLines[0].Hi()).nelem()-1]) );
1445  H21_cm_pops();
1446 
1447  hyperfine.cooling_total = HFLines[0].Coll().cool();
1448  if( PRT_DERIV )
1449  fprintf(ioQQQ,"DEBUG dCdT 5 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
1450 
1451 
1452  if( !fp_equal(phycon.te,TeEvalCS) )
1453  {
1454  /* do cross-sections for all except 21 cm */
1455  for( size_t i=1; i < HFLines.size(); i++ )
1456  {
1457  cs = HyperfineCS( i );
1458  PutCS( cs , HFLines[i] );
1459  }
1460  TeEvalCS = phycon.te;
1461  }
1462 
1463 
1464  /* now do level pops for all except 21 cm */
1465  for( size_t i=1; i < HFLines.size(); i++ )
1466  {
1467  /* remember current gas-phase abundance of this isotope */
1468  double save = dense.xIonDense[(*HFLines[i].Hi()).nelem()-1][(*HFLines[i].Hi()).IonStg()-1];
1469 
1470  /* bail if no abundance */
1471  if( save<=0. )
1472  continue;
1473 
1474  RT_line_one_escape( HFLines[i], true,0.f, GetDopplerWidth(dense.AtomicWeight[(*HFLines[i].Hi()).nelem()-1]) );
1475 
1476  /* set gas-phase abundance to total times isotope ratio */
1477  dense.xIonDense[(*HFLines[i].Hi()).nelem()-1][(*HFLines[i].Hi()).IonStg()-1] *=
1479 
1480  /* use the collision strength generated above and find pops and cooling */
1481  atom_level2( HFLines[i], true );
1482 
1483  /* put the correct gas-phase abundance back in the array */
1484  dense.xIonDense[(*HFLines[i].Hi()).nelem()-1][(*HFLines[i].Hi()).IonStg()-1] = save;
1485 
1486  /* find total cooling due to hyperfine structure lines */
1487  hyperfine.cooling_total += HFLines[i].Coll().cool();
1488  }
1489 
1490  return;
1491 }
1492 
1493 static const double EPS = 0.01;
1494 
1495 /*fndneg search cooling array to find negative values */
1496 STATIC void fndneg(void)
1497 {
1498  long int i;
1499  double trig;
1500 
1501  DEBUG_ENTRY( "fndneg()" );
1502 
1503  trig = fabs(thermal.htot*EPS);
1504  for( i=0; i < thermal.ncltot; i++ )
1505  {
1506  if( thermal.cooling[i] < 0. && fabs(thermal.cooling[i]) > trig )
1507  {
1508  fprintf( ioQQQ, " negative line=%s %.2f fraction of heating=%.3f\n",
1510  thermal.htot );
1511  }
1512 
1513  if( thermal.heatnt[i] > trig )
1514  {
1515  fprintf( ioQQQ, " heating line=%s %.2f fraction of heating=%.3f\n",
1517  thermal.htot );
1518  }
1519  }
1520  return;
1521 }
1522 
1523 /*fndstr search cooling stack to find strongest values */
1524 STATIC void fndstr(double tot,
1525  double dc)
1526 {
1527  char chStrngLab[NCOLNT_LAB_LEN+1];
1528  long int i;
1529  realnum wl;
1530  double str,
1531  strong;
1532 
1533  DEBUG_ENTRY( "fndstr()" );
1534 
1535  strong = 0.;
1536  wl = -FLT_MAX;
1537  for( i=0; i < thermal.ncltot; i++ )
1538  {
1539  if( fabs(thermal.cooling[i]) > strong )
1540  {
1541  /* this is the wavelength of the coolant, 0 for a continuum*/
1542  wl = thermal.collam[i];
1543  /* make sure labels are all valid*/
1544  /*>>chng 06 jun 06, bug fix, assert length was ==4, should be <=NCOLNT_LAB_LEN */
1545  ASSERT( strlen( thermal.chClntLab[i] ) <= NCOLNT_LAB_LEN );
1546  strcpy( chStrngLab, thermal.chClntLab[i] );
1547  strong = fabs(thermal.cooling[i]);
1548  }
1549  }
1550 
1551  str = strong;
1552 
1553  fprintf( ioQQQ,
1554  " fndstr cool: TE=%10.4e Ne=%10.4e C=%10.3e dC/dT=%10.2e ABS(%s %.1f)=%.2e nz=%ld\n",
1555  phycon.te, dense.eden, tot, dc, chStrngLab
1556  , wl, str, nzone );
1557 
1558  /* option for extensive printout of lines */
1559  if( trace.lgCoolTr )
1560  {
1561  realnum ratio;
1562 
1563  /* flag all significant coolants, first zero out the array */
1564  coolpr(ioQQQ,thermal.chClntLab[0],1,0.,"ZERO");
1565 
1566  /* push all coolants onto the stack */
1567  for( i=0; i < thermal.ncltot; i++ )
1568  {
1569  /* usually positive, although can be neg for coolants that heats,
1570  * only do positive here */
1571  ratio = (realnum)(thermal.cooling[i]/thermal.ctot);
1572  if( ratio >= EPS )
1573  {
1574  /*>>chng 99 jan 27, only cal when ratio is significant */
1575  coolpr(ioQQQ,thermal.chClntLab[i],thermal.collam[i], ratio,"DOIT");
1576  }
1577  }
1578 
1579  /* complete the printout for positive coolants */
1580  coolpr(ioQQQ,"DONE",1,0.,"DONE");
1581 
1582  /* now do cooling that was counted as a heat source if significant */
1583  if( thermal.heating(0,22)/thermal.ctot > 0.05 )
1584  {
1585  fprintf( ioQQQ,
1586  " All coolant heat greater than%6.2f%% of the total will be printed.\n",
1587  EPS*100. );
1588 
1589  coolpr(ioQQQ,"ZERO",1,0.,"ZERO");
1590  for( i=0; i < thermal.ncltot; i++ )
1591  {
1592  ratio = (realnum)(thermal.heatnt[i]/thermal.ctot);
1593  if( fabs(ratio) >=EPS )
1594  {
1596  ratio,"DOIT");
1597  }
1598  }
1599  coolpr(ioQQQ,"DONE",1,0.,"DONE");
1600  }
1601  }
1602  return;
1603 }
1604 
1605 
1606 // Nozawa et al. (2009) Table 1: Coefficients a^I_ik (taken from Itoh et al. 2002a).
1607 static const double aI[11][11] = { // a^I_ik = aI[i][k]
1608  { 3.15847e+00,-2.52430e+00, 4.04877e-01, 6.13466e-01, 6.28867e-01, 3.29441e-01,
1609  -1.71486e-01,-3.68685e-01,-7.59200e-02, 1.60187e-01, 8.37729e-02 },
1610  { 2.46819e-02, 1.03924e-01, 1.98935e-01, 2.18843e-01, 1.20482e-01,-4.82390e-02,
1611  -1.20811e-01,-4.46133e-04, 8.88749e-02, 2.50320e-02,-1.28900e-02 },
1612  { -2.11118e-02,-8.53821e-02,-1.52444e-01,-1.45660e-01,-4.63705e-02, 8.16592e-02,
1613  9.87296e-02,-3.24743e-02,-8.82637e-02,-7.52221e-03, 1.99419e-02 },
1614  { 1.24009e-02, 4.73623e-02, 7.51656e-02, 5.07201e-02,-2.25247e-02,-8.17151e-02,
1615  -4.59297e-02, 5.05096e-02, 5.58818e-02,-9.11885e-03,-1.71348e-02 },
1616  { -5.41633e-03,-1.91406e-02,-2.58034e-02,-2.23048e-03, 5.07325e-02, 5.94414e-02,
1617  -2.11247e-02,-5.05387e-02, 9.20453e-03, 1.67321e-02,-3.47663e-03 },
1618  { 1.70070e-03, 5.39773e-03, 4.13361e-03,-1.14273e-02,-3.23280e-02,-2.19399e-02,
1619  1.76310e-02, 2.23352e-02,-4.59817e-03,-8.24286e-03,-3.90032e-04 },
1620  { -3.05111e-04,-7.26681e-04, 4.67015e-03, 1.24789e-02,-1.16976e-02,-1.13488e-02,
1621  6.31446e-02, 1.33830e-02,-8.54735e-02,-6.47349e-03, 3.72266e-02 },
1622  { -1.21721e-04,-7.47266e-04,-2.20675e-03,-2.74351e-03,-1.00402e-03,-2.38863e-03,
1623  -2.28987e-03, 7.79323e-03, 7.98332e-03,-3.80435e-03,-4.25035e-03 },
1624  { 1.77611e-04, 8.73517e-04,-2.67582e-03,-4.57871e-03, 2.96622e-02, 1.89850e-02,
1625  -8.84093e-02,-2.93629e-02, 1.02966e-01, 1.38957e-02,-4.22093e-02 },
1626  { -2.05480e-05,-6.92284e-05, 2.95254e-05,-1.70374e-04,-5.43191e-04, 2.50978e-03,
1627  4.45570e-03,-2.80083e-03,-5.68093e-03, 1.10618e-03, 2.33625e-03 },
1628  { -3.58754e-05,-1.80305e-04, 1.40751e-03, 2.06757e-03,-1.23098e-02,-8.81767e-03,
1629  3.46210e-02, 1.23727e-02,-4.04801e-02,-5.68689e-03, 1.66733e-02 }
1630 };
1631 
1632 // Nozawa et al. (2009) Table 2: Coefficients b^I_i (taken from Itoh et al. 2002a).
1633 static const double bI[11] = { // b^I_i = bI[i]
1634  2.21564e+00, 1.83879e-01, -1.33575e-01, 5.89871e-02, -1.45904e-02, -7.10244e-04,
1635  2.80940e-03, -1.70485e-03, 5.26075e-04, 9.94159e-05, -1.06851e-04
1636 };
1637 
1638 // Nozawa et al. (2009) Table 3: Coefficients a^II_ij and b^II_ij.
1639 static const double aII[9][3] = { // a^II_ij = aII[j][i]
1640  { -3.7369800e+01,-9.3647000e+00, 9.2170000e-01 },
1641  { 3.8036590e+02, 9.5918600e+01,-1.3498800e+01 },
1642  { -1.4898014e+03,-3.9701720e+02, 7.6453900e+01 },
1643  { 2.8614150e+03, 8.4293760e+02,-2.1783010e+02 },
1644  { -2.3263704e+03,-9.0730760e+02, 3.2097530e+02 },
1645  { -6.9161180e+02, 3.0688020e+02,-1.8806670e+02 },
1646  { 2.8537893e+03, 2.9129830e+02,-8.2416100e+01 },
1647  { -2.0407952e+03,-2.9902530e+02, 1.6371910e+02 },
1648  { 4.9259810e+02, 7.6346100e+01,-6.0024800e+01 }
1649 };
1650 
1651 static const double bII[9][2] = { // b^II_ij = bII[j][i]
1652  { -1.1628100e+01,-8.6991000e+00 },
1653  { 1.2560660e+02, 6.3383000e+01 },
1654  { -5.3274890e+02,-1.2889390e+02 },
1655  { 1.1423873e+03,-1.3503120e+02 },
1656  { -1.1568545e+03, 9.7758380e+02 },
1657  { 7.5010200e+01,-1.6499529e+03 },
1658  { 9.9681140e+02, 1.2586812e+03 },
1659  { -8.8818950e+02,-4.0474610e+02 },
1660  { 2.5013860e+02, 2.7335400e+01 }
1661 };
1662 
1663 // Nozawa et al. (2009) Table 4: Coefficients c^II_ij.
1664 static const double cII[7][5] = { // c^II_ij = cII[j][i-2]
1665  { -5.7752000e+00, 3.0558600e+01,-5.4327200e+01, 3.6262500e+01,-8.4082000e+00 },
1666  { 4.6209700e+01,-2.4821770e+02, 4.5096760e+02,-3.1009720e+02, 7.4792500e+01 },
1667  { -1.6072800e+02, 8.7419640e+02,-1.6165987e+03, 1.1380531e+03,-2.8295400e+02 },
1668  { 3.0500700e+02,-1.6769028e+03, 3.1481061e+03,-2.2608347e+03, 5.7639300e+02 },
1669  { -3.2954200e+02, 1.8288677e+03,-3.4783930e+03, 2.5419361e+03,-6.6193900e+02 },
1670  { 1.9107700e+02,-1.0689366e+03, 2.0556693e+03,-1.5252058e+03, 4.0429300e+02 },
1671  { -4.6271800e+01, 2.6056560e+02,-5.0567890e+02, 3.8008520e+02,-1.0223300e+02 }
1672 };
1673 
1674 // Gf1[i][j-2] = Gamma(i + j/8. + 1.)
1675 static const double Gf1[3][5] = {
1676  { 9.0640248e-01, 8.8891357e-01, 8.8622693e-01, 8.9657428e-01, 9.1906253e-01 },
1677  { 1.1330031e+00, 1.2222562e+00, 1.3293404e+00, 1.4569332e+00, 1.6083594e+00 },
1678  { 2.5492570e+00, 2.9028584e+00, 3.3233510e+00, 3.8244497e+00, 4.4229884e+00 }
1679 };
1680 
1681 // Gf1[i][j-2] = Gamma(i + j/8. + 1.)/(i + j/8. + 1.)
1682 static const double Gf2[2][5] = {
1683  { 7.2512198e-01, 6.4648260e-01, 5.9081795e-01, 5.5173802e-01, 5.2517859e-01 },
1684  { 5.0355693e-01, 5.1463417e-01, 5.3173616e-01, 5.5502217e-01, 5.8485797e-01 }
1685 };
1686 
1687 // Nozawa et al. (2009) Table 5: Coefficients a^III_ij and b^III_ij.
1688 static const double aIII[9][3] = { // a^III_ij = aIII[j][i]
1689  { 5.2163300e+01, 4.9713900e+01, 6.4751200e+01 },
1690  { -2.5703130e+02,-1.8977460e+02,-2.1389560e+02 },
1691  { 4.4681610e+02, 2.7102980e+02, 1.7414320e+02 },
1692  { -2.9305850e+02,-2.6978070e+02, 1.3650880e+02 },
1693  { 0.0000000e+00, 4.2048120e+02,-2.7148990e+02 },
1694  { 7.7047400e+01,-5.7662470e+02, 8.9321000e+01 },
1695  { -2.3871800e+01, 4.3277900e+02, 5.8258400e+01 },
1696  { 0.0000000e+00,-1.6053650e+02,-4.6080700e+01 },
1697  { 1.9970000e-01, 2.3392500e+01, 8.7301000e+00 }
1698 };
1699 
1700 static const double bIII[9][2] = { // b^III_ij = bIII[j][i]
1701  { -8.5862000e+00, 3.7643220e+02 },
1702  { 3.4134800e+01,-1.2233635e+03 },
1703  { -1.1632870e+02, 6.2867870e+02 },
1704  { 2.9654510e+02, 2.2373946e+03 },
1705  { -3.9342070e+02,-3.8288387e+03 },
1706  { 2.3754970e+02, 2.1217933e+03 },
1707  { -3.0600000e+01,-5.5166700e+01 },
1708  { -2.7617000e+01,-3.4943210e+02 },
1709  { 8.8453000e+00, 9.2205900e+01 }
1710 };
1711 
1712 STATIC double eeBremsCooling(double Te, // electron temperature in K
1713  double eden) // electron density in cm^-3
1714 {
1715  DEBUG_ENTRY( "eeBremsCooling()" );
1716 
1717  // Evaluate e-e brems cooling; the fitting functions are described in
1718  // >>refer ee brems Nozawa S., Takahashi K., Kohyama Y., Itoh N., 2009, A&A 499, 661
1719  // >>refer ee brems Itoh N., Kawana Y., Nozawa S., 2002, Il Nuovo Cimento, 117B, 359
1720 
1721  double mec2 = ELECTRON_MASS*pow2(SPEEDLIGHT);
1722  double tau = BOLTZMANN*Te/mec2;
1723  // common factor, see Eq. 11 of Nozawa et al. (2009)
1724  double fac = mec2*pow2(eden)*SIGMA_THOMSON*SPEEDLIGHT*FINE_STRUCTURE*powpq(tau,3,2);
1725  double G = 0.;
1726  if( Te <= 5.754399e5 )
1727  {
1728  // powerlaw extrapolation to get reasonable behavior at low Te
1729  G = 1.680322*pow(Te/5.754399e5,0.116);
1730  }
1731  else if( Te <= 1.160451e7 )
1732  {
1733  // Eq. 28 of Nozawa et al. (2009)
1734  double Theta = (log10(tau) + 2.65)/1.35;
1735  // region I (50 eV < kTe <= 1 keV), Eq. 31 of Nozawa et al. (2009)
1736  double Thp = 1.;
1737  for( int i=0; i < 11; ++i )
1738  {
1739  G += bI[i]*Thp;
1740  Thp *= Theta;
1741  }
1742  G *= sqrt(8./(3.*PI));
1743  }
1744  else if( Te <= 3.481352e9 )
1745  {
1746  double A[3]={0.,0.,0.}, B[2]={0.,0.}, C[5]={0.,0.,0.,0.,0.};
1747  // region II (1 keV < kTe <= 300 keV), Eq. 39 of Nozawa et al. (2009)
1748  double tau8 = powpq(tau,1,8); // tau^(1/8)
1749  double t8p = 1.;
1750  for( int j=0; j < 9; ++j )
1751  {
1752  A[0] += aII[j][0]*t8p;
1753  A[1] += aII[j][1]*t8p;
1754  A[2] += aII[j][2]*t8p;
1755  B[0] += bII[j][0]*t8p;
1756  B[1] += bII[j][1]*t8p;
1757  t8p *= tau8;
1758  }
1759  double tau6 = powpq(tau,1,6); // tau^(1/6)
1760  double t6p = 1.;
1761  for( int j=0; j < 7; ++j )
1762  {
1763  C[0] += cII[j][0]*t6p;
1764  C[1] += cII[j][1]*t6p;
1765  C[2] += cII[j][2]*t6p;
1766  C[3] += cII[j][3]*t6p;
1767  C[4] += cII[j][4]*t6p;
1768  t6p *= tau6;
1769  }
1770  G = A[0] + A[1] + 2.*A[2] + B[0] + 0.5*B[1];
1771  for( int i=0; i < 3; ++i )
1772  for( int j=2; j < 7; ++j )
1773  G += Gf1[i][j-2]*A[i]*C[j-2];
1774  for( int i=0; i < 2; ++i )
1775  for( int j=2; j < 7; ++j )
1776  G += Gf2[i][j-2]*B[i]*C[j-2];
1777  }
1778  else
1779  {
1780  double A[3]={0.,0.,0.}, B[2]={0.,0.};
1781  // region III (300 keV < Te <= 7 Mev), Eq. 45 of Nozawa et al. (2009)
1782  double tau8 = powpq(tau,1,8); // tau^(1/8)
1783  double t8p = 1.;
1784  for( int j=0; j < 9; ++j )
1785  {
1786  A[0] += aIII[j][0]*t8p;
1787  A[1] += aIII[j][1]*t8p;
1788  A[2] += aIII[j][2]*t8p;
1789  B[0] += bIII[j][0]*t8p;
1790  B[1] += bIII[j][1]*t8p;
1791  t8p *= tau8;
1792  }
1793  G = A[0] + A[1] + 2.*A[2] + B[0] + 0.5*B[1];
1794  }
1795  return fac*G;
1796 }
1797 
1798 void eeBremsSpectrum(double Te, // electron temperature in K
1799  realnum jnu[], // 4pi*jnu in photons/cm^3/s/cell
1800  double knu[]) // opacity in cm^-1
1801 {
1802  DEBUG_ENTRY( "eeBremsSpectrum()" );
1803 
1804  // Calculate e-e brems spectrum; the fitting functions are described in
1805  // >>refer ee brems Nozawa S., Takahashi K., Kohyama Y., Itoh N., 2009, A&A 499, 661
1806  // >>refer ee brems Itoh N., Kawana Y., Nozawa S., 2002, Il Nuovo Cimento, 117B, 359
1807 
1808  double mec2 = ELECTRON_MASS*pow2(SPEEDLIGHT);
1809  // if Te is below 5.754399e5, use the shape for 5.754399e5 K
1810  double tau = BOLTZMANN*max(Te,5.754399e5)/mec2;
1811  // common factor, see Eq. 6 of Nozawa et al. (2009)
1812  // multiplication by ne^2 is done outside this routine to avoid needless reevaluations
1813  double fac = SIGMA_THOMSON*SPEEDLIGHT*FINE_STRUCTURE/sqrt(tau)*FR1RYD*HPLANCK/mec2;
1814  // adjust normalization if we are extrapolating below 50 eV
1815  // the extrapolation for the spectrum used below gives a Te^2 dependence,
1816  // eeBremsCooling() requires it to be Te^1.616, so here we correct the difference
1817  if( Te < 5.754399e5 )
1818  fac *= pow(Te/5.754399e5,-0.384);
1819  // fits are valid between 1e-4 <= hnu/kTe <= 10; convert these linits to Ryd here
1820  double Elo = max( 1.e-4*Te/TE1RYD, rfield.anu(0) );
1821  long klo = fp_equal(Elo,rfield.anu(0)) ? 0 : rfield.ithreshC(Elo);
1822  double Ehi = min( 10.*Te/TE1RYD, rfield.anu(rfield.nflux) );
1823  long khi = fp_equal(Ehi,rfield.anu(rfield.nflux)) ? rfield.nflux : rfield.ithreshC(Ehi);
1824  // normalization constant for Bnu
1825  double fac2 = PI8*pow3(FR1RYD)/pow2(SPEEDLIGHT);
1826  avx_ptr<double> arg(rfield.nflux), boltz(rfield.nflux);
1827  double* bfac;
1828  if( fp_equal( Te, phycon.te ) )
1829  {
1830  bfac = get_ptr(rfield.ContBoltz);
1831  }
1832  else
1833  {
1834  for( long k=klo; k < khi; ++k )
1835  arg[k] = rfield.anu(k)*TE1RYD/Te;
1836  vexp( arg.ptr0(), boltz.ptr0(), klo, khi );
1837  bfac = boltz.ptr0();
1838  }
1839 
1840  if( Te <= 1.160451e7 )
1841  {
1842  // region I (50 eV < kTe <= 1 keV), Eq. 26 of Nozawa et al. (2009)
1843  //
1844  // Below 50 eV (is 5.754399e5 K) we extrapolate the expressions assuming
1845  // that P_ee(x) is independent of Te. Test between Te = 1e6 and 1e7 K show
1846  // that this is reasonable. The normalization constant is then adjusted to
1847  // give the correct integral in accordance with what eeBremsCooling()
1848  // returns.
1849  //
1850  double Theta = (log10(tau) + 2.65)/1.35; // Eq. 28 of Nozawa et al. (2009)
1851  double Thp[11];
1852  Thp[0] = 1.;
1853  for( int i=1; i < 11; ++i )
1854  Thp[i] = Thp[i-1]*Theta;
1855 
1856  avx_ptr<double> x(klo,khi), y(klo,khi), z(klo,khi);
1857  double hokT = TE1RYD/Te;
1858  for( long k=klo; k < khi; ++k )
1859  x[k] = rfield.anu(k)*hokT;
1860  vlog10( x.ptr0(), y.ptr0(), klo, khi );
1861  vexpm1( x.ptr0(), z.ptr0(), klo, khi );
1862  for( long k=klo; k < khi; ++k )
1863  {
1864  double xx = (y[k] + 1.5)/2.5; // Eq. 29 of Nozawa et al. (2009)
1865  double xxp[11];
1866  xxp[0] = 1.;
1867  for( int i=1; i < 11; ++i )
1868  xxp[i] = xxp[i-1]*xx;
1869 
1870  double G = 0.;
1871  for( int i=0; i < 11; ++i )
1872  for( int j=0; j < 11; ++j )
1873  G += aI[i][j]*Thp[i]*xxp[j];
1874  G *= sqrt(8./(3.*PI));
1875  double j_nu = fac*bfac[k]/x[k]*G*rfield.widflx(k);
1876  // calculate opacity from Kirchhoff's law, use 4pi*Bnu in photons/cm^2/s/cell
1877  double Bnu = fac2*rfield.anu2(k)*rfield.widflx(k)/z[k];
1878  jnu[k] = realnum(j_nu);
1879  knu[k] = j_nu/Bnu;
1880  }
1881  }
1882  else if( Te <= 3.481352e9 )
1883  {
1884  // region II (1 keV < kTe <= 300 keV), Eq. 32 of Nozawa et al. (2009)
1885  double A[3]={0.,0.,0.}, B[2]={0.,0.}, C[5]={0.,0.,0.,0.,0.};
1886  double tau8 = powpq(tau,1,8); // tau^(1/8)
1887  double t8p = 1.;
1888  for( int j=0; j < 9; ++j )
1889  {
1890  A[0] += aII[j][0]*t8p;
1891  A[1] += aII[j][1]*t8p;
1892  A[2] += aII[j][2]*t8p;
1893  B[0] += bII[j][0]*t8p;
1894  B[1] += bII[j][1]*t8p;
1895  t8p *= tau8;
1896  }
1897  double tau6 = powpq(tau,1,6); // tau^(1/6)
1898  double t6p = 1.;
1899  for( int j=0; j < 7; ++j )
1900  {
1901  C[0] += cII[j][0]*t6p;
1902  C[1] += cII[j][1]*t6p;
1903  C[2] += cII[j][2]*t6p;
1904  C[3] += cII[j][3]*t6p;
1905  C[4] += cII[j][4]*t6p;
1906  t6p *= tau6;
1907  }
1908 
1909  for( long k=klo; k < khi; ++k )
1910  {
1911  double x = rfield.anu(k)*TE1RYD/Te;
1912  // use the identity -Ei(-x) = E1(x)
1913  double G = (A[2]*x + A[1])*x + A[0] + e1_scaled(x)*(B[1]*x + B[0]);
1914  double F = 1.;
1915  double x8 = powpq(x,1,8); // x^1/8
1916  double x8p = x8*x8;
1917  for( int i=2; i < 7; ++i )
1918  {
1919  F += C[i-2]*x8p;
1920  x8p *= x8;
1921  }
1922  double j_nu = fac*bfac[k]/x*G*F*rfield.widflx(k);
1923  // calculate opacity from Kirchhoff's law, use 4pi*Bnu in photons/cm^2/s/cell
1924  double Bnu = fac2*rfield.anu2(k)*rfield.widflx(k)/expm1(x);
1925  jnu[k] = realnum(j_nu);
1926  knu[k] = j_nu/Bnu;
1927  }
1928  }
1929  else
1930  {
1931  // region III (300 keV < Te <= 7 Mev), Eq. 40 of Nozawa et al. (2009)
1932  double A[3]={0.,0.,0.}, B[2]={0.,0.};
1933  double tau8 = powpq(tau,1,8); // tau^(1/8)
1934  double t8p = 1.;
1935  for( int j=0; j < 9; ++j )
1936  {
1937  A[0] += aIII[j][0]*t8p;
1938  A[1] += aIII[j][1]*t8p;
1939  A[2] += aIII[j][2]*t8p;
1940  B[0] += bIII[j][0]*t8p;
1941  B[1] += bIII[j][1]*t8p;
1942  t8p *= tau8;
1943  }
1944 
1945  for( long k=klo; k < khi; ++k )
1946  {
1947  double x = rfield.anu(k)*TE1RYD/Te;
1948  // use the identity -Ei(-x) = E1(x)
1949  double G = (A[2]*x + A[1])*x + A[0] + e1_scaled(x)*(B[1]*x + B[0]);
1950  double j_nu = fac*bfac[k]/x*G*rfield.widflx(k);
1951  // calculate opacity from Kirchhoff's law, use 4pi*Bnu in photons/cm^2/s/cell
1952  double Bnu = fac2*rfield.anu2(k)*rfield.widflx(k)/expm1(x);
1953  jnu[k] = realnum(j_nu);
1954  knu[k] = j_nu/Bnu;
1955  }
1956  }
1957 
1958  //
1959  // extrapolate the emissivity and opacity to cover the entire frequency range of Cloudy
1960  //
1961  // Itoh calculated the total cooling rate by integrating his fits from 0 to infinity.
1962  // This implies a different extrapolation than what we use below, and therefore the
1963  // integral over jnu[] is not exactly equal to what eeBremsCooling() reports above.
1964  // These two always agree within 0.03% though, with the largest discrepancy at 10 GK.
1965  //
1966  double x1l = rfield.anuln(klo);
1967  double x2l = rfield.anuln(klo+1);
1968  double y1 = jnu[klo]/rfield.widflx(klo);
1969  double y2 = jnu[klo+1]/rfield.widflx(klo+1);
1970  double z1 = knu[klo];
1971  double slope = log(y1/y2)/(x1l-x2l);
1972  // extrapolate downwards as jnu = C*nu^slope
1974  for( long k=0; k < klo; ++k )
1975  {
1976  x2l = rfield.anuln(k);
1977  arg[k] = slope*(x2l-x1l);
1978  }
1979  vexp( arg.ptr0(), efac.ptr0(), 0, klo );
1980  for( long k=0; k < klo; ++k )
1981  {
1982  jnu[k] = realnum(y1*rfield.widflx(k)*efac[k]);
1983  double af = rfield.anu(klo)/rfield.anu(k);
1984  knu[k] = z1*efac[k]*af;
1985  }
1986  x1l = rfield.anuln(khi-1);
1987  double x1 = rfield.anu(khi-1);
1988  x2l = rfield.anuln(khi-2);
1989  double x2 = rfield.anu(khi-2);
1990  y1 = jnu[khi-1]/rfield.widflx(khi-1);
1991  y2 = jnu[khi-2]/rfield.widflx(khi-2);
1992  z1 = knu[khi-1];
1993  slope = (log(y2/y1) + (x2-x1)*TE1RYD/Te)/(x2l-x1l);
1994  // extrapolate upwards as jnu = C*nu^slope*exp(-hnu/kTe)
1995  avx_ptr<double> arg2(rfield.nflux), efac2(rfield.nflux);
1996  for( long k=khi; k < rfield.nflux; ++k )
1997  {
1998  x2 = rfield.anu(k);
1999  x2l = rfield.anuln(k);
2000  arg[k] = (x1-x2)*TE1RYD/Te - slope*(x1l-x2l);
2001  arg2[k] = (slope-2.)*(x2l-x1l);
2002  }
2003  vexp( arg.ptr0(), efac.ptr0(), khi, rfield.nflux );
2004  vexp( arg2.ptr0(), efac2.ptr0(), khi, rfield.nflux );
2005  for( long k=khi; k < rfield.nflux; ++k )
2006  {
2007  jnu[k] = realnum(y1*rfield.widflx(k)*efac[k]);
2008  knu[k] = z1*efac2[k];
2009  }
2010 }
void CoolAdd(const char *chLabel, realnum lambda, double cool)
Definition: cool_etc.cpp:13
STATIC double ga08_pH2_e(double Temp)
Definition: cool_eval.cpp:1234
static const double aII[9][3]
Definition: cool_eval.cpp:1639
realnum cextpw
Definition: thermal.h:152
double HyperfineCS(size_t i)
t_mole_global mole_global
Definition: mole.cpp:7
double hmicol
Definition: hmi.h:33
bool lgBakesPAH_heat
Definition: grainvar.h:479
bool lgStancil
Definition: mole.h:337
double Radius
Definition: radius.h:31
double htot
Definition: thermal.h:169
vector< double, allocator_avx< double > > ContBoltz
Definition: rfield.h:124
t_atmdat atmdat
Definition: atmdat.cpp:6
double HMinus_photo_rate
Definition: hmi.h:53
t_thermal thermal
Definition: thermal.cpp:6
static double x2[63]
size_t size(void) const
Definition: transition.h:331
double brems_cool(long ion, double Te)
void coolpr(FILE *io, const char *chLabel, realnum lambda, double ratio, const char *chJOB)
Definition: cool_pr.cpp:10
double Cool()
Definition: dynamics.cpp:2207
STATIC double ga08_oH2_He(double Temp)
Definition: cool_eval.cpp:1135
STATIC double ga08_oH2_oH2(double Temp)
Definition: cool_eval.cpp:1075
T * ptr0()
Definition: vectorize.h:240
double exp10(double x)
Definition: cddefines.h:1368
STATIC double ga08_oH2_e_a10_b1e4(double Temp)
Definition: cool_eval.cpp:1260
STATIC double ga08_pH2_p(double Temp)
Definition: cool_eval.cpp:1165
double hmihet
Definition: hmi.h:33
long int ipEnergyBremsThin
Definition: rfield.h:226
T * get_ptr(T *v)
Definition: cddefines.h:1109
NORETURN void TotalInsanity(void)
Definition: service.cpp:971
double widflx(size_t i) const
Definition: mesh.h:156
t_opac opac
Definition: opacity.cpp:5
void vlog10(const double x[], double y[], long nlo, long nhi)
static double x1[83]
realnum ** flux
Definition: rfield.h:68
double expans
Definition: coolheavy.h:18
#define N21CM_TE
t_hyperfine hyperfine
Definition: hyperfine.cpp:5
double cooling_total
Definition: hyperfine.h:62
bool lgStoutOn
Definition: atmdat.h:402
static const double aIII[9][3]
Definition: cool_eval.cpp:1688
STATIC double ga08_pH2_H_b1000(double Temp)
Definition: cool_eval.cpp:895
const int NISO
Definition: cddefines.h:311
double H21cm_H_atom(double temp)
double GasCoolColl
Definition: grainvar.h:544
void RT_line_one_escape(const TransitionProxy &t, bool lgShield_this_zone, realnum pestrk, realnum DopplerWidth)
double HeatH2Dish_BD96
Definition: hmi.h:140
long int IonHigh[LIMELM+1]
Definition: dense.h:130
STATIC double ga08_pH2_oH2(double Temp)
Definition: cool_eval.cpp:1014
double HeatDexc_deriv
Definition: h2_priv.h:299
double den_Hepp
Definition: atmdat_gaunt.h:40
double cdsqte
Definition: dense.h:246
double den_Hep
Definition: atmdat_gaunt.h:39
t_magnetic magnetic
Definition: magnetic.cpp:17
long int ipMaxBolt
Definition: rfield.h:230
bool lgTimeDependentStatic
Definition: dynamics.h:102
t_conv conv
Definition: conv.cpp:5
TransitionList HFLines("HFLines",&AnonStates)
bool lgEvaluated
Definition: h2_priv.h:317
double brems_cool_he
Definition: coolheavy.h:36
t_phycon phycon
Definition: phycon.cpp:6
STATIC double ga08_oH2_p(double Temp)
Definition: cool_eval.cpp:1194
void CoolEvaluate(double *tot)
Definition: cool_eval.cpp:58
double cooling[NCOLNT]
Definition: thermal.h:104
T pow3(T a)
Definition: cddefines.h:988
t_CoolHeavy CoolHeavy
Definition: coolheavy.cpp:5
bool lgTemperatureConstant
Definition: thermal.h:44
double pressure
Definition: magnetic.h:41
static const bool PRT_DERIV
Definition: cool_eval.cpp:56
double char_tran_cool
Definition: thermal.h:166
FILE * ioQQQ
Definition: cddefines.cpp:7
STATIC double ga08_pH2_e_a1000_b1e4(double Temp)
Definition: cool_eval.cpp:1225
molezone * findspecieslocal(const char buf[])
STATIC double ga08_oH2_pH2_a100_b6000(double Temp)
Definition: cool_eval.cpp:1035
static const double aI[11][11]
Definition: cool_eval.cpp:1607
bool lgCExtraOn
Definition: thermal.h:151
double h2plus_heat
Definition: hmi.h:48
long int nzone
Definition: cddefines.cpp:14
double HeatH2Dexc_used
Definition: hmi.h:140
double HeatH2Dish_TH85
Definition: hmi.h:140
bool lgTalk
Definition: called.h:12
t_dynamics dynamics
Definition: dynamics.cpp:42
double eebrm
Definition: coolheavy.h:18
static const double EPS
Definition: cool_eval.cpp:1493
vector< freeBound > fb
Definition: iso.h:481
#define MIN2(a, b)
Definition: cddefines.h:803
double anu(size_t i) const
Definition: mesh.h:120
long int nSpecies
Definition: taulines.cpp:22
double dsexp(double x)
Definition: service.cpp:1038
void TempChange(double TempNew, bool lgForceUpdate)
Definition: temp_change.cpp:32
realnum deriv_HeatH2Dexc_ELWERT
Definition: hmi.h:156
double tccool
Definition: coolheavy.h:18
STATIC double ga08_oH2_p_a10_b1e4(double Temp)
Definition: cool_eval.cpp:1185
STATIC double ga08_oH2_H(double Temp)
Definition: cool_eval.cpp:850
static const double Gf2[2][5]
Definition: cool_eval.cpp:1682
bool lgNumDeriv
Definition: numderiv.h:18
double brems_cool_net
Definition: coolheavy.h:36
t_dense dense
Definition: global.cpp:15
double anuln(size_t i) const
Definition: mesh.h:132
STATIC double ga08_oH2_H_b100(double Temp)
Definition: cool_eval.cpp:770
double ** ExcitationGround
Definition: ionbal.h:124
static const double cII[7][5]
Definition: cool_eval.cpp:1664
static t_gaunt & Inst()
Definition: cddefines.h:209
STATIC void CoolHyperfine(void)
Definition: cool_eval.cpp:1367
bool lgCNegChk
Definition: thermal.h:122
t_iso_sp iso_sp[NISO][LIMELM]
Definition: iso.cpp:11
void vexp(const double x[], double y[], long nlo, long nhi)
long int nflux_with_check
Definition: rfield.h:49
double Heat()
Definition: dynamics.cpp:2193
bool lgIonStoutOn[LIMELM][LIMELM+1]
Definition: dense.h:143
STATIC double ga08_sum(double Temp, double *coeff, int ncoeff)
Definition: cool_eval.cpp:753
Wind wind
Definition: wind.cpp:5
realnum deriv_HeatH2Dexc_TH85
Definition: hmi.h:156
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
double H21cm_proton(double temp)
void CoolSum(double *total)
Definition: cool_etc.cpp:87
long int iteration
Definition: cddefines.cpp:16
static const double bIII[9][2]
Definition: cool_eval.cpp:1700
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
double ortho_density
Definition: h2_priv.h:326
double H21cm_electron(double temp)
STATIC double ga08_oH2_He_b6000(double Temp)
Definition: cool_eval.cpp:1125
STATIC double ga08_oH2_pH2(double Temp)
Definition: cool_eval.cpp:1045
STATIC double ga08_oH2_oH2_a100_b6000(double Temp)
Definition: cool_eval.cpp:1065
bool lgBallistic(void) const
Definition: wind.h:31
realnum H2PlsCool
Definition: coolheavy.h:47
bool lgIonChiantiOn[LIMELM][LIMELM+1]
Definition: dense.h:140
STATIC double eeBremsCooling(double Te, double eden)
Definition: cool_eval.cpp:1712
t_ionbal ionbal
Definition: ionbal.cpp:8
double chem_heat(void) const
void CoolZero(void)
Definition: cool_etc.cpp:50
realnum CoolExtra
Definition: thermal.h:152
long int IonLow[LIMELM+1]
Definition: dense.h:129
void atom_level2(const TransitionProxy &t, const bool lgHFS)
Definition: atom_level2.cpp:30
STATIC double ga08_pH2_H_b100(double Temp)
Definition: cool_eval.cpp:886
void CoolDima(void)
Definition: cool_dima.cpp:20
double para_density
Definition: h2_priv.h:326
STATIC double ga08_pH2_pH2(double Temp)
Definition: cool_eval.cpp:985
void iso_cool(long ipISO, long nelem)
void PutCS(double cs, const TransitionProxy &t)
Definition: transition.cpp:299
STATIC void fndneg(void)
Definition: cool_eval.cpp:1496
double HeatH2Dish_used
Definition: hmi.h:140
long int nPres2Ioniz
Definition: conv.h:145
double GasHeatTherm
Definition: grainvar.h:544
char chClntLab[NCOLNT][NCOLNT_LAB_LEN+1]
Definition: thermal.h:108
bool lgChiantiOn
Definition: atmdat.h:374
#define STATIC
Definition: cddefines.h:118
bool lgTrace
Definition: trace.h:12
realnum deriv_HeatH2Dexc_BD96
Definition: hmi.h:156
bool lgEnabled
Definition: h2_priv.h:352
double brems_cool_h
Definition: coolheavy.h:36
double teinv
Definition: phycon.h:33
t_mole_local mole
Definition: mole.cpp:8
double cmcool
Definition: rfield.h:275
STATIC double ga08_oH2_H_stitch_100(double Temp)
Definition: cool_eval.cpp:803
double EdenTrue
Definition: dense.h:232
t_pressure pressure
Definition: pressure.cpp:9
t_rfield rfield
Definition: rfield.cpp:9
void dBase_solve(void)
Definition: species2.cpp:271
double cextxx
Definition: coolheavy.h:18
static const double Gf1[3][5]
Definition: cool_eval.cpp:1675
double timestep
Definition: dynamics.h:187
float realnum
Definition: cddefines.h:124
double dHeatdT
Definition: dynamics.h:69
#define EXIT_FAILURE
Definition: cddefines.h:168
realnum AccelTotalOutward
Definition: wind.h:52
double HD_total
Definition: hmi.h:27
STATIC double ga08_oH2_H_b1000(double Temp)
Definition: cool_eval.cpp:779
vector< diatomics * > diatoms
Definition: h2.cpp:8
long max(int a, long b)
Definition: cddefines.h:817
#define cdEXIT(FAIL)
Definition: cddefines.h:482
double brems_heat_total
Definition: coolheavy.h:36
long min(int a, long b)
Definition: cddefines.h:762
void eeBremsSpectrum(double Te, realnum jnu[], double knu[])
Definition: cool_eval.cpp:1798
double anu2(size_t i) const
Definition: mesh.h:124
bool lgGrainPhysicsOn
Definition: grainvar.h:479
diatomics h2("h2", 4100.,&hmi.H2_total, Yan_H2_CS)
realnum GetDopplerWidth(realnum massAMU)
double HeatH2Dexc_TH85
Definition: hmi.h:140
t_iterations iterations
Definition: iterations.cpp:6
double HeatH2Dexc_BHT90
Definition: hmi.h:140
double PresGasCurr
Definition: pressure.h:46
bool lgCoolTr
Definition: trace.h:109
double heavfb
Definition: coolheavy.h:18
double den_ion[LIMELM+1]
Definition: atmdat_gaunt.h:42
static bool lgMustPrintHeader
Definition: save_line.cpp:270
double halfte
Definition: thermal.h:142
STATIC double ga08_oH2_H_b6000(double Temp)
Definition: cool_eval.cpp:789
t_radius radius
Definition: radius.cpp:5
double heating(long nelem, long ion)
Definition: thermal.h:186
bool lgElmtOn[LIMELM]
Definition: dense.h:160
double gauntff(long Z, double Te, double anu)
double brems_cool_hminus
Definition: coolheavy.h:36
double GasHeatPhotoEl
Definition: grainvar.h:544
void setHeating(long nelem, long ion, double heating)
Definition: thermal.h:190
double HeatDexc
Definition: h2_priv.h:297
STATIC double ga08_oH2_H_stitch_1000(double Temp)
Definition: cool_eval.cpp:827
realnum gas_phase[LIMELM]
Definition: dense.h:76
double HeatH2Dish_ELWERT
Definition: hmi.h:140
double brems_cool_metals
Definition: coolheavy.h:36
realnum AtomicWeight[LIMELM]
Definition: dense.h:80
#define ASSERT(exp)
Definition: cddefines.h:613
bool lgLastIt
Definition: iterations.h:47
double HeatDiss
Definition: h2_priv.h:296
long int ncltot
Definition: thermal.h:106
realnum collam[NCOLNT]
Definition: thermal.h:103
double HeatH2Dish_BHT90
Definition: hmi.h:140
static const double bI[11]
Definition: cool_eval.cpp:1633
double den_Hm
Definition: atmdat_gaunt.h:37
bool lgDHetOn
Definition: grainvar.h:489
char chH2_small_model_type
Definition: hmi.h:179
realnum deriv_HeatH2Dexc_BHT90
Definition: hmi.h:156
const int ipH_LIKE
Definition: iso.h:64
const int LIMELM
Definition: cddefines.h:308
T pow2(T a)
Definition: cddefines.h:981
t_NumDeriv NumDeriv
Definition: numderiv.cpp:5
double dHeatdT
Definition: thermal.h:169
STATIC double ga08_pH2_He(double Temp)
Definition: cool_eval.cpp:1105
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:723
double powpq(double x, int p, int q)
Definition: service.cpp:630
realnum scalingDensity(void)
Definition: dense.cpp:409
const int ipHELIUM
Definition: cddefines.h:350
double * FreeFreeOpacity
Definition: opacity.h:126
bool lgH2_Thermal_BigH2
Definition: hmi.h:171
STATIC double ga08_pH2_H_b6000(double Temp)
Definition: cool_eval.cpp:905
double dCooldT()
Definition: dynamics.cpp:2222
double H2_total
Definition: hmi.h:25
#define NCOLNT_LAB_LEN
Definition: thermal.h:107
double elementcool[LIMELM+1]
Definition: thermal.h:111
diatomics hd("hd", 4100.,&hmi.HD_total, Yan_H2_CS)
STATIC double ga08_oH2_e(double Temp)
Definition: cool_eval.cpp:1269
vector< species > dBaseSpecies
Definition: taulines.cpp:15
double eden
Definition: dense.h:201
double e1_scaled(double x)
size_t ithreshC(double threshold) const
Definition: mesh.h:174
#define MAX2(a, b)
Definition: cddefines.h:824
double HeatH2Dexc_ELWERT
Definition: hmi.h:140
int fprintf(const Output &stream, const char *format,...)
Definition: service.cpp:1121
STATIC double ga08_pH2_H_stitch_100(double Temp)
Definition: cool_eval.cpp:916
STATIC double ga08_pH2_e_a10_b1000(double Temp)
Definition: cool_eval.cpp:1215
realnum deriv_HeatH2Dexc_used
Definition: hmi.h:156
void dBaseTrim(void)
Definition: species2.cpp:61
sys_float SDIV(sys_float x)
Definition: cddefines.h:1002
double alogte
Definition: phycon.h:92
long int n_initial_relax
Definition: dynamics.h:132
realnum Upstream_density
Definition: dynamics.h:175
STATIC double ga08_pH2_p_a10_b1e4(double Temp)
Definition: cool_eval.cpp:1156
static const double bII[9][2]
Definition: cool_eval.cpp:1651
double sqrte
Definition: phycon.h:58
STATIC void fndstr(double tot, double dc)
Definition: cool_eval.cpp:1524
#define fixit(a)
Definition: cddefines.h:417
GrainVar gv
Definition: grainvar.cpp:5
t_hmi hmi
Definition: hmi.cpp:5
double fnzone
Definition: cddefines.cpp:15
void brems_sum_ions(t_brems_den &sum) const
double interpolate_LTE_Cooling(double Temp)
void ShowMe(void)
Definition: service.cpp:205
t_save save
Definition: save.cpp:5
double te
Definition: phycon.h:21
bool lgNoMole
Definition: mole.h:325
double dCooldT
Definition: thermal.h:139
STATIC double ga08_pH2_pH2_a100_b6000(double Temp)
Definition: cool_eval.cpp:975
double h2line
Definition: coolheavy.h:18
bool lgDustOn() const
Definition: grainvar.h:475
const int ipHYDROGEN
Definition: cddefines.h:349
void vexpm1(const double x[], double y[], long nlo, long nhi)
double heatnt[NCOLNT]
Definition: thermal.h:104
bool lgFreeOn
Definition: coolheavy.h:35
long int nflux
Definition: rfield.h:46
const int ipLITHIUM
Definition: cddefines.h:351
double cyntrn
Definition: coolheavy.h:18
double HD
Definition: coolheavy.h:18
void H21_cm_pops(void)
realnum windv
Definition: wind.h:18
double HeatH2Dexc_BD96
Definition: hmi.h:140
t_called called
Definition: called.cpp:4
vector< diatomics * >::iterator diatom_iter
Definition: h2.h:13
STATIC double ga08_pH2_oH2_a100_b6000(double Temp)
Definition: cool_eval.cpp:1004
void dBaseUpdateCollCoeffs(void)
Definition: species2.cpp:102
double den_Hp
Definition: atmdat_gaunt.h:38
STATIC double ga08_pH2_He_b6000(double Temp)
Definition: cool_eval.cpp:1096
double ** RR_rate_coef_used
Definition: ionbal.h:210
double ctot
Definition: thermal.h:130
STATIC double ga08_pH2_H(double Temp)
Definition: cool_eval.cpp:939
realnum * HFLabundance
Definition: hyperfine.h:53
double colmet
Definition: coolheavy.h:18
STATIC double CoolH2_GA08(double Temp)
Definition: cool_eval.cpp:1292