cloudy  trunk
 All Data Structures 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-2008 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 
6 #include "cddefines.h"
7 #include "phycon.h"
8 #include "dense.h"
9 #include "taulines.h"
10 #include "h2.h"
11 #include "h2_priv.h"
12 #include "mole.h"
13 
14 /* set true to print all collision rates then quit */
15 #define PRT_COLL false
16 
17 /* following are related to H2 - He collision data set
18 * H2_He_coll_init, H2_He_coll */
19 #define N_H2_HE_FIT_PAR 8
21 static bool **lgDefn_H2He_coll;
22 
23 /* compute rate coefficient for a single quenching collision */
25  /*returns collision rate coefficient, cm-3 s-1 for quenching collision
26  * from this upper state */
27  long iVibHi, long iRotHi,long iVibLo,
28  /* to this lower state */
29  long iRotLo, long ipHi , long ipLo ,
30  /* colliders are H, He, H2(ortho), H2(para), and H+ */
31  long nColl )
32 {
33  double fitted;
34  realnum rate;
35  double t3Plus1 = phycon.te/1000. + 1.;
36  double t3Plus1Squared = POW2(t3Plus1);
37  /* these are fits to the existing collision data
38  * used to create g-bar rates */
39  double gbarcoll[N_X_COLLIDER][3] =
40  {
41  {-9.9265 , -0.1048 , 0.456 },
42  {-8.281 , -0.1303 , 0.4931 },
43  {-10.0357, -0.0243 , 0.67 },
44  {-8.6213 , -0.1004 , 0.5291 },
45  {-9.2719 , -0.0001 , 1.0391 }
46  };
47 
48  DEBUG_ENTRY( "H2_CollidRateEvalOne()" );
49 
50  /* first do special cases
51  * ORNL He collision data */
52  if( nColl == 1 && mole.lgH2_He_ORNL )
53  {
54  /* ORNL in use */
55 
56  /* H2 - He collisions
57  * The H2 - He collisional data set is controlled by the
58  * atom H2 He collisions command
59  * either mole.lgH2_He_Meudon or mole.lgH2_He_ORNL must be true
60  * (this is asserted). Meudon is collider 1 and Stancil is 5. */
61  /*>>chng 07 apr 04, by GS - bugfix, had used 1 for collider for g-bar */
62  if( (fitted=H2_He_coll(ipHi , ipLo , phycon.te ))>0. )
63  {
64  /* >>chng 05 oct 02, add this collider
65  * CHECK - this is deexcitation - is that correct
66  * this is new H2 - He collision data - use it if positive
67  * not all states have collision data */
68  rate = (realnum)fitted*mole.lgColl_deexec_Calc;
69  /*if( PRT_COLL )
70  fprintf(ioQQQ,"col fit\t%li\t%li\t%.4e\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
71  ipLo,ipHi,phycon.te,nColl,
72  iVibHi,iRotHi,iVibLo,iRotLo,
73  energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
74  H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] );*/
75  }
76 
77  /* use g-bar g bar guess of rate coefficient for
78  * collision with He, this does not change ortho & para
79  * turn mole.lgColl_gbar on/off with atom h2 gbar on off */
80  else if( mole.lgColl_gbar &&
81  (H2_lgOrtho[0][iVibHi][iRotHi]-H2_lgOrtho[0][iVibLo][iRotLo]==0) )
82  {
83  /* the fit is log(K)=y_0+a*((x)^b), where K is the rate coefficient,
84  * and x is the energy in wavenumbers */
85  double ediff = energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo];
86 
87  /* do not let energy difference be smaller than 100 wn, the smallest
88  * difference for which we fit the rate coefficients */
89  ediff = MAX2(100., ediff );
90  rate = (realnum)pow(10. ,
91  gbarcoll[nColl][0] + gbarcoll[nColl][1] *
92  pow(ediff,gbarcoll[nColl][2]) )*mole.lgColl_deexec_Calc;
93 
94  /*if( PRT_COLL )
95  fprintf(ioQQQ,"col gbr\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
96  nColl+10,
97  iVibHi,iRotHi,iVibLo,iRotLo,
98  energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
99  H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] );*/
100  }
101  else
102  rate = 0;
103  }
104  else if( nColl==4 )
105  {
106  /* collisions of H2 with protons - of this group, these are only
107  * that cause ortho - para conversion */
108  /* >>refer H2 coll Hp Gerlich, D., 1990, J. Chem. Phys., 92, 2377-2388 */
109  if( CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][1] != 0 )
110  {
111  rate = CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][0] * 1e-10f *
112  /* sec fit coef was dE in milli eV */
113  (realnum)sexp( CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][1]/1000./phycon.te_eV)*mole.lgColl_deexec_Calc;
114  /*if( PRT_COLL )
115  fprintf(ioQQQ,"col fit\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
116  nColl,
117  iVibHi,iRotHi,iVibLo,iRotLo,
118  energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
119  H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] );*/
120  }
121  /* this is option to use guess of rate coefficient for ortho-para
122  * conversion by collision with protons */
123  /* turn mole.lgColl_gbar on/off with atom h2 gbar on off */
124  else if( mole.lgColl_gbar )
125  {
126  /* the fit is log(K)=y_0+a*((x)^b), where K is the rate coefficient,
127  * and x is the energy in wavenumbers */
128  double ediff = energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo];
129  ediff = MAX2(100., ediff );
130  rate = (realnum)pow(10. ,
131  gbarcoll[nColl][0] + gbarcoll[nColl][1] *
132  pow(ediff ,gbarcoll[nColl][2]) )*mole.lgColl_deexec_Calc;
133 
134  /*if( PRT_COLL )
135  fprintf(ioQQQ,"col gbr\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
136  nColl+10,
137  iVibHi,iRotHi,iVibLo,iRotLo,
138  energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
139  H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] );*/
140  }
141  else
142  rate = 0;
143  }
144  else if( CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][0]!= 0 )
145  {
146  /* these are the fits from
147  *>>refer H2 coll Le Bourlot, J., Pineau des Forets,
148  *>>refercon G., & Flower, D.R. 1999, MNRAS, 305, 802
149  * evaluate collision rates for those with real collision data */
150 
151  double r = CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][0] +
152  CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][1]/t3Plus1 +
153  CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][2]/t3Plus1Squared;
154 
155  rate = (realnum)pow(10.,r)*mole.lgColl_deexec_Calc;
156 
157  /*if( PRT_COLL )
158  fprintf(ioQQQ,"col fit\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
159  nColl,
160  iVibHi,iRotHi,iVibLo,iRotLo,
161  energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
162  H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] );*/
163  }
164  /* this is option to use guess of collision rate coefficient - but only if this is
165  * a downward transition that does not mix ortho and para */
166  /* turn mole.lgColl_gbar on/off with atom h2 gbar on off */
167  else if( mole.lgColl_gbar &&
168  (H2_lgOrtho[0][iVibHi][iRotHi]-H2_lgOrtho[0][iVibLo][iRotLo]==0) )
169  {
170  /* the fit is log(K)=y_0+a*((x)^b), where K is the rate coefficient,
171  * and x is the energy in wavenumbers */
172  double ediff = energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo];
173  /* do not let energy difference be smaller than 100 wn, the smallest
174  * difference for which we fit the rate coefficients */
175  ediff = MAX2(100., ediff );
176  rate = (realnum)pow(10. ,
177  gbarcoll[nColl][0] + gbarcoll[nColl][1] *
178  pow(ediff,gbarcoll[nColl][2]) )*mole.lgColl_deexec_Calc;
179  /* this is hack to change H2 H collision rates */
180  if( nColl == 0 && h2.lgH2_H_coll_07 )
181  H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] *= 100.;
182 
183  /*if( PRT_COLL )
184  fprintf(ioQQQ,"col gbr\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
185  nColl+10,
186  iVibHi,iRotHi,iVibLo,iRotLo,
187  energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
188  H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] );*/
189  }
190  else
191  rate = 0;
192 
193 
194  /* >>chng 05 feb 09, add option to kill ortho - para collisions */
196  (H2_lgOrtho[0][iVibHi][iRotHi]-H2_lgOrtho[0][iVibLo][iRotLo]) )
197  rate = 0.;
198 
199  {
200  enum {DEBUG_LOC=false};
201  if( DEBUG_LOC )
202  {
203  fprintf(ioQQQ,"bugcoll\tiVibHi\t%li\tiRotHi\t%li\tiVibLo\t%li\tiRotLo\t%li\tcoll\t%.2e\n",
204  iVibHi,iRotHi,iVibLo,iRotLo,
205  H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] );
206  }
207  }
208  return rate;
209 }
210 
211 /*H2_CollidRateEvalAll - set H2 collision rate coefficients */
213 {
214  long int numb_coll_trans = 0;
215  double excit;
216  long int iElecHi , iElecLo , ipHi , iVibHi , iRotHi ,
217  ipLo , iVibLo , iRotLo , nColl;
218 
219  DEBUG_ENTRY( "H2_CollidRateEvalAll()" );
220 
221  if( PRT_COLL )
222  fprintf(ioQQQ,"H2 coll deex rate coef\n"
223  "VibHi\tRotHi\tVibLo\tRotLo\trate\n");
224 
225  iElecHi = 0;
226  iElecLo = 0;
228  fprintf(ioQQQ,"H2 set collision rates\n");
229  /* loop over all possible collisional changes within X
230  * and set collision rates, which only depend on Te
231  * will go through array in energy order since coll trans do not
232  * correspond to a line
233  * collisional dissociation rate coefficient, units cm3 s-1 */
234  H2_coll_dissoc_rate_coef[0][0] = 0.;
235  H2_coll_dissoc_rate_coef_H2[0][0] = 0.;
236  for( ipHi=0; ipHi<nLevels_per_elec[0]; ++ipHi )
237  {
238  double energy;
239 
240  /* obtain the proper indices for the upper level */
241  long int ip = H2_ipX_ener_sort[ipHi];
242  iVibHi = ipVib_H2_energy_sort[ip];
243  iRotHi = ipRot_H2_energy_sort[ip];
244 
245  /* this is a guess of the collisional dissociation rate coefficient -
246  * will be multiplied by the sum of densities of all colliders
247  * except H2*/
248  energy = H2_DissocEnergies[0] - energy_wn[0][iVibHi][iRotHi];
249  ASSERT( energy > 0. );
250  /* we made this up - Boltzmann factor times rough coefficient */
251  H2_coll_dissoc_rate_coef[iVibHi][iRotHi] =
252  1e-14f * (realnum)sexp(energy/phycon.te_wn) * mole.lgColl_dissoc_coll;
253 
254  /* collisions with H2 - pre coefficient changed from 1e-8
255  * (from umist) to 1e-11 as per extensive discussion with Phillip Stancil */
256  H2_coll_dissoc_rate_coef_H2[iVibHi][iRotHi] =
257  1e-11f * (realnum)sexp(energy/phycon.te_wn) * mole.lgColl_dissoc_coll;
258 
259  /*fprintf(ioQQQ,"DEBUG coll_dissoc_rateee\t%li\t%li\t%.3e\t%.3e\n",
260  iVibHi,iRotHi,
261  H2_coll_dissoc_rate_coef[iVibHi][iRotHi],
262  H2_coll_dissoc_rate_coef_H2[iVibHi][iRotHi]);*/
263 
264  for( ipLo=0; ipLo<ipHi; ++ipLo )
265  {
266  ip = H2_ipX_ener_sort[ipLo];
267  iVibLo = ipVib_H2_energy_sort[ip];
268  iRotLo = ipRot_H2_energy_sort[ip];
269 
270  ASSERT( energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo] > 0.);
271 
272  /* in following the colliders are H, He, H2(ortho), H2(para), and H+ */
273  /* fits were read in from the following files: "H2_coll_H.dat" ,
274  * "H2_coll_He.dat" , "H2_coll_H2ortho.dat" ,"H2_coll_H2para.dat",
275  * "H2_coll_Hp.dat" */
276 
277  /* keep track of number of different collision routes */
278  ++numb_coll_trans;
279  /* this is sum over all different colliders, except last two which are special,
280  * linear rather than log formula for that one for second to last,
281  * and special fitting formula for last */
282  if( PRT_COLL && iVibHi == 1 && iRotHi==3)
283  fprintf(ioQQQ,"%li\t%li\t%li\t%li",
284  iVibHi,iRotHi,iVibLo,iRotLo);
285  for( nColl=0; nColl<N_X_COLLIDER; ++nColl )
286  {
287  H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] =
288  H2_CollidRateEvalOne( iVibHi,iRotHi,iVibLo,iRotLo,
289  ipHi , ipLo , nColl );
290  if( PRT_COLL && iVibHi == 1 && iRotHi==3)
291  fprintf(ioQQQ,"\t%.2e",H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] );
292  }
293  if( PRT_COLL && iVibHi == 1 && iRotHi==3)
294  fprintf(ioQQQ,"\n");
295  }
296  }
297  if( PRT_COLL )
298  cdEXIT( EXIT_FAILURE );
299 
300  /* at this stage the collision rates that came in from the large data files
301  * have been entered into the H2_CollRate array. Now add on three extra collision
302  * terms, the ortho para atomic H collision rates from
303  * >>>refer H2 collision Sun, Y., & Dalgarno, A., 1994, ApJ, 427, 1053-1056
304  */
305  nColl = 0;
306  iElecHi = 0;
307  iElecLo = 0;
308  iVibHi = 0;
309  iVibLo = 0;
310 
311  /* >>chng 02 nov 13, the Sun and Dalgarno rates diverge to + inf below this temp */
312  /* >>chng 05 feb 09, do not return zero when T < 100 - instead, don't let T fall below 100 */
313  double excit1;
314  double te_used = MAX2( 100. , phycon.te );
315  /* this is the J=1-0 downward collision rate */
316  iRotLo = 0;
317  iRotHi = 1;
318  excit1 = sexp( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyK/te_used);
319  excit = sexp( -(POW2(5.30-460./te_used)-21.2) )*1e-13;
320 
321  H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][0] = (realnum)(
322  excit*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->g/
323  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g /
324  /* >>chng 02 nov 13, from 2nd to first */
325  SDIV(excit1) )*mole.lgColl_deexec_Calc *
326  /* option to disable ortho-para conversion by coll with grains */
328 
329  /*if( PRT_COLL )
330  fprintf(ioQQQ,"col o-p\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
331  nColl,
332  iVibHi,iRotHi,iVibLo,iRotLo,
333  energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
334  H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] );*/
335 
336  /* this is the J=3-0 downward collision rate */
337  iRotLo = 0;
338  iRotHi = 3;
339  excit = sexp( -(POW2(6.36-373./te_used)-34.5) )*1e-13;
340  excit1 = sexp( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyK/te_used);
341  H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][0] = (realnum)(
342  excit*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->g/
343  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g /
344  SDIV(excit1) )*mole.lgColl_deexec_Calc *
345  /* option to disable ortho-para conversion by coll with grains */
347 
348  /*if( PRT_COLL )
349  fprintf(ioQQQ,"col o-p\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
350  nColl,
351  iVibHi,iRotHi,iVibLo,iRotLo,
352  energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
353  H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] );*/
354 
355  /* this is the downward J=2-1 collision rate */
356  iRotLo = 1;
357  iRotHi = 2;
358  excit = sexp( -(POW2(5.35-454./te_used)-23.1 ) )*1e-13;
359  excit1 = sexp( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyK/te_used);
360  H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][0] = (realnum)(
361  excit*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->g/
362  H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g /
363  SDIV(excit1) )*mole.lgColl_deexec_Calc *
364  /* option to disable ortho-para conversion by coll with grains */
366 
367  /* >>chng 05 nov 30, GS, rates decreases exponentially for low temperature, see Le Bourlot et al. 1999 */
368  /* Phillips mail--Apparently, the SD fit is only valid over the range of their calculations, 100-1000K.
369  * The rate should continue to fall exponentially with decreasing T, something like exp(-3900/T) for 0->1 and
370  * exp[-(3900-170.5)/T] for 1->0. It is definitely, not constant for T lower than 100 K, as far as we know.
371  * There may actually be a quantum tunneling effect which causes the rate to increase at lower T, but no
372  * one has calculated it (as far as I know) and it might happen below 1K or so.???*/
373  if( phycon.te < 100. )
374  {
375  /* first term in exp is suggested by Phillip, second temps in paren is to ensure continuity
376  * across 100K */
377  H2_CollRate[0][1][0][0][0] = (realnum)(H2_CollRate[0][0][1][0][0]*exp(-(3900-170.5)*(1./phycon.te - 1./100.)));
378  H2_CollRate[0][3][0][0][0] = (realnum)(H2_CollRate[0][0][3][0][0]*exp(-(3900-1015.1)*(1./phycon.te - 1./100.)));
379  H2_CollRate[0][2][0][1][0] = (realnum)(H2_CollRate[0][0][2][0][1]*exp(-(3900-339.3)*(1./phycon.te - 1./100.)));
380  }
381 
382  /*if( PRT_COLL )
383  fprintf(ioQQQ,"col o-p\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n",
384  nColl,
385  iVibHi,iRotHi,iVibLo,iRotLo,
386  energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo],
387  H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] );*/
388 
390  fprintf(ioQQQ,
391  " collision rates updated for new temp, number of trans is %li\n",
392  numb_coll_trans);
393 
394  /* quit it we are only printing - but do this after printing coll rates
395  if( PRT_COLL )
396  cdEXIT( EXIT_FAILURE ); */
397  return;
398 }
399 
400 /*H2_CollidRateRead read collision rates */
401 void H2_CollidRateRead( long int nColl )
402 {
403  /* the colliders are H, He, H2 ortho, H2 para, H+
404  * these are the default file names. they can be overridden with
405  * the SET ATOMIC DATA MOLECULE H2 command */
406 
407  /*NB thesee must be kept parallel with labels in chH2ColliderLabels */
408 
409  const char* cdDATAFILE[N_X_COLLIDER] =
410  {
411  /* 0 */"H2_coll_H_07.dat",
412  /* 1 */"H2_coll_He_LeBourlot.dat",
413  /* 2 */"H2_coll_H2ortho.dat",
414  /* 3 */"H2_coll_H2para.dat",
415  /* 4 */"H2_coll_Hp.dat"
416  };
417 
418  FILE *ioDATA;
419  char chLine[FILENAME_PATH_LENGTH_2];
420  const char* chFilename;
421  long int i, n1;
422  long int iVibHi , iVibLo , iRotHi , iRotLo;
423  long int magic_expect = -1;
424 
425  DEBUG_ENTRY( "H2_CollidRateRead()" );
426 
427  if( nColl == 0 )
428  {
429  /* which H2 - H data set? */
430  if( h2.lgH2_H_coll_07 )
431  {
432  magic_expect = 71106;
433  chFilename = "H2_coll_H_07.dat";
434  }
435  else
436  {
437  /* the 1999 data set */
438  magic_expect = 20429;
439  chFilename = "H2_coll_H_99.dat";
440  }
441  }
442  else if( nColl == 1 )
443  {
444  /* which H2 - He data set? */
445  if( mole.lgH2_He_ORNL )
446  {
447  /* >>chng 07 may 12, update magic number when one transition in H2 - H2 file was
448  * replaced - the fitting coefficients had terms of order -8e138 which fpe in float */
449  long int magic_found,
450  magic_expect = 70513;
451  /* special case, new data file from Oak Ridge project -
452  * call init routine and return - data is always read in when large H2 is included -
453  * but data are only used (for now, mid 2005) when command
454  * atom H2 He OLD (Meudon) NEW (Stancil) and OFF given */
455  /*H2_He_coll_init receives the name of the file that contains the fitting coefficients
456  * of all transitions and read into 3d vectors. It outputs 'test.out' to test the arrays*/
457  chFilename = "H2_coll_He_ORNL.dat";
458  if( (magic_found = H2_He_coll_init( chFilename )) != magic_expect )
459  {
460  fprintf(ioQQQ,"The H2 - He collision data file H2_coll_He_ORNL.dat does not have the correct magic number.\n");
461  fprintf(ioQQQ,"I found %li but expected %li\n", magic_found , magic_expect );
462 
463  cdEXIT(EXIT_FAILURE);
464  }
465  return;
466  }
467  else
468  {
469  magic_expect = 20429;
470  /* the Le Bourlot et al data */
471  chFilename = "H2_coll_He_LeBourlot.dat";
472  }
473  }
474  else
475  {
476  /* files with only one version */
477  magic_expect = 20429;
478  chFilename = cdDATAFILE[nColl];
479  }
480 
481  ioDATA = open_data( chFilename, "r" );
482 
483  /* read the first line and check that magic number is ok */
484  if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
485  {
486  fprintf( ioQQQ, " H2_CollidRateRead could not read first line of %s\n", chFilename );
487  cdEXIT(EXIT_FAILURE);
488  }
489 
490  /* level 1 magic number */
491  n1 = atoi( chLine );
492 
493  /* magic number
494  * the following is the set of numbers that appear at the start of level1.dat 01 08 10 */
495  if( n1 != magic_expect )
496  {
497  fprintf( ioQQQ,
498  " H2_CollidRateRead: the version of %s is not the current version.\n", chFilename );
499  fprintf( ioQQQ,
500  " I expected to find the number %li and got %li instead.\n" ,
501  magic_expect , n1 );
502  fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
503  cdEXIT(EXIT_FAILURE);
504  }
505 
506  while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
507  {
508  if( chLine[0]=='#' )
509  continue;
510  double a[3];
511  sscanf(chLine,"%li\t%li\t%li\t%li\t%le\t%le\t%le",
512  &iVibHi ,&iRotHi , &iVibLo , &iRotLo , &a[0],&a[1],&a[2] );
513  /*fprintf(ioQQQ,"DEBUG %li %li %li %li \n",
514  iVibHi , iRotHi , iVibLo , iRotLo);*/
515  ASSERT( iRotHi >= 0 && iVibHi >= 0 && iRotLo >= 0 && iVibLo >=0 );
516 
517  /* check that we actually included the levels in the model representation
518  * depending on the potential surface, the collision date set
519  * may not agree with our adopted model - skip those */
520  if( iVibHi > VIB_COLLID ||
521  iVibLo > VIB_COLLID ||
522  iRotHi > h2.nRot_hi[0][iVibHi] ||
523  iRotLo > h2.nRot_hi[0][iVibLo])
524  {
525  iVibHi = -1;
526  continue;
527  }
528 
529  /* some collision rates have the same upper and lower indices - skip them */
530  if( !( (iVibHi == iVibLo) && (iRotHi == iRotLo )) )
531  {
532  /* this is downward transition - make sure that the energy difference is positive */
533  ASSERT( (energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo] ) > 0. );
534  for( i=0; i<3; ++i )
535  {
536  CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][i] = (realnum)a[i];
537  }
538 
539  /* this prints all levels with rates
540  fprintf(ioQQQ,"no\t%li\t%li\t%li\t%li\t%.2e\t%.2e\t%.2e\n",
541  iVibHi,iRotHi,iVibLo,iRotLo,a[0],a[1],a[2]);*/
542  }
543  }
544  fclose( ioDATA );
545 
546  return;
547 }
548 /* end of H2_CollidRateRead */
549 
550 /* This code was written by Terry Yun, 2005 */
551 /* H2_He_coll_init receives the name of the file that contains the fitting
552 * coefficients of all transitions and read into 3d vectors.
553 * It returns magic number*/
554 long int H2_He_coll_init(const char FILE_NAME_IN[])
555 {
556 
557  int i, j;
558  long int magic;
559  double par[N_H2_HE_FIT_PAR];
560  char line[INPUT_LINE_LENGTH];
561  int h2_i, h2_f, he_i, he_f;/* target, projectile initial and final indices */
562  char quality, space;
563  /*'space' variable is for reading spaces between characters */
564  /* scanf can skip spaces(or tap) when it reads numbers, but not for characters. */
565  double error;
566 
567  FILE *ifp;
568 
569  DEBUG_ENTRY( "H2_He_coll_init()" );
570 
571  /* create space for two arrays - H2_He_coll_fit_par will contain the
572  * fit parameters in form [init state][final state][param]
573  * lgDefn_H2He_coll is flag saying whether data exists */
574  H2_He_coll_fit_par = (realnum ***)MALLOC(sizeof(realnum **)*(unsigned)nLevels_per_elec[0] );
575  lgDefn_H2He_coll = (bool**)MALLOC(sizeof(bool *)*(unsigned)nLevels_per_elec[0] );
576  for( i=1; i<nLevels_per_elec[0]; ++i )
577  {
578  lgDefn_H2He_coll[i] = (bool*)MALLOC(sizeof(bool)*(unsigned)i );
579  H2_He_coll_fit_par[i] = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)i );
580  for( j=0; j<i; ++j)
581  H2_He_coll_fit_par[i][j] = (realnum *)MALLOC(sizeof(realnum)*(unsigned)N_H2_HE_FIT_PAR );
582  }
583 
584  /* set a flag saying whether data exits: initially everything is '0' */
585  for( i=1; i<nLevels_per_elec[0]; i++ )
586  {
587  for( j=0; j<i; j++ )
588  {
589  lgDefn_H2He_coll[i][j] = false;
590  }
591  }
592 
593  /*open the input file and put into 3d arrays*/
594  ifp = open_data( FILE_NAME_IN, "r" );
595 
596  /*read the magic number*/
597  if( read_whole_line(line, (int)sizeof(line), ifp)==NULL )
598  {
599  printf("DISASTER H2_He_coll_init can't read first line of %s\n", FILE_NAME_IN);
600  cdEXIT(EXIT_FAILURE);
601  }
602  sscanf(line, "%li", &magic);
603 
604  /* read the file until the end of line */
605  while( read_whole_line(line, (int)sizeof(line), ifp) != NULL )
606  {
607  /* skip any line that starts with '#' */
608  if( line[0]!='#' )
609  {
610  sscanf(line, "%i%i%i%i%c%c%c%c%lf%lf%lf%lf%lf%lf%lf%lf%lf", &h2_i,
611  &h2_f, &he_i, &he_f,&space, &space, &space, &quality,
612  &error, &par[0], &par[1], &par[2], &par[3], &par[4],
613  &par[5], &par[6], &par[7]);
614 
615  /* >>chng 07 may 28, extensive testing showed rate > 1e38 for
616  * temperatures where fitting equation hit a pole. Exchange
617  * with Phillip Stancil:
618  "You are correct. The problem is the parameter d which is
619  negative in all the problem cases you found. A pole occurs
620  when d*T/1000 + 1 =0.
621 
622  I found that there are a total of 11 transitions with d<1.
623  A quick solution is to take abs(d). I checked one case and the
624  resulting function went smoothly through the pole. On either
625  side of the pole the two functions looked the same on a log-log plot.
626 
627  LeBourlot used a similar function, but had d=1. So, they never
628  got a pole. I guess I tried to make the function a little too
629  flexible.
630 
631  There is a similar term with g*T/1000 + 1. There are no fits
632  with a negative g, but you might take abs(g) to be on the space
633  side for the future."
634  */
635  if( par[3] < 0. )
636  {
637  /*fprintf(ioQQQ,"DEBUG negative [3]=%.2e\n", par[3]);*/
638  par[3] = -par[3];
639  }
640  if( par[6] < 0. )
641  {
642  /*fprintf(ioQQQ,"DEBUG negative [6]=%.2e\n", par[6]);*/
643  par[6] = -par[6];
644  }
645  /* input file is on physics or Fortran scale so lowest energy
646  * index is 1. Store on C scale */
647  --h2_f;
648  --h2_i;
649  /* set true when the indices are defined */
650  lgDefn_H2He_coll[h2_i][h2_f] = true;
651  for( i=0; i<N_H2_HE_FIT_PAR; i++ )
652  /* assigning the parameters to 3d array */
653  H2_He_coll_fit_par[h2_i][h2_f][i] = (realnum)par[i];
654  }
655  }
656  fclose(ifp);
657  return magic;
658 }
659 
660 /* This code was written by Terry Yun, 2005 */
661 /* H2_He_coll Interpolate the rate coefficients
662 * It receives initial and final transition and temperature
663 * The range of the temperature is between 2K - 1e8K */
664 double H2_He_coll(int init, int final, double temp)
665 {
666 
667  double k, t3Plus1, b2, c2, f2, t3;
668 
669  DEBUG_ENTRY( "H2_He_coll()" );
670 
671  /* Invalid entries returns '-1':the initial indices are smaller than the final indices */
672  if( temp<2 || temp > 1e8 )
673  {
674  k = -1;
675  }
676  else if( init <= final )
677  {
678  k = -1;
679  }
680  /* Invalid returns '-1': the indices are greater than 302 or smaller than 0 */
681  else if( init < 0 || init >302 || final < 0 || final > 302 )
682  {
683  k = -1;
684  }
685  /* Undefined indices returns '0' */
686  else if( !lgDefn_H2He_coll[init][final] )
687  {
688  k = -1;
689  }
690  /* defined indices */
691  else if( lgDefn_H2He_coll[init][final] )
692  {
693  /* the fitting equation we used:
694  * k_(vj,v'j') = 10^(a + b/(d*T/10^3+1) + c/t^2)
695  * + 10^(e + f/(g*T/10^3+1)**h)
696  * T in K, k in cm3/s. */
697  double sum1 , sum2;
698 
699  /* this kludge is to bypass errors in fitting formula - there
700  * is a pole possible at around 1e6 K and rates can diverge
701  * and high temperatures */
703  temp = MIN2( 1e4 , temp );
704 
705  t3 = temp/1e3;
706  /* t3Plus1 = T/1000 + 1*/
707  t3Plus1 = t3+1;
708  b2 = H2_He_coll_fit_par[init][final][1]/(H2_He_coll_fit_par[init][final][3]*t3+1);
709  c2 = H2_He_coll_fit_par[init][final][2]/(t3Plus1*t3Plus1);
710  {
711  enum {DEBUG_LOC=false};
712  if( DEBUG_LOC )
713  {
714  fprintf(ioQQQ,"bug H2 He coll\t%i %i %.3e %.3e %.3e \n",
715  init,final,
716  H2_He_coll_fit_par[init][final][5],
717  H2_He_coll_fit_par[init][final][6]*t3+1,
718  H2_He_coll_fit_par[init][final][7]
719  /*pow(H2_He_coll_fit_par[init][final][6]*t3+1,H2_He_coll_fit_par[init][final][7])*/
720  );
721  }
722  }
723  /* this is log of f2 - see whether it is within bounds */
724  sum1 = H2_He_coll_fit_par[init][final][7] * log10( H2_He_coll_fit_par[init][final][6]*t3+1. );
725  /* this protects against overflow */
726  if( fabs(sum1)< 38. )
727  {
728  /* >>chng 06 jun 07, Mitchell Martin, from following to equivalent below. This was basically a
729  * bug in the pgcc compiler that caused h2_pdr_leiden_f1 to throw an fpe
730  * the changed code is equivalent */
731  /*f2 = H2_He_coll_fit_par[init][final][5]/pow(H2_He_coll_fit_par[init][final][6]*t3+1.,H2_He_coll_fit_par[init][final][7]);*/
732  f2 = H2_He_coll_fit_par[init][final][5]/pow(10. , sum1 );
733  }
734  else
735  f2= 0.;
736  {
737  enum {DEBUG_LOC=false };
738  if( DEBUG_LOC )
739  {
740  fprintf(ioQQQ,"bug H2 He coll\t%i %i %.3e %.3e %.3e %.3e %.3e sum %.3e %.3e \n",
741  init,final,
742  H2_He_coll_fit_par[init][final][0],
743  b2,
744  c2,
745  H2_He_coll_fit_par[init][final][4],
746  f2 ,
747  H2_He_coll_fit_par[init][final][0]+b2+c2,
748  H2_He_coll_fit_par[init][final][4]+f2);
749  }
750  }
751 
752  sum1 = H2_He_coll_fit_par[init][final][0]+b2+c2;
753  sum2 = H2_He_coll_fit_par[init][final][4]+f2;
754  k = 0.;
755  if( fabs(sum1) < 38. )
756  k += pow(10., sum1 );
757  if( fabs(sum2) < 38. )
758  k += pow(10., sum2 );
759 
760  if( k>1e-6 )
761  {
762  /* rate of 1e-6 is largest possible rate according to Phillip
763  * Stancil email 07 may 27 */
764  fprintf(ioQQQ,"PROBLEM H2-He rate coefficient hit cap, upper index of %i"
765  " lower index of %i rate was %.2e resetting to 1e-6, Te=%e\n",
766  init , final , k , phycon.te );
767  k = 1e-6;
768  }
769  {
770  enum {DEBUG_LOC=false};
771  if( DEBUG_LOC )
772  {
773  fprintf(ioQQQ,"DEBUG H2 He rate %.3e \n",
774  k);
775  }
776  }
777  }
778  /*unknown invalid entries returns '99'*/
779  else
780  {
781  k = 99;
782  }
783  return k;
784 }

Generated for cloudy by doxygen 1.8.4