root/branches/newmole/source/atom_oi.cpp

Revision 2542, 14.2 kB (checked in by rjrw, 3 weeks ago)

Merged from trunk r2448:2541

  • Property svn:eol-style set to native
Line 
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_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 = Transitions[ipH_LIKE][ipHYDROGEN][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 = Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pelec_esc +
69                Transitions[ipH_LIKE][ipHYDROGEN][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                  DoppVel.doppler[0], DoppVel.doppler[7], eslb, esoi, esab );
77        }
78
79        /* find relative opacities for two lines */
80        opacoi = 2.92e-9*dense.xIonDense[ipOXYGEN][0]*0.5556/DoppVel.doppler[ipOXYGEN];
81        opaclb = 1.22e-8*dense.xIonDense[ipHYDROGEN][1]*StatesElem[ipH_LIKE][ipHYDROGEN][ipH1s].Pop/DoppVel.doppler[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(DoppVel.doppler[ipHYDROGEN],DoppVel.doppler[ipOXYGEN])/DoppVel.doppler[ipOXYGEN];
89        flb = MIN2(DoppVel.doppler[ipHYDROGEN],DoppVel.doppler[ipOXYGEN])/DoppVel.doppler[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*dense.xIonDense[ipHYDROGEN][1]*
105                                         StatesElem[ipH_LIKE][ipHYDROGEN][ipH3p].Pop*
106                                         xoi*(1. - esab)/dense.xIonDense[ipOXYGEN][0]);
107                /* net decay rate from upper level */
108                atoms.pmpo51 = (realnum)(aoi*(1. - (1. - foi)*(1. - esoi) - xoi*(1. - esab)*foi));
109        }
110        else
111        {
112                atoms.pmpo15 = 0.;
113                atoms.pmpo51 = 0.;
114        }
115
116        /* find level populations for OI */
117        oi_level_pops(dense.xIonDense[ipOXYGEN][0],coloi);
118
119        /* continuum pumping due to J=1, 0 sub states.
120         * neglect J=2 since L-Beta very optically thick */
121
122        /* lower level populations */
123        TauLines[ipT1304].Lo->Pop = atoms.popoi[0];
124        TauLines[ipTO1025].Lo->Pop = atoms.popoi[0];
125        TauLines[ipT1039].Lo->Pop = atoms.popoi[0];
126        TauLines[ipT8446].Lo->Pop = atoms.popoi[1];
127        TauLines[ipT4368].Lo->Pop = atoms.popoi[1];
128        TauLines[ipTOI13].Lo->Pop = atoms.popoi[2];
129        TauLines[ipTOI11].Lo->Pop = atoms.popoi[2];
130        TauLines[ipTOI29].Lo->Pop = atoms.popoi[3];
131        TauLines[ipTOI46].Lo->Pop = atoms.popoi[4];
132
133        /* upper level populations */
134        TauLines[ipT1304].Hi->Pop = atoms.popoi[1];
135        TauLines[ipTO1025].Hi->Pop = atoms.popoi[4];
136        TauLines[ipT1039].Hi->Pop = atoms.popoi[3];
137        TauLines[ipT8446].Hi->Pop = atoms.popoi[2];
138        TauLines[ipT4368].Hi->Pop = atoms.popoi[5];
139        TauLines[ipTOI13].Hi->Pop = atoms.popoi[3];
140        TauLines[ipTOI11].Hi->Pop = atoms.popoi[4];
141        TauLines[ipTOI29].Hi->Pop = atoms.popoi[5];
142        TauLines[ipTOI46].Hi->Pop = atoms.popoi[5];
143
144        /* opacity factor
145         * t1304(ipLnEmis.PopOpc) = (popoi(1)-popoi(2)*3.0) */
146        /** \todo       2       following needed to get badbugs/bug8.in to work */
147        TauLines[ipT1304].Emis->PopOpc = atoms.popoi[0];
148        TauLines[ipTO1025].Emis->PopOpc = (atoms.popoi[0] - atoms.popoi[4]*0.6);
149        TauLines[ipT1039].Emis->PopOpc = (atoms.popoi[0] - atoms.popoi[3]*3.0);
150
151        /* t8446(ipLnEmis.PopOpc) = (popoi(2)-popoi(3)*.33) */
152        /** \todo       2       following needed to get badbugs/bug5.in to work */
153        TauLines[ipT8446].Emis->PopOpc =
154                (MAX2(0.,atoms.popoi[1]-atoms.popoi[2]*.33));
155        TauLines[ipT4368].Emis->PopOpc = (atoms.popoi[1] - atoms.popoi[5]*.33);
156        TauLines[ipTOI13].Emis->PopOpc = (atoms.popoi[2] - atoms.popoi[3]*3.0);
157        TauLines[ipTOI11].Emis->PopOpc = (atoms.popoi[2] - atoms.popoi[4]*0.6);
158        TauLines[ipTOI29].Emis->PopOpc = (atoms.popoi[3] - atoms.popoi[5]*.33);
159        TauLines[ipTOI46].Emis->PopOpc = (atoms.popoi[4] - atoms.popoi[5]*1.67);
160        return;
161}
162
163/*oi_level_pops get OI level population with Ly-beta pumping */
164STATIC void oi_level_pops(double abundoi,
165                          double *coloi)
166{
167        bool lgNegPop;
168
169        long int i, j;
170
171        int32 ipiv[6], ner;
172
173        double a21,
174          a32,
175          a41,
176          a43,
177          a51,
178          a53,
179          a62,
180          a64,
181          a65,
182          c12,
183          c13,
184          c14,
185          c15,
186          c16,
187          c21,
188          c23,
189          c24,
190          c25,
191          c26,
192          c31,
193          c32,
194          c34,
195          c35,
196          c36,
197          c41,
198          c42,
199          c43,
200          c45,
201          c46,
202          c51,
203          c52,
204          c53,
205          c54,
206          c56,
207          c61,
208          c62,
209          c63,
210          c64,
211          c65,
212          cs,
213          deptoi[6],
214          e12,
215          e23,
216          e34,
217          e45,
218          e56,
219          simple;
220
221        double amat[6][6],
222          bvec[6],
223          zz[7][7];
224
225        static realnum g[6]={9.,3.,9.,3.,15.,9};
226
227        DEBUG_ENTRY( "oilevl()" );
228
229        /* following used for linpac matrix inversion */
230
231        /* compute emission from six level OI atom*/
232
233        /* Boltzmann factors for ContBoltz since collisions not dominant in UV tran
234         * ipoiex is array lof dEnergy for each level, set in DoPoint */
235        e12 = rfield.ContBoltz[atoms.ipoiex[0]-1];
236        e23 = rfield.ContBoltz[atoms.ipoiex[1]-1];
237        e34 = rfield.ContBoltz[atoms.ipoiex[2]-1];
238        e45 = rfield.ContBoltz[atoms.ipoiex[3]-1];
239        e56 = rfield.ContBoltz[atoms.ipoiex[4]-1];
240
241        /* total rad rates here have dest by background continuum */
242        a21 = TauLines[ipT1304].Emis->Aul*(TauLines[ipT1304].Emis->Pdest+ TauLines[ipT1304].Emis->Pesc + TauLines[ipT1304].Emis->Pelec_esc);
243        a41 = TauLines[ipT1039].Emis->Aul*(TauLines[ipT1039].Emis->Pdest+ TauLines[ipT1039].Emis->Pesc + TauLines[ipT1039].Emis->Pelec_esc);
244        a51 = TauLines[ipTO1025].Emis->Aul*(TauLines[ipTO1025].Emis->Pdest+ TauLines[ipTO1025].Emis->Pesc + TauLines[ipTO1025].Emis->Pelec_esc);
245        a51 = atoms.pmpo51;
246        a32 = TauLines[ipT8446].Emis->Aul*(TauLines[ipT8446].Emis->Pdest+ TauLines[ipT8446].Emis->Pesc + TauLines[ipT8446].Emis->Pelec_esc);
247        a62 = TauLines[ipT4368].Emis->Aul*(TauLines[ipT4368].Emis->Pdest+ TauLines[ipT4368].Emis->Pesc + TauLines[ipT4368].Emis->Pelec_esc);
248        a43 = TauLines[ipTOI13].Emis->Aul*(TauLines[ipTOI13].Emis->Pdest+ TauLines[ipTOI13].Emis->Pesc + TauLines[ipTOI13].Emis->Pelec_esc);
249        a53 = TauLines[ipTOI11].Emis->Aul*(TauLines[ipTOI11].Emis->Pdest+ TauLines[ipTOI11].Emis->Pesc + TauLines[ipTOI11].Emis->Pelec_esc);
250        a64 = TauLines[ipTOI29].Emis->Aul*(TauLines[ipTOI29].Emis->Pdest+ TauLines[ipTOI29].Emis->Pesc + TauLines[ipTOI29].Emis->Pelec_esc);
251        a65 = TauLines[ipTOI46].Emis->Aul*(TauLines[ipTOI46].Emis->Pdest+ TauLines[ipTOI46].Emis->Pesc + TauLines[ipTOI46].Emis->Pelec_esc);
252
253        /* even at density of 10^17 excited states not in lte due
254         * to fast transitions down - just need to kill a21 to get to unity at 10^17*/
255
256        /* the 2-1 transition is 1302, cs Wang and McConkey '92 Jphys B 25, 5461 */
257        cs = 2.151e-5*phycon.te/phycon.te03;
258        PutCS(cs,&TauLines[ipT1304]);
259
260        /* the 5-1 transition is 1027, cs Wang and McConkey '92 Jphys B 25, 5461 */
261        cs = 9.25e-7*phycon.te*phycon.te10/phycon.te01/phycon.te01;
262        PutCS(cs,&TauLines[ipTO1025]);
263        c21 = dense.cdsqte*TauLines[ipT1304].Coll.col_str/g[1];
264        c51 = dense.cdsqte*TauLines[ipTO1025].Coll.col_str/g[4];
265
266        /* all following are g-bar approx, g-bar = 0.2 */
267        c31 = dense.cdsqte*1.0/g[2];
268        PutCS(0.27,&TauLines[ipT1039]);
269        c41 = dense.cdsqte*TauLines[ipT1039].Coll.col_str/g[3];
270        c61 = dense.cdsqte*1./g[5];
271
272        c12 = c21*g[1]/g[0]*e12;
273        c13 = c31*g[2]/g[0]*e12*e23;
274        c14 = c41*g[3]/g[0]*e12*e23*e34;
275        c15 = c51*g[4]/g[0]*e12*e23*e34*e45;
276        c16 = c61*g[5]/g[0]*e12*e23*e34*e45*e56;
277
278        c32 = dense.cdsqte*85./g[2];
279        c42 = dense.cdsqte*85./g[3];
280        c52 = dense.cdsqte*85./g[4];
281        c62 = dense.cdsqte*85./g[5];
282
283        c23 = c32*g[2]/g[1]*e23;
284        c24 = c42*g[3]/g[1]*e23*e34;
285        c25 = c52*g[4]/g[1]*e23*e34*e45;
286        c26 = c62*g[5]/g[1]*e23*e34*e45*e56;
287
288        c43 = dense.cdsqte*70./g[3];
289        c53 = dense.cdsqte*312./g[4];
290        c63 = dense.cdsqte*1./g[5];
291
292        c34 = c43*g[3]/g[2]*e34;
293        c35 = c53*g[4]/g[2]*e34*e45;
294        c36 = c63*g[5]/g[2]*e34*e45*e56;
295
296        c54 = dense.cdsqte*50./g[4];
297        c64 = dense.cdsqte*415./g[5];
298
299        c45 = c54*g[4]/g[3]*e45;
300        c46 = c64*g[5]/g[3]*e45*e56;
301
302        c65 = dense.cdsqte*400./g[5];
303        c56 = c65*g[5]/g[4]*e56;
304
305        /** \todo       2       this must have all stimulated emission, pump by cont, etc*/
306
307        /* this is check for whether matrix inversion likely to fail */
308        simple = (c16 + atoms.pmpo15)/(c61 + c62 + c64 + a65 + a64 + a62);
309        if( simple < 1e-19 )
310        {
311                atoms.popoi[0] = (realnum)abundoi;
312                for( i=1; i < 6; i++ )
313                {
314                        atoms.popoi[i] = 0.;
315                }
316                *coloi = 0.;
317                return;
318        }
319
320        /*--------------------------------------------------------- */
321
322        for( i=0; i < 6; i++ )
323        {
324                zz[i][0] = 1.0;
325                zz[6][i] = 0.;
326        }
327
328        /* first equation is sum of populations */
329        zz[6][0] = abundoi;
330
331        /* level two, 3s 3So */
332        zz[0][1] = -c12;
333        zz[1][1] = c21 + c23 + c24 + c25 + c26 + a21;
334        zz[2][1] = -c32 - a32;
335        zz[3][1] = -c42;
336        zz[4][1] = -c52;
337        zz[5][1] = -c62 - a62;
338
339        /* level three */
340        zz[0][2] = -c13;
341        zz[1][2] = -c23;
342        zz[2][2] = c31 + c32 + c34 + c35 + c36 + a32;
343        zz[3][2] = -c43 - a43;
344        zz[4][2] = -c53 - a53;
345        zz[5][2] = -c63;
346
347        /* level four */
348        zz[0][3] = -c14;
349        zz[1][3] = -c24;
350        zz[2][3] = -c34;
351        zz[3][3] = c41 + c42 + c43 + c45 + c46 + a41 + a43;
352        zz[4][3] = -c54;
353        zz[5][3] = -c64 - a64;
354
355        /* level five */
356        zz[0][4] = -c15 - atoms.pmpo15;
357        zz[1][4] = -c25;
358        zz[2][4] = -c35;
359        zz[3][4] = -c45;
360        zz[4][4] = c51 + c52 + c53 + c54 + c56 + a51 + a53;
361        zz[5][4] = -c65 - a65;
362
363        /* level six */
364        zz[0][5] = -c16;
365        zz[1][5] = -c26;
366        zz[2][5] = -c36;
367        zz[3][5] = -c46;
368        zz[4][5] = -c56;
369        zz[5][5] = c61 + c62 + c63 + c64 + c65 + a65 + a64 + a62;
370
371        /* this one may be more robust */
372        for( j=0; j < 6; j++ )
373        {
374                for( i=0; i < 6; i++ )
375                {
376                        amat[i][j] = zz[i][j];
377                }
378                bvec[j] = zz[6][j];
379        }
380
381        ner = 0;
382
383        getrf_wrapper(6, 6, (double*)amat, 6, ipiv, &ner);
384        getrs_wrapper('N', 6, 1, (double*)amat, 6, ipiv, bvec, 6, &ner);
385
386        /*DGETRF(6,6,(double*)amat,6,ipiv,&ner);*/
387        /*DGETRS('N',6,1,(double*)amat,6,ipiv,bvec,6,&ner);*/
388        if( ner != 0 )
389        {
390                fprintf( ioQQQ, " oi_level_pops: dgetrs finds singular or ill-conditioned matrix\n" );
391                cdEXIT(EXIT_FAILURE);
392        }
393
394        /* now put results back into z so rest of code treates only
395                * one case - as if matin1 had been used */
396        for( i=0; i < 6; i++ )
397        {
398                zz[6][i] = bvec[i];
399        }
400
401        lgNegPop = false;
402        for( i=0; i < 6; i++ )
403        {
404                atoms.popoi[i] = (realnum)zz[6][i];
405                if( atoms.popoi[i] < 0. )
406                        lgNegPop = true;
407        }
408
409        /* following used to confirm that all dep coef are unity at
410         * density of 1e17, t=10,000, and all A's set to zero */
411        if( trace.lgTrace && trace.lgTr8446 )
412        {
413                deptoi[0] = 1.;
414                deptoi[1] = atoms.popoi[1]/atoms.popoi[0]/(g[1]/g[0]*
415                  e12);
416                deptoi[2] = atoms.popoi[2]/atoms.popoi[0]/(g[2]/g[0]*
417                  e12*e23);
418                deptoi[3] = atoms.popoi[3]/atoms.popoi[0]/(g[3]/g[0]*
419                  e12*e23*e34);
420                deptoi[4] = atoms.popoi[4]/atoms.popoi[0]/(g[4]/g[0]*
421                  e12*e23*e34*e45);
422                deptoi[5] = atoms.popoi[5]/atoms.popoi[0]/(g[5]/g[0]*
423                  e12*e23*e34*e45*e56);
424
425                fprintf( ioQQQ, " oilevl finds levl pop" );
426                for(i=0; i < 6; i++)
427                        fprintf( ioQQQ, "%11.3e", atoms.popoi[i] );
428                fprintf( ioQQQ, "\n" );
429
430                fprintf( ioQQQ, " oilevl finds dep coef" );
431                for(i=0; i < 6; i++)
432                        fprintf( ioQQQ, "%11.3e", deptoi[i] );
433                fprintf( ioQQQ, "\n" );
434        }
435
436        /* this happens due to numerical instability in matrix inversion routine */
437        if( lgNegPop )
438        {
439                atoms.nNegOI += 1;
440
441                fprintf( ioQQQ, " OILEVL finds negative population" );
442                for(i=0;i < 6; i++)
443                        fprintf( ioQQQ, "%10.2e", atoms.popoi[i] );
444                fprintf( ioQQQ, "\n" );
445
446                fprintf( ioQQQ, " simple 5 =%10.2e\n", simple );
447
448                atoms.popoi[5] = 0.;
449                atoms.popoi[4] = (realnum)(abundoi*(c15 + atoms.pmpo15)/(a51 + a53 +
450                  c51 + c53));
451                atoms.popoi[3] = 0.;
452                atoms.popoi[2] = (realnum)(atoms.popoi[4]*(a53 + c53)/(a32 + c32));
453                atoms.popoi[1] = (realnum)((atoms.popoi[2]*(a32 + c32) + abundoi*
454                  c12)/(a21 + c21));
455
456                atoms.popoi[0] = (realnum)abundoi;
457                /*  write(QQ,'('' OILEVL resets this to simple pop'',1P,6E10.2)')
458                 *  1   popoi */
459        }
460
461        /* this is total cooling due to model atom, can be neg (heating) */
462        *coloi =
463           (atoms.popoi[0]*c12 - atoms.popoi[1]*c21)*1.53e-11 +
464           (atoms.popoi[0]*c14 - atoms.popoi[3]*c41)*1.92e-11 +
465           (atoms.popoi[0]*c15 - atoms.popoi[4]*c51)*1.94e-11 +
466           (atoms.popoi[1]*c23 - atoms.popoi[2]*c32)*2.36e-12 +
467           (atoms.popoi[1]*c26 - atoms.popoi[5]*c62)*4.55e-12 +
468           (atoms.popoi[2]*c35 - atoms.popoi[4]*c53)*1.76e-12 +
469           (atoms.popoi[2]*c34 - atoms.popoi[3]*c43)*1.52e-12 +
470           (atoms.popoi[3]*c46 - atoms.popoi[5]*c64)*6.86e-13 +
471           (atoms.popoi[4]*c56 - atoms.popoi[5]*c65)*4.32e-13;
472
473        return;
474}
Note: See TracBrowser for help on using the browser.