cloudy trunk
Loading...
Searching...
No Matches
atom_seq_beryllium.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/*AtomSeqBeryllium compute level populations and emissivity for Be-sequence ions */
4#include "cddefines.h"
5#include "phycon.h"
6#include "transition.h"
7#include "thirdparty.h"
8#include "dense.h"
9#include "cooling.h"
10#include "thermal.h"
11#include "atoms.h"
12
13void AtomSeqBeryllium(double cs12,
14 double cs13,
15 double cs23,
16 const TransitionProxy& t,
17 double a30)
18{
19
20 char chLab[5];
21
22 bool lgNegPop;
23
24 int32 ipiv[4], nerror;
25 long int i, j;
26
27 double AbunxIon,
28 Enr01,
29 Enr10,
30 a20,
31 boltz,
32 c01,
33 c02,
34 c03,
35 c10,
36 c12,
37 c13,
38 c20,
39 c21,
40 c23,
41 c30,
42 c31,
43 c32,
44 coolng,
45 cs1u,
46 excit,
47 r01,
48 r02,
49 r03,
50 r10,
51 r12,
52 r13,
53 r20,
54 r21,
55 r23,
56 r30,
57 r31,
58 r32,
59 ri02,
60 ri20,
61 tot20;
62
63 double amat[4][4],
64 bvec[4],
65 zz[5][5];
66
67 DEBUG_ENTRY( "AtomSeqBeryllium()" );
68
69 /* function returns total emission in both components of line
70 * destruction by background opacity computed and added to otslin field
71 * here,
72 * total cooling computed and added via CoolAdd
73 *
74 * finds population of four level Be-like atom
75 *
76 * level 1 = ground, J=0
77 * levels 2,3,4 are J=0,1,2, OF 3P.
78 * levels 3-1 is the fast one, j=1 to ground j=0
79 *
80 * cs1u is coll str to all J sub-levels of 3P; C.S. goes as 2J+1
81 * when routine exits the collision strength in the fast line is 1/3 of the entry value,
82 * so that save lines data does give the right critical density */
83
84 /* AtomSeqBeryllium(cs12,cs13,cs23,(*t).t,a30,chLab)
85 * dense.xIonDense[Hi.nelem(),i] is density of ith ionization stage (cm^-3) */
86 AbunxIon = dense.xIonDense[(*t.Hi()).nelem() -1][(*t.Hi()).IonStg()-1];
87
88 /* put null terminated line label into chLab using line array*/
89 chIonLbl(chLab, t);
90 boltz = t.EnergyK()/phycon.te;
91
92 /* set the cs before if below, since we must reset to line cs in all cases */
93 cs1u = t.Coll().col_str();
94 /* >>chng 01 sep 09, reset cs coming in, to propoer cs for the ^3P_1 level only
95 * so that critical density is printed correctly with save lines data command */
96 t.Coll().col_str() /= 3.f;
97
98 /* low density cutoff to keep matrix happy */
99 if( AbunxIon <= 1e-20 || boltz > 30. )
100 {
101 /* put in pop since possible just too cool */
102 (*t.Lo()).Pop() = AbunxIon;
103 t.Emis().PopOpc() = AbunxIon;
104 (*t.Hi()).Pop() = 0.;
105 t.Emis().xIntensity() = 0.;
106 t.Coll().cool() = 0.;
107 t.Emis().phots() = 0.;
108 /*>>chng 03 oct 04, move to RT_OTS */
109 /*t.ots() = 0.;*/
110 t.Emis().ColOvTot() = 0.;
111 t.Coll().heat() = 0.;
112 /* level populations */
113 atoms.PopLevels[0] = AbunxIon;
114 atoms.PopLevels[1] = 0.;
115 atoms.PopLevels[2] = 0.;
116 atoms.PopLevels[3] = 0.;
117 atoms.DepLTELevels[0] = 1.;
118 atoms.DepLTELevels[1] = 0.;
119 atoms.DepLTELevels[2] = 0.;
120 atoms.DepLTELevels[3] = 0.;
121 CoolAdd(chLab, t.WLAng() ,0.);
122 return;
123 }
124 excit = exp(-boltz);
125
126 /* these must be the statistical weights */
127 ASSERT( (*t.Lo()).g() == 1. );
128 ASSERT( (*t.Hi()).g() == 3. );
129
130 /* collision strength must be positive */
131 ASSERT( cs1u > 0. );
132
133 /* incuded rates for fastest transition */
134 ri02 = t.Emis().pump();
135
136 /* back reaction has ratio of stat weights */
137 ri20 = t.Emis().pump()*1./3.;
138
139 /* net rate out of level 3, with destruction */
140 a20 = t.Emis().Aul()*(t.Emis().Pesc() + t.Emis().Pelec_esc() + t.Emis().Pdest());
141 tot20 = a20 + ri20;
142
143 /* rates between j=0, lowest 3P level,
144 * 1/9 is ratio of level stat weight to term stat weight */
145 c10 = cs1u*dense.cdsqte/9.;
146 c01 = c10*excit;
147 r01 = c01;
148 r10 = c10;
149
150 /* stat weights canceled out here */
151 c20 = c10;
152 c02 = c01*3.;
153 r02 = c02 + ri02;
154 r20 = c20 + tot20;
155
156 c30 = c10;
157 c03 = c01*5.;
158 r30 = c30 + a30;
159 r03 = c03;
160
161 c21 = cs12*dense.cdsqte/3.;
162 c12 = c21*3.;
163 r12 = c12;
164 r21 = c21;
165
166 c31 = cs13*dense.cdsqte/5.;
167 c13 = c31*5.;
168 r31 = c31;
169 r13 = c13;
170
171 c32 = cs23*dense.cdsqte/5.;
172 c23 = c32*1.667;
173 r32 = c32;
174 r23 = c23;
175
176 /* set up matrix */
177 for( i=0; i <= 3; i++ )
178 {
179 /* first equation will be sum to abund */
180 zz[i][0] = 1.;
181 zz[4][i] = 0.;
182 }
183
184 /* zz(0,4) = AbunxIon */
185 zz[4][0] = 1.;
186
187 /* ground level 0 is implicit
188 * level 1 balance equation */
189 zz[0][1] = -r01;
190 zz[1][1] = r10 + r12 + r13;
191 zz[2][1] = -r21;
192 zz[3][1] = -r31;
193
194 /* level 2 balance equation */
195 zz[0][2] = -r02;
196 zz[1][2] = -r12;
197 zz[2][2] = r20 + r21 + r23;
198 zz[3][2] = -r32;
199
200 /* level 3 balance equation */
201 zz[0][3] = -r03;
202 zz[1][3] = -r13;
203 zz[2][3] = -r23;
204 zz[3][3] = r30 + r31 + r32;
205
206 /* this one may be more robust */
207 for( j=0; j <= 3; j++ )
208 {
209 for( i=0; i <= 3; i++ )
210 {
211 amat[i][j] = zz[i][j];
212 }
213 bvec[j] = zz[3+1][j];
214 }
215
216 nerror = 0;
217
218 getrf_wrapper(4, 4, (double*)amat, 4, ipiv, &nerror);
219 getrs_wrapper('N', 4, 1, (double*)amat, 4, ipiv, bvec, 4, &nerror);
220
221 if( nerror != 0 )
222 {
223 fprintf( ioQQQ, " AtomSeqBeryllium: dgetrs finds singular or ill-conditioned matrix\n" );
225 }
226 /* now put results back into z so rest of code treates only
227 * one case - as if matin1 had been used */
228 for( i=0; i <= 3; i++ )
229 {
230 zz[3+1][i] = bvec[i];
231 }
232
233 lgNegPop = false;
234 for( i=0; i <= 3; i++ )
235 {
236 atoms.PopLevels[i] = zz[4][i]*AbunxIon;
237 if( atoms.PopLevels[i] < 0. )
238 lgNegPop = true;
239 }
240
241 /* check for negative level populations, this would be a major error */
242 if( lgNegPop )
243 {
244 fprintf( ioQQQ, " AtomSeqBeryllium finds non-positive pop,=" );
245 for( i=0; i <= 3; i++ )
246 {
247 fprintf( ioQQQ, "%g ", atoms.PopLevels[i] );
248 }
249 fprintf( ioQQQ, "%s \n", chLab );
250 fprintf( ioQQQ, " te=%g total abund=%g boltz=%g \n",
251 phycon.te, AbunxIon, boltz );
253 }
254
255 /* convert level populations over to departure coeficients relative to ground */
256 atoms.DepLTELevels[0] = 1.;
257 atoms.DepLTELevels[1] = (atoms.PopLevels[1]/atoms.PopLevels[0])/
258 excit;
259 atoms.DepLTELevels[2] = (atoms.PopLevels[2]/atoms.PopLevels[0])/
260 (excit*3.);
261 atoms.DepLTELevels[3] = (atoms.PopLevels[3]/atoms.PopLevels[0])/
262 (excit*5.);
263
264 /* compute ratio Aul/(Aul+Cul) */
265 t.Emis().ColOvTot() = c02/r02;
266
267 (*t.Lo()).Pop() = AbunxIon;
268
269 /* >>chng 96 sep 12, ipLnPopu had not been set before */
270 (*t.Hi()).Pop() = atoms.PopLevels[2];
271
272 t.Emis().PopOpc() = AbunxIon - atoms.PopLevels[2]*1./3.;
273
274 /* this will be escaping part of line
275 * number of escaping line photons, used elsewhere for outward beam */
276 t.Emis().phots() = atoms.PopLevels[2] * t.Emis().Aul() * ( t.Emis().Pesc() + t.Emis().Pelec_esc() );
277
278 t.Emis().xIntensity() = t.Emis().phots() * t.EnergyErg();
279
280 /* two cases - collisionally excited (usual case) or
281 * radiatively excited - in which case line can be a heat source
282 * following are correct heat exchange, they will mix to get correct deriv
283 * the sum of heat-cool will add up to EnrUL - EnrLU - this is a trick to
284 * keep stable solution by effectively dividing up heating and cooling,
285 * so that negative cooling does not occur */
286 Enr01 = atoms.PopLevels[0]*(c01 + c02 + c03)*t.EnergyErg();
287 Enr10 = (atoms.PopLevels[1]*c10 + atoms.PopLevels[2]*c20 +
288 atoms.PopLevels[3]*c30)*t.EnergyErg();
289
290 /* net cooling due to excit minus part of de-excit */
291 t.Coll().cool() = Enr01 - Enr10*t.Emis().ColOvTot();
292
293 /* net heating is remainder */
294 t.Coll().heat() = Enr10*(1. - t.Emis().ColOvTot());
295
296 /* put line cooling into cooling stack */
297 coolng = t.Coll().cool();
298 CoolAdd(chLab, t.WLAng() ,coolng);
299
300 /* derivative of cooling function */
301 thermal.dCooldT += coolng*(t.EnergyK()*thermal.tsq1 - thermal.halfte);
302
303 /* >>chhg 03 oct 04, move to RT_OTS */
304 /* destroyed part of line
305 dest = atoms.PopLevels[2]*t.Emis->Aul()*t.Pdest();
306 t.ots() = dest; */
307
308 /* now add to ots field
309 * >>chng 03 oct 03, moved to RT_OTS
310 RT_OTS_AddLine(dest , t.ipCont );*/
311 return;
312}
static double * amat
void AtomSeqBeryllium(double cs12, double cs13, double cs23, const TransitionProxy &t, double a30)
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 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 & 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
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_thermal thermal
Definition thermal.cpp:5
void getrs_wrapper(char trans, long N, long nrhs, double *A, long lda, int32 *ipiv, double *B, long ldb, int32 *info)
void getrf_wrapper(long M, long N, double *A, long lda, int32 *ipiv, int32 *info)
void chIonLbl(char *chIonLbl_v, const TransitionProxy &t)