cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cool_calc.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 /*CoolCalc compute calcium cooling */
4 #include "cddefines.h"
5 #include "taulines.h"
6 #include "doppvel.h"
7 #include "phycon.h"
8 #include "ca.h"
9 #include "dense.h"
10 #include "thermal.h"
11 #include "opacity.h"
12 #include "rfield.h"
13 #include "ligbar.h"
14 #include "lines_service.h"
15 #include "atoms.h"
16 #include "cooling.h"
17 
18 void CoolCalc(void)
19 {
20  realnum p2;
21  double a21,
22  a31,
23  a41,
24  a42,
25  a51,
26  a52,
27  a53,
28  c21,
29  Ca2pop[5] ,
30  cs,
31  cs2s2p,
32  cs2s3p ,
33  cs01,
34  cs02,
35  cs12,
36  cs14,
37  cs15,
38  d41,
39  d42,
40  d51,
41  d52,
42  d53,
43  hlgam,
44  op41,
45  op51,
46  opckh,
47  opcxyz,
48  PhotoRate2,
49  p3,
50  PhotoRate3,
51  PhotoRate4,
52  PhotoRate5,
53  r21,
54  r31,
55  r41,
56  r42,
57  r51,
58  r52,
59  r53;
60  static double gCa2[5]={2.,4.,6.,2.,4.};
61  static double exCa2[4]={13650.2,60.7,11480.6,222.9};
62  static realnum opCax = 0.f;
63  static realnum opCay = 0.f;
64  static realnum opCaz = 0.f;
65 
66  DEBUG_ENTRY( "CoolCalc()" );
67 
68  /* Ca I 4228 */
70  atom_level2(&TauLines[ipCaI4228]);
71 
72  /* photoionization of evcited levels by Ly-alpha */
73  hlgam = rfield.otslin[ Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont -1];
74  PhotoRate5 = 1.7e-18*hlgam;
75  PhotoRate4 = 8.4e-19*hlgam;
76  PhotoRate3 = 7.0e-18*hlgam;
77  PhotoRate2 = 4.8e-18*hlgam;
78 
79  /* spontaneous decays
80  * frist two trans prob from
81  * >>refer ca2 as Zeippen, C.J. 1990, A&A, 229, 248 */
82  a21 = 1.02*TauLines[ipT7324].Emis->Pesc;
83  a31 = 1.05*TauLines[ipT7291].Emis->Pesc;
84  a41 = 1.4e8*TauLines[ipT3969].Emis->Pesc;
85  a51 = 1.4e8*TauLines[ipT3934].Emis->Pesc;
86  a42 = 7.9e6*TauLines[ipT8662].Emis->Pesc;
87  a52 = 8.2e5*TauLines[ipT8498].Emis->Pesc;
88  a53 = 7.48e6*TauLines[ipT8542].Emis->Pesc;
89 
90  /* destruction of IR triplet by continuous opacities */
91  opcxyz = opac.opacity_abs[ TauLines[ipT7324].ipCont -1];
92 
93  /* opcxyz = opac(icaxyz) */
94  if( opcxyz > 0. )
95  {
96  d52 = 5.6*opcxyz/(opcxyz + opCax)*(1. - TauLines[ipT8498].Emis->Pesc);
97  d53 = 5.6*opcxyz/(opcxyz + opCay)*(1. - TauLines[ipT8542].Emis->Pesc);
98  d42 = 5.6*opcxyz/(opcxyz + opCaz)*(1. - TauLines[ipT8662].Emis->Pesc);
99  }
100  else
101  {
102  d52 = 0.;
103  d53 = 0.;
104  d42 = 0.;
105  }
106 
107  /* near UV dest of KH by background continuum */
108  opckh = opac.opacity_abs[ TauLines[ipT3969].ipCont -1];
109 
110  /* opckh = opac(icakh) */
111  if( opckh > 0. )
112  {
113  op51 = dense.xIonDense[ipCALCIUM][1]*3.89e-7/DoppVel.doppler[19];
114  d51 = 5.6*opckh/(opckh + op51);
115  op41 = dense.xIonDense[ipCALCIUM][1]*1.96e-7/DoppVel.doppler[19];
116  d41 = 5.6*opckh/(opckh + op41);
117  }
118  else
119  {
120  op51 = 0.;
121  d51 = 0.;
122  op41 = 0.;
123  d41 = 0.;
124  }
125  /* net rates */
126  r21 = PhotoRate2 + a21;
127  r31 = PhotoRate3 + a31;
128  r41 = a41 + PhotoRate4 + d41;
129  r51 = a51 + PhotoRate5 + d51;
130  r42 = a42 + d42;
131  r52 = a52 + d52;
132  r53 = a53 + d53;
133  cs14 = 0.923*phycon.te10*phycon.te10;
134  cs15 = cs14*2.;
135  TauLines[ipT3969].Coll.cs = (realnum)cs14;
136  TauLines[ipT3934].Coll.cs = (realnum)cs15;
137 
138  /* following used to correct rec contribution
139  * fcakh = a51 / ( a51 + eden*1.5e-5 / sqrte )
140  * cs 1-2 from
141  * >>refer ca2 cs Saraph, H.E. 1970, J.Phys. B, 3, 952
142  * other
143  * >>refer ca2 cs Chidichimo, M.C. 1981, J.Phys. B, 14, 4149 */
144  atom_pop5(gCa2,exCa2,5.8,8.6,cs14,cs15,20.6,22.9,9.8,3.4,44.4,1.0,
145  r21,r31,r41,r51,0.,r42,r52,0.,r53,0.,Ca2pop,dense.xIonDense[ipCALCIUM][1]);
146 
147  /* CDSQTE = 8.629E-6*EDEN/SQRTE */
148  c21 = 5.8/4.*dense.cdsqte;
149 
150  /* remember largest ratio of Ly-al removal to total */
151  if( dense.xIonDense[ipCALCIUM][1] > 0. )
152  {
153  ca.Ca2RmLya = MAX2(ca.Ca2RmLya,(realnum)(PhotoRate2/(PhotoRate2+a21+c21)));
154  }
155  ca.Cak = (realnum)(Ca2pop[4]*a51*5.06e-12);
156  ca.Cah = (realnum)(Ca2pop[3]*a41*5.01e-12);
157  ca.Cax = (realnum)(Ca2pop[4]*a52*2.34e-12);
158  ca.Cay = (realnum)(Ca2pop[4]*a53*2.33e-12);
159  ca.Caz = (realnum)(Ca2pop[3]*a42*2.30e-12);
160  ca.Caf1 = (realnum)(Ca2pop[2]*a31*2.73e-12);
161  ca.Caf2 = (realnum)(Ca2pop[1]*a21*2.72e-12);
162  ca.popca2ex = (realnum)(Ca2pop[1] + Ca2pop[2] + Ca2pop[3] + Ca2pop[4]);
163 
164  /* this is the total cooling due to the model atom */
165  ca.Cair = ca.Cax + ca.Cay + ca.Caz;
166  ca.c7306 = ca.Caf1 + ca.Caf2;
167  ca.Cakh = ca.Cak + ca.Cah;
168 
169  /* add the CaII lines to the cooling stack */
170  /* >>chng 04 feb 28, next three, deriv wrt temp had not been present
171  * extra factor of three in second term of deriv is due to extra T-1
172  * dependence of cooling on recombination coefficient. Almost always
173  * CaII cooling is most important when Ca+ is a trace stage of ionization
174  * so changes in temp change cooling due to change in Ca+ abundance */
175  CoolAdd("Ca 2",7306,ca.c7306);
177 
178  CoolAdd("Ca 2",8400,ca.Cair);
179  /* >>chng 04 feb 28, this is emperical fit to cooling deriv - does down with
180  * increasing temp */
182 
183  /* >>chng 04 feb 28 wl had been 3800, chng to 3950, closer to mean of K and R */
184  CoolAdd("Ca 2",3950,ca.Cakh);
185  /* >>chng 04 feb 28, this is emperical fit to cooling deriv - does down with
186  * increasing temp */
188 
189  /*fprintf(ioQQQ,"DEBUG ca2\t%.2f\t%.5e\t%.4e\t%.4e\n",
190  fnzone, phycon.te,ca.Cakh,dense.xIonDense[ipCALCIUM][1]);*/
191 
192  /* level populations that will be used for excited state photoionization */
193  ca.dstCala = (realnum)(Ca2pop[4]*PhotoRate5 + Ca2pop[3]*PhotoRate4);
194  ca.dCakh = (realnum)(ca.dstCala*5.03e-12);
195  ca.dCaf12 = (realnum)((Ca2pop[2]*PhotoRate3 + Ca2pop[1]*PhotoRate2)*2.31e-12);
196  opCax = (realnum)(Ca2pop[1]*1.13e-8/DoppVel.doppler[19]);
197  opCay = (realnum)(Ca2pop[2]*6.87e-8/DoppVel.doppler[19]);
198  opCaz = (realnum)(Ca2pop[1]*5.74e-8/DoppVel.doppler[19]);
199 
200  /* total rate Lalpha destroys CaII,
201  * this is only used in ioncali to increase ionization rate by
202  * adding it to the ct vector */
203  if( dense.xIonDense[ipCALCIUM][1] > 0. )
204  {
205  ca.dstCala = (realnum)(
206  (ca.dstCala + ca.dCaf12/2.31e-12)/dense.xIonDense[ipCALCIUM][1]);
207  /* take average of old and new */
208  ca.dstCala = (realnum)((ca.dstCala + ca.oldcla)/2.);
209  ca.oldcla = ca.dstCala;
210  {
211  /*@-redef@*/
212  enum {DEBUG_LOC=false};
213  /*@+redef@*/
214  if( DEBUG_LOC )
215  {
216  fprintf(ioQQQ," hlgam is %e\n", hlgam);
217  }
218  }
219  }
220  else
221  {
222  ca.dstCala = 0.;
223  }
224  ca.Ca3d = (realnum)(Ca2pop[1] + Ca2pop[2]);
225  ca.Ca4p = (realnum)(Ca2pop[3] + Ca2pop[4]);
226 
227  /* incl stimulated emission for Calcium II 5-level atom */
228  TauLines[ipT3934].Emis->PopOpc = (Ca2pop[0] - Ca2pop[4]/2.);
229  TauLines[ipT3934].Hi->Pop = Ca2pop[4];
230  TauLines[ipT3934].Lo->Pop = Ca2pop[0];
231  TauLines[ipT3969].Emis->PopOpc = (Ca2pop[0] - Ca2pop[3]);
232  TauLines[ipT3969].Hi->Pop = Ca2pop[3];
233  TauLines[ipT3969].Lo->Pop = Ca2pop[0];
234 
236  TauLines[ipT8498].Emis->PopOpc = (Ca2pop[1] - Ca2pop[4]);
237  TauLines[ipT8498].Hi->Pop = Ca2pop[4];
238  TauLines[ipT8498].Lo->Pop = Ca2pop[1];
239  TauLines[ipT8542].Emis->PopOpc = (Ca2pop[2] - Ca2pop[4]*1.5);
240  TauLines[ipT8542].Hi->Pop = Ca2pop[4];
241  TauLines[ipT8542].Lo->Pop = Ca2pop[2];
242  TauLines[ipT8662].Emis->PopOpc = (Ca2pop[1] - Ca2pop[3]*2.);
243  TauLines[ipT8662].Hi->Pop = Ca2pop[3];
244  TauLines[ipT8662].Lo->Pop = Ca2pop[1];
246  TauLines[ipT7291].Hi->Pop = 0.;
249  TauLines[ipT7324].Hi->Pop = 0.;
251 
252  /* Ca IV 3.2 micron; data from
253  * >>refer ca4 as Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103,
254  * >>refercon ed by D.R. Flower, (D. Reidel: Holland), 143
255  * Y(ik) from
256  * >>refer ca4 cs Pelan, J., & Berrington, K.A. 1995, A&A Suppl, 110, 209 */
257  if( phycon.te <= 1e5 )
258  {
259  cs = MAX2(1.0,8.854e-3*phycon.sqrte);
260  }
261  else if( phycon.te < 2.512e5 )
262  {
263  cs = 2.8;
264  }
265  else
266  {
267  cs = 641.1/(phycon.te30*phycon.te10*phycon.te02*phycon.te02/
268  phycon.te003);
269  }
270  PutCS(cs,&TauLines[ipTCa3]);
271  atom_level2(&TauLines[ipTCa3]);
272 
273  /* [Ca V] IR 4.16, 11.47 micron; A from
274  * >>refer ca5 as Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103,
275  * >>refercon ed by D.R. Flower, (D. Reidel: Holland), 143
276  * cs from
277  * >>refer ca5 cs Galavis, M.E., Mendoza, C., & Zeippen, C.J. 1995, A&AS, 111, 347
278  * >>chng 96 jul 16, big changes in cs */
279  cs = MIN2(3.3,0.392*phycon.te20/phycon.te005/phycon.te003);
280  cs = MAX2(2.2,cs);
281  PutCS(cs,&TauLines[ipTCa4]);
282 
283  /* >>chng 96 aug 02, following had error in te dep in ver 90.01 */
284  cs = MIN2(0.93,0.162*phycon.te10*phycon.te05*phycon.te003*
285  phycon.te001);
286  cs = MAX2(0.67,cs);
287  PutCS(cs,&TauLines[ipTCa12]);
288 
289  cs = MIN2(0.97,0.0894*phycon.te20*phycon.te01*phycon.te005);
290  cs = MAX2(0.60,cs);
291  PutCS(cs,&TauDummy);
292 
293  atom_level3(&TauLines[ipTCa4],&TauLines[ipTCa12],&TauDummy);
294 
295  /* Ca V lines from 1d, 1s; A from
296  * >>refer ca5 as Mendoza, C. 1982, in Planetary Nebulae, IAU Symp No. 103,
297  * >>refercon ed by D.R. Flower, (D. Reidel: Holland), 143
298  * cs from
299  * >>refer ca5 cs Galavis, M.E., Mendoza, C., & Zeippen, C.J. 1995, A&AS, 111, 347
300  * POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2) */
301  cs01 = MIN2(4.1,0.533*phycon.te20/phycon.te01);
302  cs01 = MAX2(2.8,cs01);
303  cs02 = MIN2(0.87,5.22e-03*phycon.sqrte);
304  p3 = atom_pop3(9.,5.,1.,cs01,cs02,1.35,2.326,23.2,3.73,2.57e4,3.60e4,
305  &p2,dense.xIonDense[ipCALCIUM][4],0.,0.,0.);
306 
307  ca.c3997 = p3*3.73*4.98e-12;
308  ca.c2414 = p3*23.1*8.245e-12;
309  ca.Ca6087 = p2*0.426*3.268e-12;
310  ca.c5311 = p2*1.90*3.747e-12;
311 
312  CoolAdd("Ca 5",3997,ca.c3997);
313  CoolAdd("Ca 5",2414,ca.c2414);
314  CoolAdd("Ca 5",6087,ca.Ca6087);
315  CoolAdd("Ca 5",5311,ca.c5311);
316 
317  /* Ca VII lines from 1d, 1s
318  * all cs from
319  * >>refer ca72 cs Galavis, M.E., Mendoza, C., & Zeippen, C.J. 1995, A&AS, 111, 347
320  * POP3(G1,G2,G3,O12,O13,O23,A21,A31,A32,E12,E23,P2,ABUND,GAM2) */
321  cs01 = MIN2(4.4,22.25/(phycon.te20/phycon.te02/phycon.te02));
322  cs01 = MAX2(3.5,cs01);
323 
324  cs12 = MIN2(1.20,0.303*phycon.te30*phycon.te03);
325  cs12 = MAX2(0.62,cs12);
326 
327  cs02 = MIN2(0.959,7.889/(phycon.te20*phycon.te05/phycon.te01));
328  cs02 = MAX2(0.50,cs02);
329 
330  p3 = atom_pop3(9.,5.,1.,cs01,cs02,cs12,3.124,30.4,6.81,2.91e4,3.90e4,
331  &p2,dense.xIonDense[ipCALCIUM][6],0.,0.,0.);
332 
333  ca.Ca3688 = p3*6.81*5.40e-12;
334  ca.Ca2112 = p3*30.4*9.42e-12;
335  ca.Ca5620 = p2*2.15*3.548e-12;
336  ca.Ca4941 = p2*0.974*4.037e-12;
337  CoolAdd("Ca 7",3688,ca.Ca3688);
338  CoolAdd("Ca 7",2112,ca.Ca2112);
339  CoolAdd("Ca 7",5620,ca.Ca5620);
340  CoolAdd("Ca 7",4941,ca.Ca4941);
341 
342  /* all cs from
343  * >>refer ca7 cs Galavis, M.E., Mendoza, C., & Zeippen, C.J. 1995, A&AS, 111, 347
344  * [Ca VII] 4.09, 6.15 mic 3P lines */
345  cs = MIN2(5.354,0.406*phycon.te20*phycon.te03*phycon.te01);
346  cs = MAX2(3.702,cs);
347  PutCS(cs,&TauLines[ipCa0741]);
348 
349  cs = MIN2(1.59,0.183*phycon.te20);
350  cs = MAX2(1.153,cs);
351  PutCS(cs,&TauLines[ipCa0761]);
352 
353  cs = MIN2(1.497,0.0917*phycon.te20*phycon.te05* phycon.te01);
354  cs = MAX2(1.005,cs);
355  PutCS(cs,&TauDummy);
356 
357  /* atom_level3( t10,t21,t20) */
358  atom_level3(&TauLines[ipCa0761],&TauLines[ipCa0741],&TauDummy);
359 
360  /* [Ca VIII] 2.32 microns, cs
361  * >>refer ca8 cs Saraph, H.E., & Storey, P.J. A&AS, 115, 151 */
362  cs = MIN2(6.75,22.04/(phycon.te10*phycon.te02*phycon.te005));
363  PutCS(cs,&TauLines[ipCa08232]);
364  atom_level2(&TauLines[ipCa08232]);
365 
366  /* [Ca 12] 3328.78A, cs from
367  * >>refer ca12 cs Saraph, H.E. & Tully, J.A. 1994, A&AS, 107, 29 */
368  cs = MIN2(0.172,0.0118*phycon.te20*phycon.te01);
369  cs = MAX2(0.10,cs);
370  PutCS(cs,&TauLines[ipCa12333]);
371  atom_level2(&TauLines[ipCa12333]);
372 
373  /* Li seq Ca 18 2s2p 2s3p, 2s2p as two separate lines
374  * >>refer ca18 cs Cochrane, D.M., & McWhirter, R.W.P. 1983, PhyS, 28, 25 */
375  ligbar(20,&TauLines[ipTCa302],&TauLines[ipTCa19],&cs2s2p,&cs2s3p);
376 
377  PutCS(cs2s2p,&TauLines[ipTCa302]);
378  atom_level2(&TauLines[ipTCa302]);
379 
380  /* funny factor (should have been 0.5) due to energy change */
381  PutCS(cs2s2p*0.439,&TauLines[ipTCa345]);
382  atom_level2(&TauLines[ipTCa345]);
383 
384  PutCS(cs2s3p,&TauLines[ipTCa19]);
385  atom_level2(&TauLines[ipTCa19]);
386  return;
387 }

Generated for cloudy by doxygen 1.8.1.2