cloudy trunk
Loading...
Searching...
No Matches
mole_h2_coll.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/*H2_CollidRateRead read collision rates */
4/*H2_CollidRateEvalAll - set H2 collision rates */
5#include "cddefines.h"
6#include "atmdat.h"
7#include "phycon.h"
8#include "dense.h"
9#include "taulines.h"
10#include "input.h"
11#include "h2.h"
12#include "h2_priv.h"
13
14STATIC realnum GbarRateCoeff( long nColl, double ediff );
15
16/*H2_CollidRateEvalAll - set H2 collision rate coefficients */
18{
19 DEBUG_ENTRY( "H2_CollidRateEvalAll()" );
20
21 long int numb_coll_trans = 0;
22
23 if( nTRACE >= n_trace_full )
24 fprintf(ioQQQ,"%s set collision rates\n", label.c_str());
25 /* loop over all possible collisional changes within X
26 * and set collision rates, which only depend on Te
27 * will go through array in energy order since coll trans do not
28 * correspond to a line
29 * collisional dissociation rate coefficient, units cm3 s-1 */
30 H2_coll_dissoc_rate_coef[0][0] = 0.;
32 for( long ipHi=0; ipHi<nLevels_per_elec[0]; ++ipHi )
33 {
34 /* obtain the proper indices for the upper level */
35 long iVibHi = ipVib_H2_energy_sort[ipHi];
36 long iRotHi = ipRot_H2_energy_sort[ipHi];
37
38 /* this is a guess of the collisional dissociation rate coefficient -
39 * will be multiplied by the sum of densities of all colliders
40 * except H2*/
41 double energy = H2_DissocEnergies[0] - states[ipHi].energy().WN();
42 ASSERT( energy > 0. );
43 /* we made this up - Boltzmann factor times rough coefficient */
44 H2_coll_dissoc_rate_coef[iVibHi][iRotHi] =
45 1e-14f * (realnum)sexp(energy/phycon.te_wn) * lgColl_dissoc_coll;
46
47 /* collisions with H2 - pre coefficient changed from 1e-8
48 * (from umist) to 1e-11 as per extensive discussion with Phillip Stancil */
49 H2_coll_dissoc_rate_coef_H2[iVibHi][iRotHi] =
50 1e-11f * (realnum)sexp(energy/phycon.te_wn) * lgColl_dissoc_coll;
51
52 /*fprintf(ioQQQ,"DEBUG coll_dissoc_rateee\t%li\t%li\t%.3e\t%.3e\n",
53 iVibHi,iRotHi,
54 H2_coll_dissoc_rate_coef[iVibHi][iRotHi],
55 H2_coll_dissoc_rate_coef_H2[iVibHi][iRotHi]);*/
56
57 for( long ipLo=0; ipLo<ipHi; ++ipLo )
58 {
59 long iVibLo = ipVib_H2_energy_sort[ipLo];
60 long iRotLo = ipRot_H2_energy_sort[ipLo];
61
62 ASSERT( states[ipHi].energy().WN() - states[ipLo].energy().WN() > 0. );
63
64 /* keep track of number of different collision routes */
65 ++numb_coll_trans;
66 for( long nColl=0; nColl<N_X_COLLIDER; ++nColl )
67 {
68 realnum newrate = H2_CollidRateEvalOne( iVibHi,iRotHi,iVibLo,iRotLo,
69 ipHi , ipLo , nColl, phycon.te );
70 CollRateCoeff[ipHi][ipLo][nColl] = newrate;
71 }
72 }
73 }
74
75 fixit(); // test that this does not matter with new rate table interpolation. Remove if it does not.
76#if 0
77 /* >>chng 05 nov 30, GS, rates decreases exponentially for low temperature, see Le Bourlot et al. 1999 */
78 /* Phillips mail--Apparently, the SD fit is only valid over the range of their calculations, 100-1000K.
79 * The rate should continue to fall exponentially with decreasing T, something like exp(-3900/T) for 0->1 and
80 * exp[-(3900-170.5)/T] for 1->0. It is definitely, not constant for T lower than 100 K, as far as we know.
81 * There may actually be a quantum tunneling effect which causes the rate to increase at lower T, but no
82 * one has calculated it (as far as I know) and it might happen below 1K or so.???*/
83 if( phycon.te < 100. )
84 {
85 /* first term in exp is suggested by Phillip, second temps in paren is to ensure continuity
86 * across 100K */
87 H2_CollRate[0][1][0][0][0] *= (realnum)(exp(-(3900-170.5)*(1./phycon.te - 1./100.)));
88 H2_CollRate[0][3][0][0][0] *= (realnum)(exp(-(3900-1015.1)*(1./phycon.te - 1./100.)));
89 H2_CollRate[0][2][0][1][0] *= (realnum)(exp(-(3900-339.3)*(1./phycon.te - 1./100.)));
90 }
91#endif
92
93 if( nTRACE >= n_trace_full )
94 fprintf(ioQQQ,
95 " collision rates updated for new temp, number of trans is %li\n",
96 numb_coll_trans);
97
98 return;
99}
100
101/* compute rate coefficient for a single quenching collision */
103 /*returns collision rate coefficient, cm-3 s-1 for quenching collision
104 * from this upper state */
105 long iVibHi, long iRotHi,long iVibLo,
106 /* to this lower state */
107 long iRotLo, long ipHi , long ipLo ,
108 /* colliders are H, He, H2(ortho), H2(para), and H+ */
109 long nColl,
110 double te_K )
111{
112 DEBUG_ENTRY( "H2_CollidRateEvalOne()" );
113
114 realnum rate = InterpCollRate( RateCoefTable[nColl], ipHi, ipLo, te_K );
115
116 /* this is option to use guess of collision rate coefficient */
117 if( (rate==0.f) &&
118 /* turn lgColl_gbar on/off with atom h2 gbar on off */
119 lgColl_gbar &&
120 /* but only if this does not mix ortho and para */
121 (H2_lgOrtho[0][iVibHi][iRotHi]==H2_lgOrtho[0][iVibLo][iRotLo] ) )
122 {
123 double ediff = states[ipHi].energy().WN() - states[ipLo].energy().WN();
124 rate = GbarRateCoeff( nColl, ediff );
125 }
126
127 rate *= lgColl_deexec_Calc;
128
129 /* >>chng 05 feb 09, add option to kill ortho - para collisions */
131 (H2_lgOrtho[0][iVibHi][iRotHi]-H2_lgOrtho[0][iVibLo][iRotLo]) )
132 rate = 0.;
133
134 if( lgH2_NOISE )
135 rate *= CollRateErrFac[ipHi][ipLo][nColl];
136
137 return rate;
138}
139
140STATIC realnum GbarRateCoeff( long nColl, double ediff )
141{
142 // these are fits to the existing collision data
143 // used to create g-bar rates
144 double gbarcoll[N_X_COLLIDER][3] =
145 {
146 {-9.9265 , -0.1048 , 0.456 },
147 {-8.281 , -0.1303 , 0.4931 },
148 {-10.0357, -0.0243 , 0.67 },
149 {-8.6213 , -0.1004 , 0.5291 },
150 {-9.2719 , -0.0001 , 1.0391 }
151 };
152
153 // do not let energy difference be smaller than 100 wn, the smallest
154 // difference for which we fit the rate coefficients
155 ediff = MAX2(100., ediff );
156
157 // the fit is log(K)=y_0+a*((x)^b), where K is the rate coefficient,
158 // and x is the energy in wavenumbers
159 realnum rate = (realnum)pow(10. ,
160 gbarcoll[nColl][0] + gbarcoll[nColl][1] *
161 pow(ediff,gbarcoll[nColl][2]) );
162
163 return rate;
164}
165
166void diatomics::H2_CollidRateRead( long int nColl )
167{
168 DEBUG_ENTRY( "H2_CollidRateRead()" );
169
170 if( coll_source[nColl].filename.size() == 0 && coll_source[nColl].magic == 0 )
171 return;
172
173 const char* chFilename = coll_source[nColl].filename.c_str();
174 long magic_expect = coll_source[nColl].magic;
175 char chLine[INPUT_LINE_LENGTH];
176 char chPath[FILENAME_PATH_LENGTH_2];
177 strcpy( chPath, path.c_str() );
178 strcat( chPath, input.chDelimiter );
179 strcat( chPath, chFilename );
180 FILE *ioDATA = open_data( chPath, "r" );
181
182 /* read the first line and check that magic number is ok */
183 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
184 {
185 fprintf( ioQQQ, " H2_CollidRateRead could not read first line of %s\n", chFilename );
187 }
188
189 /* magic number */
190 long n1 = atoi( chLine );
191 if( n1 != magic_expect )
192 {
193 fprintf( ioQQQ, " H2_CollidRateRead: the version of %s is not the current version.\n", chFilename );
194 fprintf( ioQQQ, " I expected to find the number %li and got %li instead.\n", magic_expect, n1 );
195 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
197 }
198
199 FunctPtr func = new FunctDiatoms( *this );
200 ReadCollisionRateTable( RateCoefTable[nColl], ioDATA, func, nLevels_per_elec[0] );
201 delete func;
202
203 fclose( ioDATA );
204
205 return;
206}
207
208void diatomics::GetIndices( long& ipHi, long& ipLo, const char* chLine, long& i ) const
209{
210 bool lgEOL;
211 long iVibHi = (long)FFmtRead( chLine, &i, strlen(chLine), &lgEOL );
212 long iRotHi = (long)FFmtRead( chLine, &i, strlen(chLine), &lgEOL );
213 long iVibLo = (long)FFmtRead( chLine, &i, strlen(chLine), &lgEOL );
214 long iRotLo = (long)FFmtRead( chLine, &i, strlen(chLine), &lgEOL );
215 ASSERT( iRotHi >= 0 && iVibHi >= 0 && iRotLo >= 0 && iVibLo >=0 );
216
217 // check that we actually included the levels in the model representation
218 // depending on the potential surface, the collision data set
219 // may not agree with our adopted model - skip those
220 if( ( iVibHi > nVib_hi[0] || iVibLo > nVib_hi[0] ) ||
221 ( iRotHi < Jlowest[0] || iRotLo < Jlowest[0] ) ||
222 ( iRotHi > nRot_hi[0][iVibHi] || iRotLo > nRot_hi[0][iVibLo] ) ||
223 /* some collision rates have the same upper and lower indices - skip them */
224 ( iVibHi == iVibLo && iRotHi == iRotLo ) )
225 {
226 ipHi = -1;
227 ipLo = -1;
228 return;
229 }
230
231 ipHi = ipEnergySort[0][iVibHi][iRotHi];
232 ipLo = ipEnergySort[0][iVibLo][iRotLo];
233 if( ipHi < ipLo )
234 swap( ipHi, ipLo );
235
236 return;
237}
238
double InterpCollRate(const CollRateCoeffArray &rate_table, const long &ipHi, const long &ipLo, const double &ftemp)
Definition atmdat.cpp:135
void ReadCollisionRateTable(CollRateCoeffArray &coll_rate_table, FILE *io, FunctPtr GetIndices, long nMolLevs, long nTemps, long nTrans)
Definition atmdat.cpp:14
Funct * FunctPtr
Definition atmdat.h:408
FILE * ioQQQ
Definition cddefines.cpp:7
const int FILENAME_PATH_LENGTH_2
Definition cddefines.h:249
#define ASSERT(exp)
Definition cddefines.h:578
sys_float sexp(sys_float x)
Definition service.cpp:914
#define STATIC
Definition cddefines.h:97
const int INPUT_LINE_LENGTH
Definition cddefines.h:254
#define EXIT_FAILURE
Definition cddefines.h:140
double FFmtRead(const char *chCard, long int *ipnt, long int last, bool *lgEOL)
Definition service.cpp:381
#define cdEXIT(FAIL)
Definition cddefines.h:434
float realnum
Definition cddefines.h:103
#define MAX2
Definition cddefines.h:782
char * read_whole_line(char *chLine, int nChar, FILE *ioIN)
Definition service.cpp:70
#define DEBUG_ENTRY(funcname)
Definition cddefines.h:684
void fixit(void)
Definition service.cpp:991
multi_arr< bool, 3 > H2_lgOrtho
Definition h2_priv.h:643
bool lgColl_gbar
Definition h2_priv.h:356
long int nVib_hi[N_ELEC]
Definition h2_priv.h:611
void H2_CollidRateEvalAll(void)
void GetIndices(long &ipHi, long &ipLo, const char *chLine, long &i) const
string path
Definition h2_priv.h:573
t_coll_source coll_source[N_X_COLLIDER]
Definition h2_priv.h:316
qList states
Definition h2_priv.h:565
multi_arr< realnum, 2 > H2_coll_dissoc_rate_coef
Definition h2_priv.h:668
string label
Definition h2_priv.h:571
long int Jlowest[N_ELEC]
Definition h2_priv.h:616
realnum H2_CollidRateEvalOne(long iVibHi, long iRotHi, long iVibLo, long iRotLo, long ipHi, long ipLo, long nColl, double temp_K)
valarray< long > ipVib_H2_energy_sort
Definition h2_priv.h:687
multi_arr< realnum, 3 > CollRateErrFac
Definition h2_priv.h:622
double H2_DissocEnergies[N_ELEC]
Definition h2_priv.h:609
vector< CollRateCoeffArray > RateCoefTable
Definition h2_priv.h:623
valarray< long > ipRot_H2_energy_sort
Definition h2_priv.h:689
multi_arr< realnum, 3 > CollRateCoeff
Definition h2_priv.h:621
bool lgColl_dissoc_coll
Definition h2_priv.h:362
bool lgH2_ortho_para_coll_on
Definition h2_priv.h:372
bool lgH2_NOISE
Definition h2_priv.h:383
bool lgColl_deexec_Calc
Definition h2_priv.h:359
int n_trace_full
Definition h2_priv.h:404
multi_arr< realnum, 2 > H2_coll_dissoc_rate_coef_H2
Definition h2_priv.h:671
valarray< long > nRot_hi[N_ELEC]
Definition h2_priv.h:613
int nTRACE
Definition h2_priv.h:399
long int nLevels_per_elec[N_ELEC]
Definition h2_priv.h:618
void H2_CollidRateRead(long int nColl)
multi_arr< long int, 3 > ipEnergySort
Definition h2_priv.h:690
void swap(count_ptr< T > &a, count_ptr< T > &b)
Definition count_ptr.h:82
FILE * open_data(const char *fname, const char *mode, access_scheme scheme)
Definition cpu.cpp:625
const int N_X_COLLIDER
Definition h2_priv.h:13
t_input input
Definition input.cpp:12
STATIC realnum GbarRateCoeff(long nColl, double ediff)
t_phycon phycon
Definition phycon.cpp:6