cloudy trunk
Loading...
Searching...
No Matches
iso_collide.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#include "cddefines.h"
4#include "atmdat.h"
5#include "conv.h"
6#include "dense.h"
7#include "heavy.h"
8#include "helike_cs.h"
9#include "hydroeinsta.h"
10#include "hydrogenic.h"
11#include "hydro_vs_rates.h"
12#include "ionbal.h"
13#include "iso.h"
14#include "opacity.h"
15#include "phycon.h"
16#include "physconst.h"
17#include "rfield.h"
18#include "secondaries.h"
19#include "trace.h"
20#include "taulines.h"
21
22/* These are masses relative to the proton mass of the electron, proton, he+, and alpha particle. */
23static double ColliderMass[4] = {ELECTRON_MASS/PROTON_MASS, 1.0, 4.0, 4.0};
24
25void iso_collisional_ionization( long ipISO, long nelem )
26{
27 ASSERT( ipISO < NISO );
28
29 DEBUG_ENTRY( "iso_collisional_ionization()" );
30
31 t_iso_sp* sp = &iso_sp[ipISO][nelem];
32
33 /* collisional ionization from ground */
34 sp->fb[0].ColIoniz = iso_ctrl.lgColl_ionize[ipISO] *
35 t_ADfA::Inst().coll_ion_wrapper( nelem, nelem-ipISO, phycon.te );
36
37 iso_put_error(ipISO,nelem,sp->numLevels_max,0,IPCOLLIS,0.20f,0.20f);
38
39 for( long ipHi=1; ipHi<sp->numLevels_max; ipHi++ )
40 {
41 if( nelem == ipISO )
42 {
43 /* use routine from Vriens and Smeets (1981). */
44 /* >>refer iso neutral col.ion. Vriens, L., & Smeets, A.H.M. 1980, Phys Rev A 22, 940 */
45 sp->fb[ipHi].ColIoniz = hydro_vs_ioniz( sp->fb[ipHi].xIsoLevNIonRyd, phycon.te );
46 }
47 else
48 {
49 /* ions */
50 /* use hydrogenic ionization rates for ions
51 * >>refer iso ions col.ion. Allen 1973, Astro. Quan. for low Te.
52 * >>refer iso ions col.ion. Sampson and Zhang 1988, ApJ, 335, 516 for High Te.
53 * */
54 sp->fb[ipHi].ColIoniz =
55 Hion_coll_ioniz_ratecoef( ipISO, nelem, N_(ipHi), sp->fb[ipHi].xIsoLevNIonRyd, phycon.te );
56 }
57
58 // iso_ctrl.lgColl_ionize is option to turn off collisions, "atom XX-like collis off" command
59 sp->fb[ipHi].ColIoniz *= iso_ctrl.lgColl_ionize[ipISO];
60
61 iso_put_error(ipISO,nelem,sp->numLevels_max,ipHi,IPCOLLIS,0.20f,0.20f);
62 }
63
64 /* Here we arbitrarily scale the highest level ionization to account for the fact
65 * that, if the atom is not full size, this level should be interacting with higher
66 * levels and not just the continuum. We did add on collisional excitation terms instead
67 * but this caused a problem at low temperatures because the collisional ionization was
68 * a sum of terms with different Boltzmann factors, while PopLTE had just one Boltzmann
69 * factor. The result was a collisional recombination that had residual exponentials of
70 * the form exp(x/kT), which blew up at small T. */
71 if( 0 && !sp->lgLevelsLowered )
72 {
73 sp->fb[sp->numLevels_max-1].ColIoniz *= 100.;
74 }
75
76 return;
77}
78
79void iso_suprathermal( long ipISO, long nelem )
80{
81 DEBUG_ENTRY( "iso_suprathermal()" );
82
83 /* check that we were called with valid parameters */
84 ASSERT( ipISO < NISO );
85 ASSERT( nelem >= ipISO );
86 ASSERT( nelem < LIMELM );
87
88 t_iso_sp* sp = &iso_sp[ipISO][nelem];
89
90 /***********************************************************************
91 * *
92 * get secondary excitation by suprathermal electrons *
93 * *
94 ***********************************************************************/
95
96 for( long i=1; i < sp->numLevels_max; i++ )
97 {
98 if( sp->trans(i,0).ipCont() > 0 )
99 {
100 /* get secondaries for all iso lines by scaling LyA
101 * excitation by ratio of cross section (oscillator strength/energy)
102 * Born approximation or plane-wave approximation based on
103 *>>refer HI excitation Shemansky, D.E., et al., 1985, ApJ, 296, 774 */
104 sp->trans(i,0).Coll().rate_lu_nontherm_set() = secondaries.x12tot *
105 (sp->trans(i,0).Emis().gf()/
106 sp->trans(i,0).EnergyWN()) /
107 (iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,0).Emis().gf()/
108 iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,0).EnergyWN()) *
109 iso_ctrl.lgColl_excite[ipISO];;
110 }
111 else
112 sp->trans(i,0).Coll().rate_lu_nontherm_set() = 0.;
113 }
114
115 return;
116}
117
118/*============================*/
119/* evaluate collisional rates */
120void iso_collide( long ipISO, long nelem )
121{
122 DEBUG_ENTRY( "iso_collide()" );
123
124 /* this array stores the last temperature at which collision data were evaluated for
125 * each species of the isoelectronic sequence. */
126 static double TeUsed[NISO][LIMELM]={
127 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0},
128 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0} };
129
130 /* check that we were called with valid parameters */
131 ASSERT( ipISO < NISO );
132 ASSERT( nelem >= ipISO );
133 ASSERT( nelem < LIMELM );
134
135 t_iso_sp* sp = &iso_sp[ipISO][nelem];
136
137 /* skip most of this routine if temperature has not changed,
138 * the check on conv.nTotalIoniz is to make sure that we redo this
139 * on the very first call in a grid calc - it is 0 on the first call */
140 if( fp_equal( TeUsed[ipISO][nelem], phycon.te ) && conv.nTotalIoniz )
141 {
142 ASSERT( sp->trans( iso_ctrl.nLyaLevel[ipISO], 0 ).Coll().ColUL( colliders ) >= 0. );
143
144 if( trace.lgTrace && (trace.lgHBug||trace.lgHeBug) )
145 {
146 fprintf( ioQQQ,
147 " iso_collide called %s nelem %li - no reeval Boltz fac, LTE dens\n",
148 iso_ctrl.chISO[ipISO], nelem );
149 }
150 }
151 else
152 {
153 TeUsed[ipISO][nelem] = phycon.te;
154
155 if( trace.lgTrace && (trace.lgHBug||trace.lgHeBug) )
156 {
157 fprintf( ioQQQ,
158 " iso_collide called %s nelem %li - will reeval Boltz fac, LTE dens\n",
159 iso_ctrl.chISO[ipISO], nelem );
160 }
161
162 /**********************************************************
163 * *
164 * Boltzmann factors for all levels, and *
165 * collisional ionization and excitation *
166 * *
167 **********************************************************/
168
169 /* HION_LTE_POP is planck^2 / (2 pi m_e k ), raised to 3/2 next */
170 double factor = HION_LTE_POP*dense.AtomicWeight[nelem]/
171 (dense.AtomicWeight[nelem]+ELECTRON_MASS/ATOMIC_MASS_UNIT);
172
173 /* term in () is stat weight of electron * ion */
174 double ConvLTEPOP = pow(factor,1.5)/(2.*iso_ctrl.stat_ion[ipISO])/phycon.te32;
175
176 sp->lgPopLTE_OK = true;
177
178 // this is the maximum value of iso.PopLTE (units cm^3) that corresponds
179 // to the minimum positive density values. A smaller density will be
180 // regarded as zero, and the product PopLTE*n_e*n_Z+ will also be zero.
181#define MAX_POP_LTE (MAX_DENSITY/dense.density_low_limit/dense.density_low_limit)
182
183 /* fully define Boltzmann factors to continuum for model levels */
184 for( long ipLo=0; ipLo<sp->numLevels_max; ipLo++ )
185 {
186 /* this Boltzmann factor is exp( +ioniz energy / Te ) */
187 sp->st[ipLo].Boltzmann() =
188 dsexp(sp->st[ipLo].energy().Ryd()/phycon.te_ryd);
189 sp->st[ipLo].ConBoltz() =
190 dsexp(sp->fb[ipLo].xIsoLevNIonRyd/phycon.te_ryd);
191
192 /***************************************
193 * *
194 * LTE abundances for all levels *
195 * *
196 ***************************************/
197
198 if( sp->st[ipLo].ConBoltz() > SMALLDOUBLE )
199 {
200 /* LTE population of given level. */
201 sp->fb[ipLo].PopLTE =
202 sp->st[ipLo].g() / sp->st[ipLo].ConBoltz() * ConvLTEPOP;
203 ASSERT( sp->fb[ipLo].PopLTE < BIGDOUBLE );
204 }
205 else
206 {
207 sp->fb[ipLo].PopLTE = 0.;
208 }
209
210 sp->fb[ipLo].PopLTE = MIN2( sp->fb[ipLo].PopLTE, MAX_POP_LTE );
211
212 /* now check for any zeros - if present then matrix cannot be used */
213 if( sp->fb[ipLo].PopLTE <= 0. )
214 {
215 sp->lgPopLTE_OK = false;
216 }
217 }
218
219 iso_collisional_ionization( ipISO, nelem );
220
221 /***********************************************************
222 * *
223 * collisional deexcitation for all lines in iso sequence *
224 * *
225 ***********************************************************/
226
227 if( iso_ctrl.lgColl_excite[ipISO] )
228 {
229 for( long ipHi=1; ipHi<sp->numLevels_max; ipHi++ )
230 {
231 for( long ipLo=0; ipLo < ipHi; ipLo++ )
232 {
233 for( long ipCollider = ipELECTRON; ipCollider <= ipALPHA; ipCollider++ )
234 {
235 double cs_temp = 0.;
236
237 if( N_(ipHi) == N_(ipLo) && !iso_ctrl.lgColl_l_mixing[ipISO] )
238 cs_temp = 0.;
239 else if( N_(ipHi)-N_(ipLo) > 2 && ipCollider > ipELECTRON )
240 cs_temp = 0.;
241 else if( ipISO == ipH_LIKE )
242 cs_temp = HydroCSInterp( nelem , ipHi , ipLo, ipCollider );
243 else if( ipISO == ipHE_LIKE )
244 cs_temp = HeCSInterp( nelem , ipHi , ipLo, ipCollider );
245 else
247
248 cs_temp *= iso_ctrl.lgColl_excite[ipISO];
249
250 if( opac.lgCaseB_HummerStorey && N_(ipHi)==N_(ipLo)+1 && abs(L_(ipHi)-L_(ipLo))==1 )
251 {
252 double Aul = HydroEinstA( N_(ipLo), N_(ipHi) );
253 cs_temp *= ( sp->trans(ipHi,ipLo).Emis().gf() / sp->st[ipLo].g() ) /
254 ( GetGF(Aul, sp->trans(ipHi,ipLo).EnergyWN(), 2.*N_(ipHi)*N_(ipHi))/(2.*N_(ipLo)*N_(ipLo)) );
255 }
256
257 // store electron collision strength in generic collision strength
258 if( ipCollider == ipELECTRON )
259 sp->trans(ipHi,ipLo).Coll().col_str() = (realnum) cs_temp;
260
261 double reduced_mass_collider_system = dense.AtomicWeight[nelem]*ColliderMass[ipCollider]/
262 (dense.AtomicWeight[nelem]+ColliderMass[ipCollider])*ATOMIC_MASS_UNIT;
263
264 if( ipCollider == ipELECTRON )
265 reduced_mass_collider_system = ELECTRON_MASS;
266
267 double rateCoef = cs_temp *
268 pow(ELECTRON_MASS/reduced_mass_collider_system, 1.5) * COLL_CONST/
269 ( phycon.sqrte * (double)sp->st[ipHi].g() );
270
271 if( !opac.lgCaseB_HummerStorey )
272 {
273 if( ipISO == ipH_LIKE )
274 {
275 if( N_(ipHi) > sp->n_HighestResolved_max &&
276 N_(ipLo) <= sp->n_HighestResolved_max )
277 {
278 rateCoef *= (8./3.)*(log(double(N_(ipHi)))+2.);
279 }
280 }
281 else if( ipISO == ipHE_LIKE )
282 {
283 fixit();
284 /* This is intended to be a trick to get the correct collisional excitation from
285 * collapsed levels to resolved levels. Is it needed or does the stat weight used
286 * above handle it automatically? If it is needed, is this correct? */
287 if( N_(ipHi) > sp->n_HighestResolved_max &&
288 N_(ipLo) <= sp->n_HighestResolved_max )
289 {
290 rateCoef *= (2./3.)*(log(double(N_(ipHi)))+2.);
291 }
292 }
293 }
294
295 sp->trans(ipHi,ipLo).Coll().rate_coef_ul_set()[ipCollider] =
296 (realnum) rateCoef;
297 }
298
299 /* check for sanity */
300 ASSERT( sp->trans(ipHi,ipLo).Coll().rate_coef_ul()[ipELECTRON] >= 0. );
301
302 if( N_(ipHi) <= 5 && N_(ipLo) <= 2 )
303 iso_put_error( ipISO, nelem, ipHi, ipLo, IPCOLLIS, 0.10f, 0.30f );
304 else
305 iso_put_error( ipISO, nelem, ipHi, ipLo, IPCOLLIS, 0.20f, 0.30f );
306 }
307 }
308 }
309
310 if( (trace.lgTrace && trace.lgIsoTraceFull[ipISO]) && (nelem == trace.ipIsoTrace[ipISO]) )
311 {
312 fprintf( ioQQQ, " iso_collide: %s Z=%li de-excitation rates coefficients\n", iso_ctrl.chISO[ipISO], nelem + 1 );
313 long upper_limit = sp->numLevels_local;
314 for( long ipHi=1; ipHi < upper_limit; ipHi++ )
315 {
316 fprintf( ioQQQ, " %li\t", ipHi );
317 for( long ipLo=0; ipLo < ipHi; ipLo++ )
318 {
319 fprintf( ioQQQ,PrintEfmt("%10.3e",
320 sp->trans(ipHi,ipLo).Coll().ColUL( colliders ) / dense.EdenHCorr ));
321 }
322 fprintf( ioQQQ, "\n" );
323 }
324
325 fprintf( ioQQQ, " iso_collide: %s Z=%li collisional ionization coefficients\n", iso_ctrl.chISO[ipISO], nelem + 1 );
326 for( long ipHi=0; ipHi < upper_limit; ipHi++ )
327 {
328 fprintf( ioQQQ,PrintEfmt("%10.3e", sp->fb[ipHi].ColIoniz ));
329 }
330 fprintf( ioQQQ, "\n" );
331
332 fprintf( ioQQQ, " iso_collide: %s Z=%li continuum boltzmann factor\n", iso_ctrl.chISO[ipISO], nelem + 1 );
333 for( long ipHi=0; ipHi < upper_limit; ipHi++ )
334 {
335 fprintf( ioQQQ,PrintEfmt("%10.3e", sp->st[ipHi].ConBoltz() ));
336 }
337 fprintf( ioQQQ, "\n" );
338
339 fprintf( ioQQQ, " iso_collide: %s Z=%li continuum boltzmann factor\n", iso_ctrl.chISO[ipISO], nelem + 1 );
340 for( long ipHi=0; ipHi < upper_limit; ipHi++ )
341 {
342 fprintf( ioQQQ,PrintEfmt("%10.3e", sp->fb[ipHi].PopLTE ));
343 }
344 fprintf( ioQQQ, "\n" );
345 }
346
347 /* the case b hummer and storey option,
348 * this kills collisional excitation and ionization from n=1 and n=2 */
349 if( opac.lgCaseB_HummerStorey )
350 {
351 for( long ipLo=0; ipLo<sp->numLevels_max-1; ipLo++ )
352 {
353 if( N_(ipLo)>=3 )
354 break;
355
356 sp->fb[ipLo].ColIoniz = 0.;
357
358 for( long ipHi=ipLo+1; ipHi<sp->numLevels_max; ipHi++ )
359 {
360 /* don't disable 2-2 collisions */
361 if( N_(ipLo)==2 && N_(ipHi)==2 )
362 continue;
363
364 sp->trans(ipHi,ipLo).Coll().col_str() = 0.;
365 for( long k=0; k<ipNCOLLIDER; ++k )
366 sp->trans(ipHi,ipLo).Coll().rate_coef_ul_set()[k] = 0.;
367 }
368 }
369 }
370 }
371
372 iso_suprathermal( ipISO, nelem );
373
374 /* this must be reevaluated every time since eden can change when Te does not */
375 /* save into main array - collisional ionization by thermal electrons */
376 ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0] =
377 sp->fb[0].ColIoniz*dense.EdenHCorr;
378
379 /* cooling due to collisional ionization, which only includes thermal electrons */
380 ionbal.CollIonRate_Ground[nelem][nelem-ipISO][1] =
381 ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0]*
382 rfield.anu[Heavy.ipHeavy[nelem][nelem-ipISO]-1]*EN1RYD;
383
384 return;
385}
FILE * ioQQQ
Definition cddefines.cpp:7
#define ASSERT(exp)
Definition cddefines.h:578
#define PrintEfmt(F, V)
Definition cddefines.h:1472
#define MIN2
Definition cddefines.h:761
const int LIMELM
Definition cddefines.h:258
const int NISO
Definition cddefines.h:261
float realnum
Definition cddefines.h:103
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
NORETURN void TotalInsanity(void)
Definition service.cpp:886
const int ipHYDROGEN
Definition cddefines.h:305
double dsexp(double x)
Definition service.cpp:953
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void fixit(void)
Definition service.cpp:991
const double * rate_coef_ul() const
Definition collision.h:176
realnum & col_str() const
Definition collision.h:167
realnum & rate_lu_nontherm_set() const
Definition collision.h:181
realnum ColUL(const ColliderList &colls) const
Definition collision.h:99
double * rate_coef_ul_set() const
Definition collision.h:172
static t_ADfA & Inst()
Definition cddefines.h:175
CollisionProxy Coll() const
Definition transition.h:424
realnum & EnergyWN() const
Definition transition.h:438
long & ipCont() const
Definition transition.h:450
double coll_ion_wrapper(long int z, long int n, double t)
long int numLevels_max
Definition iso.h:493
vector< freeBound > fb
Definition iso.h:452
TransitionProxy trans(const long ipHi, const long ipLo)
Definition iso.h:444
long int numLevels_local
Definition iso.h:498
bool lgPopLTE_OK
Definition iso.h:523
bool lgLevelsLowered
Definition iso.h:474
long int n_HighestResolved_max
Definition iso.h:505
qList st
Definition iso.h:453
ColliderList colliders
Definition collision.cpp:57
@ ipALPHA
Definition collision.h:12
@ ipNCOLLIDER
Definition collision.h:18
@ ipELECTRON
Definition collision.h:9
t_conv conv
Definition conv.cpp:5
const double BIGDOUBLE
Definition cpu.h:194
const double SMALLDOUBLE
Definition cpu.h:195
t_dense dense
Definition dense.cpp:24
t_Heavy Heavy
Definition heavy.cpp:5
realnum HeCSInterp(long int nelem, long int ipHi, long int ipLo, long int Collider)
static const double ColliderMass[4]
Definition helike_cs.cpp:87
double Hion_coll_ioniz_ratecoef(long int ipISO, long int nelem, long int n, double ionization_energy_Ryd, double Te)
double hydro_vs_ioniz(double ionization_energy_Ryd, double Te)
realnum HydroCSInterp(long int nelem, long int ipHi, long int ipLo, long int ipCollider)
double HydroEinstA(long int n1, long int n2)
t_ionbal ionbal
Definition ionbal.cpp:5
t_iso_sp iso_sp[NISO][LIMELM]
Definition iso.cpp:8
t_isoCTRL iso_ctrl
Definition iso.cpp:6
#define IPCOLLIS
Definition iso.h:87
#define N_(A_)
Definition iso.h:20
const int ipHE_LIKE
Definition iso.h:63
const int ipH2p
Definition iso.h:29
const int ipH_LIKE
Definition iso.h:62
#define L_(A_)
Definition iso.h:21
void iso_put_error(long ipISO, long nelem, long ipHi, long ipLo, long whichData, realnum errorOpt, realnum errorPess)
void iso_collide(long ipISO, long nelem)
void iso_collisional_ionization(long ipISO, long nelem)
#define MAX_POP_LTE
void iso_suprathermal(long ipISO, long nelem)
double GetGF(double trans_prob, double enercm, double gup)
t_opac opac
Definition opacity.cpp:5
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double ELECTRON_MASS
Definition physconst.h:91
UNUSED const double EN1RYD
Definition physconst.h:179
UNUSED const double PROTON_MASS
Definition physconst.h:94
UNUSED const double ATOMIC_MASS_UNIT
Definition physconst.h:88
UNUSED const double COLL_CONST
Definition physconst.h:229
UNUSED const double HION_LTE_POP
Definition physconst.h:157
t_rfield rfield
Definition rfield.cpp:8
t_secondaries secondaries
t_trace trace
Definition trace.cpp:5