cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cool_oxyg.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 /*CoolOxyg evaluate total cooling due to oxygen */
4 #include "cddefines.h"
5 #include "coolheavy.h"
6 #include "dense.h"
7 #include "taulines.h"
8 #include "h2.h"
9 #include "phycon.h"
10 #include "embesq.h"
11 #include "hmi.h"
12 #include "oxy.h"
13 #include "colden.h"
14 #include "mole.h"
15 #include "ligbar.h"
16 #include "thermal.h"
17 #include "lines_service.h"
18 #include "atoms.h"
19 #include "cooling.h"
20 /* oi3pcs make collision strengths with electrons for 3P O I ground term.
21  *>>refer o1 cs Berrington, K.A. 1988, J.Phys.B, 21, 1083 for Te > 3000K
22  *>>refer o1 cs Bell, Berrington & Thomas 1998, MNRAS 293, L83 for 50K <= Te <= 3000K
23  * numbers in cs variables refer to j value: j=0,1,2 (n=3,2,1 in 5 level atom)
24  * written by Kirk Korista, 29 may 96
25  * adapted by Peter van Hoof, 30 march 99 (to include Bell et al. data)
26  * all data above 3000K have been adjusted down by a constant factor to make
27  * them line up with Bell et al. data. */
28 STATIC void oi3Pcs(double *cs21,
29  double *cs20,
30  double *cs10);
31 
32 #if defined (__ICC) && defined(__ia64) && __INTEL_COMPILER < 910
33 #pragma optimization_level 1
34 #endif
35 void CoolOxyg(void)
36 {
37  double a21,
38  a31,
39  a32,
40  a41,
41  a42,
42  a43,
43  a51,
44  a52,
45  a53,
46  a54,
47  a6300,
48  a6363,
49  aeff,
50  cs,
51  cs2s2p,
52  cs2s3p ,
53  cs01,
54  cs02,
55  cs12,
56  cs13,
57  cs21,
58  cs23,
59  cs31,
60  cs32,
61  cs41,
62  cs42,
63  cs43,
64  cs51,
65  cs52,
66  cs53,
67  cs54,
68  Te_bounded,
69  Te_bounded_log ,
70  o3cs23,
71  p[5],
72  r12 ,
73  r13,
74  pump_rate;
75  /*double cs_ab ,cs_ac;(Commented out- Humeshkar Nemala)*/
76  /* these were added by Peter van Hoof to update the collision
77  * data within the OI ground term */
78  double cse01=-1.,cse12=-1.,cse02 =-1.,
79  csh01=-1.,cshe01=-1.,csh201=-1.,csh12=-1.,cshe12=-1.,csh212=-1.,csh02=-1.,cshe02=-1.,csh202 =-1.,
80  csh2o01=-1.,csh2o02=-1.,csh2o12=-1.,csh2p01=-1.,csh2p02=-1.,csh2p12=-1.,csp01=-1.,csp02=-1.,
81  csp12=-1.;
82  static double go2[5]={4.,6.,4.,4.,2.};
83  static double exo2[4]={26808.4,21.0,13637.5,1.5};
84  /* these will be used to update change in cs wrt temperature,
85  * and its affect on cooling derivative */
86  static double oi_cs_tsave=-1. , oi_te_tsave=-1. , oi_dcdt_tsave=-1.;
87  long int i;
88  static bool lgFirst=true;
89  static long int *ipO4Pump=NULL,
90  nO4Pump=0;
91  double rate_OH_dissoc;
92 
93  DEBUG_ENTRY( "CoolOxyg()" );
94 
95  /* following does the OI Bowen Ly-bet pumping problem */
97  CoolAdd("o 1",8446,CoolHeavy.coolOi);
98 
99  /* [OI] 6300, 6363, 5575, etc
100  * >>chng 06 oct 02, Humeshkar Nemala incorporate Barklem cs data for OI
101  * largest difference is 3P - 1D (6300+6363) which is now roughly 3x smaller */
102  /* This is the transition from 3P(J=2,1,0;the levels are reversed) to 1D2
103  * The rate coefficients were converted to CS
104  *>>refer oi cs Barklem,P.S.,2006,A&A (astroph 0609684)
105  * Data pts are avilable over 1000,3000,5000,8000,12000,20000 and 50000
106  * Fits between 1000K and 20000K are reliable;this is not the case
107  * between 20000K & 50000K*/
108  /* this was old calculation of cs12 that was replaced by Humeshkar */
109 #if 0
110  cs12 = 2.68e-5*phycon.te*(1. + 1.67e-6*phycon.te - 2.95e-10*phycon.te*
111  phycon.te);
112  cs12 = MAX2(0.1,cs12);
113 #endif
114 
115  /* [OI] electron 6300, 6363 collision strength
116  * >>chng 06 oct 02, Humeshkar Nemala code, based on Barklem paper */
117  if(phycon.te < 8E3)
118  {
119  cs12 = (4E-08)*(phycon.te*phycon.te70*phycon.te05);
120  }
121  else if(phycon.te < 2E4)
122  {
123  cs12 = (4.630155E-05)*((phycon.te/phycon.te04)*phycon.te005*phycon.te0001);
124  }
125  else
126  {
127  /* this branch T > 20,000K */
128  cs12 = (1.5394E-03)*(phycon.sqrte*phycon.te10*phycon.te01*phycon.te001*phycon.te0003);
129  }
130 
131  /* this block adds on collisional excitation by H0 */
132  {
133  /* >>chng 06 aug 18, add atomic hydrogen collisional processes using rates from
134  * >>refer O1 coll Krems, R.V., Jameson, M.J. & Dalgarno, A. 2006, ApJ, 647, 1531
135  * their equation 10 - the deecxiation rate coefficient, cm3 s-1
136  * >>referold oi cs Federman, S.R., & Shipsey, E.J. 1983, ApJ, 269, 791 */
137  double te_scale = phycon.te / 6000.;
138  double rate_H0 = (1.74*te_scale + 0.6)*1e-12*sexp(0.47*te_scale) / sqrt(te_scale ) *
140  /*fprintf(ioQQQ,"DEBUG OI H0 rate %.2e " , cs12 );*/
141  cs12 += ConvRate2CS( 5.f , (realnum)rate_H0 );
142  /*fprintf(ioQQQ,"%.2e\n" , cs12 );*/
143  }
144  /* following was old fit to 2-3 5755 transition */
145  /*cs23 = 1.05e-3*phycon.sqrte;*/
146  /*>>chng 06 oct 02, Humeshkar Nemala's new fit to Barklem paper */
147  if(phycon.te < 5E3)
148  {
150  }
151  else if(phycon.te<2E4)
152  {
153  cs23 = (1.98479e-04)*((phycon.te70/phycon.te03)*phycon.te003*phycon.te0007);
154  }
155  else
156  {
158  }
159 
160  /* following was old fit to 1,2,3 -> 5 transition */
161  /*cs13 = 3.24e-6*phycon.te;*/
162  /* following is Humeshkar's new fit to Barklem paper */
163  if(phycon.te < 2E4)
164  {
166  }
167  else
168  {
169  cs13 = (1.054E-03)*(phycon.sqrte/phycon.te04)*phycon.te003*phycon.te0005;
170  }
171 
172  /* O I 6300, 6363, A from
173  * >>refer all all Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103,
174  * >>refercon ed by D.R. Flower, (D. Reidel: Holland), 143 */
177  TauLines[ipT6300].Lo->Pop = (dense.xIonDense[ipOXYGEN][0]*5./5.);
178  TauLines[ipT6300].Hi->Pop = 0.;
179  TauLines[ipT6300].Coll.cs = (realnum)(cs12*5./9.);
181  TauLines[ipT6363].Lo->Pop = (dense.xIonDense[ipOXYGEN][0]*5./3.);
182  TauLines[ipT6363].Hi->Pop = 0.;
183  TauLines[ipT6363].Coll.cs = (realnum)(cs12*3./9.);
185 
186  a21 = a6300 + a6363 + oxy.d6300;
188 
189  PutCS(cs23,&TauLines[ipT5577]);
190 
191  /* rate of new populations of O^0 due to dissociation of OH,
192  * co.rate_OH_dissoc is rate OH -> O + H [cm-3 s-1],
193  * must make it per unit O atom, so this rate is s-1 excitations per O atom */
194  rate_OH_dissoc = CO_findrate("PHOTON,OH=>O,H");
195  r12 = rate_OH_dissoc*0.55/SDIV( dense.xIonDense[ipOXYGEN][0] );
196  r13 = rate_OH_dissoc*0.05/SDIV( dense.xIonDense[ipOXYGEN][0] );
197 
198  /* below is correction for fraction of excitations the produce emission */
199  CoolHeavy.c6300_frac_emit = (a6300+a6363)/(a6300+a6363+cs12*dense.cdsqte/5.);
200  CoolHeavy.c5577_frac_emit = (a32)/(a32+cs23*dense.cdsqte/3.);
201 
202  /* d6300 is the photoionization rate from the excited level
203  * was computed when ionization balance done */
204  CoolHeavy.c5577 = atom_pop3(9.,5.,1.,cs12,cs13,cs23,a21,7.56e-2,a32,
205  22590.,25775.7,&oxy.poiexc,dense.xIonDense[ipOXYGEN][0],0.,r12,r13)*a32*
206  3.57e-12;
207 
210  TauLines[ipT5577].Hi->Pop = 0.;
211 
212  /* convert poiexc to relative abundances */
213  /* >>chng 04 apr 20, include correction for fraction of 6300 due to OH pump */
214  CoolHeavy.c5577 *= (1.-r13/(SDIV(atoms.c13)));
215 
217  (1.-r12/(SDIV(atoms.c12)));
219  (1.-r12/(SDIV(atoms.c12)));
220  /* must introduce correction for fraction of 6300 that is photo produced */
222  /* note that atoms.c12 has all ra tes 1->2 */
223  (1.-r12/(SDIV(atoms.c12)));
224 
225  oxy.poiexc /= (realnum)MAX2(1e-20,dense.xIonDense[ipOXYGEN][0]);
226  CoolAdd("o 1",5577,CoolHeavy.c5577);
227  CoolAdd("o 1",6300,CoolHeavy.c6300);
228  CoolAdd("o 1",6363,CoolHeavy.c6363);
229 
230  /* O I fine structure lines rad data from
231  * >>refer all cs Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103,
232  * >>refercon ed by D.R. Flower, (D. Reidel: Holland), 143
233  * >>refer o1 as Berrington, K.A. 1988, J.Phys. B, 21, 1083
234  * hydrogen collisions from
235  * >>refer oi cs Tielens, A.G.G., & Hollenbach, D. 1985, ApJ, 291, 722
236  * factor in () is their rate coef
237  * assume H2 and H are same
238  * CDSQTE = 8.629E-6*EDEN/SQRTE
239  * cs01 = 9.8e-6*te + (4.2e-12*te70/te03) / cdsqte * 3. * hdcor
240  * cs12 = 2.6e-6*te + (1.5e-10*sqrte/te03/te03)/cdsqte*hdcor
241  * cs02 = 2.9e-6*te + (1.1e-12*te70*te10)/cdsqte*hdcor
242  * evaluate fits to OI electron rates, indices on var do not agree
243  * with Kirk's in sub, but are OK */
244  /*==============================================================*/
245  /* >>>chng 99 jun 01,
246  * following changed to be parallel to Peter van Hoof's changes
247  * in the Fortran C90.05 version
248  * following is collisions with electrons */
249  oi3Pcs(&cse01,&cse02,&cse12);
250  /* remember the electron part of the cs */
251  cs = cse01;
252 
253  /* rate coefficients for collisional de-excitation of O^o(3P) with H^o(2S1/2)
254  * NOTE: due to lack of data these relations are extrapolated to higher Te !
255  * >>refer o1 cs Launay & Roueff 1977, AA 56, 289
256  * the first fit is for Te <= 300K, the second for higher temps
257  * these data are valid for 50K <= Te <= 1000K*/
258  csh12 = MAX2(2.00e-11*phycon.te30*phycon.te05*phycon.te02,
259  7.67e-12*phycon.sqrte*phycon.te03);
260 
261  /* these data are valid for 100K <= Te <= 1000K */
262  csh01 = MIN2(3.12e-12*phycon.te70*phycon.te02*phycon.te02,
263  7.51e-12*phycon.sqrte*phycon.te05*phycon.te03);
264 
265  /* these data are valid for 150K <= Te <= 1000K*/
266  csh02 = MIN2(6.96e-13*phycon.te/phycon.te10/phycon.te02,
267  1.49e-12*phycon.te70*phycon.te05);
268 
269  /* rate coefficients for collisional de-excitation of O^o(3P) with He^o(1S)
270  * NOTE: due to lack of data these relations are extrapolated to higher Te !
271  * >>refer oi cs Monteiro & Flower 1987, MNRAS 228, 101
272  * the first fit is for Te <= 300K, the second for higher temps
273  * these data are valid for 100K <= Te <= 1000K */
274  cshe12 = MIN2(8.09e-16*phycon.te32*phycon.te10*phycon.te03,
275  9.72e-15*phycon.te*phycon.te20);
276 
277  cshe01 = 1.57e-12*phycon.te70/phycon.te01;
278 
279  cshe02 = MIN2(1.80e-12*phycon.te70*phycon.te05,
280  4.45e-12*phycon.te70/phycon.te10);
281 
282  if( phycon.te<=1.5e3 )
283  {
284  /* >>chng 04 mar 15, use explicit ortho-para densities */
285  double ortho_frac = h2.ortho_density/SDIV(hmi.H2_total);
286  /* rate coefficients for collisional de-excitation of O^o(3P) with H2(J=1,0)
287  * >>refer oi cs Jaquet et al. 1992, J.Phys.B 25, 285
288  * these data are valid for 40K <= Te <= 1500K
289  * the first entry is contribution from ortho H2, the second para H2.*/
290  csh2o12 = 2.21e-14*phycon.te*phycon.te10/phycon.te01;
291  csh2p12 = 9.45e-15*phycon.te*phycon.te20/phycon.te01;
292  csh212 = ortho_frac*csh2o12 + (1.-ortho_frac)*csh2p12;
293 
294  csh2o01 = 2.37e-11*phycon.te30*phycon.te10/phycon.te02;
295  csh2p01 = 3.40e-11*phycon.te30*phycon.te02;
296  csh201 = ortho_frac*csh2o01 + (1.-ortho_frac)*csh2p01;
297 
298  csh2o02 = 4.86e-11*phycon.te30*phycon.te02*phycon.te02;
299  csh2p02 = 6.82e-11*phycon.te30/phycon.te02;
300  csh202 = ortho_frac*csh2o02 + (1.-ortho_frac)*csh2p02;
301  }
302  else
303  {
304  csh212 = 0.;
305  csh201 = 0.;
306  csh202 = 0.;
307  }
308 
309  /* effective collision strength of O^o(3P) with p
310  * >>refer oi cs Pequignot, D. 1990, A&A 231, 499
311  * original data:
312  * >>refer oi cs Chambaud et al., 1980, J.Phys.B, 13, 4205 (upto 5000K)
313  * >>refer oi cs Roueff, private communication (10,000K and 20,000K)*/
314  if( phycon.te<=1000. )
315  {
316  csp01 = MAX2(2.22e-5*phycon.te/phycon.te10,
317  /* >>chng 05 jul 05, eden to dense.EdenHCorr */
318  /*2.69e-6*phycon.te*phycon.te30)*dense.xIonDense[ipHYDROGEN][1]/dense.eden;*/
320 
321  csp02 = MIN2(7.07e-8*phycon.te32*phycon.te10,2.46e-7*
322  /* >>chng 05 jul 05, eden to dense.EdenHCorr */
323  /*phycon.te32/phycon.te10)*dense.xIonDense[ipHYDROGEN][1]/dense.eden;*/
325  }
326  else
327  {
328  csp01 = MIN2(2.69e-6*phycon.te*phycon.te30,9.21e-5*phycon.te/phycon.te10/
329  /* >>chng 05 jul 05, eden to dense.EdenHCorr */
330  /*phycon.te03)*dense.xIonDense[ipHYDROGEN][1]/dense.eden;*/
332 
333  csp02 = MIN2(2.46e-7*phycon.te32/phycon.te10,5.21e-5*phycon.te/
334  /* >>chng 05 jul 05, eden to dense.EdenHCorr */
335  /*phycon.te20)*dense.xIonDense[ipHYDROGEN][1]/dense.eden;*/
337  }
338 
339  csp12 = MIN2(2.35e-6*phycon.te*phycon.te05*phycon.te01,3.98e-5*
340  /* >>chng 05 jul 05, eden to dense.EdenHCorr */
341  /*phycon.te20)*dense.xIonDense[ipHYDROGEN][1]/dense.eden;*/
342  /*phycon.te70/phycon.te01)*dense.xIonDense[ipHYDROGEN][1]/dense.eden;*/
344 
345  /*cs01 = cse01+csp01+3.*(csh01*dense.xIonDense[ipHYDROGEN][0] + cshe01*dense.xIonDense[ipHELIUM][0] + csh201*hmi.Hmolec[ipMH2g])/
346  dense.cdsqte;
347  cs12 = cse12+csp12+(csh12*dense.xIonDense[ipHYDROGEN][0] + cshe12*dense.xIonDense[ipHELIUM][0] + csh212*hmi.Hmolec[ipMH2g])/
348  dense.cdsqte;
349  cs02 = cse02+csp02+(csh02*dense.xIonDense[ipHYDROGEN][0] + cshe02*dense.xIonDense[ipHELIUM][0] + csh202*hmi.Hmolec[ipMH2g])/
350  dense.cdsqte;*/
351 
352  cs01 = cse01+csp01+3.*(csh01*dense.xIonDense[ipHYDROGEN][0] + cshe01*dense.xIonDense[ipHELIUM][0] + csh201*hmi.H2_total)/
353  dense.cdsqte;
354  cs12 = cse12+csp12+(csh12*dense.xIonDense[ipHYDROGEN][0] + cshe12*dense.xIonDense[ipHELIUM][0] + csh212*hmi.H2_total)/
355  dense.cdsqte;
356  cs02 = cse02+csp02+(csh02*dense.xIonDense[ipHYDROGEN][0] + cshe02*dense.xIonDense[ipHELIUM][0] + csh202*hmi.H2_total)/
357  dense.cdsqte;
358 
359  /* end changes */
360  /*==============================================================*/
361  PutCS(cs01,&TauLines[ipT63]);
362  PutCS(cs12,&TauLines[ipT146]);
363  PutCS(cs02,&TauDummy);
364  atom_level3(&TauLines[ipT63],&TauLines[ipT146],&TauDummy);
365 
366  /* now save pops to add col den in radinc */
367  for( i=0; i<3; ++i)
368  {
369  /* >>chng 02 oct 23, bug - had been O1Colden rather than O1Pops */
371  }
372  /*fprintf(ioQQQ,"colllll\t%li\t%.3e\t%.3e\n",
373  nzone,
374  colden.O1Pops[0],
375  colden.O1Pops[1] );*/
376 
377  /* >>>chng 99 apr 29, for orion pdr.in calc, temp oscillated due to
378  * bad dC/dT estimate, due to rapid change in cs with T. added following*/
379  /* when neutral collisions onto OI dominate the cooling the derivative
380  * produced by atom_level3 is vastly too small because of the strong temperature
381  * dependence of the neutral collision strength. Increase cooling derivative by
382  * ratio of cs due to electron (nearly constant) and total */
383  /* following was correction for this between 00 apr 29 and 02 jul 25 */
384 # if 0
385  if( cs01 > SMALLFLOAT )
386  {
387  thermal.dCooldT +=
388  (cs01-cs)/cs01*TauLines[ipT63].cool*
389  TauLines[ipT63].EnergyK*thermal.tsq1;
390  }
391 # endif
392 
393  /* >>chng 02 jul 25, keep track of change in cs for 63 micron line */
394  if( fabs(phycon.te - oi_te_tsave)/phycon.te > 1e-4 )
395  {
396  /* very first time we come through, previous values are -1 */
397  if(oi_te_tsave>0. )
398  {
399  oi_dcdt_tsave = (cs01-oi_cs_tsave) / (phycon.te-oi_te_tsave);
400  }
401  else
402  {
403  oi_dcdt_tsave = 0.;
404  }
405  oi_cs_tsave = cs01;
406  oi_te_tsave = phycon.te;
407 
408  /* >>chng 03 jan 21, this factor could become very large - it should
409  * always be positive since neutral cs's are, and not much bigger than
410  * the usual derivative */
411  /* can't be negative */
412  oi_dcdt_tsave = MAX2( 0. , oi_dcdt_tsave);
413  /* can't be bigger than several times normal dC/dT */
414  oi_dcdt_tsave = MIN2( TauLines[ipT63].EnergyK*thermal.tsq1*4. , oi_dcdt_tsave);
415  }
416  /* this is only derivative due to change in collision strength, which is capped to be
417  * less than 4x the thermal derivative just above */
418  thermal.dCooldT += TauLines[ipT63].Coll.cool*oi_dcdt_tsave;
419  /* this is a bebug print statement for the numerical cs derivative */
420  {
421  /*@-redef@*/
422  enum{DEBUG_LOC=false};
423  /*@+redef@*/
424  if( DEBUG_LOC )
425  {
426  fprintf(ioQQQ,"DEBUG OI\t%.2f\tte\t%.5e\tcool\t%.5e\tcs\t%.4e\told\t%.4e\tnew\t%.4e\n",
427  fnzone,
428  phycon.te,
429  TauLines[ipT63].Coll.cool ,
430  TauLines[ipT63].Coll.cs ,
431  TauLines[ipT63].Coll.cool*TauLines[ipT63].EnergyK*thermal.tsq1,
432  TauLines[ipT63].Coll.cool*oi_dcdt_tsave );
433  }
434  }
435  /* [O II]
436  * two sets of A's, Zeippen 1982 being the default
437  * >>chng 97 feb 19, change in last sig fig as per Bob Rubin e-mail feb 11 97
438  * A from NIST compilation
439  * >>refer o2 as Wiese, W.L., Fuhr, J.R., Deters, T.M. 1996, J Phys Chem Ref Data,
440  * >>refercon Monograph 7
441  * these are selected with the set atomic data ion oxygen 2 command */
442  if( dense.lgAsChoose[ipOXYGEN][1] )
443  {
444  /* not used unless requested with set atomic data command */
445  a21 = 3.058e-5;
446  a31 = 1.776e-4;
447  a41 = 5.22e-2;
448  a51 = 2.12e-2;
449  a32 = 1.30e-7;
450  a42 = 0.09907;
451  a52 = 0.0519;
452  a43 = 0.0534;
453  a53 = 0.08672;
454  a54 = 1.41e-10;
455  }
456  else
457  {
458  /* default - the Zeippen (1982) values
459  *>>refer O2 As Zeippen, C. J. 1982, MNRAS, 198, 111
460  * these are recommended by
461  *>>refer O2 As Wang et al. 2004 A&A, 427, 873
462  >>refer O2 As Chen, Qing & Li Phys Rev A, 2007, 76, 042507
463  */
464  a21 = 3.82e-5;
465  a31 = 1.65e-4;
466  a41 = 5.64e-2;
467  a51 = 2.3202e-2;
468  a32 = 1.20e-7;
469  a42 = 0.11678;
470  a52 = 0.0615;
471  a43 = 0.0614;
472  a53 = 0.10175;
473  a54 = 2.08e-11;
474  }
475  /* >>chng 01 may 18, update cs to
476  * >>referold o2 cs McLaughlin, B.M., & Bell, K.L. 1998, J Phys B 31, 4317
477  * >>chng 02 mar 13, go back to older values as per Seaton/Osterbrock correspondence
478  * this is a mix of the McLauthlin Bell 93 and 98 results, with the expected
479  * branching to levels within the term
480  * >>chng 04 nov 01, statistical weights were reversed, caught by Kevin Blagrave
481  * >>refer o2 cs A.K.Pradhan, M.Montenegro,S.N.Nahar & W.Eissner,2006,MNRAS,366,L6-L9
482  */
483  cs21 = (realnum)(0.8229*(phycon.te007*phycon.te0005*phycon.te0001));
484  cs31 = (realnum)(0.5976/phycon.te002);
485  cs32 = (realnum)(2.6084/(phycon.te07/(phycon.te003*phycon.te0001)));
486  cs41 = (realnum)(0.2471*phycon.te02*(phycon.te007/phycon.te0004));
487  cs42 = (realnum)(0.7075*phycon.te03*phycon.te004*phycon.te0002);
488  cs43 = (realnum)(0.415*phycon.te04*(phycon.te004/phycon.te0004));
489  cs51 = (realnum)(0.1294*(phycon.te02/(phycon.te001*phycon.te0004)));
490  cs52 = (realnum)(0.2841*phycon.te04*phycon.te0004);
491  cs53 = (realnum)(0.2737*phycon.te04*phycon.te003*phycon.te0001);
492  cs54 = (realnum)(0.2037*phycon.te04*(phycon.te002/phycon.te0004));
493  atom_pop5(go2,exo2,cs21,cs31,cs41,cs51,cs32,cs42,cs52,cs43,cs53,cs54,
494  a21,a31,a41,a51,a32,a42,a52,a43,a53,a54,p,dense.xIonDense[ipOXYGEN][1]);
495 
496  CoolHeavy.O3730 = (realnum)(p[1]*a21*5.34e-12);
497  CoolHeavy.O3726 = (realnum)(p[2]*a31*5.34e-12);
498  CoolHeavy.O2471 = (realnum)((p[3]*a41 + p[4]*a51)*8.05e-12);
499  CoolHeavy.O7323 = (realnum)((p[3]*a42 + p[4]*a52)*2.72e-12);
500  CoolHeavy.O7332 = (realnum)((p[3]*a43 + p[4]*a53)*2.71e-12);
504  CoolAdd("O 2",3727,CoolHeavy.c3727);
505  CoolAdd("O 2",7325,CoolHeavy.c7325);
506  CoolAdd("O 2",2470,CoolHeavy.c7325*0.78);
507 
508  /* remember ratio of radiative to total decays, to use for estimating
509  * recombination contribution in lines_lv1_li_ne */
510  if( (p[3] + p[4]) > SMALLFLOAT )
511  {
512  CoolHeavy.O2_A3_tot = (p[3]*(a41+a42+a43) + p[4]*(a51+a52+a53) ) /
513  ( (p[3]*(a41+a42+a43) + p[4]*(a51+a52+a53) ) +
514  ( p[3]*(cs41+cs42+cs43)/go2[3] + p[4]*(cs51+cs52+cs53)/go2[4]) *
515  dense.cdsqte );
516  }
517  else
518  {
519  CoolHeavy.O2_A3_tot = 0.;
520  }
521 
522  if( (p[1] + p[2]) > SMALLFLOAT )
523  {
524  CoolHeavy.O2_A2_tot = (p[1]*a21 + p[2]*a31 ) /
525  ( (p[1]*a21 + p[2]*a31 ) +
526  ( p[1]*cs21/go2[1] + p[2]*cs31/go2[2]) *
527  dense.cdsqte );
528  }
529  else
530  {
531  CoolHeavy.O2_A2_tot = 0.;
532  }
533 
534  /* O II 833.9A, CS
535  * >>refer o2 cs McLaughlin, B.M., & Bell, K.L. 1993, ApJ, 408, 753 */
536  /* >>chng 01 aug 10, turn this line back on - no reason given for turning it off. */
538  phycon.te001;
539  PutCS(cs,&TauLines[ipT834]);
540  atom_level2(&TauLines[ipT834]);
541 
542  /* O III 1666, crit den=2.6+10, A from
543  * >>refer all cs Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103,
544  * >>refercon ed by D.R. Flower, (D. Reidel: Holland), 143
545  * c.s.
546  * >>referold o3 cs Lennon, D.J. Burke, V.M. 1994, A&AS, 103, 273 */
547  /* >>chng 06 jun 30- Humeshkar Nemala */
548  /*Refer o3 cs Aggarwal,K.M. & Keenan,F. P.1999,ApJS,123,311*/
549  /*Data available in the temperature range 2500 K to 2E6 K*/
550  /*fits at temperatures below 30000K and at temperatures above 30000K*/
551  if(phycon.te < 30000)
552  {
554  }
555  else
556  {
557  cs = (realnum)(6.166388*(1/(phycon.te20*phycon.te01*phycon.te002)));
558  }
559  /*PutCS(0.756,&TauLines[ipT1666]);*/
560  PutCS(cs ,&TauLines[ipT1666]);
561  /*>>chng 06 jun 30- Humeshkar Nemala*/
562  /*>>Refer o3 cs Aggarwal,K.M. & Keenan,F. P.1999,ApJS,123,311*/
563  /*Data available in the temperature range 2500 K to 2E6 K*/
564  /*fits at temperatures below 30000K and at temperatures above 30000K*/
565  if(phycon.te < 30000)
566  {
568  }
569  else
570  {
571  cs=(realnum)((3.668474)*(1/(phycon.te20*phycon.te01*phycon.te001*phycon.te0002)));
572  }
573  /*PutCS(0.454,&TauLines[ipT1661]);*/
574  PutCS(cs,&TauLines[ipT1661]);
575  /*this line correspond to the transition 3P0-5S^o2*/
576  /*>>chng 06 jun 30- Humeshkar Nemala*/
577  /*>>Refer o3 cs Aggarwal,K.M. & Keenan,F. P.1999,ApJS,123,311*/
578  /*Data available in the temperature range 2500 K to 2E6 K*/
579  /*fits at temperatures below 30000K and at temperatures above 30000K*/
580  if(phycon.te < 30000)
581  {
583  }
584  else
585  {
586  cs = (realnum)(1.223633/(phycon.te20*phycon.te01*phycon.te001*phycon.te0002));
587  }
588  /*PutCS(1.3,&TauDummy);*/
589  PutCS(cs,&TauDummy);
590  atom_level3(&TauDummy,&TauLines[ipT1666],&TauLines[ipT1661]);
591  /* if( nLev3Fail.gt.0 )write(qq,'('' nlev3fail='',i3)') nLev3Fail
592  *
593  * OIII 304 not in cooling function */
597  TauLines[ipT304].Hi->Pop = 0.;
598 
599  /* o iii 5007+4959, As 96 NIST */
600  /* >>chng 01 may 04, As updated to
601  * >>refer o3 as Storey, P.J., & Zeippen, C.J., 2000, MNRAS, 312, 813-816,
602  * changed by 10 percent! */
603  aeff = 0.027205 + oxy.d5007r;
604  /* following data used in routine that deduces OIII temp from spectrum
605  *
606  * >>refer o3 cs Lennon, D.J. Burke, V.M. 1994, A&AS, 103, 273
607  * IP paper Y(ki) differ significantly from those
608  * calculated by
609  * >>refer o3 cs Burke, V.M., Lennon, D.J., & Seaton, M.J. 1989, MNRAS, 236, 353
610  * especially for 1D2-1S0.
611  * BLS89 is adopted for 1D2-1S0 and LB94 for 3P2,1-1D2.
612  * NB!! these cs's must be kept parallel with those in Peimbert analysis */
613 
614  /* >>refer o3 cs Aggarwal,K.M. & Keenan,F. P.1999,ApJS,123,311*/
615  /* >>chng 06 jul 24- Humeshkar Nemala */
616  /*cs are available over two temperature ranges: below 30000K and above 30000K*/
617  /*The old values seemed to saturate at around 100,000K.The new values are around
618  10-15% different from the old values*/
619  /*The cs of the transitions 3P0,1,2 to 1D2 are added together to give oxy.o3cs12 */
620  /*the cs of the transition 1D2-1S0 is mentioned as o3cs23*/
621  /*The cs of the transitions 3P0,1,2 to 1S0 are added together to give oxy.o3cs13*/
622  if(phycon.te < 3E4)
623  {
624  oxy.o3cs12 = (realnum)(0.9062*(phycon.te10/phycon.te002));
625  o3cs23 = (realnum)(0.0995*phycon.te10*phycon.te07*(phycon.te007/phycon.te0002));
627 
628  }
629  else if(phycon.te < 6E4)
630  {
632  oxy.o3cs13 = (realnum)(0.1963109*phycon.te05*phycon.te0007);
633  o3cs23 = (realnum)(0.781266/(phycon.te02*phycon.te003*phycon.te0001));
634  }
635  else
636  {
639  o3cs23 = (realnum)(0.781266/(phycon.te02*phycon.te003*phycon.te0001));
640  }
641 # if 0
642  /*Old code commented -Humeshkar Nemala */
643  if( phycon.te > 3981 && phycon.te <= 1e5 )
644  {
645  oxy.o3cs12 = (realnum)(3.0211144 - 101.57536/phycon.sqrte + 355.06905*
646  /*really is base e log of temperature */
647  log(phycon.te)/phycon.te);
648  o3cs23 = 0.32412181 + 79.051672/phycon.sqrte - 4374.7816/
649  phycon.te;
650  }
651  else if( phycon.te < 3981. )
652  {
653  oxy.o3cs12 = 2.15f;
654  o3cs23 = 0.4781;
655  }
656  else
657  {
658  oxy.o3cs12 = 2.6594f;
659  o3cs23 = 0.5304;
660  }
661  oxy.o3cs13 =
662  (realnum)(MIN2(0.36,0.0784*phycon.te10*phycon.te03*phycon.te01*
663  phycon.te003));
664  /*End :Old Code Commented Out-Humeshkar Nemala*/
665  /* >>chng 01 may 04, derive a21 from aeff, had been 0.02431
666  a21 = 0.02431;*/
667 # endif
668 
669  a21 = aeff - oxy.d5007r;
670  a31 = .215634;
671  a32 = 1.71;
672  oxy.o3ex23 = 32947.;
673  oxy.o3br32 = (realnum)(a32/(a31 + a32));
674  oxy.o3enro = (realnum)(4.56e-12/3.98e-12);
675  /* POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2) */
676  oxy.poiii3 =
677  (realnum)(atom_pop3(9.,5.,1.,oxy.o3cs12,oxy.o3cs13,o3cs23,
678  a21,a31,a32,28990.,oxy.o3ex23,&oxy.poiii2,dense.xIonDense[ipOXYGEN][2],
679  oxy.d5007r,0.,0.));
680  CoolHeavy.c4363 = oxy.poiii3*a32*4.56e-12;
681  CoolHeavy.c5007 = oxy.poiii2*a21*3.98e-12;
685  CoolAdd("O 3",5007,CoolHeavy.c5007);
686  CoolAdd("O 3",4363,CoolHeavy.c4363);
687  CoolAdd("O 3",2331,CoolHeavy.c4363*0.236);
688  /* >>chng 02 jun 27, prevent crash on Sun with gcc 3.1, PvH */
689  /* if( xIonDense[ipOXYGEN][2] > 0. ) */
691  {
694  }
695 
696  /* OIII IR lines, col str from iron project,
697  * >>referold o3 cs Lennon, D.J. Burke, V.M. 1994, A&AS, 103, 273 */
698  /*>>refer o3 cs Aggarwal,K.M. & Keenan,F. P.1999,ApJS,123,311*/
699  /*88 microns refers to the 3P0-3P1 transition*/
700  if(phycon.te < 3E4)
701  {
702  cs = (realnum)(0.3993*(phycon.te03/phycon.te001)*phycon.te0001);
703  }
704  else if(phycon.te < 1E5)
705  {
707  }
708  else
709  {
710  cs = (realnum)(1.310467/((phycon.te07/phycon.te001)*phycon.te0002));
711  }
712  /*PutCS(0.5590,&TauLines[ipTO88]);*/
713  PutCS(cs,&TauLines[ipTO88]);
714  /*52 microns refers to the 3P1-3P2 transition*/
715  if(phycon.te < 3E4)
716  {
717  cs = (realnum)(0.7812*(phycon.te05/phycon.te0001));
718  }
719  else if(phycon.te < 1.2E5)
720  {
721  cs = (realnum)(0.516157*phycon.te07*phycon.te02*phycon.te0001);
722  }
723  else
724  {
726  }
727  /*PutCS(1.335,&TauLines[ipT52]);*/
728  PutCS(cs,&TauLines[ipT52]);
729  /*TauDummy refers to the 3P0-3P2 transition*/
730  if(phycon.te < 3E4)
731  {
732  cs = (realnum)(0.1337*phycon.te07*phycon.te002*phycon.te0002);
733  }
734  else if(phycon.te < 1.2E5)
735  {
737  }
738  else
739  {
741  }
742  /*PutCS(0.2832,&TauDummy);*/
743  PutCS(cs,&TauDummy);
744  atom_level3(&TauLines[ipTO88],&TauLines[ipT52],&TauDummy);
745 
746  /* O III 834A, CS */
747  /* >>referold o3 cs Aggarwal, K.M., 1985 A&A 146, 149. */
748  /* >>chng 06 jul 22- Humeshkar Nemala */
749  /*>>refer o3 cs Aggarwal,K.M. & Keenan,F. P.1999,ApJS,123,311*/
750  /*the cs of OIII 834A was obtained by summing the cs of the six transitions:
751  3P0-3D^o1,3P1-3D^o1,3P2-3D^o1,3P1-3D^o2,3P2-3D^o2,3P2-3D^o3*/
752  if(phycon.te < 3E4)
753  {
754  cs = (realnum)(0.9595*phycon.te04*phycon.te0002);
755  }
756  else
757  {
759  }
760  /*PutCS(5.,&TauLines[ipT835]);*/
761  PutCS(cs,&TauLines[ipT835]);
762  atom_level2(&TauLines[ipT835]);
763  /* call DumpLine( t835 ); */
764 
765  /* these cs data only extend over a modest Temp range */
766  Te_bounded = MAX2(phycon.te,4500.);
767  Te_bounded = MIN2(Te_bounded,450000.);
768  Te_bounded_log = log(Te_bounded);
769 
770  /* O IV 789A, CS from
771  * >>refer o4 cs Zhang, H.L., Graziani, M., Pradhan, A.K. 1994, A&A, 283, 319 */
772  cs = -3.0102462 + 109.22973/Te_bounded_log - 18666.357/Te_bounded;
773  PutCS(cs,&TauLines[ipT789]);
774  atom_level2(&TauLines[ipT789]);
775 
776  /* O IV 26 micron, CS
777  * >>referold o4 cs Blum, R.D., & Pradhan, A.K. 1992, ApJS 80, 425
778  * A=
779  * >>refer o4 as Brage, T., Judge, P.G., & Brekke, P. 1996, ApJ. 464, 1030 */
780  /* O IV 26 micron */
781  /*cs = 9.2141614 - 5.1727759e-3*sqrt(Te_bounded) - 58.116447/Te_bounded_log;*/
782 
783  /* >>chng 06 jun 30- Humeshkar Nemala */
784  /*>>refer o4 cs Zhang,H.L.,Graziani,M. & Pradhan,A.K. 1994,A&A,283,319*/
785  /*cs given over the temp range :900K to 450000K
786 
787  if(phycon.te < 1.8E4)
788  {
789  cs = (realnum)(0.5363*phycon.te10*phycon.te05*(phycon.te01/phycon.te0004));
790  }
791  else if(phycon.te < 5.4E4)
792  {
793  cs = (realnum)(1.554285*phycon.te05*phycon.te001);
794  }
795  else
796  {
797  cs = (realnum)(95.373669/(phycon.te30*phycon.te02*(phycon.te007/phycon.te0002)));
798  }
799 
800  * >>chng 06 nov 08 - NPA. Update collision strength to new data from:
801  * >>refer o4 cs Tayal, S. 2006, ApJS 166, 634
802  * Equation derived by using TableCurve, and goes to zero as
803  * T => 0 and T => infinity */
804 
805  cs = (realnum)(exp(3.265723 - 0.00014686984*phycon.alnte*phycon.sqrte
806  -22.035066/phycon.alnte));
807  /* cs goes to zero at very high T, which will cause an assert in the
808  * n-level solver - don't let this happen */
809  cs = MAX2( cs , 1e-4 );
810  PutCS(cs,&TauLines[ipT26]);
811 
812  /* one time initialization if first call, and level 2 lines are on */
813  if( lgFirst && nWindLine )
814  {
815  lgFirst = false;
816  nO4Pump = 0;
817  for( i=0; i<nWindLine; ++i )
818  {
819  /* don't test on nelem==ipIRON since lines on physics, not C, scale */
820  if( TauLine2[i].Hi->nelem ==8 && TauLine2[i].Hi->IonStg==4 )
821  {
822  ++nO4Pump;
823  }
824  }
825  if( nO4Pump<0 )
826  TotalInsanity();
827  else if( nO4Pump > 0 )
828  /* create the space - can't malloc 0 bytes */
829  ipO4Pump = (long *)MALLOC((unsigned)(nO4Pump)*sizeof(long) );
830  nO4Pump = 0;
831  for( i=0; i<nWindLine; ++i )
832  {
833  /* don't test on nelem==ipIRON since lines on physics, not C, scale */
834  if( TauLine2[i].Hi->nelem ==8 && TauLine2[i].Hi->IonStg==4 )
835  {
836 # if 0
837  DumpLine( &TauLine2[i] );
838 # endif
839  ipO4Pump[nO4Pump] = i;
840  ++nO4Pump;
841  }
842  }
843  }
844  else
845  /* level 2 lines are not enabled */
846  nO4Pump = 0;
847 
848  /* now sum pump rates */
849  pump_rate = 0.;
850  for( i=0; i<nO4Pump; ++i )
851  {
852  pump_rate += TauLine2[ipO4Pump[i]].Emis->pump;
853 # if 0
854  fprintf(ioQQQ,"DEBUG C %li %.3e %.3e\n",
855  i,
856  TauLine2[ipO4Pump[i]].WLAng , TauLine2[ipO4Pump[i]].pump );
857 # endif
858  }
859 
860  /*atom_level2(&TauLines[ipT26]);*/
861  /*AtomSeqBoron compute cooling from 5-level boron sequence model atom */
862  /* >>refer o4 cs Blum, R.D., & Pradhan, A.K., 1992, ApJS 80, 425
863  * >>refer o4 cs Zhang, H.L., Graziani, M., Pradhan, A.K. 1994, A&A, 283, 319 */
864  AtomSeqBoron(&TauLines[ipT26],
865  &TauLines[ipO4_1400],
866  &TauLines[ipO4_1397],
867  &TauLines[ipO4_1407],
868  &TauLines[ipO4_1405],
869  &TauLines[ipO4_1401],
870  0.1367 , 0.1560 , 1.1564 , 0.0457 , pump_rate,"O 4");
871 
872  /* O V 630, CS from
873  * >>refer o5 cs Berrington, K.A., Burke, P.G., Dufton, P.L., Kingston, A.E. 1985,
874  * >>refercon At. Data Nucl. Data Tables, 33, 195 */
875  cs = MIN2(4.0,1.317*phycon.te10/phycon.te03);
876  PutCS(cs,&TauLines[ipT630]);
877  atom_level2(&TauLines[ipT630]);
878 
879  /* O V 1218; coll data from
880  * >>refer o4 cs Berrington, K.A., Burke, P.G., Dufton, P.L., Kingston, A.E. 1985,
881  * >>refercon At. Data Nucl. Data Tables, 33, 195 */
882  if( phycon.te <= 3.16e4 )
883  {
884  cs = 3.224/(phycon.te10*phycon.te03*phycon.te03*phycon.te003);
885  }
886  else
887  {
888  cs = 10.549/(phycon.te10*phycon.te10*phycon.te10/phycon.te03);
889  }
890  /* >>chng 01 sep 09, AtomSeqBeryllium will reset this to 1/3 so critical density correct */
891  PutCS(cs,&TauLines[ipT1214]);
892  /* c1214 = AtomSeqBeryllium( .87,1.05,3.32, t1214, .0216) * 1.64E-11
893  * AtomSeqBeryllium(CS23,CS24,CS34,tarray,A41)
894  * A's
895  * >>refer o5 cs Fleming, J., Bell, K.L, Hibbert, A., Vaeck, N., Godefroid, M.R.
896  * >>refercon 1996, MNRAS, 279, 1289 */
897  AtomSeqBeryllium(.87,0.602,2.86,&TauLines[ipT1214],.02198);
898  embesq.em1218 = (realnum)(atoms.PopLevels[3]*0.0216*1.64e-11);
899 
900  /* O VI 1035 li seq
901  * generate collision strengths, then stuff them in
902  * >>refer o6 vs Cochrane, D.M., & McWhirter, R.W.P. 1983, PhyS, 28, 25 */
903  ligbar(8,&TauLines[ipT1032],&TauLines[ipT150],&cs2s2p,&cs2s3p);
904  PutCS(cs2s2p,&TauLines[ipT1032]);
905  PutCS(cs2s2p*0.5,&TauLines[ipT1037]);
906  /* no data for the 2-3 transition */
907  PutCS(1.0,&TauDummy);
908  /* solve the 3 level atom */
909  atom_level3(&TauLines[ipT1037],&TauDummy,&TauLines[ipT1032]);
910 
911  PutCS(cs2s3p,&TauLines[ipT150]);
912  atom_level2(&TauLines[ipT150]);
913  return;
914 }
915 
916 /*"oi3pcs.h" make collision strengths with electrons for 3P O I ground term.
917  *>>refer o1 cs Berrington, K.A. 1988, J.Phys.B, 21, 1083 for Te > 3000K
918  *>>refer o1 cs Bell, Berrington & Thomas 1998, MNRAS 293, L83 for 50K <= Te <= 3000K
919  * numbers in cs variables refer to j value: j=0,1,2 (n=3,2,1 in 5 level atom)
920  * written by Kirk Korista, 29 may 96
921  * adapted by Peter van Hoof, 30 march 99 (to include Bell et al. data)
922  * all data above 3000K have been adjusted down by a constant factor to make
923  * them line up with Bell et al. data. */
924 STATIC void oi3Pcs(double *cs21,
925  double *cs20,
926  double *cs10)
927 {
928  double a,
929  b,
930  c,
931  d;
932 
933  DEBUG_ENTRY( "oi3Pcs()" );
934  /* local variables */
935 
936  /* 3P2 - 3P1 */
937  if( phycon.te <= 3.0e3 )
938  {
939  *cs21 = 1.49e-4*phycon.sqrte/phycon.te02/phycon.te02;
940  }
941  else if( phycon.te <= 1.0e4 )
942  {
943  a = -5.5634127e-04;
944  b = 8.3458102e-08;
945  c = 2.3068232e-04;
946  *cs21 = 0.228f*(a + b*phycon.te32 + c*phycon.sqrte);
947  }
948  else
949  {
950  *cs21 = 0.228*MIN2(0.222,5.563e-06*phycon.te*phycon.te05*
951  phycon.te02);
952  }
953 
954  /* 3P2 - 3P0 */
955  if( phycon.te <= 3.0e3 )
956  {
957  *cs20 = 4.98e-5*phycon.sqrte;
958  }
959  else if( phycon.te <= 1.0e4 )
960  {
961  a = -3.7178028e-04;
962  b = 2.0569267e-08;
963  c = 1.1898539e-04;
964  *cs20 = 0.288*(a + b*phycon.te32 + c*phycon.sqrte);
965  }
966  else
967  {
968  *cs20 = 0.288*MIN2(0.0589,1.015e-05*phycon.te/phycon.te10/
970  }
971 
972  /* 3P1 - 3P0 */
973  if( phycon.te <= 3.0e3 )
974  {
975  *cs10 = 1.83e-9*phycon.te*phycon.te30*phycon.te05;
976  }
977  else if( phycon.te <= 1.0e4 )
978  {
979  a = -5.9364373e-04;
980  b = 0.02946867;
981  c = 10768.675;
982  d = 3871.9826;
983  *cs10 = 0.0269*(a + b*exp(-0.5*POW2((phycon.te-c)/d)));
984  }
985  else
986  {
987  *cs10 = 0.0269*MIN2(0.074,7.794e-08*phycon.te32/phycon.te10/
988  phycon.te01);
989  }
990 
991  return;
992 }

Generated for cloudy by doxygen 1.8.4