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