cloudy trunk
Loading...
Searching...
No Matches
mole_dissociate.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 "physconst.h"
5#include "h2.h"
6#include "h2_priv.h"
7#include "hmi.h"
8#include "rfield.h"
9#include "thirdparty.h"
10#include "ipoint.h"
11#include "mole.h"
12
14
16{
17 DEBUG_ENTRY( "diatomics::Read_Mol_Diss_cross_sections()" );
18
19 FILE *ioDATA;
20
21 char chPath[FILENAME_PATH_LENGTH_2],
22 chDirectory[FILENAME_PATH_LENGTH_2],
24
25 /* these can be generalized later for other molecules */
26 const int ipNUM_FILES = 1;
27
28 char chFileNames[ipNUM_FILES][17] =
29 {"cont_diss.dat"} ;
30
31 // these cross sections are not meaningful without big H2 enabled.
33
34 // for the summed-over-final-state rates, allocate space for all states in model.
35 // Will be zero-filled and overwritten with whatever data is available.
37 for( int iElecLo=0; iElecLo<n_elec_states; ++iElecLo )
38 {
39 Cont_Dissoc_Rate.reserve( iElecLo, nVib_hi[iElecLo]+1 );
40 for( int iVibLo=0; iVibLo<=nVib_hi[iElecLo]; ++iVibLo )
41 {
42 Cont_Dissoc_Rate.reserve( iElecLo, iVibLo, nRot_hi[iElecLo][iVibLo]+1 );
43 }
44 }
45 Cont_Dissoc_Rate.alloc();
46
47 /* drill down to the directory where H2_cont_diss.dat lives */
48# ifdef _WIN32
49 strcpy( chDirectory, "h2\\" );
50# else
51 strcpy( chDirectory, "h2/" );
52# endif
53
54 /* now open the data file */
55 strcpy( chPath, chDirectory );
56 strcat( chPath, chFileNames[0] );
57 ioDATA = open_data( chPath, "r" );
58
59 Diss_Trans.clear();
60
61 /* track number of datasets in file */
62 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
63 {
64 static bool skipData = false;
65 long ini=0, inf, iv, ij;
66
67 /* ignore all lines that begin with a '#', unless they are telling us the next quantum numbers. If so,
68 * then read the quantum numbers and start a new dataset */
69 if(sizeof(chLine) >= 2 && chLine[0] == '#' && chLine[1] == '!' && chLine[4] == 'n' && chLine[5] =='e' )
70 {
71 sscanf(chLine,"#! nei=%li, nef=%li, vi= %li, ji= %li",&ini, &inf, &iv, &ij);
72 // skip data that is for levels beyond our model
73 if( ini > n_elec_states )
74 skipData = true;
75 else if( iv > nVib_hi[ini] )
76 skipData = true;
77 else if( ij > nRot_hi[ini][iv] )
78 skipData = true;
79 else
80 {
81 skipData = false;
82 diss_level Leveli = {ini, iv, ij};
83 diss_level Levelf = {inf, -1, -1};
84 diss_tran thisTran( Leveli, Levelf );
85 Diss_Trans.push_back( thisTran );
86 }
87 }
88
89 if( skipData )
90 continue;
91
92 /* if line does not begin with a '#', then they are numbers in the datset. Read in energy in
93 * wavenumbers and cross section in Angstroms ^2 */
94 if(chLine[0] != '#')
95 {
96 double energy, xsection;
97 const double AngSquared = 1e-16;
98 sscanf(chLine,"%lf,%lf",&energy, &xsection);
99
100 /* convert energy in wavenumbers to Ryd */
101 Diss_Trans.back().energies.push_back( energy*WAVNRYD );
102
103 /* convert cross section from Angstrom^2 to cm^2 by multiplying by 1e-16 */
104 Diss_Trans.back().xsections.push_back( xsection*AngSquared );
105 }
106 }
107
108 fclose(ioDATA);
109 return;
110}
111
112double diatomics::MolDissocOpacity( const diss_tran& tran, const double& Mol_Ene )
113{
114 DEBUG_ENTRY( "diatomics::MolDissocOpacity()" );
115
116 double cross_section = MolDissocCrossSection( tran, Mol_Ene ) *
117 states[ ipEnergySort[ tran.initial.n ][ tran.initial.v ][ tran.initial.j ] ].Pop();
118 return cross_section;
119}
120
121double MolDissocCrossSection( const diss_tran& tran, const double& Mol_Ene )
122{
123 DEBUG_ENTRY( "MolDissocCrossSection()" );
124
125 double photodiss_cs = 0.;
126
127 /* if energy is less than threshold, then set cross section is zero */
128 if (Mol_Ene < tran.energies[0])
129 photodiss_cs = 0.;
130 /* if energy is greater than the highest energy defined in Philip's data, set
131 * cross section to zero for now. This may need to be changed later */
132 else if(Mol_Ene > tran.energies.back())
133 photodiss_cs = tran.xsections.back()/sqrt(powi(Mol_Ene/tran.energies.back(),7));
134 /* If energy is in between high and low energy limits, do linear interpolation
135 * on Philip's data to compute cross section */
136 else
137 {
138 ASSERT( Mol_Ene > tran.energies[0] && Mol_Ene < tran.energies.back() );
139 photodiss_cs = linint(&tran.energies[0],
140 &tran.xsections[0],
141 tran.xsections.size(),
142 Mol_Ene);
143 }
144
145 return photodiss_cs;
146}
147
149{
150 DEBUG_ENTRY( "diatomics::Mol_Photo_Diss_Rates()" );
151
152 /* Start Calculation of H2 continuum dissociation rate here */
153
154 // these cross sections are not meaningful without big H2 enabled.
155 ASSERT( lgEnabled && mole_global.lgStancil );
156
157 /* Zero out dissociation rates */
158 Cont_Dissoc_Rate.zero();
161
162 for_each( Diss_Trans.begin(), Diss_Trans.end(), GetDissociationRateCoeff );
163 for( vector< diss_tran >::const_iterator dt = Diss_Trans.begin(); dt != Diss_Trans.end(); ++dt )
164 {
165 double rate = (*this).GetDissociationRate( *dt );
166 // this is summed over final electron state, tran.rate_coeff is not (because it's used to conserve energy)
167 Cont_Dissoc_Rate[dt->initial.n][dt->initial.v][dt->initial.j] += dt->rate_coeff;
168 if( states[ ipEnergySort[dt->initial.n][dt->initial.v][dt->initial.j] ].energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
169 {
170 /* Add up total rate due to levels which are above arbitrarily defined
171 * H2* threshold. Call this dissociation rate of H2* */
172 Cont_Dissoc_Rate_H2s += rate;
173 }
174 else
175 {
176 /* Add up total rate due to levels which are below arbitrarily defined
177 * H2* threshold. Call this dissociation rate of H2g (H2 ground) */
178 Cont_Dissoc_Rate_H2g += rate;
179 }
180 }
181
182 /* This has current units of cm-3 s-1. We want just the total rate coefficient for
183 * insertion into h_step. Therefore, divide rate by density of H2* and H2g to get rate
184 * coefficient for reaction */
187
188 return;
189}
190
192{
193 DEBUG_ENTRY( "GetDissociationRateCoeff()" );
194
195 long index_min = ipoint( tran.energies[0] );
196 long index_max = ipoint( tran.energies.back() );
197 index_max = MIN2( index_max, rfield.nflux-1 );
198
199 tran.rate_coeff = 0.;
200 for(long i = index_min; i <= index_max; ++i)
201 {
202 /* Find cross section at given n, v, j, and photon energy */
203 double photodiss_cs = MolDissocCrossSection( tran, rfield.anu[i] );
204
205 /* Compute integral of cross section and photon energy -- rate coefficient */
206 tran.rate_coeff += photodiss_cs*( rfield.flux[0][i] + rfield.ConInterOut[i]+
207 rfield.outlin[0][i]+ rfield.outlin_noplot[i] );
208 }
209}
210
212{
213 DEBUG_ENTRY( "diatomics::GetDissociationRate()" );
214
215 long iElecLo = tran.initial.n;
216 long iVibLo = tran.initial.v;
217 long iRotLo = tran.initial.j;
218
219 /* Compute rate, product of rate coefficient times density */
220 return states[ ipEnergySort[iElecLo][iVibLo][iRotLo] ].Pop() * tran.rate_coeff;
221}
222
224{
225 DEBUG_ENTRY( "diatomics::Cont_Diss_Heat_Rate()" );
226
227 /* 29 Jul 10 - NPA - this tells code to not compute detailed H2 dissociation heating rate
228 * if Stancil rate is not used. However, we also do not want to compute this rate if Stancil
229 * rate is on, but small H2 molecule is used. Insert logic hear to not use rate if
230 * small H2 atom is on */
231
232 if( !mole_global.lgStancil || !lgEnabled )
233 return 0.;
234
236
237 double Cont_Dissoc_Heat_Rate = 0.0;
238 for( vector< diss_tran >::const_iterator dt = Diss_Trans.begin(); dt != Diss_Trans.end(); ++dt )
239 Cont_Dissoc_Heat_Rate += (*this).GetHeatRate( *dt );
240
241 return Cont_Dissoc_Heat_Rate;
242}
243
245{
246 DEBUG_ENTRY( "diatomics::GetHeatRate()" );
247
248 long index_min = ipoint( tran.energies[0] );
249 long index_max = ipoint( tran.energies.back() );
250 index_max = MIN2( index_max, rfield.nflux-1 );
251
252 long iElecLo = tran.initial.n;
253 long iVibLo = tran.initial.v;
254 long iRotLo = tran.initial.j;
255 double rate = 0.;
256
257 for( long i = index_min; i<= index_max; ++i )
258 {
259 /* Photon energy minus threshold (Cloudy convention) */
260 double energy = MAX2( 0., rfield.anu[i] - tran.energies[0] );
261
262 /* The density of the current H2(v,J) level */
263 double density = states[ ipEnergySort[iElecLo][iVibLo][iRotLo] ].Pop();
264
265 double photodiss_cs = MolDissocCrossSection( tran, rfield.anu[i] );
266 double Rate_Coeff = photodiss_cs*( rfield.flux[0][i] + rfield.ConInterOut[i]+
267 rfield.outlin[0][i]+ rfield.outlin_noplot[i] );
268
269 /* The heating rate. This is the product of the H2(v,J) density, the dissociation rate,
270 * and the energy of the photon (minus threshold). Units are erg cm-3 s-1.
271 * Sum over all possible H2 levels and dissociation energies. */
272 rate += EN1RYD * energy * Rate_Coeff * density;
273 }
274
275 return rate;
276}
277
const int FILENAME_PATH_LENGTH_2
Definition cddefines.h:249
#define ASSERT(exp)
Definition cddefines.h:578
#define STATIC
Definition cddefines.h:97
#define MIN2
Definition cddefines.h:761
#define MAX2
Definition cddefines.h:782
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition service.cpp:70
sys_float SDIV(sys_float x)
Definition cddefines.h:952
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
double powi(double, long int)
Definition service.cpp:604
long int n_elec_states
Definition h2_priv.h:409
double Cont_Dissoc_Rate_H2g
Definition h2_priv.h:278
void Mol_Photo_Diss_Rates(void)
long int nVib_hi[N_ELEC]
Definition h2_priv.h:611
bool lgEnabled
Definition h2_priv.h:345
double Cont_Diss_Heat_Rate(void)
qList states
Definition h2_priv.h:565
double GetDissociationRate(const diss_tran &tran)
double MolDissocOpacity(const diss_tran &tran, const double &Mol_Ene)
double H2_den_s
Definition h2_priv.h:682
vector< diss_tran > Diss_Trans
Definition h2_priv.h:568
const double ENERGY_H2_STAR
Definition h2_priv.h:585
void Read_Mol_Diss_cross_sections(void)
double GetHeatRate(const diss_tran &tran)
double Cont_Dissoc_Rate_H2s
Definition h2_priv.h:277
double H2_den_g
Definition h2_priv.h:682
valarray< long > nRot_hi[N_ELEC]
Definition h2_priv.h:613
multi_arr< long int, 3 > ipEnergySort
Definition h2_priv.h:690
multi_arr< double, 3 > Cont_Dissoc_Rate
Definition h2_priv.h:279
vector< double > xsections
Definition h2_priv.h:50
diss_level initial
Definition h2_priv.h:48
double rate_coeff
Definition h2_priv.h:51
vector< double > energies
Definition h2_priv.h:49
long ipoint(double energy_ryd)
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
double MolDissocCrossSection(const diss_tran &tran, const double &Mol_Ene)
STATIC double cross_section(double EgammaRyd, double EthRyd, long nelem, long n, long l, long s)
t_hmi hmi
Definition hmi.cpp:5
t_mole_global mole_global
Definition mole.cpp:6
STATIC void GetDissociationRateCoeff(diss_tran &tran)
double MolDissocCrossSection(const diss_tran &tran, const double &Mol_Ene)
UNUSED const double EN1RYD
Definition physconst.h:179
UNUSED const double WAVNRYD
Definition physconst.h:173
t_rfield rfield
Definition rfield.cpp:8
long v
Definition h2_priv.h:34
long j
Definition h2_priv.h:34
long n
Definition h2_priv.h:34
double linint(const double x[], const double y[], long n, double xval)