cloudy trunk
Loading...
Searching...
No Matches
atom_level3.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/*atom_level3 compute three level atom, 10, 21, and 20 are line */
4#include "cddefines.h"
5#include "phycon.h"
6#include "physconst.h"
7#include "dense.h"
8#include "transition.h"
9#include "thermal.h"
10#include "cooling.h"
11/* NB - use the nLev3Fail failures mode indicators when debugging this routine */
12#include "atoms.h"
13#include "rfield.h"
14
16 const TransitionProxy& t21,
17 const TransitionProxy& t20)
18{
19 char chLab[5],
20 chLab10[11];
21 double AbunxIon,
22 a,
23 a10,
24 a20,
25 a21,
26 b,
27 beta,
28 bolt01,
29 bolt02,
30 bolt12,
31 c,
32 ener10,
33 ener20,
34 ener21,
35 g0,
36 g010,
37 g020,
38 g1,
39 g110,
40 g121,
41 g2,
42 g220,
43 g221,
44 o10,
45 o20,
46 o21,
47 p0,
48 p1,
49 p2,
50 pump01,
51 pump02,
52 pump12;
53
54 int nLev3Fail;
55
56 double TotCool,
57 TotHeat,
58 TotInten,
59 alpha,
60 alpha1,
61 alpha2,
62 c01,
63 c02,
64 c10,
65 c12,
66 c20,
67 c21,
68 cnet01,
69 cnet02,
70 cnet12,
71 cool01,
72 cool02,
73 cool12,
74 heat10,
75 heat20,
76 heat21,
77 hnet01,
78 hnet02,
79 hnet12,
80 pump10,
81 pump20,
82 pump21,
83 r01,
84 r02,
85 r10,
86 r12,
87 r20,
88 r21,
89 temp01,
90 temp02,
91 temp12;
92
93 DEBUG_ENTRY( "atom_level3()" );
94
95 /* compute three level atom, 10, 21, and 20 are line
96 * arrays for 10, 21, and 20 transitions.
97 * one can be a dummy */
98 /* >>chng 96 dec 06, to double precision due to round off problems below */
99
100 /* generalized three level atom for any ion
101 * sum of three levels normalized to total abundance
102 *
103 * stat weights of all three lines
104 * sanity check will confirm ok */
105 g010 = (*t10.Lo()).g();
106 g110 = (*t10.Hi()).g();
107
108 g121 = (*t21.Lo()).g();
109 g221 = (*t21.Hi()).g();
110
111 g020 = (*t20.Lo()).g();
112 g220 = (*t20.Hi()).g();
113
114 /* these are statistical weights */
115 if( g010 > 0. )
116 {
117 g0 = g010;
118 }
119
120 else if( g020 > 0. )
121 {
122 g0 = g020;
123 }
124
125 else
126 {
127 g0 = -1.;
128 strcpy( chLab10, chLineLbl(t10) );
129 fprintf( ioQQQ, " PROBLEM atom_level3: insane stat weights g0 :%10.2e%10.2e %10.10s\n",
130 g010, g020, chLab10 );
132 }
133
134 if( g110 > 0. )
135 {
136 g1 = g110;
137 }
138
139 else if( g121 > 0. )
140 {
141 g1 = g121;
142 }
143
144 else
145 {
146 g1 = -1.;
147 strcpy( chLab10, chLineLbl(t10) );
148 fprintf( ioQQQ, " PROBLEM atom_level3: insane stat weights g1 :%10.2e%10.2e %10.10s\n",
149 g010, g020, chLab );
151 }
152
153 if( g220 > 0. )
154 {
155 g2 = g220;
156 }
157 else if( g221 > 0. )
158 {
159 g2 = g221;
160 }
161
162 else
163 {
164 g2 = -1.;
165 strcpy( chLab10, chLineLbl(t20) );
166 fprintf( ioQQQ, " PROBLEM atom_level3: insane stat weights g2 :%10.2e%10.2e %10.10s\n",
167 g010, g020, chLab10 );
169 }
170
171 /* abundances from the elements grid
172 * one of these must be a true line */
173 if( g010 > 0. )
174 {
175 /* put null terminated line label into chLab using line array*/
176 chIonLbl(chLab, t10);
177 AbunxIon = dense.xIonDense[ (*t10.Hi()).nelem() -1][(*t10.Hi()).IonStg()-1];
178 }
179
180 else if( g121 > 0. )
181 {
182 /* put null terminated line label into chLab using line array*/
183 chIonLbl(chLab, t21);
184 AbunxIon = dense.xIonDense[(*t21.Hi()).nelem() -1][(*t21.Hi()).IonStg()-1];
185 }
186
187 else
188 /* this cannot possibly happen */
189 {
190 chLab[0] = ' ';
191 AbunxIon = 0.;
192 fprintf( ioQQQ, " PROBLEM atom_level3: insanity at g010 g121 branch \n" );
194 }
195
196 a = t10.EnergyK()*phycon.teinv;
197 b = t21.EnergyK()*phycon.teinv;
198 c = t20.EnergyK()*phycon.teinv;
199
200 if( c == 0. )
201 {
202 c = a + b;
203 }
204
205 /* if still neg at end, then success!, so possible to
206 * to check why zero returned */
207 nLev3Fail = -1;
208
209 /* if two of the lines are below the plasma frequency, hand the third in atom_level2 */
210 if( ( t10.EnergyErg()/EN1RYD < rfield.plsfrq && t20.EnergyErg()/EN1RYD < rfield.plsfrq ) )
211 {
212 atom_level2( t21 );
213 return;
214 }
215 else if( ( t10.EnergyErg()/EN1RYD < rfield.plsfrq && t21.EnergyErg()/EN1RYD < rfield.plsfrq ) )
216 {
217 atom_level2( t20 );
218 return;
219 }
220 else if( ( t20.EnergyErg()/EN1RYD < rfield.plsfrq && t21.EnergyErg()/EN1RYD < rfield.plsfrq ) )
221 {
222 atom_level2( t10 );
223 return;
224 }
225
228 if( AbunxIon <= 1e-30 || c > 60. )
229 {
230 nLev3Fail = 0;
231
232 /* all populations are zero */
233 atoms.PopLevels[0] = AbunxIon;
234 atoms.PopLevels[1] = 0.;
235 atoms.PopLevels[2] = 0.;
236
238 atoms.DepLTELevels[0] = 1.;
239 atoms.DepLTELevels[1] = 0.;
240 atoms.DepLTELevels[2] = 0.;
241
242 /* level populations */
243 t21.Emis().PopOpc() = 0.;
244 t10.Emis().PopOpc() = AbunxIon;
245 t20.Emis().PopOpc() = AbunxIon;
246 (*t21.Lo()).Pop() = 0.;
247 (*t10.Lo()).Pop() = AbunxIon;
248 (*t20.Lo()).Pop() = AbunxIon;
249 (*t21.Hi()).Pop() = 0.;
250 (*t10.Hi()).Pop() = 0.;
251 (*t20.Hi()).Pop() = 0.;
252
253 /* line heating */
254 t20.Coll().heat() = 0.;
255 t21.Coll().heat() = 0.;
256 t10.Coll().heat() = 0.;
257
258 /* intensity of line */
259 t21.Emis().xIntensity() = 0.;
260 t10.Emis().xIntensity() = 0.;
261 t20.Emis().xIntensity() = 0.;
262
263 /* line cooling */
264 t20.Coll().cool() = 0.;
265 t21.Coll().cool() = 0.;
266 t10.Coll().cool() = 0.;
267
268 /* local ots rates */
269# if 0
270 /* >>chng 03 oct 04, move to RT_OTS */
271 t20.ots() = 0.;
272 t21.ots() = 0.;
273 t10.ots() = 0.;
274# endif
275
276 /* number of photons in line zero */
277 t21.Emis().phots() = 0.;
278 t10.Emis().phots() = 0.;
279 t20.Emis().phots() = 0.;
280
281 /* ratio coll over total excitation */
282 t21.Emis().ColOvTot() = 0.;
283 t10.Emis().ColOvTot() = 0.;
284 t20.Emis().ColOvTot() = 0.;
285
286 /* add zero to cooling */
287 CoolAdd(chLab, t21.WLAng() ,0.);
288 CoolAdd(chLab, t10.WLAng() ,0.);
289 CoolAdd(chLab, t20.WLAng() ,0.);
290 return;
291 }
292
293 /* collision strengths */
294 o10 = t10.Coll().col_str();
295 o21 = t21.Coll().col_str();
296 o20 = t20.Coll().col_str();
297
298 /* begin sanity checks, check statistic weights,
299 * first check is protection against dummy lines */
300 ASSERT( g010*g020 == 0. || fp_equal( g010, g020 ) );
301
302 ASSERT( g110*g121 == 0. || fp_equal( g110, g121 ) );
303
304 ASSERT( g221*g220 == 0. || fp_equal( g221, g220 ) );
305
306 /* both abundances must be same, equal abundance
307 * dense.xIonDense(nelem,i) is density of ith ionization stage (cm^-3) */
308 ASSERT( ((*t10.Hi()).IonStg()*(*t21.Hi()).IonStg() == 0) || ((*t10.Hi()).IonStg() == (*t21.Hi()).IonStg()));
309
310 ASSERT( ((*t20.Hi()).IonStg()*(*t21.Hi()).IonStg() == 0) || ((*t20.Hi()).IonStg() == (*t21.Hi()).IonStg() ) );
311
312 ASSERT( ((*t10.Hi()).nelem() * (*t21.Hi()).nelem() == 0) || ((*t10.Hi()).nelem() == (*t21.Hi()).nelem()) );
313
314 ASSERT( ((*t20.Hi()).nelem() * (*t21.Hi()).nelem() == 0) || ((*t20.Hi()).nelem() == (*t21.Hi()).nelem()) );
315
316 ASSERT( o10 > 0. && o21 > 0. && o20 > 0. );
317
318 /*end sanity checks */
319
320 /* net loss of line escape and destruction */
321 a21 = t21.Emis().Aul() * (t21.Emis().Pesc()+ t21.Emis().Pelec_esc() + t21.Emis().Pdest());
322 a10 = t10.Emis().Aul() * (t10.Emis().Pesc()+ t10.Emis().Pelec_esc() + t10.Emis().Pdest());
323 a20 = t20.Emis().Aul() * (t20.Emis().Pesc()+ t20.Emis().Pelec_esc() + t20.Emis().Pdest());
324
325 /* find energies of all transitions - one line could be a dummy
326 * also find Boltzmann factors */
327 if( t10.Emis().Aul() == 0. )
328 {
329 ener20 = t20.EnergyErg();
330 ener21 = t21.EnergyErg();
331 ener10 = ener20 - ener21;
332 bolt12 = exp(-t21.EnergyK()/phycon.te);
333 bolt02 = exp(-t20.EnergyK()/phycon.te);
334 bolt01 = bolt02/bolt12;
335 temp12 = t21.EnergyK();
336 temp02 = t20.EnergyK();
337 temp01 = temp02 - temp12;
338 }
339
340 else if( t21.Emis().Aul() == 0. )
341 {
342 ener10 = t10.EnergyErg();
343 ener20 = t20.EnergyErg();
344 ener21 = ener20 - ener10;
345 bolt01 = exp(-t10.EnergyK()/phycon.te);
346 bolt02 = exp(-t20.EnergyK()/phycon.te);
347 bolt12 = bolt02/bolt01;
348 temp02 = t20.EnergyK();
349 temp01 = t10.EnergyK();
350 temp12 = temp02 - temp01;
351 }
352
353 else if( t20.Emis().Aul() == 0. )
354 {
355 ener10 = t10.EnergyErg();
356 ener21 = t21.EnergyErg();
357 ener20 = ener21 + ener10;
358 bolt01 = exp(-t10.EnergyK()/phycon.te);
359 bolt12 = exp(-t21.EnergyK()/phycon.te);
360 bolt02 = bolt01*bolt12;
361 temp01 = t10.EnergyK();
362 temp12 = t21.EnergyK();
363 temp02 = temp01 + temp12;
364 }
365
366 else
367 {
368 /* all lines are ok */
369 ener10 = t10.EnergyErg();
370 ener20 = t20.EnergyErg();
371 ener21 = t21.EnergyErg();
372 bolt01 = exp(-t10.EnergyK()/phycon.te);
373 bolt12 = exp(-t21.EnergyK()/phycon.te);
374 bolt02 = bolt01*bolt12;
375 temp02 = t20.EnergyK();
376 temp01 = t10.EnergyK();
377 temp12 = t21.EnergyK();
378 }
379
380 /* check all energies positive */
381 ASSERT( ener10 > 0. && ener20 > 0. && ener21 > 0. );
382
383 /* check if energy order is ok */
384 ASSERT( ener10 < ener20 && ener21 < ener20 );
385
386 /* check if energy scale is ok */
387 ASSERT( fabs((ener10+ener21)/ener20-1.) < 1e-4 );
388
389 pump01 = t10.Emis().pump();
390 pump10 = pump01*g0/g1;
391 pump12 = t21.Emis().pump();
392 pump21 = pump12*g1/g2;
393 pump02 = t20.Emis().pump();
394 pump20 = pump02*g0/g2;
395
396 /* cdsqte is 8.629E-6 / sqrte * eden */
397 c01 = o10*bolt01*dense.cdsqte/g0;
398 r01 = c01 + pump01;
399 c10 = o10*dense.cdsqte/g1;
400 r10 = c10 + a10 + pump10;
401 c20 = o20*dense.cdsqte/g2;
402 r20 = c20 + a20 + pump20;
403 c02 = o20*bolt02*dense.cdsqte/g0;
404 r02 = c02 + pump02;
405 c12 = o21*bolt12*dense.cdsqte/g1;
406 r12 = c12 + pump12;
407 c21 = o21*dense.cdsqte/g2;
408 r21 = c21 + a21 + pump21;
409
410 alpha1 = (double)(AbunxIon)*(r01+r02)/(r10+r01+r02);
411 alpha2 = (double)(AbunxIon)*(r01)/(r10+r12+r01);
412 alpha = alpha1 - alpha2;
413
414 /* 1( DBLE(r01+r02)/DBLE(r10+r01+r02) - DBLE(r01)/DBLE(r10+r12+r01) )
415 * beta is factor with n2 */
416 beta = (r21 - r01)/(r10 + r12 + r01) + (r20 + r01 + r02)/(r10 +
417 r01 + r02);
418
419 if( alpha/MAX2(alpha1,alpha2) < 1e-11 )
420 {
421 /* this catches both negative and round off */
422 p2 = 0.;
423 alpha = 0.;
424 nLev3Fail = 1;
425 }
426
427 else
428 {
429 p2 = alpha/beta;
430 }
431 atoms.PopLevels[2] = p2;
432
433 if( alpha < 0. || beta < 0. )
434 {
435 fprintf( ioQQQ, " atom_level3: insane n2 pop alf, bet, p2=%10.2e%10.2e%10.2e %10.10s t=%10.2e\n",
436 alpha, beta, p2, chLab, phycon.te );
437 fprintf( ioQQQ, " gs are%5.1f%5.1f%5.1f\n", g0, g1,
438 g2 );
439 fprintf( ioQQQ, " Bolts are%10.2e%10.2e%10.2e\n",
440 bolt01, bolt12, bolt02 );
441 fprintf( ioQQQ, " As are%10.2e%10.2e%10.2e\n", a10,
442 a21, a20 );
443 fprintf( ioQQQ, " Energies are%10.2e%10.2e%10.2e\n",
444 ener10, ener21, ener20 );
445 fprintf( ioQQQ, " 2 terms, dif of alpha are%15.6e%15.6e\n",
446 (r01 + r02)/(r10 + r01 + r02), r01/(r10 + r12 + r01) );
447 ShowMe();
449 }
450
451 alpha = (double)(AbunxIon)*(r01+r02) - (double)(p2)*(r20+r01+r02);
452 /* guard against roundoff - this should really have been zero
453 * >>chng 96 nov 14, protection against round-off to zero
454 * >>chng 96 dec 03, made r01, etc, double, and changed limit to 1e-9 */
455 if( fabs(alpha)/(MAX2(AbunxIon*(r01+r02),p2*(r20+r01+r02))) < 1e-9 )
456 {
457 alpha = 0.;
458 nLev3Fail = 2;
459 }
460
461 beta = r10 + r01 + r02;
462 p1 = alpha/beta;
463 atoms.PopLevels[1] = p1;
464
465 if( p1 < 0. )
466 {
467 if( p1 > -1e-37 )
468 {
469 /* slightly negative solution, probably just round-off, zero it */
470 p1 = 0.;
471 atoms.PopLevels[1] = p1;
472 nLev3Fail = 3;
473 }
474
475 else
476 {
477 /* very negative solution, better punt */
478 fprintf( ioQQQ, " atom_level3: insane n1 pop alf, bet, p1=%10.2e%10.2e%10.2e %10.10s%5f\n",
479 alpha, beta, p1, chLab, t10.WLAng() );
480 fprintf( ioQQQ, " local electron density and temperature were%10.2e%10.2e\n",
481 dense.eden, phycon.te );
482 ShowMe();
484 }
485 }
486
487 p0 = AbunxIon - p2 - p1;
488
489 /* population of lowest level */
490 atoms.PopLevels[0] = p0;
491 if( p0 <= 0. )
492 {
493 fprintf( ioQQQ, " atom_level3: insane n0 pop p1, 2, abun=%10.2e%10.2e%10.2e \n",
494 p1, p2, AbunxIon );
495 ShowMe();
497 }
498
499 /* level populations for line opacities */
500 (*t21.Lo()).Pop() = p1;
501 (*t10.Lo()).Pop() = p0;
502 (*t20.Lo()).Pop() = p0;
503 t21.Emis().PopOpc() = (p1 - p2*g1/g2);
504 t10.Emis().PopOpc() = (p0 - p1*g0/g1);
505 t20.Emis().PopOpc() = (p0 - p2*g0/g2);
506 (*t21.Hi()).Pop() = p2;
507 (*t10.Hi()).Pop() = p1;
508 (*t20.Hi()).Pop() = p2;
509
510 /* line emission - net emission in line */
511 t21.Emis().phots() = t21.Emis().Aul() * (t21.Emis().Pesc() + t21.Emis().Pelec_esc())*p2;
512 t21.Emis().xIntensity() = t21.Emis().phots() * t21.EnergyErg();
513
514 t20.Emis().phots() = t20.Emis().Aul() * (t20.Emis().Pesc() + t20.Emis().Pelec_esc())*p2;
515 t20.Emis().xIntensity() = t20.Emis().phots() * t20.EnergyErg();
516
517 t10.Emis().phots() = t10.Emis().Aul() * (t10.Emis().Pesc() + t10.Emis().Pelec_esc())*p1;
518 t10.Emis().xIntensity() = t10.Emis().phots() * t10.EnergyErg();
519
520# if 0
521 /* >>chng 03 oct 04, move to RT_OTS */
522 t21.ots() = p2 * t21.Emis().Aul() * t21.Emis().Pdest();
523 t20.ots() = p2 * t20.Emis().Aul() * t20.Emis().Pdest();
524 t10.ots() = p2 * t10.Emis().Aul() * t10.Emis().Pdest();
525 /* now add thess lines to ots field, routine works on f not c scale */
526 RT_OTS_AddLine( t21.ots() , t21.ipCont);
527 RT_OTS_AddLine( t20.ots() , t20.ipCont);
528 RT_OTS_AddLine( t10.ots() , t10.ipCont);
529# endif
530
531 /* total intensity used to divide line up - one may be 0 */
532 /* >>chng 99 nov 30, rewrite algebra so double prec throughout,
533 * for very low density models */
534 /*TotInten = t21.Emis().xIntensity() + t20.xIntensity() + t10.xIntensity();*/
535 TotInten = t21.Emis().phots() * (double)t21.EnergyErg()
536 + t20.Emis().phots() * (double)t20.EnergyErg() +
537 t10.Emis().phots() * (double)t10.EnergyErg();
538
539 /* fraction that was collisionally excited */
540 if( r12 > 0. )
541 {
542 t21.Emis().ColOvTot() = c12/r12;
543 }
544 else
545 {
546 t21.Emis().ColOvTot() = 0.;
547 }
548
549 if( r01 > 0. )
550 {
551 t10.Emis().ColOvTot() = c01/r01;
552 }
553 else
554 {
555 t10.Emis().ColOvTot() = 0.;
556 }
557
558 if( r02 > 0. )
559 {
560 t20.Emis().ColOvTot() = c02/r02;
561 }
562 else
563 {
564 t20.Emis().ColOvTot() = 0.;
565 }
566
567 /* heating or cooling due to each line */
568 heat20 = p2*c20*ener20;
569 cool02 = p0*c02*ener20;
570 heat21 = p2*c21*ener21;
571 cool12 = p1*c12*ener21;
572 heat10 = p1*c10*ener10;
573 cool01 = p0*c01*ener10;
574
575 /* two cases - collisionally excited (usual case) or
576 * radiatively excited - in which case line can be a heat source
577 * following are correct heat exchange, they will mix to get correct deriv
578 * the sum of heat-cool will add up to EnrUL - EnrLU - this is a trick to
579 * keep stable solution by effectively dividing up heating and cooling,
580 * so that negative cooling does not occur */
581
582 /* now get net heating or cooling */
583 cnet02 = cool02 - heat20*t20.Emis().ColOvTot();
584 hnet02 = heat20 *(1. - t20.Emis().ColOvTot());
585 cnet12 = cool12 - heat21*t21.Emis().ColOvTot();
586 hnet12 = heat21 *(1. - t21.Emis().ColOvTot());
587 cnet01 = cool01 - heat10*t10.Emis().ColOvTot();
588 hnet01 = heat10 *(1. - t10.Emis().ColOvTot());
589
590 /*TotalCooling = p0*(c01*ener10 + c02*ener20) + p1*c12*ener21 -
591 (p2*(c21*ener21 + c20*ener20) + p1*c10*ener10);*/
592 /* >>chng 96 nov 22, very dense cool models, roundoff error
593 *could cause [OI] 63 mic to be dominant heating, cooling, or
594 *just zero
595 * >>chng 96 dec 06, above change caused o iii 1666 cooling to
596 * be zeroed when important in a model in kirk's grid. was at 1e-6,
597 * set to 1e-7
598 * >>chng 96 dec 17, from 1e-7 to 1e-8, caused temp fail */
599 /* >>chng 99 nov 29, had been 1e-30 to prevent div by zero (?),
600 * change to dble max since 1e-30 was very large compared to
601 * cooling for 1e-10 cm-3 den models */
602 /*if( fabs(cnet01/MAX2(1e-30,cool01)) < 1e-8 )*/
603
604 /* >>chng 02 jan 28, min from 1e-8 to 1e-10, conserve.in had massive
605 * temp failures when no molecules turned on, due to this tripping */
606 /*if( fabs(cnet01/MAX2(DBL_MIN,cool01)) < 1e-8 )*/
607 if( fabs(cnet01/MAX2(DBL_MIN,cool01)) < 1e-10 )
608 {
609 nLev3Fail = 4;
610 cnet02 = 0.;
611 hnet02 = 0.;
612 cnet12 = 0.;
613 hnet12 = 0.;
614 cnet01 = 0.;
615 hnet01 = 0.;
616 }
617
618 TotCool = cnet02 + cnet12 + cnet01;
619 TotHeat = hnet02 + hnet12 + hnet01;
620
621
622 if( TotInten > 0. )
623 {
624 cool02 = TotCool * t20.Emis().phots() * (double)t20.EnergyErg()/TotInten;
625 cool12 = TotCool * t21.Emis().phots() * (double)t21.EnergyErg()/TotInten;
626 cool01 = TotCool * t10.Emis().phots() * (double)t10.EnergyErg()/TotInten;
627 heat20 = TotHeat * t20.Emis().phots() * (double)t20.EnergyErg()/TotInten;
628 heat21 = TotHeat * t21.Emis().phots() * (double)t21.EnergyErg()/TotInten;
629 heat10 = TotHeat * t10.Emis().phots() * (double)t10.EnergyErg()/TotInten;
630 t20.Coll().cool() = cool02;
631 t21.Coll().cool() = cool12;
632 t10.Coll().cool() = cool01;
633 t20.Coll().heat() = heat20;
634 t21.Coll().heat() = heat21;
635 t10.Coll().heat() = heat10;
636 }
637 else
638 {
639 nLev3Fail = 5;
640 cool02 = 0.;
641 cool12 = 0.;
642 cool01 = 0.;
643 heat20 = 0.;
644 heat21 = 0.;
645 heat10 = 0.;
646 t20.Coll().cool() = 0.;
647 t21.Coll().cool() = 0.;
648 t10.Coll().cool() = 0.;
649 t20.Coll().heat() = 0.;
650 t21.Coll().heat() = 0.;
651 t10.Coll().heat() = 0.;
652 }
653
654 /* add cooling due to each line,
655 * heating broken out above, will be added to thermal.heating[0][22] in CoolEvaluate*/
656 /* >>chng 99 nov 30, rewrite algebra to keep precision on very low density models*/
657 CoolAdd(chLab, t21.WLAng() ,cool12);
658 CoolAdd(chLab, t10.WLAng() ,cool01);
659 CoolAdd(chLab, t20.WLAng() ,cool02);
660
661 /* derivative of cooling function
662 * dC/dT = Cooling * ( T(excitation)/T_e^2 - 1/(2T) )
663 * in following I assume that a 1-2 exciation will have the 0-2
664 * exponential in the dcdt term = NOT A TYPO */
665
666 thermal.dCooldT += t10.Coll().cool()*(temp01*thermal.tsq1 - thermal.halfte) +
667 (t20.Coll().cool() + t21.Coll().cool())*(temp02*thermal.tsq1 - thermal.halfte);
668 /* two t20.t's above are not a typo!
669 * */
670 {
671 enum{DEBUG_LOC=false};
672 if( DEBUG_LOC )
673 {
674 fprintf(ioQQQ,"atom_level3 nLev3Fail %i\n",
675 nLev3Fail );
676 }
677 }
678 return;
679}
void atom_level2(const TransitionProxy &t)
void atom_level3(const TransitionProxy &t10, const TransitionProxy &t21, const TransitionProxy &t20)
t_atoms atoms
Definition atoms.cpp:5
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
#define MAX2
Definition cddefines.h:782
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
NORETURN void TotalInsanity(void)
Definition service.cpp:886
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void ShowMe(void)
Definition service.cpp:181
double & cool() const
Definition collision.h:190
realnum & col_str() const
Definition collision.h:167
double & heat() const
Definition collision.h:194
double & ColOvTot() const
Definition emission.h:573
double & PopOpc() const
Definition emission.h:603
double & xIntensity() const
Definition emission.h:483
realnum & Pelec_esc() const
Definition emission.h:533
realnum & Aul() const
Definition emission.h:613
double & pump() const
Definition emission.h:473
realnum & Pdest() const
Definition emission.h:543
double & phots() const
Definition emission.h:503
CollisionProxy Coll() const
Definition transition.h:424
realnum & WLAng() const
Definition transition.h:429
realnum EnergyErg() const
Definition transition.h:78
qList::iterator Lo() const
Definition transition.h:392
long & ipCont() const
Definition transition.h:450
realnum EnergyK() const
Definition transition.h:73
qList::iterator Hi() const
Definition transition.h:396
EmissionList::reference Emis() const
Definition transition.h:408
void CoolAdd(const char *chLabel, realnum lambda, double cool)
Definition cool_etc.cpp:13
t_dense dense
Definition dense.cpp:24
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double EN1RYD
Definition physconst.h:179
t_rfield rfield
Definition rfield.cpp:8
void RT_OTS_AddLine(double ots, long int ip)
Definition rt_ots.cpp:402
t_thermal thermal
Definition thermal.cpp:5
void chIonLbl(char *chIonLbl_v, const TransitionProxy &t)
char * chLineLbl(const TransitionProxy &t)