cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
atom_seq_beryllium.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 /*AtomSeqBeryllium compute level populations and emissivity for Be-sequence ions */
4 #include "cddefines.h"
5 #include "phycon.h"
6 #include "lines_service.h"
7 #include "thirdparty.h"
8 #include "dense.h"
9 #include "cooling.h"
10 #include "thermal.h"
11 #include "atoms.h"
12 
13 void AtomSeqBeryllium(double cs12,
14  double cs13,
15  double cs23,
16  transition * 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 punch 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.cs;
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 punch lines data command */
96  t->Coll.cs /= 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" );
224  cdEXIT(EXIT_FAILURE);
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 );
252  cdEXIT(EXIT_FAILURE);
253  }
254 
255  /* convert level populations over to departure coeficients relative to ground */
256  atoms.DepLTELevels[0] = 1.;
258  excit;
260  (excit*3.);
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 }

Generated for cloudy by doxygen 1.8.1.2