cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
atom_level3.cpp
Go to the documentation of this file.
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 
15 void 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 
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.cs;
295  o21 = t21->Coll.cs;
296  o20 = t20->Coll.cs;
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 }

Generated for cloudy by doxygen 1.8.3.1