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