cloudy trunk
Loading...
Searching...
No Matches
mole_h2_etc.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/* mole_H2_LTE sets Boltzmann factors and LTE unit population of large H2 molecular */
4/* H2_zero_pops_too_low - zero out some H2 variables if we decide not to compute
5 * the full sim, called by H2_LevelPops*/
6/* H2_Solomon_rate find rates between H2s and H2g and other levels,
7 * for eventual use in the chemistry */
8/* gs_rate evaluate rates between ground and star states */
9/* H2_He_coll_init initialize H2 - He collision data set
10 * H2_He_coll interpolate on h2 - He collision data set to return rate at temp*/
11#include "cddefines.h"
12#include "phycon.h"
13#include "hmi.h"
14#include "taulines.h"
15#include "h2.h"
16#include "h2_priv.h"
17#include "rfield.h"
18#include "thermal.h"
19#include "gammas.h"
20#include "ionbal.h"
21
22/*H2_Solomon_rate find rates between H2s and H2g and other levels,
23 * for eventual use in the chemistry */
25{
26 DEBUG_ENTRY( "H2_Solomon_rate()" );
27
28 /* iElecLo will always be X in this routine */
29
30 /* find rate (s-1) h2 dissociation into X continuum by Solomon process and
31 * assign to the TH85 g and s states
32 * these will go back into the chemistry network */
33
34 /* rates [s-1] for dissociation from s or g, into electronic excited states
35 * followed by dissociation */
38
39 /* these are used in a print statement - are they needed? */
42
43 /* at this point we have already evaluated the sum of the radiative rates out
44 * of the electronic excited states - this is H2_rad_rate_out[electronic][vib][rot]
45 * and this includes both decays into the continuum and bound states of X */
46
47 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
48 {
49 qList::iterator Hi = (*tr).Hi() ;
50 if( (*Hi).n() < 1 ) continue;
51 double factor = (double)H2_dissprob[(*Hi).n()][(*Hi).v()][(*Hi).J()]/H2_rad_rate_out[(*Hi).n()][(*Hi).v()][(*Hi).J()];
52 /* this is the rate [cm-3 s-1] that mole goes from
53 * lower level into electronic excited states then
54 * into continuum */
55 double rate_up_cont = (*(*tr).Lo()).Pop() * (*tr).Emis().pump() * factor;
56
57 /* rate electronic state decays into H2g */
58 double elec_decay = (*(*tr).Hi()).Pop() * (*tr).Emis().Aul() *
59 ((*tr).Emis().Pesc() + (*tr).Emis().Pdest() + (*tr).Emis().Pelec_esc());
60 if( (*(*tr).Lo()).energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
61 {
62 /* this is H2g up to excited then to continuum -
63 * cm-3 s-1 at this point */
64 Solomon_dissoc_rate_s += rate_up_cont;
65 /* rate electronic state decays into H2g */
66 Solomon_elec_decay_s += elec_decay;
67 }
68 else
69 {
70 /* this is H2g up to excited then to continuum -
71 * cm-3 s-1 at this point */
72 Solomon_dissoc_rate_g += rate_up_cont;
73 /* rate electronic state decays into H2g */
74 Solomon_elec_decay_g += elec_decay;
75 }
76 }
77
78 double H2_sum_excit_elec_den = GetExcitedElecDensity();
79
80 /* at this point units of Solomon_elec_decay_g, H2s are cm-3 s-1
81 * since populations are included -
82 * div by pops to get actual dissocation rate, s-1 */
83 if( *dense_total > SMALLFLOAT )
84 {
85 Solomon_elec_decay_g /= SDIV( H2_sum_excit_elec_den );
86 Solomon_elec_decay_s /= SDIV( H2_sum_excit_elec_den );
87
88 /* will be used for H2s-> H + H */
90
91 /* will be used for H2g-> H + H */
93
94 }
95 else
96 {
99
100 }
101 /*fprintf(ioQQQ,"DEBUG H2 new %.2e %.2e %.2e %.2e \n",
102 Solomon_elec_decay_g ,
103 Solomon_elec_decay_s ,
104 Solomon_dissoc_rate_s,
105 Solomon_dissoc_rate_g );*/
106
107 return;
108}
109
110/* gs_rate evaluate rate between ground and star states */
111double diatomics::gs_rate( void )
112{
113 DEBUG_ENTRY( "diatomics::gs_rate()" );
114
115 /* rate g goes to s */
116 double ground_to_star_rate = 0.;
117
118 /* loop over all levels in H2g */
119 for( long ipLoX=0; ipLoX < nEner_H2_ground; ++ipLoX )
120 {
121 /* now find all pumps up to electronic excited states */
122 /* sum over all electronic states, finding dissociation rate */
123 for( long iElecHi=1; iElecHi<n_elec_states; ++iElecHi )
124 {
125 for( long iVibHi=0; iVibHi<=nVib_hi[iElecHi]; ++iVibHi )
126 {
127 for( long iRotHi=Jlowest[iElecHi]; iRotHi<=nRot_hi[iElecHi][iVibHi]; ++iRotHi )
128 {
129 long ipHi = ipEnergySort[iElecHi][iVibHi][iRotHi];
130 if( lgH2_radiative[ipHi][ipLoX] )
131 {
132 /* this is the rate [cm-3 s-1] that mole goes from
133 * lower level into electronic excited states then
134 * into continuum */
135 double rate_up_cont =
136 states[ipLoX].Pop() *
137 trans[ ipTransitionSort[ipHi][ipLoX] ].Emis().pump();
138 double decay_star = H2_rad_rate_out[iElecHi][iVibHi][iRotHi] - H2_dissprob[iElecHi][iVibHi][iRotHi];
139 /* loop over all other levels in H2g, subtracting
140 * their rate - remainder is rate into star, this is
141 * usually only a few levels */
142 for( long ipOther=0; ipOther < nEner_H2_ground; ++ipOther )
143 {
144 if( lgH2_radiative[ipHi][ipOther] )
145 {
146 EmissionProxy em = trans[ ipTransitionSort[ipHi][ipOther] ].Emis();
147 decay_star -= em.Aul() * ( em.Pesc() + em.Pdest() + em.Pelec_esc() );
148 }
149 }
150 /* MAX because may underflow to negative numbers is rates very large
151 * this is fraction that returns to H2s */
152 decay_star = MAX2(0., decay_star)/SDIV(H2_rad_rate_out[iElecHi][iVibHi][iRotHi]);
153 ground_to_star_rate += rate_up_cont*decay_star;
154
155 }/* end if line exists */
156 }/* end loop rot electronic excited */
157 }/* end loop vib electronic excited */
158 }/* end loop electronic electronic excited */
159 }
160
161 /* at this point units are cm-3 s-1 - convert to rate s-1 */
162 ground_to_star_rate /= SDIV( H2_den_g );
163 return ground_to_star_rate;
164}
165
166long diatomics::OpacityCreate( double *stack )
167{
168 DEBUG_ENTRY( "diatomics::OpacityCreate()" );
169
170 ASSERT( photoion_opacity_fun != NULL );
171
172 for( long i=ip_photo_opac_thresh-1; i < rfield.nupper; i++ )
173 {
174 /* >>chng 05 nov 24, add H2 photoionization cross section */
175 //opac.OpacStack[i-ip_photo_opac_thresh + ip_photo_opac_offset] =
177 photoion_opacity_fun( rfield.AnuOrg[i] );
178 //Yan_H2_CS(rfield.AnuOrg[i]);
179 }
180
181 //opac.nOpacTot += rfield.nupper - ip_photo_opac_thresh + 1;
182 return rfield.nupper - ip_photo_opac_thresh + 1;
183}
184
185/* H2_zero_pops_too_low - zero out some H2 variables if we decide not to compute
186 * the full sim, called by H2_LevelPops*/
188{
189 DEBUG_ENTRY( "H2_zero_pops_too_low()" );
190
191 // zero populations
192 fill_n( pops_per_elec, N_ELEC, 0. );
193 pops_per_vib.zero();
194 for( qList::iterator st = states.begin(); st != states.end(); ++st )
195 {
196 double pop = H2_populations_LTE[ (*st).n() ][ (*st).v() ][ (*st).J() ] * (*dense_total);
197 H2_old_populations[ (*st).n() ][ (*st).v() ][ (*st).J() ] = pop;
198 (*st).Pop() = pop;
199 }
200
201
202 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
203 {
204 /* population of lower level with correction for stim emission */
205 (*tr).Emis().PopOpc() = (*(*tr).Lo()).Pop() - (*(*tr).Hi()).Pop() * (*(*tr).Lo()).g() / (*(*tr).Hi()).g();
206
207 /* following two heat exchange excitation, deexcitation */
208 (*tr).Coll().cool() = 0.;
209 (*tr).Coll().heat() = 0.;
210
211 /* intensity of line */
212 (*tr).Emis().xIntensity() = 0.;
213 (*tr).Emis().phots() = 0.;
214 (*tr).Emis().ots() = 0.;
215 }
216
219 HeatDiss = 0.;
220 HeatDexc = 0.;
221 HeatDexc_deriv = 0.;
224 return;
225}
226
227/*mole_H2_LTE sets Boltzmann factors and LTE unit population of large H2 molecular */
229{
230 DEBUG_ENTRY( "mole_H2_LTE()" );
231
232 /* do we need to update the Boltzmann factors and unit LTE populations? */
233 if( ! fp_equal( phycon.te, TeUsedBoltz ) )
234 {
235 double part_fun = 0.;
236 TeUsedBoltz = phycon.te;
237 /* loop over all levels setting H2_Boltzmann and deriving partition function */
238 for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
239 {
240 long iElec = (*st).n();
241 long iVib = (*st).v();
242 long iRot = (*st).J();
243 H2_Boltzmann[iElec][iVib][iRot] =
244 /* energy is relative to lowest level in the molecule, v=0, J=0,
245 * so Boltzmann factor is relative to this level */
246 dsexp( (*st).energy().K() / phycon.te );
247 /* sum the partition function - Boltzmann factor times statistical weight */
248 part_fun += H2_Boltzmann[iElec][iVib][iRot] * (*st).g();
249 ASSERT( part_fun > 0 );
250 }
251 /* have partition function, set H2_populations_LTE (populations for unit H2 density) */
252 for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
253 {
254 long iElec = (*st).n();
255 long iVib = (*st).v();
256 long iRot = (*st).J();
257 /* these are the H2 LTE populations for a unit H2 density -
258 * these populations will sum up to unity */
259 H2_populations_LTE[iElec][iVib][iRot] =
260 H2_Boltzmann[iElec][iVib][iRot] *
261 (*st).g() / part_fun;
262 }
263 if( nTRACE >= n_trace_full )
264 fprintf(ioQQQ,
265 "mole_H2_LTE set H2_Boltzmann factors, T=%.2f, partition function is %.2f\n",
266 phycon.te,
267 part_fun);
268 }
269
270 return;
271}
272
273/*H2_Reset called by IterRestart to reset variables that are needed after an iteration */
275{
276
277 DEBUG_ENTRY( "H2_Reset()" );
278
279 if(nTRACE)
280 fprintf(ioQQQ,
281 "\n***************H2_Reset called, resetting nCall_this_iteration, zone %.2f iteration %li\n",
282 fnzone,
283 iteration );
284
285 /* number of times large molecules evaluated in this iteration,
286 * is false if never evaluated, on next evaluation will start with LTE populations */
288 nCall_this_zone = 0;
289
290 /* these remember the largest and smallest factors needed to
291 * renormalize the H2 chemistry */
292 renorm_max = 1.;
293 renorm_min = 1.;
294
295 /* counters used by H2_itrzn to find number of calls of h2 per zone */
296 nH2_pops = 0;
297 nH2_zone = 0;
298 /* this is used to establish zone number for evaluation of number of levels in matrix */
300
301 nzoneAsEval = -1;
302 iterationAsEval = -1;
303
304 TeUsedColl = -1;
305 TeUsedBoltz = -1;
306
307 lgEvaluated = false;
308
309 /* zero out array used to save emission line intensities */
310 H2_SaveLine.zero();
311
312 /* this was option to print all electronic states in the main printout - but
313 * number of electronic states was not known at initialization so set to -1,
314 * will set properly now */
315 if( nElecLevelOutput < 1 )
317
318 return;
319}
320
321/*Yan_H2_CS - cross sections for the photoionization of H2
322 * may be represented analytically in the paper
323 *>>refer H2 photo cs Yan, M.; Sadeghpour, H. R.; Dalgarno, A. 1998, ApJ, 496, 1044
324 * This is revised version given in ERRATUM 2001, ApJ, 559, 1194
325 * return value is photo cs in cm-2 */
326double Yan_H2_CS( double energy_ryd /* photon energy in ryd */)
327{
328
329 double energy_keV; /* keV */
330 double cross_section; /* barns */
331 double x; /* x = E/15.4 */
332 double xsqrt , x15 , x2;
333 double energy = energy_ryd * EVRYD;
334
335 DEBUG_ENTRY( "Yan_H2_CS()" );
336
337 /* energy relative to threshold */
338 x = energy / 15.4;
339 energy_keV = energy/1000.0;
340 xsqrt = sqrt(x);
341 x15 = x * xsqrt;
342 x2 = x*x;
343
344 if( energy < 15.4 )
345 {
346 cross_section = 0.;
347 }
348
349 else if(energy >= 15.4 && energy < 18.)
350 {
351 cross_section = 1e7 * (1 - 197.448/xsqrt + 438.823/x - 260.481/x15 + 17.915/x2);
352 /* this equation is obviously negative for x = 1 */
354 }
355
356 else if(energy >= 18. && energy <= 30.)
357 {
358 cross_section = (-145.528 +351.394*xsqrt - 274.294*x + 74.320*x15)/pow(energy_keV,3.5);
359 }
360
361 else if(energy > 30. && energy <= 85.)
362 {
363 cross_section = (65.304 - 91.762*xsqrt + 51.778*x - 9.364*x15)/pow(energy_keV,3.5);
364 }
365
366 /* the high-energy tail */
367 else
368 {
369 cross_section = 45.57*(1 - 2.003/xsqrt - 4.806/x + 50.577/x15 - 171.044/x2 + 231.608/(xsqrt*x2) - 81.885/(x*x2))/pow(energy_keV,3.5);
370 }
371
372 return( cross_section * 1e-24 );
373}
374
376{
377 DEBUG_ENTRY( "diatomics::CalcPhotoionizationRate()" );
378
379 /* must reevaluate During search phase */
381 {
382 t_phoHeat photoHeat;
383 /* generally not important, do one time per zone */
386 rfield.nupper,
388 &photoHeat)*
389 ionbal.lgPhotoIoniz_On +
390 /* Compton recoil ionization - we include this in the H2 photoionization
391 * rate but not the heating rate - factor of two since assume 2H
392 * is same as two H0 at such high energies */
393 2.*ionbal.CompRecoilIonRate[ipHYDROGEN][0];
394
395 /* photo heating - this has units s-1 - needs H2 density
396 * to become vol heat rate */
397 photo_heat_soft = photoHeat.HeatLowEnr * ionbal.lgPhotoIoniz_On;
398 photo_heat_hard = photoHeat.HeatHiEnr * ionbal.lgPhotoIoniz_On;
401 }
402
403 return;
404}
405
static double x2[63]
long int nzone
Definition cddefines.cpp:14
FILE * ioQQQ
Definition cddefines.cpp:7
long int iteration
Definition cddefines.cpp:16
double fnzone
Definition cddefines.cpp:15
#define ASSERT(exp)
Definition cddefines.h:578
#define MAX2
Definition cddefines.h:782
bool fp_equal(sys_float x, sys_float y, int n=3)
Definition cddefines.h:812
sys_float SDIV(sys_float x)
Definition cddefines.h:952
const int ipHYDROGEN
Definition cddefines.h:305
double dsexp(double x)
Definition service.cpp:953
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
realnum & Pesc() const
Definition emission.h:523
realnum & Pelec_esc() const
Definition emission.h:533
realnum & Aul() const
Definition emission.h:613
realnum & Pdest() const
Definition emission.h:543
TransitionProxy::iterator iterator
Definition transition.h:280
long int iterationAsEval
Definition h2_priv.h:646
long int n_elec_states
Definition h2_priv.h:409
long int nzone_eval
Definition h2_priv.h:706
double Solomon_elec_decay_s
Definition h2_priv.h:269
double photodissoc_BigH2_H2s
Definition h2_priv.h:257
long int nVib_hi[N_ELEC]
Definition h2_priv.h:611
double pops_per_elec[N_ELEC]
Definition h2_priv.h:620
long int nzoneAsEval
Definition h2_priv.h:646
long ip_photo_opac_offset
Definition h2_priv.h:314
double Solomon_elec_decay_g
Definition h2_priv.h:268
multi_arr< realnum, 3 > H2_dissprob
Definition h2_priv.h:632
double TeUsedColl
Definition h2_priv.h:416
void H2_Solomon_rate(void)
double photo_heat_soft
Definition h2_priv.h:255
qList states
Definition h2_priv.h:565
long int Jlowest[N_ELEC]
Definition h2_priv.h:616
double TeUsedBoltz
Definition h2_priv.h:415
multi_arr< long int, 2 > ipTransitionSort
Definition h2_priv.h:691
int nElecLevelOutput
Definition h2_priv.h:349
long int nH2_zone
Definition h2_priv.h:718
double(* photoion_opacity_fun)(double energy)
Definition h2_priv.h:74
double HeatDiss
Definition h2_priv.h:289
const double *const dense_total
Definition h2_priv.h:589
void CalcPhotoionizationRate(void)
double HeatDexc_deriv
Definition h2_priv.h:292
double Solomon_dissoc_rate_g
Definition h2_priv.h:264
multi_arr< realnum, 6 > H2_SaveLine
Definition h2_priv.h:710
bool lgEvaluated
Definition h2_priv.h:310
multi_arr< double, 3 > H2_populations_LTE
Definition h2_priv.h:639
double gs_rate(void)
double H2_den_s
Definition h2_priv.h:682
long int iteration_evaluated
Definition h2_priv.h:707
double photodissoc_BigH2_H2g
Definition h2_priv.h:258
long OpacityCreate(double *stack)
double photoionize_rate
Definition h2_priv.h:254
double GetExcitedElecDensity(void)
Definition mole_h2.cpp:2538
multi_arr< double, 3 > H2_Boltzmann
Definition h2_priv.h:638
const double ENERGY_H2_STAR
Definition h2_priv.h:585
void H2_zero_pops_too_low(void)
double renorm_max
Definition h2_priv.h:336
double photo_heat_hard
Definition h2_priv.h:256
int n_trace_full
Definition h2_priv.h:404
void mole_H2_LTE(void)
long ip_photo_opac_thresh
Definition h2_priv.h:313
long int nCall_this_iteration
Definition h2_priv.h:726
multi_arr< bool, 2 > lgH2_radiative
Definition h2_priv.h:714
multi_arr< double, 2 > pops_per_vib
Definition h2_priv.h:600
double HeatDexc
Definition h2_priv.h:290
double H2_den_g
Definition h2_priv.h:682
TransitionList trans
Definition h2_priv.h:566
double Solomon_dissoc_rate_s
Definition h2_priv.h:265
long int nCall_this_zone
Definition h2_priv.h:341
valarray< long > nRot_hi[N_ELEC]
Definition h2_priv.h:613
double renorm_min
Definition h2_priv.h:337
multi_arr< double, 3 > H2_old_populations
Definition h2_priv.h:637
long int nH2_pops
Definition h2_priv.h:717
void H2_Reset(void)
int nTRACE
Definition h2_priv.h:399
multi_arr< double, 3 > H2_rad_rate_out
Definition h2_priv.h:634
TransitionList::iterator rad_end
Definition h2_priv.h:567
long int nzone_nlevel_set
Definition h2_priv.h:721
multi_arr< long int, 3 > ipEnergySort
Definition h2_priv.h:690
long int nEner_H2_ground
Definition h2_priv.h:597
ProxyIterator< qStateConstProxy, qStateConstProxy > const_iterator
ProxyIterator< qStateProxy, qStateConstProxy > iterator
double GammaK(long int ipLoEnr, long int ipHiEnr, long int ipOpac, double yield1, t_phoHeat *photoHeat)
const realnum SMALLFLOAT
Definition cpu.h:191
const int N_ELEC
Definition h2_priv.h:27
STATIC double cross_section(double EgammaRyd, double EthRyd, long nelem, long n, long l, long s)
t_hmi hmi
Definition hmi.cpp:5
t_ionbal ionbal
Definition ionbal.cpp:5
double Yan_H2_CS(double energy_ryd)
t_phycon phycon
Definition phycon.cpp:6
UNUSED const double EVRYD
Definition physconst.h:189
t_rfield rfield
Definition rfield.cpp:8
double HeatHiEnr
Definition thermal.h:176
double HeatLowEnr
Definition thermal.h:174