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