cloudy trunk
Loading...
Searching...
No Matches
cool_eval.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/*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 "rt.h"
14#include "doppvel.h"
15#include "opacity.h"
16#include "ionbal.h"
17#include "dense.h"
18#include "trace.h"
19#include "dynamics.h"
20#include "rfield.h"
21#include "grainvar.h"
22#include "atmdat.h"
23#include "atoms.h"
24#include "called.h"
25#include "mole.h"
26#include "hmi.h"
27#include "numderiv.h"
28#include "magnetic.h"
29#include "phycon.h"
30#include "lines_service.h"
31#include "hyperfine.h"
32#include "iso.h"
33#include "thermal.h"
34#include "cooling.h"
35#include "pressure.h"
36/*fndneg search cooling array to find negative values */
37STATIC void fndneg(void);
38/*fndstr search cooling stack to find strongest values */
39STATIC void fndstr(double tot,
40 double dc);
41
42/* set true to debug derivative of heating and cooling */
43static const bool PRT_DERIV = false;
44
45void CoolEvaluate(double *tot)
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 long int coolnum, coolcal;
71
72 DEBUG_ENTRY( "CoolEvaluate()" );
73
74 /* returns tot, the total cooling,
75 * and dc, the derivative of the cooling */
76
77 /* routine atom_level2( t10 )
78 * routine atom_level3( abund , t10,t21,t20)
79 * tsq1 = 1. / (te**2)
80 * POPEXC( O12,g1,g2,A21,excit,abund); result already*a21
81 * POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2)
82 * AtomSeqBeryllium(cs23,cs24,cs34,tarray,a41)
83 * FIVEL( G(1-5) , ex(wn,1-5), cs12,cs13,14,15,23,24,25,34,35,45,
84 * A21,31,41,51,32,42,52,43,53,54, pop(1-5), abund) */
85
86 if( trace.lgTrace )
87 fprintf( ioQQQ, " COOLR TE:%.4e zone %li %li Cool:%.4e Heat:%.4e eden:%.4e edenTrue:%.4e\n",
88 phycon.te,
89 nzone, conv.nPres2Ioniz ,
90 thermal.ctot , thermal.htot,dense.eden,dense.EdenTrue );
91
92 /* must call TempChange since ionization has changed, there are some
93 * terms that affect collision rates (H0 term in electron collision) */
94 TempChange(phycon.te , false);
95
96 /* now zero out the cooling stack */
97 CoolZero();
98 if( PRT_DERIV )
99 fprintf(ioQQQ,"DEBUG dCdT 0 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
100 if( gv.lgGrainPhysicsOn )
101 {
102 /* grain heating and cooling */
103 /* grain recombination cooling, evaluated elsewhere
104 * can either heat or cool the gas, do cooling here */
105 CoolAdd("dust",0,MAX2(0.,gv.GasCoolColl));
106
107 /* grain cooling proportional to temperature ^3/2 */
108 thermal.dCooldT += MAX2(0.,gv.GasCoolColl)*3./(2.*phycon.te);
109
110 /* these are the various heat agents from grains */
111 /* options to force gas heating or cooling by grains to zero - for tests only ! */
112 if( gv.lgDustOn() && gv.lgDHetOn )
113 {
114 /* rate dust heats gas by photoelectric effect */
115 thermal.heating[0][13] = gv.GasHeatPhotoEl;
116
117 /* if grains hotter than gas then collisions with gas act
118 * to heat the gas, add this in here
119 * a symmetric statement appears in COOLR, where cooling is added on */
120 thermal.heating[0][14] = MAX2(0.,-gv.GasCoolColl);
121
122 /* this is gas heating due to thermionic emissions */
123 thermal.heating[0][25] = gv.GasHeatTherm;
124 }
125 else
126 {
127 thermal.heating[0][13] = 0.;
128 thermal.heating[0][14] = 0.;
129 thermal.heating[0][25] = 0.;
130 }
131 }
132 else if( gv.lgBakesPAH_heat )
133 {
134 /* >>chng 06 jul 21, option to include Bakes PAH hack with grain physics off,
135 * needed to test dynamics models */
136 thermal.heating[0][13] = gv.GasHeatPhotoEl;
137 }
138
139 if( PRT_DERIV )
140 fprintf(ioQQQ,"DEBUG dCdT 1 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
141
142 /* molecular molecules molecule cooling */
143 if( mole_global.lgNoMole )
144 {
145 /* this branch - do not include molecules */
146 hmi.hmicol = 0.;
147 CoolHeavy.brems_cool_hminus = 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 hmi.deriv_HeatH2Dexc_used = 0.;
163 }
164
165 else
166 {
167 /* save various molecular heating/cooling agent */
168 thermal.heating[0][15] = hmi.hmihet;
169 thermal.heating[0][16] = hmi.h2plus_heat;
170 /* now get heating from H2 molecule, either simple or from big one */
171 if( h2.lgEnabled && hmi.lgH2_Thermal_BigH2 )
172 {
173 if( h2.lgEvaluated )
174 {
175 /* these are explicitly from big H2 molecule,
176 * first is heating due to radiative pump of excited states, followed by
177 * radiative decay into continuum of X, followed by dissociation of molecule
178 * with kinetic energy, typically 0.25 - 0.5 eV per event */
179 hmi.HeatH2Dish_used = h2.HeatDiss;
180 hmi.HeatH2Dexc_used = h2.HeatDexc;
181 if (0)
182 fprintf(ioQQQ,"DEBUG big %.2f\t%.5e\t%.2e\t%.2e\t%.2e\n",
183 fnzone , phycon.te, hmi.HeatH2Dexc_used,
184 hmi.H2_total, dense.gas_phase[ipHYDROGEN] );
185 /* negative sign because right term is really deriv of heating,
186 * but will be used below as deriv of cooling */
187 hmi.deriv_HeatH2Dexc_used = -h2.HeatDexc_deriv;
188 }
189 else
190 {
191 hmi.HeatH2Dish_used = 0;
192 hmi.HeatH2Dexc_used = 0;
193 hmi.deriv_HeatH2Dexc_used = 0;
194 }
195 }
196
197 else if( hmi.chH2_small_model_type == 'T' )
198 {
199 /* TH85 dissociation heating */
200 /* these come from approximations in TH85, see comments above */
201 hmi.HeatH2Dish_used = hmi.HeatH2Dish_TH85;
202 hmi.HeatH2Dexc_used = hmi.HeatH2Dexc_TH85;
203 hmi.deriv_HeatH2Dexc_used = hmi.deriv_HeatH2Dexc_TH85;
204 }
205 else if( hmi.chH2_small_model_type == 'H' )
206 {
207 /* Burton et al. 1990 */
208 hmi.HeatH2Dish_used = hmi.HeatH2Dish_BHT90;
209 hmi.HeatH2Dexc_used = hmi.HeatH2Dexc_BHT90;
210 hmi.deriv_HeatH2Dexc_used = hmi.deriv_HeatH2Dexc_BHT90;
211 }
212 else if( hmi.chH2_small_model_type == 'B')
213 {
214 /* Bertoldi & Draine */
215 hmi.HeatH2Dish_used = hmi.HeatH2Dish_BD96;
216 hmi.HeatH2Dexc_used = hmi.HeatH2Dexc_BD96;
217 hmi.deriv_HeatH2Dexc_used = hmi.deriv_HeatH2Dexc_BD96;
218 }
219 else if(hmi.chH2_small_model_type == 'E')
220 {
221 /* this is the default when small H2 used */
222 hmi.HeatH2Dish_used = hmi.HeatH2Dish_ELWERT;
223 hmi.HeatH2Dexc_used = hmi.HeatH2Dexc_ELWERT;
224 hmi.deriv_HeatH2Dexc_used = hmi.deriv_HeatH2Dexc_ELWERT;
225 }
226 else
228
229 /* heating due to photodissociation heating */
230 thermal.heating[0][17] = hmi.HeatH2Dish_used;
231
232 /* heating due to continuum photodissociation */
233 thermal.heating[0][28] = 0.;
234 {
235 double heat = 0.;
236 for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
237 {
238 if( (*diatom)->lgEnabled && mole_global.lgStancil )
239 {
240 heat += (*diatom)->Cont_Diss_Heat_Rate();
241 }
242 }
243 thermal.heating[0][28] += MAX2( 0., heat );
244 CoolAdd("H2cD",0,MAX2(0.,-heat));
245 }
246
247 /* heating (usually cooling in big H2) due to collisions within X */
248 /* add to heating is net heating is positive */
249 thermal.heating[0][8] = MAX2(0.,hmi.HeatH2Dexc_used);
250
251 /* add to cooling if net heating is negative */
252 CoolAdd("H2cX",0,MAX2(0.,-hmi.HeatH2Dexc_used));
253 /*fprintf(ioQQQ,"DEBUG coolh2\t%.2f\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
254 fnzone, phycon.te, dense.eden, hmi.H2_total, thermal.ctot, -hmi.HeatH2Dexc_used );*/
255 /* add to net derivative */
256 /*thermal.dCooldT += MAX2(0.,-hmi.HeatH2Dexc_used)* ( 30172. * thermal.tsq1 - thermal.halfte );*/
257 /* >>chng 04 jan 25, check sign to prevent cooling from entering here,
258 * also enter neg sign since going into cooling stack (bug), in heatsum
259 * same term adds to deriv of heating */
260 if( hmi.HeatH2Dexc_used < 0. )
261 thermal.dCooldT -= hmi.deriv_HeatH2Dexc_used;
262
263 /* H + H+ => H2+ cooling */
264 CoolHeavy.H2PlsCool = (realnum)(MAX2((2.325*phycon.te-1875.)*1e-20,0.)*
265 dense.xIonDense[ipHYDROGEN][0]*dense.xIonDense[ipHYDROGEN][1]*1.66e-11);
266
267 if( h2.lgEnabled )
268 {
269 /* this is simplified approximation to H2 rotation cooling,
270 * big molecule goes this far better */
271 CoolHeavy.h2line = 0.;
272 }
273 else
274 {
275 /* rate for rotation lines from
276 * >>refer h2 cool Lepp, S., & Shull, J.M. 1983, ApJ, 270, 578 */
277 x = phycon.alogte - 4.;
278 if( phycon.te > 1087. )
279 {
280 rothi = 3.90e-19*sexp(6118./phycon.te);
281 }
282 else
283 {
284 rothi = pow(10.,-19.24 + 0.474*x - 1.247*x*x);
285 }
286
287 /* low density rotation cooling */
288 /*&qn = pow(MAX2(findspecieslocal("H2")->den,1e-37),0.77) + 1.2*pow(MAX2(dense.xIonDense[ipHYDROGEN][0],1e-37),0.77);*/
289 qn = pow(MAX2(hmi.H2_total,1e-37),0.77) + 1.2*pow(MAX2(dense.xIonDense[ipHYDROGEN][0],1e-37),0.77);
290 /* these are equations 11 from LS83 */
291 if( phycon.te > 4031. )
292 {
293 rotlow = 1.38e-22*sexp(9243./phycon.te)*qn;
294 }
295 else
296 {
297 rotlow = pow(10.,-22.90 - 0.553*x - 1.148*x*x)*qn;
298 }
299
300 CoolHeavy.h2line = 0.;
301 if( rotlow > 0. )
302 CoolHeavy.h2line += hmi.H2_total*rothi/(1. + rothi/rotlow);
303 /* \todo 1 add this from LS83 or (better yet) update to another form. See Galli & Palla 1998, A5-7. */
304 //if( viblow > 0. )
305 // CoolHeavy.h2line += hmi.H2_total*vibhi/(1. + vibhi/viblow);
306 }
307
308 {
309 enum {DEBUG_LOC=false};
310 if( DEBUG_LOC && nzone>187&& iteration > 1/**/)
311 {
312 fprintf(ioQQQ,"h2coolbug\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
313 phycon.te,
314 CoolHeavy.h2line,
315 hmi.H2_total,
316 findspecieslocal("H-")->den,
317 hmi.HMinus_photo_rate,
318 rothi,
319 rotlow );
320 }
321 }
322
323 if( hd.lgEnabled )
324 {
325 CoolHeavy.HD = 0.;
326 }
327 else
328 {
329 /* >>chng 02 mar 07, add DH cooling using rates (eqn 6) from
330 * >>refer HD cooling Puy, D., Grenacher, L, & Jetzer, P., 1999, A&A, 345, 723 */
331 factor = sexp(128.6/phycon.te);
332 CoolHeavy.HD = 2.66e-21 * hydro.D2H_ratio * POW2((double)hmi.H2_total) * phycon.sqrte *
333 factor/(1416.+phycon.sqrte*hmi.H2_total * (1. + 3.*factor));
334 }
335 }
336
337 fixit(); // test and enable this by default
338#if 0
339 double chemical_heating = mole.chem_heat();
340 thermal.heating[0][29] = MAX2(0.,chemical_heating);
341 /* add to cooling if net heating is negative */
342 CoolAdd("Chem",0,MAX2(0.,-chemical_heating));
343#endif
344
345 /* cooling due to charge transfer ionization / recombination */
346 CoolAdd("CT C" , 0. , thermal.char_tran_cool );
347
348 /* H- FB; H + e -> H- + hnu */
349 /* H- FF is in with H ff */
350 CoolAdd("H-fb",0,hmi.hmicol);
351
352 /* >>chng 96 nov 15, fac of 2 in deriv to help convergence in very dense
353 * models where H- is important, this takes change in eden into
354 * partial account */
355 thermal.dCooldT += 2.*hmi.hmicol*phycon.teinv;
356 if( PRT_DERIV )
357 fprintf(ioQQQ,"DEBUG dCdT 2 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
358
359 CoolAdd("H2ln",0,CoolHeavy.h2line);
360 /* >>chng 00 oct 21, added coef of 3.5, sign had been wrong */
361 /*thermal.dCooldT += CoolHeavy.h2line*phycon.teinv;*/
362 /* >>chng 03 mar 17, change 3.5 to 15 as per behavior in primal.in */
363 /*thermal.dCooldT += 3.5*CoolHeavy.h2line*phycon.teinv;*/
364 /* >>chng 03 may 18, from 15 to 30 as per behavior in primal.in - overshoots happen */
365 /*thermal.dCooldT += 15.*CoolHeavy.h2line*phycon.teinv;*/
366 /*>>chng 03 oct 03, from 30 to 3, no more overshoots in primalin */
367 /*thermal.dCooldT += 30.*CoolHeavy.h2line*phycon.teinv;*/
368 thermal.dCooldT += 3.0*CoolHeavy.h2line*phycon.teinv;
369
370 {
371 /* problems with H2 cooling */
372 enum {DEBUG_LOC=false};
373 if( DEBUG_LOC /*&& nzone>300 && iteration > 1*/)
374 {
375 fprintf(ioQQQ,"CoolEvaluate debuggg\t%.2f\t%.5e\t%.5e\t%.5e\t%.5e\t%.5e\n",
376 fnzone,
377 phycon.te,
378 hmi.H2_total ,
379 CoolHeavy.h2line,
380 findspecieslocal("H-")->den ,
381 dense.eden);
382 }
383 }
384
385 CoolAdd("HDro",0,CoolHeavy.HD);
386 thermal.dCooldT += CoolHeavy.HD*phycon.teinv;
387
388 CoolAdd("H2+ ",0,CoolHeavy.H2PlsCool);
389 thermal.dCooldT += CoolHeavy.H2PlsCool*phycon.teinv;
390
391 /* heating due to three-body, will be incremented in iso_cool*/
392 thermal.heating[0][3] = 0.;
393 /* heating due to hydrogen lines */
394 thermal.heating[0][23] = 0.;
395 /* heating due to photoionization of all excited states of hydrogen species */
396 thermal.heating[0][1] = 0.;
397
398 /* isoelectronic species cooling, mainly lines, and level ionization */
399 for( long int ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
400 {
401 for( long int nelem=ipISO; nelem < LIMELM; nelem++ )
402 {
403 /* must always call iso_cool since we must zero those variables
404 * that would have been set had the species been present */
405 iso_cool( ipISO , nelem );
406 }
407 }
408
409 /* >>chng 02 jun 18, don't reevaluate needlessly */
410 /* >>chng 03 nov 28, even faster - special logic for when ff is pretty
411 * negligible - eval of ff is pretty slow */
412 /* >>chng 04 feb 19, must not test on temp not changing, since ionization
413 * can change at constant temperature
414 * >>chng 04 sep 14, above introduced bug since brems never reevaluated
415 * now test is zone or temp has changed */
416 if( fabs(CoolHeavy.brems_cool_net)/SDIV(thermal.ctot) > conv.HeatCoolRelErrorAllowed/5. ||
417 conv.lgSearch ||
418 !fp_equal(phycon.te,TeUsedBrems) ||
419 nzone != nzoneUsedBrems )
420 {
421 double BremsThisEnergy;
422 /*double OpacityThisIon;*/
423 long int limit;
424 /* free-free free free brems emission for all ions */
425
426 TeUsedBrems = phycon.te;
427 nzoneUsedBrems = nzone;
428 /* highest frequency where we have non-zero Boltzmann factors */
429 limit = MIN2( rfield.ipMaxBolt , rfield.nflux );
430
431 CoolHeavy.brems_cool_hminus = 0.;
432 CoolHeavy.brems_cool_h = 0.;
433 CoolHeavy.brems_cool_metals = 0.;
434 CoolHeavy.brems_cool_he = 0.;
435 CoolHeavy.brems_heat_total = 0.;
436
437 {
438 double bhfac, bhMinusfac;
439 realnum sumion[LIMELM+1];
440 long int ion_lo , ion_hi;
441
442 ASSERT(rfield.ipEnergyBremsThin < rfield.nupper);
443 ASSERT(limit < rfield.nupper);
444
445 /* ipEnergyBremsThin is index to energy where gas becomes optically thin to brems,
446 * so this loop is over optically thin frequencies
447 * do not include optically thick part as net emission since self absorbed */
448
449 /* do hydrogen first, before main loop since want to break out as separate
450 * coolant, and what to add on H- brems */
451 CoolHeavy.brems_cool_h = 0.;
452 CoolHeavy.brems_cool_hminus = 0.;
453 /* this is all done in opacity_addtotal - why do here too? */
454 for( long int i=rfield.ipEnergyBremsThin; i < limit; i++ )
455 {
456 long int ion = 1;
457
458 /* in all following CoolHeavy.lgFreeOn is flag set with 'no free-free' to
459 * turn off brems heating and cooling */
460 BremsThisEnergy = rfield.gff[ion][i]*rfield.widflx[i]*rfield.ContBoltz[i];
461 /*ASSERT( BremsThisEnergy >= 0. );*/
462 CoolHeavy.brems_cool_h += BremsThisEnergy;
463
464 /* for case of hydrogen, add H- brems - OpacStack contains the ratio
465 * of the H- to H brems cross section - multiply by this and H(1s) population */
466 CoolHeavy.brems_cool_hminus += BremsThisEnergy * opac.OpacStack[i-1+opac.iphmra];
467 }
468 bhfac = (dense.xIonDense[ipHYDROGEN][1]+findspecieslocal("H2+")->den+findspecieslocal("H3+")->den)*CoolHeavy.lgFreeOn* dense.eden*1.032e-11/phycon.sqrte*EN1RYD;
469 bhMinusfac = iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop()*CoolHeavy.lgFreeOn* dense.eden*1.032e-11/phycon.sqrte*EN1RYD;
470 CoolHeavy.brems_cool_h *= bhfac;
471 CoolHeavy.brems_cool_hminus *= bhMinusfac;
472
473 /* now do helium, both He+ and He++ */
474 CoolHeavy.brems_cool_he = 0.;
475 for( long int i=rfield.ipEnergyBremsThin; i < limit; i++ )
476 {
477 long int nelem = ipHELIUM;
478 /* eff. charge is ion, so first rfield.gff argument must be "ion". */
479 BremsThisEnergy =
480 (dense.xIonDense[nelem][1]*rfield.gff[1][i] + 4.*dense.xIonDense[nelem][2]*rfield.gff[2][i])*
481 rfield.widflx[i]*rfield.ContBoltz[i];
482 CoolHeavy.brems_cool_he += BremsThisEnergy;
483 }
484 CoolHeavy.brems_cool_he *=
485 CoolHeavy.lgFreeOn * dense.eden*1.032e-11/phycon.sqrte*EN1RYD;
486
487 /* >>chng 05 jul 13, rewrite this for speed */
488 /* gaunt factors depend only on photon energy and ion charge, so do
489 * sum of ions here before entering into loop over photon energy */
490 sumion[0] = 0.;
491 for( long int ion=1; ion<=LIMELM; ++ion )
492 {
493 sumion[ion] = 0.;
494 for( long int nelem=ipLITHIUM; nelem < LIMELM; ++nelem )
495 {
496 if( dense.lgElmtOn[nelem] && ion<=nelem+1 )
497 {
498 sumion[ion] += dense.xIonDense[nelem][ion];
499 }
500 }
501 /* now include the charge, density, and temperature */
502 sumion[ion] *= POW2((realnum)ion);
503 }
504
505 /* add molecular ions */
506 for( long ipMol = 0; ipMol<mole_global.num_calc; ipMol++ )
507 {
508 ASSERT( (mole_global.list[ipMol]->n_nuclei() != 1) ==
509 (!mole_global.list[ipMol]->isMonatomic()));
510
511 if( !mole_global.list[ipMol]->isMonatomic() &&
512 mole_global.list[ipMol]->parentLabel.empty() &&
513 mole_global.list[ipMol]->charge > 0 &&
514 mole_global.list[ipMol]->label != "H2+" &&
515 mole_global.list[ipMol]->label != "H3+" )
516 {
517 ASSERT( mole_global.list[ipMol]->charge < LIMELM+1 );
518 sumion[mole_global.list[ipMol]->charge] += (realnum)mole.species[ipMol].den * POW2((realnum)mole_global.list[ipMol]->charge);
519 }
520 }
521
522 /* now find lowest and highest ion we need to consider - following loop is over
523 * full continuum and eats time
524 * >>chng 05 oct 19, bounds check had been on ion, rather than ion_lo and ion_hi, so
525 * array bounds were exceeded */
526 ion_lo = 1;
527 while( sumion[ion_lo]==0 && ion_lo<LIMELM-1 )
528 ++ion_lo;
529 ion_hi = LIMELM;
530 while( sumion[ion_hi]==0 && ion_hi>0 )
531 --ion_hi;
532
533 /* heavy elements */
534 CoolHeavy.brems_cool_metals = 0.;
535 CoolHeavy.brems_heat_total = 0.;
536 for( long int i=rfield.ipEnergyBremsThin; i < limit; i++ )
537 {
538 BremsThisEnergy = 0.;
539 for(long int ion=ion_lo; ion<=ion_hi; ++ion )
540 BremsThisEnergy += sumion[ion]*rfield.gff[ion][i];
541
542 CoolHeavy.brems_cool_metals += BremsThisEnergy*rfield.widflx[i]*rfield.ContBoltz[i];
543 /* the total heating due to bremsstrahlung */
544 CoolHeavy.brems_heat_total += opac.FreeFreeOpacity[i]*rfield.flux[0][i]*rfield.anu[i]*EN1RYD*CoolHeavy.lgFreeOn;
545 }
546 CoolHeavy.brems_cool_metals *=
547 CoolHeavy.lgFreeOn * dense.eden*1.032e-11/phycon.sqrte*EN1RYD;
548
549 {
550 enum {DEBUG_LOC=false};
551 if( DEBUG_LOC && nzone>60 /*&& iteration > 1*/)
552 {
553 double sumfield = 0., sumtot=0., sum1=0., sum2=0.;
554 for( long int i=rfield.ipEnergyBremsThin; i<limit; i++ )
555 {
556 sumtot += opac.FreeFreeOpacity[i]*rfield.flux[0][i]*rfield.anu[i];
557 sumfield += rfield.flux[0][i]*rfield.anu[i];
558 sum1 += opac.FreeFreeOpacity[i]*rfield.flux[0][i]*rfield.anu[i];
559 sum2 += opac.FreeFreeOpacity[i]*rfield.flux[0][i];
560 }
561 fprintf(ioQQQ,"DEBUG brems heat\t%.2f\t%.3e\t%.3e\t%.3e\t%e\t%.3e\t%.3e\n",
562 fnzone,
563 CoolHeavy.brems_heat_total,
564 sumtot/SDIV(sumfield) ,
565 sum1/SDIV(sum2),
566 phycon.te ,
567 rfield.gff[1][1218],
568 opac.FreeFreeOpacity[1218]);
569 }
570 }
571 }
572 }
573
574 /* these two terms are both large, nearly canceling, near lte */
575 CoolHeavy.brems_cool_net =
576 CoolHeavy.brems_cool_h +
577 CoolHeavy.brems_cool_he +
578 CoolHeavy.brems_cool_hminus +
579 CoolHeavy.brems_cool_metals -
580 CoolHeavy.brems_heat_total;
581 /*fprintf(ioQQQ,"DEBUG brems\t%.2f\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\t%.4e\n",
582 fnzone,
583 phycon.te,
584 CoolHeavy.brems_cool_net,
585 CoolHeavy.brems_cool_h ,
586 CoolHeavy.brems_cool_he ,
587 CoolHeavy.brems_cool_hminus,
588 CoolHeavy.brems_cool_metals ,
589 CoolHeavy.brems_heat_total);*/
590
591 /* net free free brems cooling, count as cooling if positive */
592 CoolAdd( "FF c" , 0, MAX2(0.,CoolHeavy.brems_cool_net) );
593
594 /* now stuff into heating array if negative */
595 thermal.heating[0][11] = MAX2(0.,-CoolHeavy.brems_cool_net);
596
597 /* >>chng 96 oct 30, from HFFNet to just FreeFreeCool,
598 * since HeatSum picks up CoolHeavy.brems_heat_total */
599 thermal.dCooldT += CoolHeavy.brems_cool_h*thermal.halfte;
600 if( PRT_DERIV )
601 fprintf(ioQQQ,"DEBUG dCdT 3 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
602
603 /* >>chng 02 jun 21, net cooling already includes this */
604 /* end of brems cooling */
605
606 /* heavy element recombination cooling, do not count hydrogenic since
607 * already done above, also helium singlets have been done */
608 /* >>chng 02 jul 21, put in charge dependent rec term */
609 CoolHeavy.heavfb = 0.;
610 for( long int nelem=ipLITHIUM; nelem < LIMELM; nelem++ )
611 {
612 if( dense.lgElmtOn[nelem] )
613 {
614 /* note that detailed iso seq atoms are done in iso_cool */
615 long limit_lo = MAX2( 1 , dense.IonLow[nelem] );
616 long limit_hi = MIN2( nelem-NISO+1, dense.IonHigh[nelem] );
617 for( long int ion=limit_lo; ion<=limit_hi; ++ion )
618 {
619 /* factor of 0.9 is roughly correct for nebular conditions, see
620 * >>refer H rec cooling LaMothe, J., & Ferland, G.J., 2001, PASP, 113, 165 */
621 /* note that ionbal.RR_rate_coef_used is the rate coef, cm3 s-1, needs eden */
622 /* >>chng 02 nov 07, move rec arrays around, this now has ONLY rad rec,
623 * previously had included grain rec and three body */
624 /* recombination cooling for iso-seq atoms are done in iso_cool */
625 double one = dense.xIonDense[nelem][ion] * ionbal.RR_rate_coef_used[nelem][ion-1]*
626 dense.eden * phycon.te * BOLTZMANN;
627 /*fprintf(ioQQQ,"debugggfb\t%li\t%li\t%.3e\t%.3e\t%.3e\n", nelem, ion, one,
628 dense.xIonDense[nelem][ion] , ionbal.RR_rate_coef_used[nelem][ion]);*/
629 CoolHeavy.heavfb += one;
630 }
631 }
632 }
633
634 /*fprintf(ioQQQ,"debuggg hvFB\t%i\t%.2f\t%.2e\t%.2e\n",iteration, fnzone,CoolHeavy.heavfb, dense.eden);*/
635
636 CoolAdd("hvFB",0,CoolHeavy.heavfb);
637 thermal.dCooldT += CoolHeavy.heavfb*.113*phycon.teinv;
638 if( PRT_DERIV )
639 fprintf(ioQQQ,"DEBUG dCdT 4 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
640
641 /* electron-electron brems, approx form from
642 * >>refer ee brems Stepney and Guilbert, MNRAS 204, 1269 (1983)
643 * ok for T<10**9 */
644 CoolHeavy.eebrm = POW2(dense.eden*phycon.te*1.84e-21);
645
646 /* >>chng 97 mar 12, added deriv */
647 thermal.dCooldT += CoolHeavy.eebrm*thermal.halfte;
648 CoolAdd("eeff",0,CoolHeavy.eebrm);
649
650 /* add advective heating and cooling */
651 /* this is cooling due to loss of matter from this region */
652 CoolAdd("adve",0,dynamics.Cool() );
653 /* >>chng02 dec 04, rm factor of 8 in front of dCooldT */
654 thermal.dCooldT += dynamics.dCooldT();
655 /* local heating due to matter moving into this location */
656 thermal.heating[1][5] = dynamics.Heat();
657 thermal.dHeatdT += dynamics.dHeatdT;
658
659 /* total Compton cooling */
660 CoolHeavy.tccool = rfield.cmcool*phycon.te;
661 CoolAdd("Comp",0,CoolHeavy.tccool);
662 thermal.dCooldT += rfield.cmcool;
663
664 /* option for "extra" cooling, expressed as power-law in temperature, these
665 * are set with the CEXTRA command */
666 if( thermal.lgCExtraOn )
667 {
668 CoolHeavy.cextxx =
669 (realnum)(thermal.CoolExtra*pow(phycon.te/1e4,(double)thermal.cextpw));
670 }
671 else
672 {
673 CoolHeavy.cextxx = 0.;
674 }
675 CoolAdd("Extr",0,CoolHeavy.cextxx);
676
677 realnum dDensityDT;
678
679 /* cooling due to wind expansion, only for winds expansion cooling */
680 if( wind.lgBallistic() )
681 {
682 dDensityDT = -(realnum)(wind.AccelTotalOutward/wind.windv + 2.*wind.windv/
683 radius.Radius);
684 CoolHeavy.expans = -2.5*pressure.PresGasCurr*dDensityDT;
685 }
686 else if( dynamics.lgTimeDependentStatic &&
687 iteration > dynamics.n_initial_relax)
688 {
689 realnum dens = scalingDensity();
690 dDensityDT =
691 (realnum)((dens-dynamics.Upstream_density)/
692 (dynamics.timestep*0.5*(dens+dynamics.Upstream_density)));
693 // pdV work term
694 CoolHeavy.expans = -pressure.PresGasCurr*dDensityDT;
695 }
696 else
697 {
698 dDensityDT = 0.;
699 CoolHeavy.expans = 0.;
700 }
701 CoolAdd("Expn",0,CoolHeavy.expans);
702 thermal.dCooldT += CoolHeavy.expans/phycon.te;
703
704 /* cyclotron cooling */
705 /* coef is 4/3 /8pi /c * vtherm(elec) */
706 CoolHeavy.cyntrn = 4.5433e-25f*magnetic.pressure*PI8*dense.eden*phycon.te;
707 CoolAdd("Cycl",0,CoolHeavy.cyntrn);
708 thermal.dCooldT += CoolHeavy.cyntrn/phycon.te;
709
710 /* heavy element collisional ionization
711 * derivative should be zero since increased coll ion rate
712 * decreases neutral fraction by proportional amount */
713 CoolAdd("Hvin",0,CoolHeavy.colmet);
714
715 /* evaluate H 21 cm spin changing collisions */
716 coolnum = thermal.ncltot;
717 if( !fp_equal(phycon.te,TeEvalCS_21cm) )
718 {
719 {
720 /* this prints table of rates at points given in original data paper */
721 enum {DEBUG_LOC=false};
722 if( DEBUG_LOC )
723 {
724# define N21CM_TE 16
725 int n;
726 double teval[N21CM_TE]={2.,5.,10.,20.,50.,100.,200.,500.,1000.,
727 2000.,3000.,5000.,7000.,10000.,15000.,20000.};
728 for( n = 0; n<N21CM_TE; ++n )
729 {
730 fprintf(
731 ioQQQ,"DEBUG 21 cm deex Te=\t%.2e\tH0=\t%.2e\tp=\t%.2e\te=\t%.2e\n",
732 teval[n],
733 H21cm_H_atom( teval[n] ),
734 H21cm_proton( teval[n] ),
735 H21cm_electron( teval[n] ) );
736 }
738# undef N21CM_TE
739 }
740 }
741 /*only evaluate T dependent part when Te changes, but update
742 * density part below since densities may constantly change */
743 atomic_rate_21cm = H21cm_H_atom( phycon.te );
744 proton_rate_21cm = H21cm_proton( phycon.te );
745 electron_rate_21cm = H21cm_electron( phycon.te );
746 TeEvalCS_21cm = phycon.te;
747 }
748 /* H 21 cm emission/population,
749 * cs will be sum of e cs and H cs converted from rate */
750 cs = (electron_rate_21cm * dense.eden +
751 atomic_rate_21cm * dense.xIonDense[ipHYDROGEN][0] +
752 proton_rate_21cm * dense.xIonDense[ipHYDROGEN][1] ) *
753 3./ dense.cdsqte;
754 PutCS( cs , HFLines[0] );
755
756 /* fine structure lines */
757 if( !fp_equal(phycon.te,TeEvalCS) )
758 {
759 /* H 21 cm done above, now loop over remaining lines to get collision strengths */
760 for( long int i=1; i < nHFLines; i++ )
761 {
762 cs = HyperfineCS( i );
763 /* now generate the collision strength and put into the line array */
764 PutCS( cs , HFLines[i] );
765 }
766 TeEvalCS = phycon.te;
767 }
768
769 /* do level pops for H 21 cm which is a special case since Lya pumping in included */
770 RT_line_one( HFLines[0], true,0.f, GetDopplerWidth(dense.AtomicWeight[(*HFLines[0].Hi()).nelem()-1]) );
771 H21_cm_pops();
772 if( PRT_DERIV )
773 fprintf(ioQQQ,"DEBUG dCdT 5 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
774
775 /* find total cooling due to hyperfine structure lines */
776 hyperfine.cooling_total = HFLines[0].Coll().cool();
777
778 /* now do level pops for all except 21 cm */
779 for( long int i=1; i < nHFLines; i++ )
780 {
781 /* remember current gas-phase abundance of this isotope */
782 realnum save = dense.xIonDense[(*HFLines[i].Hi()).nelem()-1][(*HFLines[i].Hi()).IonStg()-1];
783
784 /* bail if no abundance */
785 if( save<=0. )
786 continue;
787
788 RT_line_one( HFLines[i], true,0.f, GetDopplerWidth(dense.AtomicWeight[(*HFLines[i].Hi()).nelem()-1]) );
789
790 /* set gas-phase abundance to total times isotope ratio */
791 dense.xIonDense[(*HFLines[i].Hi()).nelem()-1][(*HFLines[i].Hi()).IonStg()-1] *=
792 hyperfine.HFLabundance[i];
793
794 /* use the collision strength generated above and find pops and cooling */
795 atom_level2( HFLines[i] );
796
797 /* put the correct gas-phase abundance back in the array */
798 dense.xIonDense[(*HFLines[i].Hi()).nelem()-1][(*HFLines[i].Hi()).IonStg()-1] = save;
799
800 /* find total cooling due to hyperfine structure lines */
801 hyperfine.cooling_total += HFLines[i].Coll().cool();
802 }
803 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
804 thermal.elementcool[ipHYDROGEN] += thermal.cooling[coolcal];
805
806 if( PRT_DERIV )
807 fprintf(ioQQQ,"DEBUG dCdT 6 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
808
809 double xIonDenseSave[LIMELM][LIMELM+1];
810 if( atmdat.lgChiantiOn ||atmdat.lgStoutOn)
811 {
812 for( int nelem=0; nelem < LIMELM; nelem++ )
813 {
814 for( int ion=0; ion<=nelem+1; ++ion )
815 {
816 xIonDenseSave[nelem][ion] = dense.xIonDense[nelem][ion];
817 // zero abundance of species if we are using Chianti for this ion
818 if( dense.lgIonChiantiOn[nelem][ion] || dense.lgIonStoutOn[nelem][ion] )
819 dense.xIonDense[nelem][ion] = 0.;
820 }
821 }
822 }
823
824 /* Carbon cooling */
825 coolnum = thermal.ncltot;
826 CoolCarb();
827 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
828 thermal.elementcool[ipCARBON] += thermal.cooling[coolcal];
829 if( PRT_DERIV )
830 fprintf(ioQQQ,"DEBUG dCdT C %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
831
832 /* Nitrogen cooling */
833 coolnum = thermal.ncltot;
834 CoolNitr();
835 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
836 thermal.elementcool[ipNITROGEN] += thermal.cooling[coolcal];
837 if( PRT_DERIV )
838 fprintf(ioQQQ,"DEBUG dCdT N %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
839
840 /* Oxygen cooling */
841 coolnum = thermal.ncltot;
842 CoolOxyg();
843 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
844 thermal.elementcool[ipOXYGEN] += thermal.cooling[coolcal];
845 if( PRT_DERIV )
846 fprintf(ioQQQ,"DEBUG dCdT 7 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
847
848 /* Neon cooling */
849 coolnum = thermal.ncltot;
850 CoolNeon();
851 if( PRT_DERIV )
852 fprintf(ioQQQ,"DEBUG dCdT Ne %.3e dHdT %.3e\n",thermal.dCooldT
853 , thermal.dHeatdT);
854 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
855 thermal.elementcool[ipNEON] += thermal.cooling[coolcal];
856
857 /* Magnesium cooling */
858 coolnum = thermal.ncltot;
859 CoolMagn();
860 if( PRT_DERIV )
861 fprintf(ioQQQ,"DEBUG dCdT 8 %.3e dHdT %.3e\n",thermal.dCooldT
862 , thermal.dHeatdT);
863 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
864 thermal.elementcool[ipMAGNESIUM] += thermal.cooling[coolcal];
865
866 /* Sodium cooling */
867 coolnum = thermal.ncltot;
868 CoolSodi();
869 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
870 thermal.elementcool[ipSODIUM] += thermal.cooling[coolcal];
871 if( PRT_DERIV )
872 fprintf(ioQQQ,"DEBUG dCdT Na %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
873
874 /* Aluminum cooling */
875 coolnum = thermal.ncltot;
876 CoolAlum();
877 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
878 thermal.elementcool[ipALUMINIUM] += thermal.cooling[coolcal];
879 if( PRT_DERIV )
880 fprintf(ioQQQ,"DEBUG dCdT Al %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
881
882 /* Silicon cooling */
883 coolnum = thermal.ncltot;
884 CoolSili();
885 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
886 thermal.elementcool[ipSILICON] += thermal.cooling[coolcal];
887 if( PRT_DERIV )
888 fprintf(ioQQQ,"DEBUG dCdT 9 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
889
890 /* Phosphorus */
891 coolnum = thermal.ncltot;
892 CoolPhos();
893 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
894 thermal.elementcool[ipPHOSPHORUS] += thermal.cooling[coolcal];
895
896 /* Sulphur cooling */
897 coolnum = thermal.ncltot;
898 CoolSulf();
899 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
900 thermal.elementcool[ipSULPHUR] += thermal.cooling[coolcal];
901
902 /* Chlorine cooling */
903 coolnum = thermal.ncltot;
904 CoolChlo();
905 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
906 thermal.elementcool[ipCHLORINE] += thermal.cooling[coolcal];
907
908 /* Argon cooling */
909 coolnum = thermal.ncltot;
910 CoolArgo();
911 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
912 thermal.elementcool[ipARGON] += thermal.cooling[coolcal];
913 if( PRT_DERIV )
914 fprintf(ioQQQ,"DEBUG dCdT 10 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
915
916 /* Potasium cooling */
917 coolnum = thermal.ncltot;
918 CoolPota();
919 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
920 thermal.elementcool[ipPOTASSIUM] += thermal.cooling[coolcal];
921
922 /* Calcium cooling */
923 coolnum = thermal.ncltot;
924 CoolCalc();
925 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
926 thermal.elementcool[ipCALCIUM] += thermal.cooling[coolcal];
927
928 /* Scandium cooling */
929 coolnum = thermal.ncltot;
930 CoolScan();
931 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
932 thermal.elementcool[ipSCANDIUM] += thermal.cooling[coolcal];
933
934 /* Chromium cooling */
935 coolnum = thermal.ncltot;
936 CoolChro();
937 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
938 thermal.elementcool[ipCHROMIUM] += thermal.cooling[coolcal];
939
940
941 /* Iron cooling */
942 coolnum = thermal.ncltot;
943 CoolIron();
944 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
945 thermal.elementcool[ipIRON] += thermal.cooling[coolcal];
946 if( PRT_DERIV )
947 fprintf(ioQQQ,"DEBUG dCdT 12 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
948
949 /* Cobalt cooling */
950 coolnum = thermal.ncltot;
951 CoolCoba();
952 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
953 thermal.elementcool[ipCOBALT] += thermal.cooling[coolcal];
954
955 /* Nickel cooling */
956 coolnum = thermal.ncltot;
957 CoolNick();
958 for( coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
959 thermal.elementcool[ipNICKEL] += thermal.cooling[coolcal];
960
961 coolnum = thermal.ncltot;
962
963 // reset abundances to original values, may have been set zero to protect against old cloudy lines
964 if( atmdat.lgChiantiOn || atmdat.lgStoutOn)
965 {
966 // this clause, first reset abundances set to zero when Chianti included
967 for( int nelem=0; nelem < LIMELM; nelem++ )
968 {
969 for( int ion=0; ion<=nelem+1; ++ion )
970 {
971 dense.xIonDense[nelem][ion] = xIonDenseSave[nelem][ion];
972 }
973 }
974 }
975
976 /* opacity project lines Dima Verner added with g-bar approximation */
977 CoolDima();
978
979 for( int coolcal = coolnum; coolcal < thermal.ncltot; coolcal++ )
980 thermal.dima += thermal.cooling[coolcal];
981
982 /* do external database lines */
983 dBase_solve();
984
985 /* Print number of levels for each species */
986 {
987 enum {DEBUG_LOC=false};
988 if( DEBUG_LOC )
989 {
990 static bool lgMustPrintHeader=true;
992 {
993 lgMustPrintHeader = false;
994 printf("DEBUG Levels\t%s",dBaseSpecies[0].chLabel );
995 for( long ipSpecies=1; ipSpecies<nSpecies; ipSpecies++ )
996 {
997 printf("\t%s",dBaseSpecies[ipSpecies].chLabel );
998 }
999 printf("\n" );
1000 printf("DEBUG Max\t%li" ,dBaseSpecies[0].numLevels_max );
1001 for( long ipSpecies=1; ipSpecies<nSpecies; ipSpecies++ )
1002 {
1003 printf( "\t%li" ,dBaseSpecies[ipSpecies].numLevels_max );
1004 }
1005 printf("\n");
1006 }
1007 printf("DEBUG Local\t%li" ,dBaseSpecies[0].numLevels_local );
1008 for( long ipSpecies=1; ipSpecies<nSpecies; ipSpecies++ )
1009 {
1010 printf("\t%li" ,dBaseSpecies[ipSpecies].numLevels_local );
1011 }
1012 printf("\n");
1013 }
1014 }
1015
1016 /* now add up all the coolants */
1017 CoolSum(tot);
1018 if( PRT_DERIV )
1019 fprintf(ioQQQ,"DEBUG dCdT 14 %.3e dHdT %.3e\n",thermal.dCooldT , thermal.dHeatdT);
1020
1021 /* negative cooling */
1022 if( *tot <= 0. )
1023 {
1024 fprintf( ioQQQ, " COOLR; cooling is <=0, this is impossible.\n" );
1025 ShowMe();
1027 }
1028
1029 /* bad derivative */
1030 if( thermal.dCooldT == 0. )
1031 {
1032 fprintf( ioQQQ, " COOLR; cooling slope <=0, this is impossible.\n" );
1033 if( *tot > 0. && dense.gas_phase[ipHYDROGEN] < 1e-4 )
1034 {
1035 fprintf( ioQQQ, " Probably due to very low density.\n" );
1036 }
1037 ShowMe();
1039 }
1040
1041 if( trace.lgTrace )
1042 {
1043 fndstr(*tot,thermal.dCooldT);
1044 }
1045
1046 /* lgTSetOn true for constant temperature model */
1047 if( (((((!thermal.lgTemperatureConstant) && *tot < 0.) && called.lgTalk) &&
1048 !conv.lgSearch) && thermal.lgCNegChk) && nzone > 0 )
1049 {
1050 fprintf( ioQQQ,
1051 " NOTE Negative cooling, zone %4ld, =%10.2e coola=%10.2e CHION=%10.2e Te=%10.2e\n",
1052 nzone,
1053 *tot,
1054 iso_sp[ipH_LIKE][ipHYDROGEN].cLya_cool,
1055 iso_sp[ipH_LIKE][ipHYDROGEN].coll_ion,
1056 phycon.te );
1057 fndneg();
1058 }
1059
1060 /* possibility of getting empirical cooling derivative
1061 * normally false, set true with 'set numerical derivatives' command */
1062 if( NumDeriv.lgNumDeriv )
1063 {
1064 if( ((nzone > 2 && nzone == nzSave) && ! fp_equal( oldtemp, phycon.te )) && nhit > 4 )
1065 {
1066 /* hnit is number of tries on this zone - use to stop numerical problems
1067 * do not evaluate numerical deriv until well into solution */
1068 deriv = (oltcool - *tot)/(oldtemp - phycon.te);
1069 thermal.dCooldT = deriv;
1070 }
1071 else
1072 {
1073 deriv = thermal.dCooldT;
1074 }
1075 if( nzone != nzSave )
1076 nhit = 0;
1077
1078 nzSave = nzone;
1079 nhit += 1;
1080 oltcool = *tot;
1081 oldtemp = phycon.te;
1082 }
1083 return;
1084}
1085
1086/* */
1087#ifdef EPS
1088# undef EPS
1089#endif
1090#define EPS 0.01
1091
1092/*fndneg search cooling array to find negative values */
1093STATIC void fndneg(void)
1094{
1095 long int i;
1096 double trig;
1097
1098 DEBUG_ENTRY( "fndneg()" );
1099
1100 trig = fabs(thermal.htot*EPS);
1101 for( i=0; i < thermal.ncltot; i++ )
1102 {
1103 if( thermal.cooling[i] < 0. && fabs(thermal.cooling[i]) > trig )
1104 {
1105 fprintf( ioQQQ, " negative line=%s %.2f fraction of heating=%.3f\n",
1106 thermal.chClntLab[i], thermal.collam[i], thermal.cooling[i]/
1107 thermal.htot );
1108 }
1109
1110 if( thermal.heatnt[i] > trig )
1111 {
1112 fprintf( ioQQQ, " heating line=%s %.2f fraction of heating=%.3f\n",
1113 thermal.chClntLab[i], thermal.collam[i], thermal.heatnt[i]/
1114 thermal.htot );
1115 }
1116 }
1117 return;
1118}
1119
1120/*fndstr search cooling stack to find strongest values */
1121STATIC void fndstr(double tot,
1122 double dc)
1123{
1124 char chStrngLab[NCOLNT_LAB_LEN+1];
1125 long int i;
1126 realnum wl;
1127 double str,
1128 strong;
1129
1130 DEBUG_ENTRY( "fndstr()" );
1131
1132 strong = 0.;
1133 wl = -FLT_MAX;
1134 for( i=0; i < thermal.ncltot; i++ )
1135 {
1136 if( fabs(thermal.cooling[i]) > strong )
1137 {
1138 /* this is the wavelength of the coolant, 0 for a continuum*/
1139 wl = thermal.collam[i];
1140 /* make sure labels are all valid*/
1141 /*>>chng 06 jun 06, bug fix, assert length was ==4, should be <=NCOLNT_LAB_LEN */
1142 ASSERT( strlen( thermal.chClntLab[i] ) <= NCOLNT_LAB_LEN );
1143 strcpy( chStrngLab, thermal.chClntLab[i] );
1144 strong = fabs(thermal.cooling[i]);
1145 }
1146 }
1147
1148 str = strong;
1149
1150 fprintf( ioQQQ,
1151 " fndstr cool: TE=%10.4e Ne=%10.4e C=%10.3e dC/dT=%10.2e ABS(%s %.1f)=%.2e nz=%ld\n",
1152 phycon.te, dense.eden, tot, dc, chStrngLab
1153 , wl, str, nzone );
1154
1155 /* option for extensive printout of lines */
1156 if( trace.lgCoolTr )
1157 {
1158 realnum ratio;
1159
1160 /* flag all significant coolants, first zero out the array */
1161 coolpr(ioQQQ,(char*)thermal.chClntLab[0],1,0.,"ZERO");
1162
1163 /* push all coolants onto the stack */
1164 for( i=0; i < thermal.ncltot; i++ )
1165 {
1166 /* usually positive, although can be neg for coolants that heats,
1167 * only do positive here */
1168 ratio = (realnum)(thermal.cooling[i]/thermal.ctot);
1169 if( ratio >= EPS )
1170 {
1171 /*>>chng 99 jan 27, only cal when ratio is significant */
1172 coolpr(ioQQQ,(char*)thermal.chClntLab[i],thermal.collam[i], ratio,"DOIT");
1173 }
1174 }
1175
1176 /* complete the printout for positive coolants */
1177 coolpr(ioQQQ,"DONE",1,0.,"DONE");
1178
1179 /* now do cooling that was counted as a heat source if significant */
1180 if( thermal.heating[0][22]/thermal.ctot > 0.05 )
1181 {
1182 fprintf( ioQQQ,
1183 " All coolant heat greater than%6.2f%% of the total will be printed.\n",
1184 EPS*100. );
1185
1186 coolpr(ioQQQ,"ZERO",1,0.,"ZERO");
1187 for( i=0; i < thermal.ncltot; i++ )
1188 {
1189 ratio = (realnum)(thermal.heatnt[i]/thermal.ctot);
1190 if( fabs(ratio) >=EPS )
1191 {
1192 coolpr(ioQQQ,(char*)thermal.chClntLab[i],thermal.collam[i],
1193 ratio,"DOIT");
1194 }
1195 }
1196 coolpr(ioQQQ,"DONE",1,0.,"DONE");
1197 }
1198 }
1199 return;
1200}
t_atmdat atmdat
Definition atmdat.cpp:6
double H21cm_electron(double temp)
double H21cm_H_atom(double temp)
double H21cm_proton(double temp)
double HyperfineCS(long i)
void H21_cm_pops(void)
void atom_level2(const TransitionProxy &t)
t_called called
Definition called.cpp:5
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
long int iteration
Definition cddefines.cpp:16
double fnzone
Definition cddefines.cpp:15
#define ASSERT(exp)
Definition cddefines.h:578
sys_float sexp(sys_float x)
Definition service.cpp:914
const int ipSILICON
Definition cddefines.h:318
#define STATIC
Definition cddefines.h:97
const int ipCHROMIUM
Definition cddefines.h:328
const int ipNITROGEN
Definition cddefines.h:311
const int ipIRON
Definition cddefines.h:330
const int ipCARBON
Definition cddefines.h:310
const int ipSCANDIUM
Definition cddefines.h:325
const int ipALUMINIUM
Definition cddefines.h:317
#define MIN2
Definition cddefines.h:761
const int ipOXYGEN
Definition cddefines.h:312
const int LIMELM
Definition cddefines.h:258
const int ipPHOSPHORUS
Definition cddefines.h:319
#define POW2
Definition cddefines.h:929
const int ipLITHIUM
Definition cddefines.h:307
const int ipMAGNESIUM
Definition cddefines.h:316
#define EXIT_FAILURE
Definition cddefines.h:140
const int ipSODIUM
Definition cddefines.h:315
const int ipCOBALT
Definition cddefines.h:331
#define cdEXIT(FAIL)
Definition cddefines.h:434
const int NISO
Definition cddefines.h:261
const int ipHELIUM
Definition cddefines.h:306
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
const int ipNICKEL
Definition cddefines.h:332
const int ipSULPHUR
Definition cddefines.h:320
const int ipCALCIUM
Definition cddefines.h:324
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
const int ipPOTASSIUM
Definition cddefines.h:323
NORETURN void TotalInsanity(void)
Definition service.cpp:886
sys_float SDIV(sys_float x)
Definition cddefines.h:952
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
const int ipNEON
Definition cddefines.h:314
const int ipCHLORINE
Definition cddefines.h:321
const int ipARGON
Definition cddefines.h:322
void ShowMe(void)
Definition service.cpp:181
void fixit(void)
Definition service.cpp:991
double den
Definition mole.h:358
t_conv conv
Definition conv.cpp:5
void CoolAlum(void)
Definition cool_alum.cpp:15
void CoolArgo(void)
Definition cool_argo.cpp:15
void CoolCalc(void)
Definition cool_calc.cpp:19
void CoolCarb(void)
Definition cool_carb.cpp:22
void CoolChlo(void)
Definition cool_chlo.cpp:14
void CoolChro(void)
Definition cool_chro.cpp:14
void CoolCoba(void)
Definition cool_coba.cpp:10
void CoolDima(void)
Definition cool_dima.cpp:25
void CoolSum(double *total)
Definition cool_etc.cpp:76
void CoolAdd(const char *chLabel, realnum lambda, double cool)
Definition cool_etc.cpp:13
void CoolZero(void)
Definition cool_etc.cpp:50
void CoolEvaluate(double *tot)
Definition cool_eval.cpp:45
#define N21CM_TE
static const bool PRT_DERIV
Definition cool_eval.cpp:43
#define EPS
STATIC void fndneg(void)
STATIC void fndstr(double tot, double dc)
void CoolIron(void)
void CoolMagn(void)
Definition cool_magn.cpp:14
void CoolNeon(void)
Definition cool_neon.cpp:17
void CoolNick(void)
Definition cool_nick.cpp:12
void CoolNitr(void)
void CoolOxyg(void)
Definition cool_oxyg.cpp:24
void CoolPhos(void)
Definition cool_phos.cpp:13
void CoolPota(void)
Definition cool_pota.cpp:11
void coolpr(FILE *io, const char *chLabel, realnum lambda, double ratio, const char *chJOB)
Definition cool_pr.cpp:9
void CoolScan(void)
Definition cool_scan.cpp:12
void CoolSili(void)
Definition cool_sili.cpp:16
void CoolSodi(void)
Definition cool_sodi.cpp:13
void CoolSulf(void)
Definition cool_sulf.cpp:20
t_CoolHeavy CoolHeavy
Definition coolheavy.cpp:5
const realnum SMALLFLOAT
Definition cpu.h:191
t_dense dense
Definition dense.cpp:24
realnum scalingDensity(void)
Definition dense.cpp:378
realnum GetDopplerWidth(realnum massAMU)
t_dynamics dynamics
Definition dynamics.cpp:44
GrainVar gv
Definition grainvar.cpp:5
diatomics hd("hd", 4100., &hmi.HD_total, Yan_H2_CS)
vector< diatomics * > diatoms
Definition h2.cpp:8
diatomics h2("h2", 4100., &hmi.H2_total, Yan_H2_CS)
vector< diatomics * >::iterator diatom_iter
Definition h2.h:13
t_hmi hmi
Definition hmi.cpp:5
t_hydro hydro
Definition hydrogenic.cpp:5
t_hyperfine hyperfine
Definition hyperfine.cpp:5
t_ionbal ionbal
Definition ionbal.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipH1s
Definition iso.h:27
void iso_cool(long ipISO, long nelem)
const int ipH_LIKE
Definition iso.h:62
t_magnetic magnetic
Definition magnetic.cpp:17
t_mole_global mole_global
Definition mole.cpp:6
t_mole_local mole
Definition mole.cpp:7
molezone * findspecieslocal(const char buf[])
t_NumDeriv NumDeriv
Definition numderiv.cpp:5
t_opac opac
Definition opacity.cpp:5
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 PI8
Definition physconst.h:38
t_pressure pressure
Definition pressure.cpp:5
t_radius radius
Definition radius.cpp:5
t_rfield rfield
Definition rfield.cpp:8
void RT_line_one(const TransitionProxy &t, bool lgShield_this_zone, realnum pestrk, realnum DopplerWidth)
t_save save
Definition save.cpp:5
static bool lgMustPrintHeader
void dBase_solve(void)
Definition species2.cpp:33
long int nSpecies
Definition taulines.cpp:21
TransitionList HFLines("HFLines", &AnonStates)
long int nHFLines
Definition taulines.cpp:31
species * dBaseSpecies
Definition taulines.cpp:14
void TempChange(double TempNew, bool lgForceUpdate)
t_thermal thermal
Definition thermal.cpp:5
#define NCOLNT_LAB_LEN
Definition thermal.h:91
t_trace trace
Definition trace.cpp:5
void PutCS(double cs, const TransitionProxy &t)
Wind wind
Definition wind.cpp:5