cloudy trunk
Loading...
Searching...
No Matches
atom_oi.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_oi drive the solution of OI level populations, Ly-beta pumping */
4/*oi_level_pops get OI level population with Ly-beta pumping */
5#include "cddefines.h"
6#include "taulines.h"
7#include "doppvel.h"
8#include "iso.h"
9#include "trace.h"
10#include "dense.h"
11#include "rt.h"
12#include "rfield.h"
13#include "phycon.h"
14#include "lines_service.h"
15#include "thirdparty.h"
16#include "atoms.h"
17
18/*oi_level_pops get OI level population with Ly-beta pumping */
19STATIC void oi_level_pops(double abundoi,
20 double *coloi);
21
22/*atom_oi drive the solution of OI level populations, Ly-beta pumping */
23void atom_oi_calc(double *coloi)
24{
25 double esab,
26 eslb,
27 esoi,
28 flb,
29 foi,
30 opaclb,
31 opacoi,
32 xlb,
33 xoi;
34 double aoi = TauLines[ipTO1025].Emis().Aul();
35 double alb = iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH3p,ipH1s).Emis().Aul();
36
37 DEBUG_ENTRY( "atom_oi_calc()" );
38
39 fixit(); // ticket #78 refers
40 // The code below should be calculating the O I 1025 pumping by H Ly beta, as well as
41 // the inverse process (this can become important in hydrogen-deficient environments).
42 // It now uses the Elitzur & Netzer (1985, ApJ, 291, 464) theory, which is no longer
43 // valid since the line overlap code prevents us from getting at the escape probability
44 // of individual lines.
45
46 /* A's from Pradhan; OI pump line; Ly beta, 8446 */
47
48 /* called by LINES before calc really start, protect here
49 * also used for cases where OI not present */
50 if( dense.xIonDense[ipOXYGEN][0] <= 0. )
51 {
52 for( int i=0; i < 6; i++ )
53 {
54 atoms.popoi[i] = 0.;
55 }
56 *coloi = 0.;
57 atoms.pmpo15 = 0.;
58 atoms.pmpo51 = 0.;
59 return;
60 }
61
62 // line overlap code makes this the escape probability of the combined lines
63 esab = TauLines[ipTO1025].Emis().Pelec_esc() + TauLines[ipTO1025].Emis().Pesc();
64
65 // these two are no longer correct, the line overlap code makes it impossible
66 // to get at the escape probabilities of the individual lines...
67 esoi = TauLines[ipTO1025].Emis().Pelec_esc() + TauLines[ipTO1025].Emis().Pesc();
68 eslb = iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH3p,ipH1s).Emis().Pelec_esc() +
69 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH3p,ipH1s).Emis().Pesc();
70
71 /* all trace output turned on with "trace ly beta command' */
72 if( trace.lgTr8446 && trace.lgTrace )
73 {
74 fprintf( ioQQQ,
75 " P8446 finds Lbeta, OI widths=%10.2e%10.2e and esc prob=%10.2e%10.2e esAB=%10.2e\n",
76 GetDopplerWidth(dense.AtomicWeight[ipHYDROGEN]), GetDopplerWidth(dense.AtomicWeight[ipOXYGEN]), eslb, esoi, esab );
77 }
78
79 /* find relative opacities for two lines */
80 opacoi = 2.92e-9*dense.xIonDense[ipOXYGEN][0]*0.5556/GetDopplerWidth(dense.AtomicWeight[ipOXYGEN]);
81 opaclb = 1.22e-8*iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH1s].Pop()/GetDopplerWidth(dense.AtomicWeight[ipHYDROGEN]);
82
83 /* these are x sub a (OI) and x sub b (ly beta) defined in Elitz+Netz */
84 xoi = opacoi/(opacoi + opaclb);
85 xlb = opaclb/(opacoi + opaclb);
86
87 /* find relative line-widths, assume same rest freq */
88 foi = MIN2(GetDopplerWidth(dense.AtomicWeight[ipHYDROGEN]),GetDopplerWidth(dense.AtomicWeight[ipOXYGEN]))/GetDopplerWidth(dense.AtomicWeight[ipOXYGEN]);
89 flb = MIN2(GetDopplerWidth(dense.AtomicWeight[ipHYDROGEN]),GetDopplerWidth(dense.AtomicWeight[ipOXYGEN]))/GetDopplerWidth(dense.AtomicWeight[ipHYDROGEN])*
90 MAX2(0.,1.- TauLines[ipTO1025].Emis().Pelec_esc() - TauLines[ipTO1025].Emis().Pesc());
91
92 if( trace.lgTr8446 && trace.lgTrace )
93 {
94 fprintf( ioQQQ,
95 " P8446 opac Lb, OI=%10.2e%10.2e X Lb, OI=%10.2e%10.2e FLb, OI=%10.2e%10.2e\n",
96 opaclb, opacoi, xlb, xoi, flb, foi );
97 }
98
99 /* pumping of OI by L-beta - this goes into OI matrix as 1-5 rate
100 * lgInducProcess set false with no induced command, usually true */
101 /* Notation: pmpo15 net rate oxygen level 5 populated from 1 */
102 if( rfield.lgInducProcess )
103 {
104 atoms.pmpo15 = (realnum)((flb*alb*iso_sp[ipH_LIKE][ipHYDROGEN].st[ipH3p].Pop()*
105 xoi*(1. - esab)/dense.xIonDense[ipOXYGEN][0]));
106 /* net decay rate from upper level */
107 atoms.pmpo51 = (realnum)(aoi*(1. - (1. - foi)*(1. - esoi) - xoi*(1. - esab)*foi));
108 }
109 else
110 {
111 atoms.pmpo15 = 0.;
112 atoms.pmpo51 = 0.;
113 }
114
115 /* find level populations for OI */
116 oi_level_pops(dense.xIonDense[ipOXYGEN][0],coloi);
117
118 /* continuum pumping due to J=1, 0 sub states.
119 * neglect J=2 since L-Beta very optically thick */
120
121 /* lower level populations */
122 (*TauLines[ipT1304].Lo()).Pop() = atoms.popoi[0];
123 (*TauLines[ipTO1025].Lo()).Pop() = atoms.popoi[0];
124 (*TauLines[ipT1039].Lo()).Pop() = atoms.popoi[0];
125 (*TauLines[ipT8446].Lo()).Pop() = atoms.popoi[1];
126 (*TauLines[ipT4368].Lo()).Pop() = atoms.popoi[1];
127 (*TauLines[ipTOI13].Lo()).Pop() = atoms.popoi[2];
128 (*TauLines[ipTOI11].Lo()).Pop() = atoms.popoi[2];
129 (*TauLines[ipTOI29].Lo()).Pop() = atoms.popoi[3];
130 (*TauLines[ipTOI46].Lo()).Pop() = atoms.popoi[4];
131
132 /* upper level populations */
133 (*TauLines[ipT1304].Hi()).Pop() = atoms.popoi[1];
134 (*TauLines[ipTO1025].Hi()).Pop() = atoms.popoi[4];
135 (*TauLines[ipT1039].Hi()).Pop() = atoms.popoi[3];
136 (*TauLines[ipT8446].Hi()).Pop() = atoms.popoi[2];
137 (*TauLines[ipT4368].Hi()).Pop() = atoms.popoi[5];
138 (*TauLines[ipTOI13].Hi()).Pop() = atoms.popoi[3];
139 (*TauLines[ipTOI11].Hi()).Pop() = atoms.popoi[4];
140 (*TauLines[ipTOI29].Hi()).Pop() = atoms.popoi[5];
141 (*TauLines[ipTOI46].Hi()).Pop() = atoms.popoi[5];
142
143 TauLines[ipT1304].Emis().PopOpc() = atoms.popoi[0] - atoms.popoi[1]*3.0;
144 TauLines[ipTO1025].Emis().PopOpc() = (atoms.popoi[0] - atoms.popoi[4]*0.6);
145 TauLines[ipT1039].Emis().PopOpc() = (atoms.popoi[0] - atoms.popoi[3]*3.0);
146
147 /* t8446(ipLnEmis.PopOpc()) = (popoi(2)-popoi(3)*.33) */
149 TauLines[ipT8446].Emis().PopOpc() =
150 (MAX2(0.,atoms.popoi[1]-atoms.popoi[2]*.33));
151 TauLines[ipT4368].Emis().PopOpc() = (atoms.popoi[1] - atoms.popoi[5]*.33);
152 TauLines[ipTOI13].Emis().PopOpc() = (atoms.popoi[2] - atoms.popoi[3]*3.0);
153 TauLines[ipTOI11].Emis().PopOpc() = (atoms.popoi[2] - atoms.popoi[4]*0.6);
154 TauLines[ipTOI29].Emis().PopOpc() = (atoms.popoi[3] - atoms.popoi[5]*.33);
155 TauLines[ipTOI46].Emis().PopOpc() = (atoms.popoi[4] - atoms.popoi[5]*1.67);
156 return;
157}
158
159/*oi_level_pops get OI level population with Ly-beta pumping */
160STATIC void oi_level_pops(double abundoi,
161 double *coloi)
162{
163 bool lgNegPop;
164
165 long int i, j;
166
167 int32 ipiv[6], ner;
168
169 double a21,
170 a32,
171 a41,
172 a43,
173 a51,
174 a53,
175 a62,
176 a64,
177 a65,
178 c12,
179 c13,
180 c14,
181 c15,
182 c16,
183 c21,
184 c23,
185 c24,
186 c25,
187 c26,
188 c31,
189 c32,
190 c34,
191 c35,
192 c36,
193 c41,
194 c42,
195 c43,
196 c45,
197 c46,
198 c51,
199 c52,
200 c53,
201 c54,
202 c56,
203 c61,
204 c62,
205 c63,
206 c64,
207 c65,
208 cs,
209 deptoi[6],
210 e12,
211 e23,
212 e34,
213 e45,
214 e56,
215 simple;
216
217 double amat[6][6],
218 bvec[6],
219 zz[7][7];
220
221 static double g[6]={9.,3.,9.,3.,15.,9};
222
223 DEBUG_ENTRY( "oilevl()" );
224
225 /* following used for linpac matrix inversion */
226
227 /* compute emission from six level OI atom*/
228
229 /* Boltzmann factors for ContBoltz since collisions not dominant in UV tran
230 * ipoiex is array lof dEnergy for each level, set in DoPoint */
231 e12 = rfield.ContBoltz[atoms.ipoiex[0]-1];
232 e23 = rfield.ContBoltz[atoms.ipoiex[1]-1];
233 e34 = rfield.ContBoltz[atoms.ipoiex[2]-1];
234 e45 = rfield.ContBoltz[atoms.ipoiex[3]-1];
235 e56 = rfield.ContBoltz[atoms.ipoiex[4]-1];
236
237 /* total rad rates here have dest by background continuum */
238 a21 = TauLines[ipT1304].Emis().Aul()*(TauLines[ipT1304].Emis().Pdest()+ TauLines[ipT1304].Emis().Pesc() + TauLines[ipT1304].Emis().Pelec_esc());
239 a41 = TauLines[ipT1039].Emis().Aul()*(TauLines[ipT1039].Emis().Pdest()+ TauLines[ipT1039].Emis().Pesc() + TauLines[ipT1039].Emis().Pelec_esc());
240 a51 = TauLines[ipTO1025].Emis().Aul()*(TauLines[ipTO1025].Emis().Pdest()+ TauLines[ipTO1025].Emis().Pesc() + TauLines[ipTO1025].Emis().Pelec_esc());
241 a51 = atoms.pmpo51;
242 a32 = TauLines[ipT8446].Emis().Aul()*(TauLines[ipT8446].Emis().Pdest()+ TauLines[ipT8446].Emis().Pesc() + TauLines[ipT8446].Emis().Pelec_esc());
243 a62 = TauLines[ipT4368].Emis().Aul()*(TauLines[ipT4368].Emis().Pdest()+ TauLines[ipT4368].Emis().Pesc() + TauLines[ipT4368].Emis().Pelec_esc());
244 a43 = TauLines[ipTOI13].Emis().Aul()*(TauLines[ipTOI13].Emis().Pdest()+ TauLines[ipTOI13].Emis().Pesc() + TauLines[ipTOI13].Emis().Pelec_esc());
245 a53 = TauLines[ipTOI11].Emis().Aul()*(TauLines[ipTOI11].Emis().Pdest()+ TauLines[ipTOI11].Emis().Pesc() + TauLines[ipTOI11].Emis().Pelec_esc());
246 a64 = TauLines[ipTOI29].Emis().Aul()*(TauLines[ipTOI29].Emis().Pdest()+ TauLines[ipTOI29].Emis().Pesc() + TauLines[ipTOI29].Emis().Pelec_esc());
247 a65 = TauLines[ipTOI46].Emis().Aul()*(TauLines[ipTOI46].Emis().Pdest()+ TauLines[ipTOI46].Emis().Pesc() + TauLines[ipTOI46].Emis().Pelec_esc());
248
249 /* even at density of 10^17 excited states not in lte due
250 * to fast transitions down - just need to kill a21 to get to unity at 10^17*/
251
252 /* the 2-1 transition is 1302, cs Wang and McConkey '92 Jphys B 25, 5461 */
253 cs = 2.151e-5*phycon.te/phycon.te03;
255
256 /* the 5-1 transition is 1027, cs Wang and McConkey '92 Jphys B 25, 5461 */
257 cs = 9.25e-7*phycon.te*phycon.te10/phycon.te01/phycon.te01;
259 c21 = dense.cdsqte*TauLines[ipT1304].Coll().col_str()/g[1];
260 c51 = dense.cdsqte*TauLines[ipTO1025].Coll().col_str()/g[4];
261
262 /* all following are g-bar approx, g-bar = 0.2 */
263 c31 = dense.cdsqte*1.0/g[2];
264 PutCS(0.27,TauLines[ipT1039]);
265 c41 = dense.cdsqte*TauLines[ipT1039].Coll().col_str()/g[3];
266 c61 = dense.cdsqte*1./g[5];
267
268 c12 = c21*g[1]/g[0]*e12;
269 c13 = c31*g[2]/g[0]*e12*e23;
270 c14 = c41*g[3]/g[0]*e12*e23*e34;
271 c15 = c51*g[4]/g[0]*e12*e23*e34*e45;
272 c16 = c61*g[5]/g[0]*e12*e23*e34*e45*e56;
273
274 c32 = dense.cdsqte*85./g[2];
275 c42 = dense.cdsqte*85./g[3];
276 c52 = dense.cdsqte*85./g[4];
277 c62 = dense.cdsqte*85./g[5];
278
279 c23 = c32*g[2]/g[1]*e23;
280 c24 = c42*g[3]/g[1]*e23*e34;
281 c25 = c52*g[4]/g[1]*e23*e34*e45;
282 c26 = c62*g[5]/g[1]*e23*e34*e45*e56;
283
284 c43 = dense.cdsqte*70./g[3];
285 c53 = dense.cdsqte*312./g[4];
286 c63 = dense.cdsqte*1./g[5];
287
288 c34 = c43*g[3]/g[2]*e34;
289 c35 = c53*g[4]/g[2]*e34*e45;
290 c36 = c63*g[5]/g[2]*e34*e45*e56;
291
292 c54 = dense.cdsqte*50./g[4];
293 c64 = dense.cdsqte*415./g[5];
294
295 c45 = c54*g[4]/g[3]*e45;
296 c46 = c64*g[5]/g[3]*e45*e56;
297
298 c65 = dense.cdsqte*400./g[5];
299 c56 = c65*g[5]/g[4]*e56;
300
302
303 /* this is check for whether matrix inversion likely to fail */
304 simple = (c16 + atoms.pmpo15)/(c61 + c62 + c64 + a65 + a64 + a62);
305 if( simple < 1e-19 )
306 {
307 atoms.popoi[0] = abundoi;
308 for( i=1; i < 6; i++ )
309 {
310 atoms.popoi[i] = 0.;
311 }
312 *coloi = 0.;
313 return;
314 }
315
316 /*--------------------------------------------------------- */
317
318 for( i=0; i < 6; i++ )
319 {
320 zz[i][0] = 1.0;
321 zz[6][i] = 0.;
322 }
323
324 /* first equation is sum of populations */
325 zz[6][0] = abundoi;
326
327 /* level two, 3s 3So */
328 zz[0][1] = -c12;
329 zz[1][1] = c21 + c23 + c24 + c25 + c26 + a21;
330 zz[2][1] = -c32 - a32;
331 zz[3][1] = -c42;
332 zz[4][1] = -c52;
333 zz[5][1] = -c62 - a62;
334
335 /* level three */
336 zz[0][2] = -c13;
337 zz[1][2] = -c23;
338 zz[2][2] = c31 + c32 + c34 + c35 + c36 + a32;
339 zz[3][2] = -c43 - a43;
340 zz[4][2] = -c53 - a53;
341 zz[5][2] = -c63;
342
343 /* level four */
344 zz[0][3] = -c14;
345 zz[1][3] = -c24;
346 zz[2][3] = -c34;
347 zz[3][3] = c41 + c42 + c43 + c45 + c46 + a41 + a43;
348 zz[4][3] = -c54;
349 zz[5][3] = -c64 - a64;
350
351 /* level five */
352 zz[0][4] = -c15 - atoms.pmpo15;
353 zz[1][4] = -c25;
354 zz[2][4] = -c35;
355 zz[3][4] = -c45;
356 zz[4][4] = c51 + c52 + c53 + c54 + c56 + a51 + a53;
357 zz[5][4] = -c65 - a65;
358
359 /* level six */
360 zz[0][5] = -c16;
361 zz[1][5] = -c26;
362 zz[2][5] = -c36;
363 zz[3][5] = -c46;
364 zz[4][5] = -c56;
365 zz[5][5] = c61 + c62 + c63 + c64 + c65 + a65 + a64 + a62;
366
367 /* this one may be more robust */
368 for( j=0; j < 6; j++ )
369 {
370 for( i=0; i < 6; i++ )
371 {
372 amat[i][j] = zz[i][j];
373 }
374 bvec[j] = zz[6][j];
375 }
376
377 ner = 0;
378
379 getrf_wrapper(6, 6, (double*)amat, 6, ipiv, &ner);
380 getrs_wrapper('N', 6, 1, (double*)amat, 6, ipiv, bvec, 6, &ner);
381
382 /*DGETRF(6,6,(double*)amat,6,ipiv,&ner);*/
383 /*DGETRS('N',6,1,(double*)amat,6,ipiv,bvec,6,&ner);*/
384 if( ner != 0 )
385 {
386 fprintf( ioQQQ, " oi_level_pops: dgetrs finds singular or ill-conditioned matrix\n" );
388 }
389
390 /* now put results back into z so rest of code treates only
391 * one case - as if matin1 had been used */
392 for( i=0; i < 6; i++ )
393 {
394 zz[6][i] = bvec[i];
395 }
396
397 lgNegPop = false;
398 for( i=0; i < 6; i++ )
399 {
400 atoms.popoi[i] = zz[6][i];
401 if( atoms.popoi[i] < 0. )
402 lgNegPop = true;
403 }
404
405 /* following used to confirm that all dep coef are unity at
406 * density of 1e17, t=10,000, and all A's set to zero */
407 if( trace.lgTrace && trace.lgTr8446 )
408 {
409 deptoi[0] = 1.;
410 deptoi[1] = atoms.popoi[1]/atoms.popoi[0]/(g[1]/g[0]*
411 e12);
412 deptoi[2] = atoms.popoi[2]/atoms.popoi[0]/(g[2]/g[0]*
413 e12*e23);
414 deptoi[3] = atoms.popoi[3]/atoms.popoi[0]/(g[3]/g[0]*
415 e12*e23*e34);
416 deptoi[4] = atoms.popoi[4]/atoms.popoi[0]/(g[4]/g[0]*
417 e12*e23*e34*e45);
418 deptoi[5] = atoms.popoi[5]/atoms.popoi[0]/(g[5]/g[0]*
419 e12*e23*e34*e45*e56);
420
421 fprintf( ioQQQ, " oilevl finds levl pop" );
422 for(i=0; i < 6; i++)
423 fprintf( ioQQQ, "%11.3e", atoms.popoi[i] );
424 fprintf( ioQQQ, "\n" );
425
426 fprintf( ioQQQ, " oilevl finds dep coef" );
427 for(i=0; i < 6; i++)
428 fprintf( ioQQQ, "%11.3e", deptoi[i] );
429 fprintf( ioQQQ, "\n" );
430 }
431
432 /* this happens due to numerical instability in matrix inversion routine */
433 if( lgNegPop )
434 {
435 atoms.nNegOI += 1;
436
437 fprintf( ioQQQ, " OILEVL finds negative population" );
438 for(i=0;i < 6; i++)
439 fprintf( ioQQQ, "%10.2e", atoms.popoi[i] );
440 fprintf( ioQQQ, "\n" );
441
442 fprintf( ioQQQ, " simple 5 =%10.2e\n", simple );
443
444 atoms.popoi[5] = 0.;
445 atoms.popoi[4] = abundoi*(c15 + atoms.pmpo15)/(a51 + a53 + c51 + c53);
446 atoms.popoi[3] = 0.;
447 atoms.popoi[2] = (atoms.popoi[4]*(a53 + c53)) / (a32 + c32);
448 atoms.popoi[1] = (atoms.popoi[4]*(a53 + c53) + abundoi*c12) / (a21 + c21);
449 atoms.popoi[0] = abundoi;
450 /* write(QQ,'('' OILEVL resets this to simple pop'',1P,6E10.2)')
451 * 1 popoi */
452 }
453
454 /* this is total cooling due to model atom, can be neg (heating) */
455 *coloi =
456 (atoms.popoi[0]*c12 - atoms.popoi[1]*c21)*1.53e-11 +
457 (atoms.popoi[0]*c14 - atoms.popoi[3]*c41)*1.92e-11 +
458 (atoms.popoi[0]*c15 - atoms.popoi[4]*c51)*1.94e-11 +
459 (atoms.popoi[1]*c23 - atoms.popoi[2]*c32)*2.36e-12 +
460 (atoms.popoi[1]*c26 - atoms.popoi[5]*c62)*4.55e-12 +
461 (atoms.popoi[2]*c35 - atoms.popoi[4]*c53)*1.76e-12 +
462 (atoms.popoi[2]*c34 - atoms.popoi[3]*c43)*1.52e-12 +
463 (atoms.popoi[3]*c46 - atoms.popoi[5]*c64)*6.86e-13 +
464 (atoms.popoi[4]*c56 - atoms.popoi[5]*c65)*4.32e-13;
465
466 return;
467}
long ipTOI46
long ipTO1025
long ipTOI11
long ipTOI13
long ipT8446
long ipT4368
long ipT1304
long ipTOI29
long ipT1039
static double * amat
STATIC void oi_level_pops(double abundoi, double *coloi)
Definition atom_oi.cpp:160
void atom_oi_calc(double *coloi)
Definition atom_oi.cpp:23
t_atoms atoms
Definition atoms.cpp:5
FILE * ioQQQ
Definition cddefines.cpp:7
#define STATIC
Definition cddefines.h:97
#define MIN2
Definition cddefines.h:761
const int ipOXYGEN
Definition cddefines.h:312
#define EXIT_FAILURE
Definition cddefines.h:140
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
const int ipHYDROGEN
Definition cddefines.h:305
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void fixit(void)
Definition service.cpp:991
t_dense dense
Definition dense.cpp:24
realnum GetDopplerWidth(realnum massAMU)
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
const int ipH3p
Definition iso.h:31
const int ipH1s
Definition iso.h:27
const int ipH_LIKE
Definition iso.h:62
t_phycon phycon
Definition phycon.cpp:6
t_rfield rfield
Definition rfield.cpp:8
static double * g
Definition species2.cpp:28
TransitionList TauLines("TauLines", &AnonStates)
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)
t_trace trace
Definition trace.cpp:5
void PutCS(double cs, const TransitionProxy &t)