cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
atom_oi.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 /*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 */
19 STATIC void oi_level_pops(double abundoi,
20  double *coloi);
21 
22 /*atom_oi drive the solution of OI level populations, Ly-beta pumping */
23 void atom_oi_calc(double *coloi)
24 {
25  long int i;
26  double esin,
27  eslb,
28  esoi,
29  flb,
30  foi,
31  opaclb,
32  opacoi,
33  tin,
34  tout,
35  xlb,
36  xoi;
37  static double esab;
38  static double aoi = TauLines[ipTO1025].Emis->Aul;
39  static double alb = Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Aul;
40 
41  DEBUG_ENTRY( "atom_oi_calc()" );
42 
43  /* A's from Pradhan; OI pump line; Ly beta, 8446 */
44 
45  /* Notation;
46  * PMPH31 net rate hydrogen level 3 depopulated to 1
47  * PMPO15 and PMPO51 are similar, but for oxygen
48  * ESCH31 is emission in 1025 transition */
49 
50  /* called by LINES before calc really start, protect here
51  * also used for cases where OI not present */
52  if( dense.xIonDense[ipOXYGEN][0] <= 0. )
53  {
54  for( i=0; i < 6; i++ )
55  {
56  atoms.popoi[i] = 0.;
57  }
58  /* return zero */
59  *coloi = 0.;
60  atoms.pmpo15 = 0.;
61  atoms.pmpo51 = 0.;
63  (Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pesc +
64  Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pelec_esc);
65 
67  (Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pesc +
68  Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pelec_esc);
69 
70  /* all trace output turned on with "trace ly beta command' */
71  if( trace.lgTr8446 && trace.lgTrace )
72  {
73  fprintf( ioQQQ,
74  " P8446 called for first time, finds A*escape prob 3-1 =%10.2e\n",
75  atoms.esch31 );
76  }
77  return;
78  }
79 
80  if( iteration > 1 )
81  {
82  /* two sided escape prob */
83  tin = TauLines[ipTO1025].Emis->TauIn;
84  esin = esc_CRDwing_1side(tin,1e-4);
85  tout = (TauLines[ipTO1025].Emis->TauTot)*0.9 -
86  TauLines[ipTO1025].Emis->TauIn;
87 
88  if( trace.lgTr8446 && trace.lgTrace )
89  {
90  fprintf( ioQQQ, " P8446 tin, tout=%10.2e%10.2e\n", tin, tout );
91  }
92 
93  if( tout > 0. )
94  {
96 
97  /* do not update esab if we overran optical depth scale */
98  esab = 0.5*(esin + esc_CRDwing_1side(tout,1e-4));
99  }
100  }
101  else
102  {
103  /* one-sided escape probability */
104  esab = esc_CRDwing_1side(TauLines[ipTO1025].Emis->TauIn,1e-4);
105  }
106 
108  eslb =Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pelec_esc +
109  Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pesc;
110 
111  /* all trace output turned on with "trace ly beta command' */
112  if( trace.lgTr8446 && trace.lgTrace )
113  {
114  fprintf( ioQQQ,
115  " P8446 finds Lbeta, OI widths=%10.2e%10.2e and esc prob=%10.2e%10.2e esAB=%10.2e\n",
116  DoppVel.doppler[0], DoppVel.doppler[7], eslb, esoi, esab );
117  }
118 
119  /* find relative opacities for two lines */
120  opacoi = 2.92e-9*dense.xIonDense[ipOXYGEN][0]*0.5556/DoppVel.doppler[7];
122 
123  /* these are x sub a (OI) and x sub b (ly beta) defined in Elitz+Netz */
124  xoi = opacoi/(opacoi + opaclb);
125  xlb = opaclb/(opacoi + opaclb);
126 
127  /* find relative line-widths, assume same rest freq */
130  MAX2(0.,1.- TauLines[ipTO1025].Emis->Pelec_esc - TauLines[ipTO1025].Emis->Pesc);
131 
132  if( trace.lgTr8446 && trace.lgTrace )
133  {
134  fprintf( ioQQQ,
135  " P8446 opac Lb, OI=%10.2e%10.2e X Lb, OI=%10.2e%10.2e FLb, OI=%10.2e%10.2e\n",
136  opaclb, opacoi, xlb, xoi, flb, foi );
137  }
138 
139  /* pumping of OI by L-beta - this goes into OI matrix as 1-5 rate
140  * lgInducProcess set false with no induced command, usually true */
141  if( rfield.lgInducProcess )
142  {
143  atoms.pmpo15 = (realnum)(flb*alb*(StatesElem[ipH_LIKE][ipHYDROGEN][ipH3p].Pop*dense.xIonDense[ipHYDROGEN][1])*xoi*(1. - esab)/
144  dense.xIonDense[ipOXYGEN][0]);
145  /* net decay rate from upper level */
146  atoms.pmpo51 = (realnum)(aoi*(1. - (1. - foi)*(1. - esoi) - xoi*(1. -
147  esab)*foi));
148  }
149  else
150  {
151  atoms.pmpo15 = 0.;
152  atoms.pmpo51 = 0.;
153  }
154 
155  /* find level populations for OI */
157 
158  /* escape term from n=3 of H; followed by pump term to n=3 */
160  //atoms.pmph31 = alb*(1. - (1. - flb)*(1. - eslb) - xlb*(1. - esab)*flb);
161 
162  /*pmph13 = xlb*(1. - esab)*foi*aoi*atoms.popoi[4];*/
164  (Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pelec_esc + Transitions[ipH_LIKE][ipHYDROGEN][ipH3p][ipH1s].Emis->Pesc);
165 
166  /* following is actual emission rate, used to predict Ly-beta in lines.for */
167  atoms.esch31 = alb*(eslb*(1. - flb) + esab*flb);
168 
169  /* all trace output turned on with "trace ly beta command' */
170  if( trace.lgTr8446 && trace.lgTrace )
171  {
172  fprintf( ioQQQ, " P8446 PMPH31=%10.2e EscP3-1%10.2e ESCH31=%10.2e\n",
173  atoms.pmph31, atoms.pmph31/alb, atoms.esch31/alb );
174  }
175 
176  /* continuum pumping due to J=1, 0 sub states.
177  * neglect J=2 since L-Beta very optically thick */
178 
179  /* lower level populations */
180  TauLines[ipT1304].Lo->Pop = atoms.popoi[0];
182  TauLines[ipT1039].Lo->Pop = atoms.popoi[0];
183  TauLines[ipT8446].Lo->Pop = atoms.popoi[1];
184  TauLines[ipT4368].Lo->Pop = atoms.popoi[1];
185  TauLines[ipTOI13].Lo->Pop = atoms.popoi[2];
186  TauLines[ipTOI11].Lo->Pop = atoms.popoi[2];
187  TauLines[ipTOI29].Lo->Pop = atoms.popoi[3];
188  TauLines[ipTOI46].Lo->Pop = atoms.popoi[4];
189 
190  /* upper level populations */
191  TauLines[ipT1304].Hi->Pop = atoms.popoi[1];
193  TauLines[ipT1039].Hi->Pop = atoms.popoi[3];
194  TauLines[ipT8446].Hi->Pop = atoms.popoi[2];
195  TauLines[ipT4368].Hi->Pop = atoms.popoi[5];
196  TauLines[ipTOI13].Hi->Pop = atoms.popoi[3];
197  TauLines[ipTOI11].Hi->Pop = atoms.popoi[4];
198  TauLines[ipTOI29].Hi->Pop = atoms.popoi[5];
199  TauLines[ipTOI46].Hi->Pop = atoms.popoi[5];
200 
201  /* opacity factor
202  * t1304(ipLnEmis.PopOpc) = (popoi(1)-popoi(2)*3.0) */
205  TauLines[ipTO1025].Emis->PopOpc = (atoms.popoi[0] - atoms.popoi[4]*0.6);
206  TauLines[ipT1039].Emis->PopOpc = (atoms.popoi[0] - atoms.popoi[3]*3.0);
207 
208  /* t8446(ipLnEmis.PopOpc) = (popoi(2)-popoi(3)*.33) */
211  (MAX2(0.,atoms.popoi[1]-atoms.popoi[2]*.33));
212  TauLines[ipT4368].Emis->PopOpc = (atoms.popoi[1] - atoms.popoi[5]*.33);
213  TauLines[ipTOI13].Emis->PopOpc = (atoms.popoi[2] - atoms.popoi[3]*3.0);
214  TauLines[ipTOI11].Emis->PopOpc = (atoms.popoi[2] - atoms.popoi[4]*0.6);
215  TauLines[ipTOI29].Emis->PopOpc = (atoms.popoi[3] - atoms.popoi[5]*.33);
216  TauLines[ipTOI46].Emis->PopOpc = (atoms.popoi[4] - atoms.popoi[5]*1.67);
217  return;
218 }
219 
220 /*oi_level_pops get OI level population with Ly-beta pumping */
221 STATIC void oi_level_pops(double abundoi,
222  double *coloi)
223 {
224  bool lgNegPop;
225 
226  long int i, j;
227 
228  int32 ipiv[6], ner;
229 
230  double a21,
231  a32,
232  a41,
233  a43,
234  a51,
235  a53,
236  a62,
237  a64,
238  a65,
239  c12,
240  c13,
241  c14,
242  c15,
243  c16,
244  c21,
245  c23,
246  c24,
247  c25,
248  c26,
249  c31,
250  c32,
251  c34,
252  c35,
253  c36,
254  c41,
255  c42,
256  c43,
257  c45,
258  c46,
259  c51,
260  c52,
261  c53,
262  c54,
263  c56,
264  c61,
265  c62,
266  c63,
267  c64,
268  c65,
269  cs,
270  deptoi[6],
271  e12,
272  e23,
273  e34,
274  e45,
275  e56,
276  simple;
277 
278  double amat[6][6],
279  bvec[6],
280  zz[7][7];
281 
282  static realnum g[6]={9.,3.,9.,3.,15.,9};
283 
284  DEBUG_ENTRY( "oilevl()" );
285 
286  /* following used for linpac matrix inversion */
287 
288  /* compute emission from six level OI atom*/
289 
290  /* Boltzmann factors for ContBoltz since collisions not dominant in UV tran
291  * ipoiex is array lof dEnergy for each level, set in DoPoint */
292  e12 = rfield.ContBoltz[atoms.ipoiex[0]-1];
293  e23 = rfield.ContBoltz[atoms.ipoiex[1]-1];
294  e34 = rfield.ContBoltz[atoms.ipoiex[2]-1];
295  e45 = rfield.ContBoltz[atoms.ipoiex[3]-1];
296  e56 = rfield.ContBoltz[atoms.ipoiex[4]-1];
297 
298  /* total rad rates here have dest by background continuum */
302  a51 = atoms.pmpo51;
309 
310  /* even at density of 10^17 excited states not in lte due
311  * to fast transitions down - just need to kill a21 to get to unity at 10^17*/
312 
313  /* the 2-1 transition is 1302, cs Wang and McConkey '92 Jphys B 25, 5461 */
314  cs = 2.151e-5*phycon.te/phycon.te03;
315  PutCS(cs,&TauLines[ipT1304]);
316 
317  /* the 5-1 transition is 1027, cs Wang and McConkey '92 Jphys B 25, 5461 */
318  cs = 9.25e-7*phycon.te*phycon.te10/phycon.te01/phycon.te01;
319  PutCS(cs,&TauLines[ipTO1025]);
320  c21 = dense.cdsqte*TauLines[ipT1304].Coll.cs/g[1];
321  c51 = dense.cdsqte*TauLines[ipTO1025].Coll.cs/g[4];
322 
323  /* all following are g-bar approx, g-bar = 0.2 */
324  c31 = dense.cdsqte*1.0/g[2];
325  PutCS(0.27,&TauLines[ipT1039]);
326  c41 = dense.cdsqte*TauLines[ipT1039].Coll.cs/g[3];
327  c61 = dense.cdsqte*1./g[5];
328 
329  c12 = c21*g[1]/g[0]*e12;
330  c13 = c31*g[2]/g[0]*e12*e23;
331  c14 = c41*g[3]/g[0]*e12*e23*e34;
332  c15 = c51*g[4]/g[0]*e12*e23*e34*e45;
333  c16 = c61*g[5]/g[0]*e12*e23*e34*e45*e56;
334 
335  c32 = dense.cdsqte*85./g[2];
336  c42 = dense.cdsqte*85./g[3];
337  c52 = dense.cdsqte*85./g[4];
338  c62 = dense.cdsqte*85./g[5];
339 
340  c23 = c32*g[2]/g[1]*e23;
341  c24 = c42*g[3]/g[1]*e23*e34;
342  c25 = c52*g[4]/g[1]*e23*e34*e45;
343  c26 = c62*g[5]/g[1]*e23*e34*e45*e56;
344 
345  c43 = dense.cdsqte*70./g[3];
346  c53 = dense.cdsqte*312./g[4];
347  c63 = dense.cdsqte*1./g[5];
348 
349  c34 = c43*g[3]/g[2]*e34;
350  c35 = c53*g[4]/g[2]*e34*e45;
351  c36 = c63*g[5]/g[2]*e34*e45*e56;
352 
353  c54 = dense.cdsqte*50./g[4];
354  c64 = dense.cdsqte*415./g[5];
355 
356  c45 = c54*g[4]/g[3]*e45;
357  c46 = c64*g[5]/g[3]*e45*e56;
358 
359  c65 = dense.cdsqte*400./g[5];
360  c56 = c65*g[5]/g[4]*e56;
361 
364  /* this is check for whether matrix inversion likely to fail */
365  simple = (c16 + atoms.pmpo15)/(c61 + c62 + c64 + a65 + a64 + a62);
366  if( simple < 1e-19 )
367  {
368  atoms.popoi[0] = (realnum)abundoi;
369  for( i=1; i < 6; i++ )
370  {
371  atoms.popoi[i] = 0.;
372  }
373  *coloi = 0.;
374  return;
375  }
376 
377  /*--------------------------------------------------------- */
378 
379  for( i=0; i < 6; i++ )
380  {
381  zz[i][0] = 1.0;
382  zz[6][i] = 0.;
383  }
384 
385  /* first equation is sum of populations */
386  zz[6][0] = abundoi;
387 
388  /* level two, 3s 3So */
389  zz[0][1] = -c12;
390  zz[1][1] = c21 + c23 + c24 + c25 + c26 + a21;
391  zz[2][1] = -c32 - a32;
392  zz[3][1] = -c42;
393  zz[4][1] = -c52;
394  zz[5][1] = -c62 - a62;
395 
396  /* level three */
397  zz[0][2] = -c13;
398  zz[1][2] = -c23;
399  zz[2][2] = c31 + c32 + c34 + c35 + c36 + a32;
400  zz[3][2] = -c43 - a43;
401  zz[4][2] = -c53 - a53;
402  zz[5][2] = -c63;
403 
404  /* level four */
405  zz[0][3] = -c14;
406  zz[1][3] = -c24;
407  zz[2][3] = -c34;
408  zz[3][3] = c41 + c42 + c43 + c45 + c46 + a41 + a43;
409  zz[4][3] = -c54;
410  zz[5][3] = -c64 - a64;
411 
412  /* level five */
413  zz[0][4] = -c15 - atoms.pmpo15;
414  zz[1][4] = -c25;
415  zz[2][4] = -c35;
416  zz[3][4] = -c45;
417  zz[4][4] = c51 + c52 + c53 + c54 + c56 + a51 + a53;
418  zz[5][4] = -c65 - a65;
419 
420  /* level six */
421  zz[0][5] = -c16;
422  zz[1][5] = -c26;
423  zz[2][5] = -c36;
424  zz[3][5] = -c46;
425  zz[4][5] = -c56;
426  zz[5][5] = c61 + c62 + c63 + c64 + c65 + a65 + a64 + a62;
427 
428  /* this one may be more robust */
429  for( j=0; j < 6; j++ )
430  {
431  for( i=0; i < 6; i++ )
432  {
433  amat[i][j] = zz[i][j];
434  }
435  bvec[j] = zz[6][j];
436  }
437 
438  ner = 0;
439 
440  getrf_wrapper(6, 6, (double*)amat, 6, ipiv, &ner);
441  getrs_wrapper('N', 6, 1, (double*)amat, 6, ipiv, bvec, 6, &ner);
442 
443  /*DGETRF(6,6,(double*)amat,6,ipiv,&ner);*/
444  /*DGETRS('N',6,1,(double*)amat,6,ipiv,bvec,6,&ner);*/
445  if( ner != 0 )
446  {
447  fprintf( ioQQQ, " oi_level_pops: dgetrs finds singular or ill-conditioned matrix\n" );
448  cdEXIT(EXIT_FAILURE);
449  }
450 
451  /* now put results back into z so rest of code treates only
452  * one case - as if matin1 had been used */
453  for( i=0; i < 6; i++ )
454  {
455  zz[6][i] = bvec[i];
456  }
457 
458  lgNegPop = false;
459  for( i=0; i < 6; i++ )
460  {
461  atoms.popoi[i] = (realnum)zz[6][i];
462  if( atoms.popoi[i] < 0. )
463  lgNegPop = true;
464  }
465 
466  /* following used to confirm that all dep coef are unity at
467  * density of 1e17, t=10,000, and all A's set to zero */
468  if( trace.lgTrace && trace.lgTr8446 )
469  {
470  deptoi[0] = 1.;
471  deptoi[1] = atoms.popoi[1]/atoms.popoi[0]/(g[1]/g[0]*
472  e12);
473  deptoi[2] = atoms.popoi[2]/atoms.popoi[0]/(g[2]/g[0]*
474  e12*e23);
475  deptoi[3] = atoms.popoi[3]/atoms.popoi[0]/(g[3]/g[0]*
476  e12*e23*e34);
477  deptoi[4] = atoms.popoi[4]/atoms.popoi[0]/(g[4]/g[0]*
478  e12*e23*e34*e45);
479  deptoi[5] = atoms.popoi[5]/atoms.popoi[0]/(g[5]/g[0]*
480  e12*e23*e34*e45*e56);
481 
482  fprintf( ioQQQ, " oilevl finds levl pop" );
483  for(i=0; i < 6; i++)
484  fprintf( ioQQQ, "%11.3e", atoms.popoi[i] );
485  fprintf( ioQQQ, "\n" );
486 
487  fprintf( ioQQQ, " oilevl finds dep coef" );
488  for(i=0; i < 6; i++)
489  fprintf( ioQQQ, "%11.3e", deptoi[i] );
490  fprintf( ioQQQ, "\n" );
491  }
492 
493  /* this happens due to numerical instability in matrix inversion routine */
494  if( lgNegPop )
495  {
496  atoms.nNegOI += 1;
497 
498  fprintf( ioQQQ, " OILEVL finds negative population" );
499  for(i=0;i < 6; i++)
500  fprintf( ioQQQ, "%10.2e", atoms.popoi[i] );
501  fprintf( ioQQQ, "\n" );
502 
503  fprintf( ioQQQ, " simple 5 =%10.2e\n", simple );
504 
505  atoms.popoi[5] = 0.;
506  atoms.popoi[4] = (realnum)(abundoi*(c15 + atoms.pmpo15)/(a51 + a53 +
507  c51 + c53));
508  atoms.popoi[3] = 0.;
509  atoms.popoi[2] = (realnum)(atoms.popoi[4]*(a53 + c53)/(a32 + c32));
510  atoms.popoi[1] = (realnum)((atoms.popoi[2]*(a32 + c32) + abundoi*
511  c12)/(a21 + c21));
512 
513  atoms.popoi[0] = (realnum)abundoi;
514  /* write(QQ,'('' OILEVL resets this to simple pop'',1P,6E10.2)')
515  * 1 popoi */
516  }
517 
518  /* this is total cooling due to model atom, can be neg (heating) */
519  *coloi =
520  (atoms.popoi[0]*c12 - atoms.popoi[1]*c21)*1.53e-11 +
521  (atoms.popoi[0]*c14 - atoms.popoi[3]*c41)*1.92e-11 +
522  (atoms.popoi[0]*c15 - atoms.popoi[4]*c51)*1.94e-11 +
523  (atoms.popoi[1]*c23 - atoms.popoi[2]*c32)*2.36e-12 +
524  (atoms.popoi[1]*c26 - atoms.popoi[5]*c62)*4.55e-12 +
525  (atoms.popoi[2]*c35 - atoms.popoi[4]*c53)*1.76e-12 +
526  (atoms.popoi[2]*c34 - atoms.popoi[3]*c43)*1.52e-12 +
527  (atoms.popoi[3]*c46 - atoms.popoi[5]*c64)*6.86e-13 +
528  (atoms.popoi[4]*c56 - atoms.popoi[5]*c65)*4.32e-13;
529 
530  return;
531 }

Generated for cloudy by doxygen 1.8.3.1