cloudy  trunk
 All Data Structures 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-2008 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 "physconst.h"
6 #include "hydrogenic.h"
7 #include "taulines.h"
8 #include "wind.h"
9 #include "coolheavy.h"
10 #include "radius.h"
11 #include "conv.h"
12 #include "h2.h"
13 #include "opacity.h"
14 #include "ionbal.h"
15 #include "dense.h"
16 #include "trace.h"
17 #include "dynamics.h"
18 #include "rfield.h"
19 #include "grainvar.h"
20 #include "atoms.h"
21 #include "called.h"
22 #include "hmi.h"
23 #include "numderiv.h"
24 #include "magnetic.h"
25 #include "phycon.h"
26 #include "lines_service.h"
27 #include "hyperfine.h"
28 #include "iso.h"
29 #include "thermal.h"
30 #include "cooling.h"
31 /*fndneg search cooling array to find negative values */
32 STATIC void fndneg(void);
33 /*fndstr search cooling stack to find strongest values */
34 STATIC void fndstr(double tot,
35  double dc);
36 
37 /* set true to debug derivative of heating and cooling */
38 static const bool PRT_DERIV = false;
39 
40 void CoolEvaluate(double *tot)
41 {
42  long int ion,
43  i,
44  ipISO,
45  nelem;
46 
47  static long int nhit = 0,
48  nzSave=0;
49 
50  static double TeEvalCS = 0., TeEvalCS_21cm=0.;
51  static double TeUsedBrems=-1.f;
52  static int nzoneUsedBrems=-1;
53 
54  static double electron_rate_21cm,
55  atomic_rate_21cm,
56  proton_rate_21cm;
57 
58  double
59  cs ,
60  deriv,
61  factor,
62  qn,
63  rothi=-SMALLFLOAT,
64  rotlow=-SMALLFLOAT,
65  x;
66 
67  static double oltcool=0.,
68  oldtemp=0.;
69 
70  DEBUG_ENTRY( "CoolEvaluate()" );
71 
72  /* returns tot, the total cooling,
73  * and dc, the derivative of the cooling */
74 
75  /* routine atom_level2( t10 )
76  * routine atom_level3( abund , t10,t21,t20)
77  * tsq1 = 1. / (te**2)
78  * POPEXC( O12,g1,g2,A21,excit,abund); result already*a21
79  * POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2)
80  * AtomSeqBeryllium(cs23,cs24,cs34,tarray,a41)
81  * FIVEL( G(1-5) , ex(wn,1-5), cs12,cs13,14,15,23,24,25,34,35,45,
82  * A21,31,41,51,32,42,52,43,53,54, pop(1-5), abund) */
83 
84  if( trace.lgTrace )
85  fprintf( ioQQQ, " COOLR TE:%.4e zone %li %li Cool:%.4e Heat:%.4e eden:%.4e edenTrue:%.4e\n",
86  phycon.te,
89 
90  /* must call TempChange since ionization has changed, there are some
91  * terms that affect collision rates (H0 term in electron collision) */
92  TempChange(phycon.te , false);
93 
94  /* now zero out the cooling stack */
95  CoolZero();
96  if( PRT_DERIV )
97  fprintf(ioQQQ,"DEBUG dcdt 0 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);
98  if( gv.lgGrainPhysicsOn )
99  {
100  /* grain heating and cooling */
101  /* grain recombination cooling, evaluated elsewhere
102  * can either heat or cool the gas, do cooling here */
103  CoolAdd("dust",0,MAX2(0.,gv.GasCoolColl));
104 
105  /* grain cooling proportional to temperature ^3/2 */
106  thermal.dCooldT += MAX2(0.,gv.GasCoolColl)*3./(2.*phycon.te);
107 
108  /* these are the various heat agents from grains */
109  /* options to force gas heating or cooling by grains to zero - for tests only ! */
110  if( gv.lgDustOn && gv.lgDHetOn )
111  {
112  /* rate dust heats gas by photoelectric effect */
113  thermal.heating[0][13] = gv.GasHeatPhotoEl;
114 
115  /* if grains hotter than gas then collisions with gas act
116  * to heat the gas, add this in here
117  * a symmetric statement appears in COOLR, where cooling is added on */
118  thermal.heating[0][14] = MAX2(0.,-gv.GasCoolColl);
119 
120  /* this is gas heating due to thermionic emissions */
121  thermal.heating[0][25] = gv.GasHeatTherm;
122  }
123  else
124  {
125  thermal.heating[0][13] = 0.;
126  thermal.heating[0][14] = 0.;
127  thermal.heating[0][25] = 0.;
128  }
129  }
130  else if( gv.lgBakesPAH_heat )
131  {
132  /* >>chng 06 jul 21, option to include Bakes PAH hack with grain physics off,
133  * needed to test dynamics models */
134  thermal.heating[0][13] = gv.GasHeatPhotoEl;
135  }
136 
137  /* find net heating - cooling due to large H2 molecule */
138  H2_Cooling("CoolEvaluate");
139  if( PRT_DERIV )
140  fprintf(ioQQQ,"DEBUG dcdt 1 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);
141 
142  /* molecular molecules molecule cooling */
143  if( hmi.lgNoH2Mole )
144  {
145  /* this branch - do not include molecules */
146  hmi.hmicol = 0.;
148  /* line cooling within simple H2 molecule - zero when big used */
149  CoolHeavy.h2line = 0.;
150  /* H + H+ => H2+ cooling */
151  CoolHeavy.H2PlsCool = 0.;
152  CoolHeavy.HD = 0.;
153 
154  /* thermal.heating[0][8] is heating due to collisions within X of H2 */
155  thermal.heating[0][8] = 0.;
156  /* thermal.heating[0][15] is H minus heating*/
157  thermal.heating[0][15] = 0.;
158  /* thermal.heating[0][16] is H2+ heating */
159  thermal.heating[0][16] = 0.;
160  hmi.HeatH2Dish_used = 0.;
161  hmi.HeatH2Dexc_used = 0.;
162  }
163 
164  else
165  {
166  /* save various molecular heating/cooling agent */
167  thermal.heating[0][15] = hmi.hmihet;
168  thermal.heating[0][16] = hmi.h2plus_heat;
169  /* now get heating from H2 molecule, either simple or from big one */
171  {
172  if( hmi.lgBigH2_evaluated )
173  {
174  /* these are explicitly from big H2 molecule,
175  * first is heating due to radiative pump of excited states, followed by
176  * radiative decay into continuum of X, followed by dissociation of molecule
177  * with kinetic energy, typically 0.25 - 0.5 eV per event */
180  /*fprintf(ioQQQ,"DEBUG big %.2f\t%.5e\t%.2e\t%.2e\t%.2e\n",
181  fnzone , phycon.te, hmi.HeatH2Dexc_used,
182  hmi.H2_total , hmi.H2_star_BigH2);*/
183  /* negative sign because right term is really deriv of heating,
184  * but will be used below as deriv of cooling */
186  }
187  else
188  {
189  hmi.HeatH2Dish_used = 0;
190  hmi.HeatH2Dexc_used = 0;
192  }
193  }
194 
195  else if( hmi.chH2_small_model_type == 'T' )
196  {
197  /* TH85 dissociation heating */
198  /* these come from approximations in TH85, see comments above */
202  }
203  else if( hmi.chH2_small_model_type == 'H' )
204  {
205  /* Burton et al. 1990 */
209  }
210  else if( hmi.chH2_small_model_type == 'B')
211  {
212  /* Bertoldi & Draine */
216  }
217  else if(hmi.chH2_small_model_type == 'E')
218  {
219  /* this is the default when small H2 used */
223  }
224  else
225  TotalInsanity();
226 
227  /* heating due to photodissociation heating */
229 
230  /* heating (usually cooling in big H2) due to collisions within X */
231  /* add to heating is net heating is positive */
232  thermal.heating[0][8] = MAX2(0.,hmi.HeatH2Dexc_used);
233 
234  /* add to cooling if net heating is negative */
235  CoolAdd("H2cX",0,MAX2(0.,-hmi.HeatH2Dexc_used));
236  /*fprintf(ioQQQ,"DEBUG coolh2\t%.2f\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
237  fnzone, phycon.te, dense.eden, hmi.H2_total, thermal.ctot, -hmi.HeatH2Dexc_used );*/
238  /* add to net derivative */
239  /*thermal.dCooldT += MAX2(0.,-hmi.HeatH2Dexc_used)* ( 30172. * thermal.tsq1 - thermal.halfte );*/
240  /* >>chng 04 jan 25, check sign to prevent cooling from entering here,
241  * also enter neg sign since going into cooling stack (bug), in heatsum
242  * same term adds to deriv of heating */
243  if( hmi.HeatH2Dexc_used < 0. )
245 
246  /* H + H+ => H2+ cooling */
247  CoolHeavy.H2PlsCool = (realnum)(MAX2((2.325*phycon.te-1875.)*1e-20,0.)*
249 
250  if( h2.lgH2ON )
251  {
252  /* this is simplified approximation to H2 rotation cooling,
253  * big molecule goes this far better */
254  CoolHeavy.h2line = 0.;
255  }
256  else
257  {
258  /* rate for rotation lines from
259  * >>refer h2 cool Lepp, S., & Shull, J.M. 1983, ApJ, 270, 578 */
260  x = phycon.alogte - 4.;
261  if( phycon.te > 1087. )
262  {
263  rothi = 3.90e-19*sexp(6118./phycon.te);
264  }
265  else
266  {
267  rothi = pow(10.,-19.24 + 0.474*x - 1.247*x*x);
268  }
269 
270  /* low density rotation cooling */
271  /*&qn = pow(MAX2(hmi.Hmolec[ipMH2g],1e-37),0.77) + 1.2*pow(MAX2(dense.xIonDense[ipHYDROGEN][0],1e-37),0.77);*/
272  qn = pow(MAX2(hmi.H2_total,1e-37),0.77) + 1.2*pow(MAX2(dense.xIonDense[ipHYDROGEN][0],1e-37),0.77);
273  /* these are equations 11 from LS83 */
274  if( phycon.te > 4031. )
275  {
276  rotlow = 1.38e-22*sexp(9243./phycon.te)*qn;
277  }
278  else
279  {
280  rotlow = pow(10.,-22.90 - 0.553*x - 1.148*x*x)*qn;
281  }
282 
283  if( rotlow > 0. )
284  {
285  /*CoolHeavy.h2line = hmi.Hmolec[ipMH2g]*rothi/(1. + rothi/rotlow);*/
286  CoolHeavy.h2line = hmi.H2_total*rothi/(1. + rothi/rotlow);
287  }
288  else
289  {
290  CoolHeavy.h2line = 0.;
291  }
292  }
293 
294  {
295  enum {DEBUG_LOC=false};
296  if( DEBUG_LOC && nzone>187&& iteration > 1)
297  {
298  fprintf(ioQQQ,"h2coolbug\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
299  phycon.te,
300  CoolHeavy.h2line,
301  hmi.H2_total,
302  hmi.Hmolec[ipMHm],
304  rothi,
305  rotlow );
306  }
307  }
308 
309  /* >>chng 02 mar 07, add DH cooling using rates (eqn 6) from
310  * >>refer HD cooling Puy, D., Grenacher, L, & Jetzer, P., 1999, A&A, 345, 723 */
311  factor = sexp(128.6/phycon.te);
312  /*CoolHeavy.HD = 2.66e-21 * hydro.D2H_ratio * POW2(hmi.Hmolec[ipMH2g]) * phycon.sqrte *
313  factor/(1416.+phycon.sqrte*hmi.Hmolec[ipMH2g] * (1. + 3.*factor));*/
314  CoolHeavy.HD = 2.66e-21 * hydro.D2H_ratio * POW2((double)hmi.H2_total) * phycon.sqrte *
315  /*factor/(1416.+phycon.sqrte*hmi.Hmolec[ipMH2g] * (1. + 3.*factor));*/
316  factor/(1416.+phycon.sqrte*hmi.H2_total * (1. + 3.*factor));
317  }
318 
319  /* cooling due to charge transfer ionization / recombination */
320  CoolAdd("CT C" , 0. , thermal.char_tran_cool );
321 
322  /* H- FB; H + e -> H- + hnu */
323  /* H- FF is in with H ff */
324  CoolAdd("H-fb",0,hmi.hmicol);
325 
326  /* >>chng 96 nov 15, fac of 2 in deriv to help convergence in very dense
327  * models where H- is important, this takes change in eden into
328  * partial account */
330  if( PRT_DERIV )
331  fprintf(ioQQQ,"DEBUG dcdt 2 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);
332 
333  CoolAdd("H2ln",0,CoolHeavy.h2line);
334  /* >>chng 00 oct 21, added coef of 3.5, sign had been wrong */
335  /*thermal.dCooldT += CoolHeavy.h2line*phycon.teinv;*/
336  /* >>chng 03 mar 17, change 3.5 to 15 as per behavior in primal.in */
337  /*thermal.dCooldT += 3.5*CoolHeavy.h2line*phycon.teinv;*/
338  /* >>chng 03 may 18, from 15 to 30 as per behavior in primal.in - overshoots happen */
339  /*thermal.dCooldT += 15.*CoolHeavy.h2line*phycon.teinv;*/
340  /*>>chng 03 oct 03, from 30 to 3, no more overshoots in primalin */
341  /*thermal.dCooldT += 30.*CoolHeavy.h2line*phycon.teinv;*/
343 
344  {
345  /* problems with H2 cooling */
346  enum {DEBUG_LOC=false};
347  if( DEBUG_LOC /*&& nzone>300 && iteration > 1*/)
348  {
349  fprintf(ioQQQ,"CoolEvaluate debuggg\t%.2f\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\n",
350  fnzone,
351  phycon.te,
352  hmi.H2_total ,
354  hmi.Hmolec[ipMHm] ,
355  dense.eden);
356  }
357  }
358 
359  CoolAdd("HDro",0,CoolHeavy.HD);
361 
362  CoolAdd("H2+ ",0,CoolHeavy.H2PlsCool);
364 
365  /* cooling due to C12O16 */
366  CoolAdd("CO C",12,CoolHeavy.C12O16Rot);
368  /* >>chng 00 oct 25, add C13O16 cooling */
369  /* cooling due to C13O16 */
370  CoolAdd("CO C",13,CoolHeavy.C13O16Rot);
372 
373  /* heating due to three-body, will be incremented in iso_cool*/
374  thermal.heating[0][3] = 0.;
375  /* heating due to hydrogen lines */
376  thermal.heating[0][23] = 0.;
377  /* heating due to photoionization of all excited states of hydrogen species */
378  thermal.heating[0][1] = 0.;
379 
380  /* isoelectronic species cooling, mainly lines, and level ionization */
381  for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
382  {
383  for( nelem=ipISO; nelem < LIMELM; nelem++ )
384  {
385  /* must always call iso_cool since we must zero those variables
386  * that would have been set had the species been present */
387  iso_cool( ipISO , nelem );
388  }
389  }
390 
391  /* >>chng 02 jun 18, don't reevaluate needlessly */
392  /* >>chng 03 nov 28, even faster - special logic for when ff is pretty
393  * negligible - eval of ff is pretty slow */
394  /* >>chng 04 feb 19, must not test on temp not changing, since ionization
395  * can change at constant temperature
396  * >>chng 04 sep 14, above introduced bug since brems never reevaluated
397  * now test is zone or temp has changed */
398  if( (fabs(CoolHeavy.brems_cool_net)/SDIV(thermal.ctot)>0.001 ) ||
399  conv.lgSearch ||
400  fabs((phycon.te-TeUsedBrems)/phycon.te-1.)>0.01 ||
401  nzone!=nzoneUsedBrems )
402  {
403  double BremsThisEnergy;
404  /*double OpacityThisIon;*/
405  long int limit;
406  /* free-free free free brems emission for all ions */
407 
408  TeUsedBrems = phycon.te;
409  nzoneUsedBrems = nzone;
410  /* highest frequency where we have non-zero Boltzmann factors */
411  limit = MIN2( rfield.ipMaxBolt , rfield.nflux );
412 
414  CoolHeavy.brems_cool_h = 0.;
418 
419  {
420  double bhfac;
421  realnum sumion[LIMELM+1];
422  long int ion_lo , ion_hi;
423 
425  ASSERT(limit < rfield.nupper);
426 
427  /* ipEnergyBremsThin is index to energy where gas becomes optically thin to brems,
428  * so this loop is over optically thin frequencies
429  * do not include optically thick part as net emission since self absorbed */
430 
431  /* do hydrogen first, before main loop since want to break out as separate
432  * coolant, and what to add on H- brems */
433  nelem = ipHYDROGEN;
434  ion = 1;
435  CoolHeavy.brems_cool_h = 0.;
437  /* this is all done in opacity_addtotal - why do here too? */
438  for( i=rfield.ipEnergyBremsThin; i < limit; i++ )
439  {
440 
441  /* in all following CoolHeavy.lgFreeOn is flag set with 'no free-free' to
442  * turn off brems heating and cooling */
443  BremsThisEnergy = rfield.gff[ion][i]*rfield.widflx[i]*rfield.ContBoltz[i];
444  /*ASSERT( BremsThisEnergy >= 0. );*/
445  CoolHeavy.brems_cool_h += BremsThisEnergy;
446 
447  /* for case of hydrogen, do H- brems - OpacStack contains the ratio
448  * of the H- to H brems cross section - multiply by this and the fraction in ground */
450 
451  }
452  bhfac = dense.xIonDense[nelem][ion]*CoolHeavy.lgFreeOn* dense.eden*1.032e-11/phycon.sqrte*EN1RYD;
453  CoolHeavy.brems_cool_h *= bhfac;
455 
456  /* now do helium, both He+ and He++ */
457  nelem = ipHELIUM;
459  for( i=rfield.ipEnergyBremsThin; i < limit; i++ )
460  {
461  /* eff. charge is ion, so first rfield.gff argument must be "ion". */
462  BremsThisEnergy =
463  (dense.xIonDense[nelem][1]*rfield.gff[1][i] + 4.*dense.xIonDense[nelem][2]*rfield.gff[2][i])*
465  CoolHeavy.brems_cool_he += BremsThisEnergy;
466  }
469 
470  /* >>chng 05 jul 13, rewrite this for speed */
471  /* gaunt factors depend only on photon energy and ion charge, so do
472  * sum of ions here before entering into loop over photon energy */
473  sumion[0] = 0.;
474  for(ion=1; ion<=LIMELM; ++ion )
475  {
476  sumion[ion] = 0.;
477  for( nelem=ipLITHIUM; nelem < LIMELM; ++nelem )
478  {
479  if( dense.lgElmtOn[nelem] && ion<=nelem+1 )
480  {
481  sumion[ion] += dense.xIonDense[nelem][ion];
482  }
483  }
484  /* now include the charge, density, and temperature */
485  sumion[ion] *= POW2((realnum)ion);
486  }
487  /* now find lowest and highest ion we need to consider - following loop is over
488  * full continuum and eats time
489  * >>chng 05 oct 19, bounds check had been on ion, rather than ion_lo and ion_hi, so
490  * array bounds were exceeded */
491  ion_lo = 1;
492  while( sumion[ion_lo]==0 && ion_lo<LIMELM-1 )
493  ++ion_lo;
494  ion_hi = LIMELM;
495  while( sumion[ion_hi]==0 && ion_hi>0 )
496  --ion_hi;
497 
498  /* heavy elements */
501  for( i=rfield.ipEnergyBremsThin; i < limit; i++ )
502  {
503  BremsThisEnergy = 0.;
504  for(ion=ion_lo; ion<=ion_hi; ++ion )
505  {
506  BremsThisEnergy += sumion[ion]*rfield.gff[ion][i];
507  }
508  CoolHeavy.brems_cool_metals += BremsThisEnergy*rfield.widflx[i]*rfield.ContBoltz[i];
509  /* the total heating due to bremsstrahlung */
511  }
514  {
515  enum {DEBUG_LOC=false};
516  if( DEBUG_LOC && nzone>60 /*&& iteration > 1*/)
517  {
518  double sumfield = 0., sumtot=0., sum1=0., sum2=0.;
519  for( i=rfield.ipEnergyBremsThin; i<limit; i++ )
520  {
521  sumtot += opac.FreeFreeOpacity[i]*rfield.flux[i]*rfield.anu[i];
522  sumfield += rfield.flux[i]*rfield.anu[i];
523  sum1 += opac.FreeFreeOpacity[i]*rfield.flux[i]*rfield.anu[i];
524  sum2 += opac.FreeFreeOpacity[i]*rfield.flux[i];
525  }
526  fprintf(ioQQQ,"DEBUG brems heat\t%.2f\t%.3e\t%.3e\t%.3e\t%e\t%.3e\t%.3e\n",
527  fnzone,
529  sumtot/SDIV(sumfield) ,
530  sum1/SDIV(sum2),
531  phycon.te ,
532  rfield.gff[1][1218],
533  opac.FreeFreeOpacity[1218]);
534  }
535  }
536  }
537  }
538 
539 
540  /* these two terms are both large, nearly canceling, near lte */
547  /*fprintf(ioQQQ,"DEBUG brems\t%.2f\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
548  fnzone,
549  phycon.te,
550  CoolHeavy.brems_cool_net,
551  CoolHeavy.brems_cool_h ,
552  CoolHeavy.brems_cool_he ,
553  CoolHeavy.brems_cool_hminus,
554  CoolHeavy.brems_cool_metals ,
555  CoolHeavy.brems_heat_total);*/
556 
557  /* net free free brems cooling, count as cooling if positive */
558  CoolAdd( "FF c" , 0, MAX2(0.,CoolHeavy.brems_cool_net) );
559 
560  /* now stuff into heating array if negative */
562 
563  /* >>chng 96 oct 30, from HFFNet to just FreeFreeCool,
564  * since HeatSum picks up CoolHeavy.brems_heat_total */
566  if( PRT_DERIV )
567  fprintf(ioQQQ,"DEBUG dcdt 3 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);
568 
569  /* >>chng 02 jun 21, net cooling already includes this */
570  /* end of brems cooling */
571 
572  /* heavy element recombination cooling, do not count hydrogenic since
573  * already done above, also helium singlets have been done */
574  /* >>chng 02 jul 21, put in charge dependent rec term */
575  CoolHeavy.heavfb = 0.;
576  for( nelem=ipLITHIUM; nelem < LIMELM; nelem++ )
577  {
578  if( dense.lgElmtOn[nelem] )
579  {
580  /*>>chng 02 jul 21, upper bound now uses NISO,
581  * note that detailed iso seq atoms are done in iso_cool */
582  /* >>chng 03 sep 03, do not do zeroed stages of ionization */
583  long limit_lo = MAX2( 1 , dense.IonLow[nelem] );
584  /* there is a -1 on the NISO since loop will be up to <= not < */
585  long limit_hi = MIN2( nelem-NISO+1-1 , dense.IonHigh[nelem] );
586  for( ion=limit_lo; ion<=limit_hi; ++ion )
587  /*for( ion=1; ion < nelem-NISO+1; ion++ )*/
588  {
589  /* factor of 0.9 is roughly correct for nebular conditions, see
590  * >>refer H rec cooling LaMothe, J., & Ferland, G.J., 2001, PASP, 113, 165 */
591  /* note that ionbal.RR_rate_coef_used is the rate coef, cm3 s-1, needs eden */
592  /* >>chng 02 nov 07, move rec arrays around, this now has ONLY rad rec,
593  * previously had included grain rec and three body */
594  /* recombination cooling for iso-seq atoms are done in iso_cool */
595  double one = dense.xIonDense[nelem][ion] * ionbal.RR_rate_coef_used[nelem][ion-1]*
597  /*fprintf(ioQQQ,"debugggfb\t%li\t%li\t%.3e\t%.3e\t%.3e\n", nelem, ion, one,
598  dense.xIonDense[nelem][ion] , ionbal.RR_rate_coef_used[nelem][ion]);*/
599  CoolHeavy.heavfb += one;
600  }
601  }
602  }
603 
604  /*fprintf(ioQQQ,"debuggg hvFB\t%i\t%.2f\t%.2e\t%.2e\n",iteration, fnzone,CoolHeavy.heavfb, dense.eden);*/
605 
606  CoolAdd("hvFB",0,CoolHeavy.heavfb);
608  if( PRT_DERIV )
609  fprintf(ioQQQ,"DEBUG dcdt 4 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);
610 
611  /* electron-electron brems, approx form from
612  * >>refer ee brems Stepney and Guilbert, MNRAS 204, 1269 (1983)
613  * ok for T<10**9 */
614  CoolHeavy.eebrm = POW2(dense.eden*phycon.te*1.84e-21);
615 
616  /* >>chng 97 mar 12, added deriv */
618  CoolAdd("eeff",0,CoolHeavy.eebrm);
619 
620  /* add advective heating and cooling */
621  /* this is cooling due to loss of matter from this region */
622  CoolAdd("adve",0,dynamics.Cool );
623  /* >>chng02 dec 04, rm factor of 8 in front of dCooldT */
625  /* local heating due to matter moving into this location */
626  thermal.heating[1][5] = dynamics.Heat;
628 
629  /* total Compton cooling */
631  CoolAdd("Comp",0,CoolHeavy.tccool);
633 
634  /* option for "extra" cooling, expressed as power-law in temperature, these
635  * are set with the CEXTRA command */
636  if( thermal.lgCExtraOn )
637  {
638  CoolHeavy.cextxx =
639  (realnum)(thermal.CoolExtra*pow(phycon.te/1e4,(double)thermal.cextpw));
640  }
641  else
642  {
643  CoolHeavy.cextxx = 0.;
644  }
645  CoolAdd("Extr",0,CoolHeavy.cextxx);
646 
647  /* cooling due to wind expansion, only for winds expansion cooling */
648  if( wind.windv > 0. )
649  {
651  radius.Radius);
653  }
654  else
655  {
656  dynamics.dDensityDT = 0.;
657  CoolHeavy.expans = 0.;
658  }
659  CoolAdd("Expn",0,CoolHeavy.expans);
661 
662  /* cyclotron cooling */
663  /* coef is 4/3 /8pi /c * vtherm(elec) */
665  CoolAdd("Cycl",0,CoolHeavy.cyntrn);
667 
668  /* heavy element collisional ionization
669  * derivative should be zero since increased coll ion rate
670  * decreases neutral fraction by proportional amount */
671  CoolAdd("Hvin",0,CoolHeavy.colmet);
672 
673  /* evaluate H 21 cm spin changing collisions */
674  if( fabs(phycon.te-TeEvalCS_21cm)/phycon.te > 0.001 )
675  {
676  {
677  /* this prints table of rates at points given in original data paper */
678  enum {DEBUG_LOC=false};
679  if( DEBUG_LOC )
680  {
681 # define N21CM_TE 16
682  int n;
683  double teval[N21CM_TE]={2.,5.,10.,20.,50.,100.,200.,500.,1000.,
684  2000.,3000.,5000.,7000.,10000.,15000.,20000.};
685  for( n = 0; n<N21CM_TE; ++n )
686  {
687  fprintf(
688  ioQQQ,"DEBUG 21 cm deex Te=\t%.2e\tH0=\t%.2e\tp=\t%.2e\te=\t%.2e\n",
689  teval[n],
690  H21cm_H_atom( teval[n] ),
691  H21cm_proton( teval[n] ),
692  H21cm_electron( teval[n] ) );
693  }
694  cdEXIT(EXIT_FAILURE);
695 # undef N21CM_TE
696  }
697  }
698  /*only evaluate T dependent part when Te changes, but update
699  * density part below since densities may constantly change */
700  atomic_rate_21cm = H21cm_H_atom( phycon.te );
701  proton_rate_21cm = H21cm_proton( phycon.te );
702  electron_rate_21cm = H21cm_electron( phycon.te );
703  TeEvalCS_21cm = phycon.te;
704  }
705  /* H 21 cm emission/population,
706  * cs will be sum of e cs and H cs converted from rate */
707  cs = (electron_rate_21cm * dense.eden +
708  atomic_rate_21cm * dense.xIonDense[ipHYDROGEN][0] +
709  proton_rate_21cm * dense.xIonDense[ipHYDROGEN][1] ) *
710  3./ dense.cdsqte;
711  PutCS( cs , &HFLines[0] );
712 
713  /* fine structure lines */
714  if( fabs(phycon.te-TeEvalCS)/phycon.te > 0.05 )
715  {
716  /* H 21 cm done above, now loop over remaining lines to get collision strengths */
717  for( i=1; i < nHFLines; i++ )
718  {
719  cs = HyperfineCS( i );
720  /* now generate the collision strength and put into the line array */
721  PutCS( cs , &HFLines[i] );
722  }
723  TeEvalCS = phycon.te;
724  }
725 
726  /* do level pops for H 21 cm which is a special case since Lya pumping in included */
727  H21_cm_pops();
728  if( PRT_DERIV )
729  fprintf(ioQQQ,"DEBUG dcdt 5 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);
730 
731  /* find total cooling due to hyperfine structure lines */
733 
734  /* now do level pops for all except 21 cm */
735  for( i=1; i < nHFLines; i++ )
736  {
737  /* remember current gas-phase abundances */
738  realnum save = dense.xIonDense[HFLines[i].Hi->nelem-1][HFLines[i].Hi->IonStg-1];
739 
740  /* bail if no abundance */
741  if( save<=0. ) continue;
742 
743  /* set gas-phase abundance to total times isotope ratio */
744  dense.xIonDense[HFLines[i].Hi->nelem-1][HFLines[i].Hi->IonStg-1] *=
746 
747  /* use the collision strength generated above and find pops and cooling */
748  atom_level2( &HFLines[i] );
749 
750  /* put the correct gas-phase abundance back in the array */
751  dense.xIonDense[HFLines[i].Hi->nelem-1][HFLines[i].Hi->IonStg-1] = save;
752 
753  /* find total cooling due to hyperfine structure lines */
754  hyperfine.cooling_total += HFLines[i].Coll.cool;
755  }
756  if( PRT_DERIV )
757  fprintf(ioQQQ,"DEBUG dcdt 6 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);
758 
759  /* Carbon cooling */
760  CoolCarb();
761  if( PRT_DERIV )
762  fprintf(ioQQQ,"DEBUG dcdt C %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);
763 
764  /* Nitrogen cooling */
765  CoolNitr();
766  if( PRT_DERIV )
767  fprintf(ioQQQ,"DEBUG dcdt N %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);
768 
769  /* Oxygen cooling */
770  CoolOxyg();
771  if( PRT_DERIV )
772  fprintf(ioQQQ,"DEBUG dcdt 7 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);
773 
774  /* Fluorine cooling */
775  CoolFluo();
776 
777  /* Neon cooling */
778  CoolNeon();
779  if( PRT_DERIV )
780  fprintf(ioQQQ,"DEBUG dcdt Ne %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);
781 
782  /* Magnesium cooling */
783  CoolMagn();
784  if( PRT_DERIV )
785  fprintf(ioQQQ,"DEBUG dcdt 8 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);
786 
787  /* Sodium cooling */
788  CoolSodi();
789  if( PRT_DERIV )
790  fprintf(ioQQQ,"DEBUG dcdt Na %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);
791 
792  /* Aluminum cooling */
793  CoolAlum();
794  if( PRT_DERIV )
795  fprintf(ioQQQ,"DEBUG dcdt Al %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);
796 
797  /* Silicon cooling */
798  CoolSili();
799  if( PRT_DERIV )
800  fprintf(ioQQQ,"DEBUG dcdt 9 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);
801 
802  /* Phosphorus */
803  CoolPhos();
804 
805  /* Sulphur cooling */
806  CoolSulf();
807 
808  /* Chlorine cooling */
809  CoolChlo();
810 
811  /* Argon cooling */
812  CoolArgo();
813  if( PRT_DERIV )
814  fprintf(ioQQQ,"DEBUG dcdt 10 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);
815 
816  /* Potasium cooling */
817  CoolPota();
818 
819  /* Calcium cooling */
820  CoolCalc();
821 
822  /* Scandium cooling */
823  CoolScan();
824 
825  /* Titanium cooling */
826  CoolTita();
827  if( PRT_DERIV )
828  fprintf(ioQQQ,"DEBUG dcdt 11 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);
829 
830  /* Vanadium cooling */
831  CoolVana();
832 
833  /* Chromium cooling */
834  CoolChro();
835 
836  /* Iron cooling */
837  CoolIron();
838  if( PRT_DERIV )
839  fprintf(ioQQQ,"DEBUG dcdt 12 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);
840 
841  /* Manganese cooling */
842  CoolMang();
843 
844  /* Cobalt cooling */
845  CoolCoba();
846 
847  /* Nickel cooling */
848  CoolNick();
849 
850  /* Zinc cooling */
851  CoolZinc();
852  if( PRT_DERIV )
853  fprintf(ioQQQ,"DEBUG dcdt 13 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);
854 
855  /* do all the thousands of lines Dima added with g-bar approximation */
856  CoolDima();
857 
858  /* do Chianti & Leiden database atoms and molecules here */
859  atmol_popsolve();
860 
861  /* now add up all the coolants */
862  CoolSum(tot);
863  if( PRT_DERIV )
864  fprintf(ioQQQ,"DEBUG dcdt 14 %.3e %.3e\n",thermal.dCooldT , thermal.dHeatdT);
865 
866  /* negative cooling */
867  if( *tot <= 0. )
868  {
869  fprintf( ioQQQ, " COOLR; cooling is <=0, this is impossible.\n" );
870  ShowMe();
871  cdEXIT(EXIT_FAILURE);
872  }
873 
874  /* bad derivative */
875  if( thermal.dCooldT == 0. )
876  {
877  fprintf( ioQQQ, " COOLR; cooling slope <=0, this is impossible.\n" );
878  if( *tot > 0. && dense.gas_phase[ipHYDROGEN] < 1e-4 )
879  {
880  fprintf( ioQQQ, " Probably due to very low density.\n" );
881  }
882  ShowMe();
883  cdEXIT(EXIT_FAILURE);
884  }
885 
886  if( trace.lgTrace )
887  {
888  fndstr(*tot,thermal.dCooldT);
889  }
890 
891  /* lgTSetOn true for constant temperature model */
892  if( (((((!thermal.lgTemperatureConstant) && *tot < 0.) && called.lgTalk) &&
893  !conv.lgSearch) && thermal.lgCNegChk) && nzone > 0 )
894  {
895  fprintf( ioQQQ,
896  " NOTE Negative cooling, zone %4ld, =%10.2e coola=%10.2e CHION=%10.2e Te=%10.2e\n",
897  nzone,
898  *tot,
900  iso.coll_ion[ipH_LIKE][ipHYDROGEN],
901  phycon.te );
902  fndneg();
903  }
904 
905  /* possibility of getting empirical cooling derivative
906  * normally false, set true with 'set numerical derivatives' command */
907  if( NumDeriv.lgNumDeriv )
908  {
909  if( ((nzone > 2 && nzone == nzSave) && ! fp_equal( oldtemp, phycon.te )) && nhit > 4 )
910  {
911  /* hnit is number of tries on this zone - use to stop numerical problems
912  * do not evaluate numerical deriv until well into solution */
913  deriv = (oltcool - *tot)/(oldtemp - phycon.te);
914  thermal.dCooldT = deriv;
915  }
916  else
917  {
918  deriv = thermal.dCooldT;
919  }
920  if( nzone != nzSave )
921  nhit = 0;
922 
923  nzSave = nzone;
924  nhit += 1;
925  oltcool = *tot;
926  oldtemp = phycon.te;
927  }
928  return;
929 }
930 
931 /* */
932 #ifdef EPS
933 # undef EPS
934 #endif
935 #define EPS 0.01
936 
937 /*fndneg search cooling array to find negative values */
938 STATIC void fndneg(void)
939 {
940  long int i;
941  double trig;
942 
943  DEBUG_ENTRY( "fndneg()" );
944 
945  trig = fabs(thermal.htot*EPS);
946  for( i=0; i < thermal.ncltot; i++ )
947  {
948  if( thermal.cooling[i] < 0. && fabs(thermal.cooling[i]) > trig )
949  {
950  fprintf( ioQQQ, " negative line=%s %.2f fraction of heating=%.3f\n",
952  thermal.htot );
953  }
954 
955  if( thermal.heatnt[i] > trig )
956  {
957  fprintf( ioQQQ, " heating line=%s %.2f fraction of heating=%.3f\n",
959  thermal.htot );
960  }
961  }
962  return;
963 }
964 
965 /*fndstr search cooling stack to find strongest values */
966 STATIC void fndstr(double tot,
967  double dc)
968 {
969  char chStrngLab[NCOLNT_LAB_LEN+1];
970  long int i;
971  realnum wl;
972  double str,
973  strong;
974 
975  DEBUG_ENTRY( "fndstr()" );
976 
977  strong = 0.;
978  wl = -FLT_MAX;
979  for( i=0; i < thermal.ncltot; i++ )
980  {
981  if( fabs(thermal.cooling[i]) > strong )
982  {
983  /* this is the wavelength of the coolant, 0 for a continuum*/
984  wl = thermal.collam[i];
985  /* make sure labels are all valid*/
986  /*>>chng 06 jun 06, bug fix, assert length was ==4, should be <=NCOLNT_LAB_LEN */
987  ASSERT( strlen( thermal.chClntLab[i] ) <= NCOLNT_LAB_LEN );
988  strcpy( chStrngLab, thermal.chClntLab[i] );
989  strong = fabs(thermal.cooling[i]);
990  }
991  }
992 
993  str = strong;
994 
995  fprintf( ioQQQ,
996  " fndstr cool: TE=%10.4e Ne=%10.4e C=%10.3e dC/dT=%10.2e ABS(%s %.1f)=%.2e nz=%ld\n",
997  phycon.te, dense.eden, tot, dc, chStrngLab
998  , wl, str, nzone );
999 
1000  /* option for extensive printout of lines */
1001  if( trace.lgCoolTr )
1002  {
1003  realnum ratio;
1004 
1005  /* flag all significant coolants, first zero out the array */
1006  coolpr(ioQQQ,(char*)thermal.chClntLab[0],1,0.,"ZERO");
1007 
1008  /* push all coolants onto the stack */
1009  for( i=0; i < thermal.ncltot; i++ )
1010  {
1011  /* usually positive, although can be neg for coolants that heats,
1012  * only do positive here */
1013  ratio = (realnum)(thermal.cooling[i]/thermal.ctot);
1014  if( ratio >= EPS )
1015  {
1016  /*>>chng 99 jan 27, only cal when ratio is significant */
1017  coolpr(ioQQQ,(char*)thermal.chClntLab[i],thermal.collam[i], ratio,"DOIT");
1018  }
1019  }
1020 
1021  /* complete the printout for positive coolants */
1022  coolpr(ioQQQ,"DONE",1,0.,"DONE");
1023 
1024  /* now do cooling that was counted as a heat source if significant */
1025  if( thermal.heating[0][22]/thermal.ctot > 0.05 )
1026  {
1027  fprintf( ioQQQ,
1028  " All coolant heat greater than%6.2f%% of the total will be printed.\n",
1029  EPS*100. );
1030 
1031  coolpr(ioQQQ,"ZERO",1,0.,"ZERO");
1032  for( i=0; i < thermal.ncltot; i++ )
1033  {
1034  ratio = (realnum)(thermal.heatnt[i]/thermal.ctot);
1035  if( fabs(ratio) >=EPS )
1036  {
1037  coolpr(ioQQQ,(char*)thermal.chClntLab[i],thermal.collam[i],
1038  ratio,"DOIT");
1039  }
1040  }
1041  coolpr(ioQQQ,"DONE",1,0.,"DONE");
1042  }
1043  }
1044  return;
1045 }

Generated for cloudy by doxygen 1.8.1.1