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