root/branches/newmole/source/atom_level3.cpp

Revision 2139, 16.6 kB (checked in by rjrw, 6 months ago)

Merged from trunk r2124:2137

  • Property svn:eol-style set to native
Line 
1/* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
2 * others.  For conditions of distribution and use see copyright notice in license.txt */
3/*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 "lines_service.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
15void atom_level3(transition * t10,
16  transition * t21,
17  transition * 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 );
131                TotalInsanity();
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 );
150                TotalInsanity();
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 );
168                TotalInsanity();
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" );
193                TotalInsanity();
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
226        /** \todo       2       test on c checks whether collisions are possible at this temperature,
227         * should add photo excitation */
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
237                /** \todo       2       these pops ARE NOT defined below */
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( a10 == 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( a21 == 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( a20 == 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();
448                cdEXIT(EXIT_FAILURE);
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();
483                        cdEXIT(EXIT_FAILURE);
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();
496                cdEXIT(EXIT_FAILURE);
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                /*@-redef@*/
672                enum{DEBUG_LOC=false};
673                /*@+redef@*/
674                if( DEBUG_LOC )
675                {
676                        fprintf(ioQQQ,"atom_level3 nLev3Fail %i\n",
677                                nLev3Fail );
678                }
679        }
680        return;
681}
Note: See TracBrowser for help on using the browser.