cloudy  trunk
 All Data Structures Namespaces 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-2017 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 "taulines.h"
7 #include "dense.h"
8 #include "rt.h"
9 #include "doppvel.h"
10 #include "phycon.h"
11 #include "mewecoef.h"
12 #include "atoms.h"
13 #include "atmdat.h"
14 #include "thermal.h"
15 #include "cooling.h"
16 
17 /*ColStrGBar generate g-bar collision strengths for level 2 line2 */
18 STATIC double ColStrGBar(const TransitionProxy::iterator& t , realnum cs1 );
19 
20 void CoolDima(void)
21 {
22  long int i,
23  ion,
24  nelem;
25  double cs;
26 
27  DEBUG_ENTRY( "CoolDima()" );
28 
29  /* no level2 command sets nWindLine to -1 */
30  if( nWindLine<0 )
31  return;
32 
33  static vector<realnum> DopplerWidth(LIMELM);
34  for(nelem = ipHYDROGEN; nelem < LIMELM; ++nelem)
35  DopplerWidth[nelem] = GetDopplerWidth(dense.AtomicWeight[nelem]);
36 
37  for( i=0; i < nWindLine; i++ )
38  {
39  ion = (*TauLine2[i].Hi()).IonStg();
40  nelem = (*TauLine2[i].Hi()).nelem();
41 
42  if( (dense.lgIonChiantiOn[nelem-1][ion-1] && !atmdat.lgChiantiHybrid) ||
43  (dense.lgIonStoutOn[nelem-1][ion-1] && !atmdat.lgStoutHybrid) )
44  {
45  /* If a species uses Chianti or Stout and hybrid is off, skip the level 2 lines */
46  continue;
47  }
48  // iso sequence ions are done elsewhere
49  if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO &&
50  // dense.maxWN is positive if hybrid chianti is turned on and this element is included
51  // in CloudyChianti.ini - zero otherwise
52  TauLine2[i].EnergyWN() >= dense.maxWN[nelem-1][ion-1])
53  {
54  /* only evaluate cs if positive abundance */
55  if( dense.xIonDense[nelem-1][ion-1] > 0. )
56  {
57  /* now generate the collision strength */
58  cs = ColStrGBar(TauLine2.begin()+i , cs1_flag_lev2[i] );
59  }
60  else
61  {
62  cs = 1.;
63  }
64  /* now put the cs into the line array */
65  PutCS(cs,TauLine2[i] );
66  RT_line_one_escape( TauLine2[i], true,0.f, DopplerWidth[(*TauLine2[i].Hi()).nelem()-1] );
67  atom_level2(TauLine2[i] );
68  thermal.dima += TauLine2[i].Coll().cool();
69  }
70  }
71 
72  return;
73 }
74 
75 /*ColStrGBar generate g-bar collision strengths for level 2 line2 */
77 {
78  long int i,
79  j;
80  double ColStrGBar_v,
81  a,
82  b,
83  c,
84  d,
85  e1,
86  gb,
87  x,
88  y;
89  double xx,
90  yy;
91 
92  DEBUG_ENTRY( "ColStrGBar()" );
93 
94  /* Calculation of the collision strengths of multiplets.
95  * Neutrals are recalculated from
96  * >>refer cs gbar Fisher et al. (1996)
97  * >>refer cs gbar Gaetz & Salpeter (1983, ApJS 52, 155) and
98  * >>refer cs gbar Mewe (1972, A&A 20, 215)
99  * fits for ions. */
100 
101  /* routine to implement g-bar data taken from
102  *>>refer cs gbar Mewe, R.; Gronenschild, E. H. B. M.; van den Oord, G. H. J., 1985,
103  *>>refercon A&AS, 62, 197 */
104 
105  /* zero hydrogenic lines since they are done by iso-sequence */
106  if( (*(*t).Hi()).nelem() == (*(*t).Hi()).IonStg() )
107  {
108  ColStrGBar_v = 0.0;
109  return( ColStrGBar_v );
110  }
111 
112  /*was the block data linked in? */
113  ASSERT( MeweCoef.g[1][0] != 0.);
114 
115  /* which type of transition is this? cs1 is the flag */
116 
117  /* >>chng 01 may 30 - cs1 < 0 means a forced collision strength */
118  if( cs1 < 0. )
119  {
120  ColStrGBar_v = -cs1;
121  return( ColStrGBar_v );
122  }
123 
124  /* >>chng 99 feb 27, change to assert */
125  ASSERT( cs1 >= 0.05 );
126 
127  /* excitation energy over kT */
128  y = (*t).EnergyK()/phycon.te;
129  if( cs1 < 1.5 )
130  {
131  xx = -log10(y);
132 
133  if( cs1 < 0.5 )
134  {
135  yy = (1.398813573838321 + xx*(0.02943050869177121 + xx*
136  (-0.4439783893114510 + xx*(0.2316073358577902 + xx*(0.001870493481643103 +
137  xx*(-0.008227246351067403))))))/(1.0 + xx*(-0.6064792600526370 +
138  xx*(0.1958559534507252 + xx*(-0.02110452007196644 +
139  xx*(0.01348743933722316 + xx*(-0.0001944731334371711))))));
140  }
141 
142  else
143  {
144  yy = (1.359675968512206 + xx*(0.04636500015069853 + xx*
145  (-0.4491620298246676 + xx*(0.2498199231048967 + xx*(0.005053803073345794 +
146  xx*(-0.01015647880244268))))))/(1.0 + xx*(-0.5904799485819767 +
147  xx*(0.1877833737815317 + xx*(-0.01536634911179847 +
148  xx*(0.01530712091180953 + xx*(-0.0001909176790831023))))));
149  }
150 
151  ColStrGBar_v = exp10(yy)*(*t).Emis().gf()/((*t).EnergyRyd() * 13.6);
152  }
153  else
154  {
155  i = (long int)cs1;
156 
157  if( i < 26 )
158  {
159  e1 = log(1.0+1.0/y) - 0.4/POW2(y + 1.0);
160  a = MeweCoef.g[i-1][0];
161  b = MeweCoef.g[i-1][1];
162  c = MeweCoef.g[i-1][2];
163  d = MeweCoef.g[i-1][3];
164  x = (double)(*(*t).Hi()).nelem() - 3.0;
165 
166  if( i == 14 )
167  {
168  a *= 1.0 - 0.5/x;
169  b = 1.0 - 0.8/x;
170  c *= 1.0 - 1.0/x;
171  }
172 
173  else if( i == 16 )
174  {
175  a *= 1.0 - 0.9/x;
176  b *= 1.0 - 1.7/x;
177  c *= 1.0 - 2.1/x;
178  }
179 
180  else if( i == 18 )
181  {
182  a *= 1.0 + 2.0/x;
183  b *= 1.0 - 0.7/x;
184  }
185 
186  gb = a + (b*y - c*y*y + d)*e1 + c*y;
187 
188  /* ipLnRyd is exciation energy in Rydbergs */
189  ColStrGBar_v = 14.510395*(*t).Emis().gf()*gb/((*t).EnergyRyd() );
190  /* following i>=26 */
191  }
192 
193  else
194  {
195  /* 210 is the dimem of g, so [209] is largest val */
196  if( i < 210 )
197  {
198  j = (long)(MeweCoef.g[i-1][3]);
199  if( j == 1 )
200  {
201  ColStrGBar_v = (*(*t).Lo()).g()*MeweCoef.g[i-1][0]*
202  pow(phycon.te/exp10((double)MeweCoef.g[i-1][2]),(double)MeweCoef.g[i-1][1]);
203  }
204  else
205  {
206  ColStrGBar_v = (*(*t).Lo()).g()*MeweCoef.g[i-1][0]*
207  sexp(MeweCoef.g[i-1][1]*(exp10((double)MeweCoef.g[i-1][2])/
208  phycon.te));
209  }
210  }
211 
212  else
213  {
214  /* This is for AlII 1670 line only!
215  * ColStrGBar=0.0125*te**0.603 */
216  /* 98 dec 27, this is still in use */
217  ColStrGBar_v = 0.0125*phycon.sqrte*phycon.te10*
218  phycon.te003;
219  }
220  }
221  }
222 
223  /* following to make sure that negative values not returned */
224  ColStrGBar_v = MAX2(ColStrGBar_v,1e-10);
225  return( ColStrGBar_v );
226 }
t_atmdat atmdat
Definition: atmdat.cpp:6
t_thermal thermal
Definition: thermal.cpp:6
bool lgStoutHybrid
Definition: atmdat.h:404
double exp10(double x)
Definition: cddefines.h:1383
const int NISO
Definition: cddefines.h:310
double e1(double x)
void RT_line_one_escape(const TransitionProxy &t, bool lgShield_this_zone, realnum pestrk, realnum DopplerWidth)
t_phycon phycon
Definition: phycon.cpp:6
iterator begin(void)
Definition: transition.h:339
sys_float sexp(sys_float x)
Definition: service.cpp:1095
realnum g[210][4]
Definition: mewecoef.h:13
TransitionList TauLine2("TauLine2",&AnonStates)
t_dense dense
Definition: global.cpp:15
double te003
Definition: phycon.h:58
bool lgIonStoutOn[LIMELM][LIMELM+1]
Definition: dense.h:143
double xIonDense[LIMELM][LIMELM+1]
Definition: dense.h:135
bool lgIonChiantiOn[LIMELM][LIMELM+1]
Definition: dense.h:140
bool lgChiantiHybrid
Definition: atmdat.h:376
void atom_level2(const TransitionProxy &t, const bool lgHFS)
Definition: atom_level2.cpp:30
void CoolDima(void)
Definition: cool_dima.cpp:20
#define POW2
Definition: cddefines.h:983
void PutCS(double cs, const TransitionProxy &t)
Definition: transition.cpp:296
#define STATIC
Definition: cddefines.h:118
float realnum
Definition: cddefines.h:124
double dima
Definition: thermal.h:118
realnum GetDopplerWidth(realnum massAMU)
long nWindLine
Definition: cdinit.cpp:19
realnum * cs1_flag_lev2
Definition: taulines.cpp:37
realnum AtomicWeight[LIMELM]
Definition: dense.h:80
#define ASSERT(exp)
Definition: cddefines.h:617
const int LIMELM
Definition: cddefines.h:307
#define DEBUG_ENTRY(funcname)
Definition: cddefines.h:729
double te10
Definition: phycon.h:58
#define MAX2(a, b)
Definition: cddefines.h:828
double pow(double x, int i)
Definition: cddefines.h:782
double maxWN[LIMELM][LIMELM+1]
Definition: dense.h:146
STATIC double ColStrGBar(const TransitionProxy::iterator &t, realnum cs1)
Definition: cool_dima.cpp:76
double sqrte
Definition: phycon.h:58
double te
Definition: phycon.h:21
t_MeweCoef MeweCoef
Definition: mewecoef.cpp:5
const int ipHYDROGEN
Definition: cddefines.h:348