cloudy  trunk
 All Data Structures Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cool_dima.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 /*CoolDima compute cooling due to level 2 lines */
4 /*ColStrGBar generate g-bar collision strengths for level 2 line2 */
5 #include "cddefines.h"
6 #include "physconst.h"
7 #include "taulines.h"
8 #include "dense.h"
9 #include "phycon.h"
10 #include "conv.h"
11 #include "thermal.h"
12 #include "opacity.h"
13 #include "lines_service.h"
14 #include "rfield.h"
15 #include "mewecoef.h"
16 #include "atoms.h"
17 #include "cooling.h"
18 
19 /*ColStrGBar generate g-bar collision strengths for level 2 line2 */
20 STATIC double ColStrGBar(transition * t , realnum cs1 );
21 
22 void CoolDima(void)
23 {
24  long int i,
25  ion,
26  nelem;
27  double cs;
28  double OTSLevel2,
29  sumots,
30  CoolLevel2;
31  static double TeEvalCS = -1.;
32  static long int nzoneEval=0;
33 
34  DEBUG_ENTRY( "CoolDima()" );
35 
36  /* >>chng 03 nov 29, add this option to skip rest of routine */
37  /* no level2 command sets nWindLine to -1 */
38  if( nWindLine<0 )
39  return;
40 
41  /* only force evaluation of cooling when significant, or first call
42  * in this zone, or not into first zone */
43  if( nzone != nzoneEval || conv.lgLevel2_Cool_Imp || !nzone )
44  {
45  nzoneEval = nzone;
46 
47  /* check whether we need to reevaluate the collision strenghts.
48  * Do so if large change in temperature or any stage of ionization has been lowered. */
50  fabs(phycon.te-TeEvalCS)/phycon.te > 0.05) )
51  {
52  for( i=0; i < nWindLine; i++ )
53  {
54  if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
55  {
56  /* only evaluate cs if positive abundance */
57  ion = TauLine2[i].Hi->IonStg;
58  nelem = TauLine2[i].Hi->nelem;
59 
60  if( dense.xIonDense[nelem-1][ion-1] > 0. )
61  {
62  /* now generate the collision strength */
63  cs = ColStrGBar(&TauLine2[i] , cs1_flag_lev2[i] );
64  }
65  else
66  {
67  cs = 1.;
68  }
69  /* now put the cs into the line array */
70  PutCS(cs,&TauLine2[i] );
71  }
72  }
73  TeEvalCS = phycon.te;
74  }
75 
76  /* now always update the cooling since this adds ots flux */
77  for( i=0; i < nWindLine; i++ )
78  {
79  /* we need to call this even if nothing present since then sets zero
80  * ignore all lines with negative atomic number
81  * >>chng 97 aug 30, get rid of test for non-pos ipLnNelem
82  * pointer tells atom_level2 to use existing WineEtc labels */
83  /* only call non-hydrogenic, non-he-like lines */
84  /* >>chng 02 sug 11, do not include he-like in sum */
85  if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
86  {
87  atom_level2(&TauLine2[i] );
88  }
89  }
90  }
91  else
92  {
93  /* now use old cooling and ots flux */
94  for( i=0; i < nWindLine; i++ )
95  {
96  if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
97  {
98  /* recall that these should be trivial cooling, ots rates */
99  CoolAdd( " ", TauLine2[i].WLAng , TauLine2[i].Coll.cool);
100  /*>>chng 03 oct 04, move to RT_OTS */
101  /*RT_OTS_AddLine( TauLine2[i].ots , TauLine2[i].ipCont);*/
102  }
103  }
104  }
105 
106  /* find ratio of level 1 to level 2 ots to see how important level 2 is.
107  * if level 2 ots not important then will update dest probs only first time
108  * in each zone. If significant then continuously update */
109  OTSLevel2 = 0.;
110  CoolLevel2 = 0.;
111  if( nzone > 1 )
112  {
113  for( i=0; i < nWindLine; i++ )
114  {
115  long int ip = TauLine2[i].ipCont-1;
116  CoolLevel2 += TauLine2[i].Coll.cool;
117  if( opac.opacity_abs[ip] > SMALLFLOAT )
118  {
119  OTSLevel2 += TauLine2[i].Emis->ots / ( (realnum)opac.opacity_abs[ip]);
120  ASSERT( fabs(OTSLevel2) < 1e37 );
121  }
122  }
123 
124  /* now sum all ots rates */
125  sumots = 0.;
126  for( i=0; i < rfield.nflux; i++ )
127  {
128  sumots += rfield.otslin[i];
129  }
130  /* is the ots contribution significant ? */
131  if( OTSLevel2/SDIV(sumots) > 1e-4 )
132  {
133  /* true if level 2 lines were contributors to the ots rates, set in dimacool */
134  conv.lgLevel2_OTS_Imp = true;
135  }
136  else
137  {
138  /* true if level 2 lines were contributors to the ots rates, set in dimacool */
139  conv.lgLevel2_OTS_Imp = false;
140  }
141 
142  /* set flags saying if important */
143  if( CoolLevel2/SDIV(thermal.ctot) > 1e-4 )
144  {
145  /* true if level 2 lines were contributors to the cooling, set in dimacool */
146  conv.lgLevel2_Cool_Imp = true;
147  }
148  else
149  {
150  /* true if level 2 lines were contributors to the cooling, set in dimacool */
151  conv.lgLevel2_Cool_Imp = false;
152  }
153 
154  {
155  /* debugging code for Lya problems */
156  /*@-redef@*/
157  enum {DEBUG_LOC=false};
158  /*@+redef@*/
159  if( DEBUG_LOC )
160  {
161  fprintf(ioQQQ,"%li\t%.2e\t%.2e\n",
162  nzone, OTSLevel2/SDIV(sumots) , CoolLevel2/SDIV(thermal.ctot) );
163  }
164  }
165  }
166  else
167  {
168  /* true if level 2 lines were contributors to the ots rates, set in dimacool */
169  conv.lgLevel2_OTS_Imp = true;
170  /* true if level 2 lines were contributors to the cooling, set in dimacool */
171  conv.lgLevel2_Cool_Imp = true;
172  }
173  return;
174 }
175 
176 /*ColStrGBar generate g-bar collision strengths for level 2 line2 */
178 {
179  long int i,
180  j;
181  double ColStrGBar_v,
182  a,
183  b,
184  c,
185  d,
186  e1,
187  gb,
188  x,
189  y;
190  double xx,
191  yy;
192 
193  DEBUG_ENTRY( "ColStrGBar()" );
194 
195  /* Calculation of the collision strengths of multiplets.
196  * Neutrals are recalculated from
197  * >>refer cs gbar Fisher et al. (1996)
198  * >>refer cs gbar Gaetz & Salpeter (1983, ApJS 52, 155) and
199  * >>refer cs gbar Mewe (1972, A&A 20, 215)
200  * fits for ions. */
201 
202  /* routine to implement g-bar data taken from
203  *>>refer cs gbar Mewe, R.; Gronenschild, E. H. B. M.; van den Oord, G. H. J., 1985,
204  *>>refercon A&AS, 62, 197 */
205 
206  /* zero hydrogenic lines since they are done by iso-sequence */
207  if( t->Hi->nelem == t->Hi->IonStg )
208  {
209  ColStrGBar_v = 0.0;
210  return( ColStrGBar_v );
211  }
212 
213  /*was the block data linked in? */
214  ASSERT( MeweCoef.g[1][0] != 0.);
215 
216  /* which type of transition is this? cs1 is the flag */
217 
218  /* >>chng 01 may 30 - cs1 < 0 means a forced collision strength */
219  if( cs1 < 0. )
220  {
221  ColStrGBar_v = -cs1;
222  return( ColStrGBar_v );
223  }
224 
225  /* >>chng 99 feb 27, change to assert */
226  ASSERT( cs1 >= 0.05 );
227 
228  /* excitation energy over kT */
229  y = t->EnergyK/phycon.te;
230  if( cs1 < 1.5 )
231  {
232  xx = -log10(y);
233 
234  if( cs1 < 0.5 )
235  {
236  yy = (1.398813573838321 + xx*(0.02943050869177121 + xx*
237  (-0.4439783893114510 + xx*(0.2316073358577902 + xx*(0.001870493481643103 +
238  xx*(-0.008227246351067403))))))/(1.0 + xx*(-0.6064792600526370 +
239  xx*(0.1958559534507252 + xx*(-0.02110452007196644 +
240  xx*(0.01348743933722316 + xx*(-0.0001944731334371711))))));
241  }
242 
243  else
244  {
245  yy = (1.359675968512206 + xx*(0.04636500015069853 + xx*
246  (-0.4491620298246676 + xx*(0.2498199231048967 + xx*(0.005053803073345794 +
247  xx*(-0.01015647880244268))))))/(1.0 + xx*(-0.5904799485819767 +
248  xx*(0.1877833737815317 + xx*(-0.01536634911179847 +
249  xx*(0.01530712091180953 + xx*(-0.0001909176790831023))))));
250  }
251 
252  ColStrGBar_v = pow((double)10,yy)*t->Emis->gf/(t->EnergyWN * WAVNRYD* 13.6);
253  }
254  else
255  {
256  i = (long int)cs1;
257 
258  if( i < 26 )
259  {
260  e1 = log(1.0+1.0/y) - 0.4/POW2(y + 1.0);
261  a = MeweCoef.g[i-1][0];
262  b = MeweCoef.g[i-1][1];
263  c = MeweCoef.g[i-1][2];
264  d = MeweCoef.g[i-1][3];
265  x = (double)t->Hi->nelem - 3.0;
266 
267  if( i == 14 )
268  {
269  a *= 1.0 - 0.5/x;
270  b = 1.0 - 0.8/x;
271  c *= 1.0 - 1.0/x;
272  }
273 
274  else if( i == 16 )
275  {
276  a *= 1.0 - 0.9/x;
277  b *= 1.0 - 1.7/x;
278  c *= 1.0 - 2.1/x;
279  }
280 
281  else if( i == 18 )
282  {
283  a *= 1.0 + 2.0/x;
284  b *= 1.0 - 0.7/x;
285  }
286 
287  gb = a + (b*y - c*y*y + d)*e1 + c*y;
288 
289  /* ipLnRyd is exciation energy in Rydbergs */
290  ColStrGBar_v = 14.510395*t->Emis->gf*gb/(t->EnergyWN * WAVNRYD);
291  /* following i>=26 */
292  }
293 
294  else
295  {
296  /* 210 is the dimem of g, so [209] is largest val */
297  if( i < 210 )
298  {
299  j = (long)(MeweCoef.g[i-1][3]);
300  if( j == 1 )
301  {
302  ColStrGBar_v = t->Lo->g*MeweCoef.g[i-1][0]*
303  pow(phycon.te/pow(10.,(double)MeweCoef.g[i-1][2]),(double)MeweCoef.g[i-1][1]);
304  }
305  else
306  {
307  ColStrGBar_v = t->Lo->g*MeweCoef.g[i-1][0]*
308  sexp(MeweCoef.g[i-1][1]*(pow(10.,(double)MeweCoef.g[i-1][2])/
309  phycon.te));
310  }
311  }
312 
313  else
314  {
315  /* This is for AlII 1670 line only!
316  * ColStrGBar=0.0125*te**0.603 */
317  /* 98 dec 27, this is still in use */
318  ColStrGBar_v = 0.0125*phycon.sqrte*phycon.te10*
319  phycon.te003;
320  }
321  }
322  }
323 
324  /* following to make sure that negative values not returned */
325  ColStrGBar_v = MAX2(ColStrGBar_v,1e-10);
326  return( ColStrGBar_v );
327 }

Generated for cloudy by doxygen 1.8.1.1