cloudy trunk
Loading...
Searching...
No Matches
cool_dima.cpp
Go to the documentation of this file.
1/* This file is part of Cloudy and is copyright (C)1978-2013 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 "rt.h"
10#include "doppvel.h"
11#include "phycon.h"
12#include "conv.h"
13#include "thermal.h"
14#include "opacity.h"
15#include "lines_service.h"
16#include "rfield.h"
17#include "mewecoef.h"
18#include "atoms.h"
19#include "cooling.h"
20#include "atmdat.h"
21
22/*ColStrGBar generate g-bar collision strengths for level 2 line2 */
24
25void CoolDima(void)
26{
27 long int i,
28 ion,
29 nelem;
30 double cs;
31
32 DEBUG_ENTRY( "CoolDima()" );
33
34 /* no level2 command sets nWindLine to -1 */
35 if( nWindLine<0 )
36 return;
37
38 for( i=0; i < nWindLine; i++ )
39 {
40 ion = (*TauLine2[i].Hi()).IonStg();
41 nelem = (*TauLine2[i].Hi()).nelem();
42
43 if( (dense.lgIonChiantiOn[nelem-1][ion-1] && !atmdat.lgChiantiHybrid) ||
44 (dense.lgIonStoutOn[nelem-1][ion-1] && !atmdat.lgStoutHybrid) )
45 {
46 /* If a species uses Chianti or Stout and hybrid is off, skip the level 2 lines */
47 continue;
48 }
49 // iso sequence ions are done elsewhere
50 if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO &&
51 // dense.maxWN is positive if hybrid chianti is turned on and this element is included
52 // in CloudyChianti.ini - zero otherwise
53 TauLine2[i].EnergyWN() >= dense.maxWN[nelem-1][ion-1])
54 {
55 /* only evaluate cs if positive abundance */
56 if( dense.xIonDense[nelem-1][ion-1] > 0. )
57 {
58 /* now generate the collision strength */
59 cs = ColStrGBar(TauLine2.begin()+i , cs1_flag_lev2[i] );
60 }
61 else
62 {
63 cs = 1.;
64 }
65 /* now put the cs into the line array */
66 PutCS(cs,TauLine2[i] );
67 RT_line_one( TauLine2[i], true,0.f, GetDopplerWidth(dense.AtomicWeight[(*TauLine2[i].Hi()).nelem()-1]) );
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 = pow((double)10,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/pow(10.,(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]*(pow(10.,(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
void atom_level2(const TransitionProxy &t)
#define ASSERT(exp)
Definition cddefines.h:578
sys_float sexp(sys_float x)
Definition service.cpp:914
#define STATIC
Definition cddefines.h:97
#define POW2
Definition cddefines.h:929
const int NISO
Definition cddefines.h:261
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
long nWindLine
Definition cdinit.cpp:19
TransitionProxy::iterator iterator
Definition transition.h:280
ProxyIterator< TransitionProxy, TransitionConstProxy > iterator
Definition transition.h:27
void CoolDima(void)
Definition cool_dima.cpp:25
STATIC double ColStrGBar(const TransitionProxy::iterator &t, realnum cs1)
Definition cool_dima.cpp:76
t_dense dense
Definition dense.cpp:24
realnum GetDopplerWidth(realnum massAMU)
t_MeweCoef MeweCoef
Definition mewecoef.cpp:5
t_phycon phycon
Definition phycon.cpp:6
void RT_line_one(const TransitionProxy &t, bool lgShield_this_zone, realnum pestrk, realnum DopplerWidth)
TransitionList TauLine2("TauLine2", &AnonStates)
realnum * cs1_flag_lev2
Definition taulines.cpp:40
void PutCS(double cs, const TransitionProxy &t)