cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
atom_level2.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_level2 do level population and cooling for two level atom,
4  * side effects:
5  * set elements of transition struc
6  * cooling via CoolAdd( chLab, (long)t->WLAng , t->cool);
7  * cooling derivative */
8 #include "cddefines.h"
9 #include "phycon.h"
10 #include "lines_service.h"
11 #include "dense.h"
12 #include "rfield.h"
13 #include "thermal.h"
14 #include "cooling.h"
15 #include "atoms.h"
16 
18 {
19  char chLab[5];
20  long int ion,
21  ip,
22  nelem;
23  double AbunxIon,
24  a21,
25  boltz,
26  col12,
27  col21,
28  coolng,
29  g1,
30  g2,
31  omega,
32  pfs1,
33  pfs2,
34  plower,
35  r,
36  rate12,
37  ri21;
38 
39  DEBUG_ENTRY( "atom_level2()" );
40 
41  /* result is density (cm-3) of excited state times A21
42  * result normalized to N1+N2=ABUND
43  * routine increments dCooldT, call CoolAdd
44  * CDSQTE is EDEN / SQRTE * 8.629E-6
45  */
46 
47  ion = t->Hi->IonStg;
48  nelem = t->Hi->nelem;
49 
50  /* dense.xIonDense[nelem][i] is density of ith ionization stage (cm^-3) */
51  AbunxIon = dense.xIonDense[nelem-1][ion-1];
52 
53  /* continuum pointer */
54  ip = t->ipCont;
55 
56  /* approximate Boltzmann factor to see if results zero */
57  boltz = rfield.ContBoltz[ip-1];
58 
59  /* collision strength for this transition, omega is zero for hydrogenic
60  * species which are done in special hydro routines */
61  omega = t->Coll.cs;
62 
63  /* ROUGH check whether upper level has significant population,*/
64  r = (boltz*dense.cdsqte + t->Emis->pump)/(dense.cdsqte + t->Emis->Aul);
65 
66  /* following first test needed for 32 bit cpu on search phase
67  * >>chng 96 jul 02, was 1e-30 crashed on dec, change to 1e-25
68  * if( AbunxIon.lt.1e-25 .or. boltz.gt.30. ) then
69  * >>chng 96 jul 11, to below, since can be strong pumping when
70  * Boltzmann factor essentially zero */
71  /* omega in following is zero for hydrogenic species, since done
72  * in hydro routines, so this should cause us to quit on this test */
73  /* >>chng 99 nov 29, from 1e-25 to 1e-30 to keep same result for
74  * very low density models, where AbunxIon is very small but still significant*/
75  /*if( omega*AbunxIon < 1e-25 || r < 1e-25 )*/
76  if( omega*AbunxIon < 1e-30 || r < 1e-25 )
77  {
78  /* put in pop since possible just too cool */
79  t->Lo->Pop = AbunxIon;
80  t->Emis->PopOpc = AbunxIon;
81  t->Hi->Pop = 0.;
82  t->Emis->xIntensity = 0.;
83  t->Coll.cool = 0.;
84  t->Emis->ots = 0.;
85  t->Emis->phots = 0.;
86  t->Emis->ColOvTot = 0.;
87  t->Coll.heat = 0.;
88  /* level populations */
89  atoms.PopLevels[0] = AbunxIon;
90  atoms.PopLevels[1] = 0.;
91  atoms.DepLTELevels[0] = 1.;
92  atoms.DepLTELevels[1] = 0.;
93  return;
94  }
95 
96  /* net rate down A21*(escape + destruction) */
97  a21 = t->Emis->Aul*(t->Emis->Pesc+ t->Emis->Pdest + t->Emis->Pelec_esc);
98 
99  /* put null terminated line label into chLab using line array*/
100  chIonLbl(chLab,t);
101 
102  /* statistical weights of upper and lower levels */
103  g1 = t->Lo->g;
104  g2 = t->Hi->g;
105 
106  /* now get real Boltzmann factor */
107  boltz = t->EnergyK/phycon.te;
108 
109  ASSERT( boltz > 0. );
110  boltz = sexp(boltz);
111 
112  ASSERT( g1 > 0. && g2 > 0. );
113 
114  /* this lacks the upper statistical weight */
115  col21 = dense.cdsqte*omega;
116  /* upward coll rate s-1 */
117  col12 = col21/g1*boltz;
118  /* downward coll rate s-1 */
119  col21 /= g2;
120 
121  /* rate 1 to 2 is both collisions and pumping */
122  /* the total excitation rate from lower to upper, collisional and pumping */
123  rate12 = col12 + t->Emis->pump;
124 
125  /* induced emissions down */
126  ri21 = t->Emis->pump*g1/g2;
127 
128  /* this is the ratio of lower to upper level populations */
129  r = (a21 + col21 + ri21)/rate12;
130 
131  /* upper level pop */
132  pfs2 = AbunxIon/(r + 1.);
133  atoms.PopLevels[1] = pfs2;
134  t->Hi->Pop = pfs2;
135 
136  /* pop of ground */
137  pfs1 = pfs2*r;
138  atoms.PopLevels[0] = pfs1;
139 
140  /* compute ratio Aul/(Aul+Cul) */
141  /* level population with correction for stimulated emission */
142  t->Lo->Pop = atoms.PopLevels[0];
143 
144  t->Emis->PopOpc = (atoms.PopLevels[0] - atoms.PopLevels[1]*g1/g2 );
145 
146  /* departure coef of excited state rel to ground */
147  atoms.DepLTELevels[0] = 1.;
148  if( boltz > 1e-20 && atoms.PopLevels[1] > 1e-20 )
149  {
150  /* this line could have zero boltz factor but radiatively excited
151  * dec alpha does not obey () in fast mode?? */
153  (boltz*g2/g1);
154  }
155  else
156  {
157  atoms.DepLTELevels[1] = 0.;
158  }
159 
160  /* number of escaping line photons, used elsewhere for outward beam */
161  t->Emis->phots = t->Emis->Aul*(t->Emis->Pesc + t->Emis->Pelec_esc)*pfs2;
162 
163  /* intensity of line */
164  t->Emis->xIntensity = t->Emis->phots*t->EnergyErg;
165  plower = AbunxIon - pfs2;
166 
167  /* ratio of collisional to total (collisional + pumped) excitation */
168  t->Emis->ColOvTot = col12/rate12;
169 
170  /* two cases - collisionally excited (usual case) or
171  * radiatively excited - in which case line can be a heat source
172  * following are correct heat exchange, they will mix to get correct deriv
173  * the sum of heat-cool will add up to EnrUL - EnrLU - this is a trick to
174  * keep stable solution by effectively dividing up heating and cooling,
175  * so that negative cooling does not occur */
176 
177  //double Enr12 = plower*col12*t->EnergyErg;
178  //double Enr21 = pfs2*col21*t->EnergyErg;
179 
180  /* energy exchange due to this level
181  * net cooling due to excit minus part of de-excit -
182  * note that ColOvTot cancels out in the sum heat - cool */
183  //coolng = Enr12 - Enr21*t->Emis->ColOvTot;
184 
185  /* this form of coolng is guaranteed to be non-negative */
186  coolng = t->EnergyErg*AbunxIon*col12*(a21 + ri21)/(a21 + col21 + ri21 + rate12);
187  ASSERT( coolng >= 0. );
188 
189  t->Coll.cool = coolng;
190 
191  /* net heating is remainder */
192  //t->Coll.heat = Enr21*(1. - t->Emis->ColOvTot);
193  t->Coll.heat = t->EnergyErg*AbunxIon*col21*(t->Emis->pump)/(a21 + col21 + ri21 + rate12);
194 
195 
196  /* expression pre jul 3 95, changed for case where line heating dominates
197  * coolng = (plower*col12 - pfs2*col21)*t->t(ipLnEnrErg)
198  * t->t(ipLnCool) = cooling */
199 
200  /* add to cooling - heating part was taken out above,
201  * and is not added in here - it will be added to thermal.heating[0][22]
202  * in CoolSum */
203  CoolAdd( chLab, t->WLAng , t->Coll.cool);
204 
205  /* derivative of cooling function */
206  thermal.dCooldT += coolng * (t->EnergyK * thermal.tsq1 - thermal.halfte );
207  return;
208 }

Generated for cloudy by doxygen 1.8.3.1