cloudy trunk
|
00001 /* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*H2_CollidRateRead read collision rates */ 00004 /*H2_CollidRateEvalAll - set H2 collision rates */ 00005 00006 #include "cddefines.h" 00007 #include "phycon.h" 00008 #include "dense.h" 00009 #include "taulines.h" 00010 #include "h2.h" 00011 #include "h2_priv.h" 00012 #include "mole.h" 00013 00014 /* set true to print all collision rates then quit */ 00015 #define PRT_COLL false 00016 00017 /* following are related to H2 - He collision data set 00018 * H2_He_coll_init, H2_He_coll */ 00019 #define N_H2_HE_FIT_PAR 8 00020 static realnum ***H2_He_coll_fit_par; 00021 static bool **lgDefn_H2He_coll; 00022 00023 /* compute rate coefficient for a single quenching collision */ 00024 STATIC realnum H2_CollidRateEvalOne( 00025 /*returns collision rate coefficient, cm-3 s-1 for quenching collision 00026 * from this upper state */ 00027 long iVibHi, long iRotHi,long iVibLo, 00028 /* to this lower state */ 00029 long iRotLo, long ipHi , long ipLo , 00030 /* colliders are H, He, H2(ortho), H2(para), and H+ */ 00031 long nColl ) 00032 { 00033 double fitted; 00034 realnum rate; 00035 double t3Plus1 = phycon.te/1000. + 1.; 00036 double t3Plus1Squared = POW2(t3Plus1); 00037 /* these are fits to the existing collision data 00038 * used to create g-bar rates */ 00039 double gbarcoll[N_X_COLLIDER][3] = 00040 { 00041 {-9.9265 , -0.1048 , 0.456 }, 00042 {-8.281 , -0.1303 , 0.4931 }, 00043 {-10.0357, -0.0243 , 0.67 }, 00044 {-8.6213 , -0.1004 , 0.5291 }, 00045 {-9.2719 , -0.0001 , 1.0391 } 00046 }; 00047 00048 DEBUG_ENTRY( "H2_CollidRateEvalOne()" ); 00049 00050 /* first do special cases 00051 * ORNL He collision data */ 00052 if( nColl == 1 && mole.lgH2_He_ORNL ) 00053 { 00054 /* ORNL in use */ 00055 00056 /* H2 - He collisions 00057 * The H2 - He collisional data set is controlled by the 00058 * atom H2 He collisions command 00059 * either mole.lgH2_He_Meudon or mole.lgH2_He_ORNL must be true 00060 * (this is asserted). Meudon is collider 1 and Stancil is 5. */ 00061 /*>>chng 07 apr 04, by GS - bugfix, had used 1 for collider for g-bar */ 00062 if( (fitted=H2_He_coll(ipHi , ipLo , phycon.te ))>0. ) 00063 { 00064 /* >>chng 05 oct 02, add this collider 00065 * CHECK - this is deexcitation - is that correct 00066 * this is new H2 - He collision data - use it if positive 00067 * not all states have collision data */ 00068 rate = (realnum)fitted*mole.lgColl_deexec_Calc; 00069 /*if( PRT_COLL ) 00070 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", 00071 ipLo,ipHi,phycon.te,nColl, 00072 iVibHi,iRotHi,iVibLo,iRotLo, 00073 energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo], 00074 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] );*/ 00075 } 00076 00077 /* use g-bar g bar guess of rate coefficient for 00078 * collision with He, this does not change ortho & para 00079 * turn mole.lgColl_gbar on/off with atom h2 gbar on off */ 00080 else if( mole.lgColl_gbar && 00081 (H2_lgOrtho[0][iVibHi][iRotHi]-H2_lgOrtho[0][iVibLo][iRotLo]==0) ) 00082 { 00083 /* the fit is log(K)=y_0+a*((x)^b), where K is the rate coefficient, 00084 * and x is the energy in wavenumbers */ 00085 double ediff = energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo]; 00086 00087 /* do not let energy difference be smaller than 100 wn, the smallest 00088 * difference for which we fit the rate coefficients */ 00089 ediff = MAX2(100., ediff ); 00090 rate = (realnum)pow(10. , 00091 gbarcoll[nColl][0] + gbarcoll[nColl][1] * 00092 pow(ediff,gbarcoll[nColl][2]) )*mole.lgColl_deexec_Calc; 00093 00094 /*if( PRT_COLL ) 00095 fprintf(ioQQQ,"col gbr\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n", 00096 nColl+10, 00097 iVibHi,iRotHi,iVibLo,iRotLo, 00098 energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo], 00099 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] );*/ 00100 } 00101 else 00102 rate = 0; 00103 } 00104 else if( nColl==4 ) 00105 { 00106 /* collisions of H2 with protons - of this group, these are only 00107 * that cause ortho - para conversion */ 00108 /* >>refer H2 coll Hp Gerlich, D., 1990, J. Chem. Phys., 92, 2377-2388 */ 00109 if( CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][1] != 0 ) 00110 { 00111 rate = CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][0] * 1e-10f * 00112 /* sec fit coef was dE in milli eV */ 00113 (realnum)sexp( CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][1]/1000./phycon.te_eV)*mole.lgColl_deexec_Calc; 00114 /*if( PRT_COLL ) 00115 fprintf(ioQQQ,"col fit\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n", 00116 nColl, 00117 iVibHi,iRotHi,iVibLo,iRotLo, 00118 energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo], 00119 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] );*/ 00120 } 00121 /* this is option to use guess of rate coefficient for ortho-para 00122 * conversion by collision with protons */ 00123 /* turn mole.lgColl_gbar on/off with atom h2 gbar on off */ 00124 else if( mole.lgColl_gbar ) 00125 { 00126 /* the fit is log(K)=y_0+a*((x)^b), where K is the rate coefficient, 00127 * and x is the energy in wavenumbers */ 00128 double ediff = energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo]; 00129 ediff = MAX2(100., ediff ); 00130 rate = (realnum)pow(10. , 00131 gbarcoll[nColl][0] + gbarcoll[nColl][1] * 00132 pow(ediff ,gbarcoll[nColl][2]) )*mole.lgColl_deexec_Calc; 00133 00134 /*if( PRT_COLL ) 00135 fprintf(ioQQQ,"col gbr\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n", 00136 nColl+10, 00137 iVibHi,iRotHi,iVibLo,iRotLo, 00138 energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo], 00139 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] );*/ 00140 } 00141 else 00142 rate = 0; 00143 } 00144 else if( CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][0]!= 0 ) 00145 { 00146 /* these are the fits from 00147 *>>refer H2 coll Le Bourlot, J., Pineau des Forets, 00148 *>>refercon G., & Flower, D.R. 1999, MNRAS, 305, 802 00149 * evaluate collision rates for those with real collision data */ 00150 00151 double r = CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][0] + 00152 CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][1]/t3Plus1 + 00153 CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][2]/t3Plus1Squared; 00154 00155 rate = (realnum)pow(10.,r)*mole.lgColl_deexec_Calc; 00156 00157 /*if( PRT_COLL ) 00158 fprintf(ioQQQ,"col fit\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n", 00159 nColl, 00160 iVibHi,iRotHi,iVibLo,iRotLo, 00161 energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo], 00162 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] );*/ 00163 } 00164 /* this is option to use guess of collision rate coefficient - but only if this is 00165 * a downward transition that does not mix ortho and para */ 00166 /* turn mole.lgColl_gbar on/off with atom h2 gbar on off */ 00167 else if( mole.lgColl_gbar && 00168 (H2_lgOrtho[0][iVibHi][iRotHi]-H2_lgOrtho[0][iVibLo][iRotLo]==0) ) 00169 { 00170 /* the fit is log(K)=y_0+a*((x)^b), where K is the rate coefficient, 00171 * and x is the energy in wavenumbers */ 00172 double ediff = energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo]; 00173 /* do not let energy difference be smaller than 100 wn, the smallest 00174 * difference for which we fit the rate coefficients */ 00175 ediff = MAX2(100., ediff ); 00176 rate = (realnum)pow(10. , 00177 gbarcoll[nColl][0] + gbarcoll[nColl][1] * 00178 pow(ediff,gbarcoll[nColl][2]) )*mole.lgColl_deexec_Calc; 00179 /* this is hack to change H2 H collision rates */ 00180 if( nColl == 0 && h2.lgH2_H_coll_07 ) 00181 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] *= 100.; 00182 00183 /*if( PRT_COLL ) 00184 fprintf(ioQQQ,"col gbr\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n", 00185 nColl+10, 00186 iVibHi,iRotHi,iVibLo,iRotLo, 00187 energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo], 00188 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] );*/ 00189 } 00190 else 00191 rate = 0; 00192 00193 00194 /* >>chng 05 feb 09, add option to kill ortho - para collisions */ 00195 if( !mole.lgH2_ortho_para_coll_on && 00196 (H2_lgOrtho[0][iVibHi][iRotHi]-H2_lgOrtho[0][iVibLo][iRotLo]) ) 00197 rate = 0.; 00198 00199 { 00200 enum {DEBUG_LOC=false}; 00201 if( DEBUG_LOC ) 00202 { 00203 fprintf(ioQQQ,"bugcoll\tiVibHi\t%li\tiRotHi\t%li\tiVibLo\t%li\tiRotLo\t%li\tcoll\t%.2e\n", 00204 iVibHi,iRotHi,iVibLo,iRotLo, 00205 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] ); 00206 } 00207 } 00208 return rate; 00209 } 00210 00211 /*H2_CollidRateEvalAll - set H2 collision rate coefficients */ 00212 void H2_CollidRateEvalAll( void ) 00213 { 00214 long int numb_coll_trans = 0; 00215 double excit; 00216 long int iElecHi , iElecLo , ipHi , iVibHi , iRotHi , 00217 ipLo , iVibLo , iRotLo , nColl; 00218 00219 DEBUG_ENTRY( "H2_CollidRateEvalAll()" ); 00220 00221 if( PRT_COLL ) 00222 fprintf(ioQQQ,"H2 coll deex rate coef\n" 00223 "VibHi\tRotHi\tVibLo\tRotLo\trate\n"); 00224 00225 iElecHi = 0; 00226 iElecLo = 0; 00227 if(mole.nH2_TRACE >= mole.nH2_trace_full) 00228 fprintf(ioQQQ,"H2 set collision rates\n"); 00229 /* loop over all possible collisional changes within X 00230 * and set collision rates, which only depend on Te 00231 * will go through array in energy order since coll trans do not 00232 * correspond to a line 00233 * collisional dissociation rate coefficient, units cm3 s-1 */ 00234 H2_coll_dissoc_rate_coef[0][0] = 0.; 00235 H2_coll_dissoc_rate_coef_H2[0][0] = 0.; 00236 for( ipHi=0; ipHi<nLevels_per_elec[0]; ++ipHi ) 00237 { 00238 double energy; 00239 00240 /* obtain the proper indices for the upper level */ 00241 long int ip = H2_ipX_ener_sort[ipHi]; 00242 iVibHi = ipVib_H2_energy_sort[ip]; 00243 iRotHi = ipRot_H2_energy_sort[ip]; 00244 00245 /* this is a guess of the collisional dissociation rate coefficient - 00246 * will be multiplied by the sum of densities of all colliders 00247 * except H2*/ 00248 energy = H2_DissocEnergies[0] - energy_wn[0][iVibHi][iRotHi]; 00249 ASSERT( energy > 0. ); 00250 /* we made this up - Boltzmann factor times rough coefficient */ 00251 H2_coll_dissoc_rate_coef[iVibHi][iRotHi] = 00252 1e-14f * (realnum)sexp(energy/phycon.te_wn) * mole.lgColl_dissoc_coll; 00253 00254 /* collisions with H2 - pre coefficient changed from 1e-8 00255 * (from umist) to 1e-11 as per extensive discussion with Phillip Stancil */ 00256 H2_coll_dissoc_rate_coef_H2[iVibHi][iRotHi] = 00257 1e-11f * (realnum)sexp(energy/phycon.te_wn) * mole.lgColl_dissoc_coll; 00258 00259 /*fprintf(ioQQQ,"DEBUG coll_dissoc_rateee\t%li\t%li\t%.3e\t%.3e\n", 00260 iVibHi,iRotHi, 00261 H2_coll_dissoc_rate_coef[iVibHi][iRotHi], 00262 H2_coll_dissoc_rate_coef_H2[iVibHi][iRotHi]);*/ 00263 00264 for( ipLo=0; ipLo<ipHi; ++ipLo ) 00265 { 00266 ip = H2_ipX_ener_sort[ipLo]; 00267 iVibLo = ipVib_H2_energy_sort[ip]; 00268 iRotLo = ipRot_H2_energy_sort[ip]; 00269 00270 ASSERT( energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo] > 0.); 00271 00272 /* in following the colliders are H, He, H2(ortho), H2(para), and H+ */ 00273 /* fits were read in from the following files: "H2_coll_H.dat" , 00274 * "H2_coll_He.dat" , "H2_coll_H2ortho.dat" ,"H2_coll_H2para.dat", 00275 * "H2_coll_Hp.dat" */ 00276 00277 /* keep track of number of different collision routes */ 00278 ++numb_coll_trans; 00279 /* this is sum over all different colliders, except last two which are special, 00280 * linear rather than log formula for that one for second to last, 00281 * and special fitting formula for last */ 00282 if( PRT_COLL && iVibHi == 1 && iRotHi==3) 00283 fprintf(ioQQQ,"%li\t%li\t%li\t%li", 00284 iVibHi,iRotHi,iVibLo,iRotLo); 00285 for( nColl=0; nColl<N_X_COLLIDER; ++nColl ) 00286 { 00287 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] = 00288 H2_CollidRateEvalOne( iVibHi,iRotHi,iVibLo,iRotLo, 00289 ipHi , ipLo , nColl ); 00290 if( PRT_COLL && iVibHi == 1 && iRotHi==3) 00291 fprintf(ioQQQ,"\t%.2e",H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] ); 00292 } 00293 if( PRT_COLL && iVibHi == 1 && iRotHi==3) 00294 fprintf(ioQQQ,"\n"); 00295 } 00296 } 00297 if( PRT_COLL ) 00298 cdEXIT( EXIT_FAILURE ); 00299 00300 /* at this stage the collision rates that came in from the large data files 00301 * have been entered into the H2_CollRate array. Now add on three extra collision 00302 * terms, the ortho para atomic H collision rates from 00303 * >>>refer H2 collision Sun, Y., & Dalgarno, A., 1994, ApJ, 427, 1053-1056 00304 */ 00305 nColl = 0; 00306 iElecHi = 0; 00307 iElecLo = 0; 00308 iVibHi = 0; 00309 iVibLo = 0; 00310 00311 /* >>chng 02 nov 13, the Sun and Dalgarno rates diverge to + inf below this temp */ 00312 /* >>chng 05 feb 09, do not return zero when T < 100 - instead, don't let T fall below 100 */ 00313 double excit1; 00314 double te_used = MAX2( 100. , phycon.te ); 00315 /* this is the J=1-0 downward collision rate */ 00316 iRotLo = 0; 00317 iRotHi = 1; 00318 excit1 = sexp( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyK/te_used); 00319 excit = sexp( -(POW2(5.30-460./te_used)-21.2) )*1e-13; 00320 00321 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][0] = (realnum)( 00322 excit*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->g/ 00323 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g / 00324 /* >>chng 02 nov 13, from 2nd to first */ 00325 SDIV(excit1) )*mole.lgColl_deexec_Calc * 00326 /* option to disable ortho-para conversion by coll with grains */ 00327 mole.lgH2_ortho_para_coll_on; 00328 00329 /*if( PRT_COLL ) 00330 fprintf(ioQQQ,"col o-p\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n", 00331 nColl, 00332 iVibHi,iRotHi,iVibLo,iRotLo, 00333 energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo], 00334 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] );*/ 00335 00336 /* this is the J=3-0 downward collision rate */ 00337 iRotLo = 0; 00338 iRotHi = 3; 00339 excit = sexp( -(POW2(6.36-373./te_used)-34.5) )*1e-13; 00340 excit1 = sexp( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyK/te_used); 00341 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][0] = (realnum)( 00342 excit*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->g/ 00343 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g / 00344 SDIV(excit1) )*mole.lgColl_deexec_Calc * 00345 /* option to disable ortho-para conversion by coll with grains */ 00346 mole.lgH2_ortho_para_coll_on; 00347 00348 /*if( PRT_COLL ) 00349 fprintf(ioQQQ,"col o-p\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n", 00350 nColl, 00351 iVibHi,iRotHi,iVibLo,iRotLo, 00352 energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo], 00353 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] );*/ 00354 00355 /* this is the downward J=2-1 collision rate */ 00356 iRotLo = 1; 00357 iRotHi = 2; 00358 excit = sexp( -(POW2(5.35-454./te_used)-23.1 ) )*1e-13; 00359 excit1 = sexp( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyK/te_used); 00360 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][0] = (realnum)( 00361 excit*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->g/ 00362 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g / 00363 SDIV(excit1) )*mole.lgColl_deexec_Calc * 00364 /* option to disable ortho-para conversion by coll with grains */ 00365 mole.lgH2_ortho_para_coll_on; 00366 00367 /* >>chng 05 nov 30, GS, rates decreases exponentially for low temperature, see Le Bourlot et al. 1999 */ 00368 /* Phillips mail--Apparently, the SD fit is only valid over the range of their calculations, 100-1000K. 00369 * The rate should continue to fall exponentially with decreasing T, something like exp(-3900/T) for 0->1 and 00370 * exp[-(3900-170.5)/T] for 1->0. It is definitely, not constant for T lower than 100 K, as far as we know. 00371 * There may actually be a quantum tunneling effect which causes the rate to increase at lower T, but no 00372 * one has calculated it (as far as I know) and it might happen below 1K or so.???*/ 00373 if( phycon.te < 100. ) 00374 { 00375 /* first term in exp is suggested by Phillip, second temps in paren is to ensure continuity 00376 * across 100K */ 00377 H2_CollRate[0][1][0][0][0] = (realnum)(H2_CollRate[0][0][1][0][0]*exp(-(3900-170.5)*(1./phycon.te - 1./100.))); 00378 H2_CollRate[0][3][0][0][0] = (realnum)(H2_CollRate[0][0][3][0][0]*exp(-(3900-1015.1)*(1./phycon.te - 1./100.))); 00379 H2_CollRate[0][2][0][1][0] = (realnum)(H2_CollRate[0][0][2][0][1]*exp(-(3900-339.3)*(1./phycon.te - 1./100.))); 00380 } 00381 00382 /*if( PRT_COLL ) 00383 fprintf(ioQQQ,"col o-p\t%li\t%li\t%li\t%li\t%li\t%.4e\t%.4e\n", 00384 nColl, 00385 iVibHi,iRotHi,iVibLo,iRotLo, 00386 energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo], 00387 H2_CollRate[iVibHi][iRotHi][iVibLo][iRotLo][nColl] );*/ 00388 00389 if( mole.nH2_TRACE >= mole.nH2_trace_full ) 00390 fprintf(ioQQQ, 00391 " collision rates updated for new temp, number of trans is %li\n", 00392 numb_coll_trans); 00393 00394 /* quit it we are only printing - but do this after printing coll rates 00395 if( PRT_COLL ) 00396 cdEXIT( EXIT_FAILURE ); */ 00397 return; 00398 } 00399 00400 /*H2_CollidRateRead read collision rates */ 00401 void H2_CollidRateRead( long int nColl ) 00402 { 00403 /* the colliders are H, He, H2 ortho, H2 para, H+ 00404 * these are the default file names. they can be overridden with 00405 * the SET ATOMIC DATA MOLECULE H2 command */ 00406 00407 /*NB thesee must be kept parallel with labels in chH2ColliderLabels */ 00408 00409 const char* cdDATAFILE[N_X_COLLIDER] = 00410 { 00411 /* 0 */"H2_coll_H_07.dat", 00412 /* 1 */"H2_coll_He_LeBourlot.dat", 00413 /* 2 */"H2_coll_H2ortho.dat", 00414 /* 3 */"H2_coll_H2para.dat", 00415 /* 4 */"H2_coll_Hp.dat" 00416 }; 00417 00418 FILE *ioDATA; 00419 char chLine[FILENAME_PATH_LENGTH_2]; 00420 const char* chFilename; 00421 long int i, n1; 00422 long int iVibHi , iVibLo , iRotHi , iRotLo; 00423 long int magic_expect = -1; 00424 00425 DEBUG_ENTRY( "H2_CollidRateRead()" ); 00426 00427 if( nColl == 0 ) 00428 { 00429 /* which H2 - H data set? */ 00430 if( h2.lgH2_H_coll_07 ) 00431 { 00432 magic_expect = 71106; 00433 chFilename = "H2_coll_H_07.dat"; 00434 } 00435 else 00436 { 00437 /* the 1999 data set */ 00438 magic_expect = 20429; 00439 chFilename = "H2_coll_H_99.dat"; 00440 } 00441 } 00442 else if( nColl == 1 ) 00443 { 00444 /* which H2 - He data set? */ 00445 if( mole.lgH2_He_ORNL ) 00446 { 00447 /* >>chng 07 may 12, update magic number when one transition in H2 - H2 file was 00448 * replaced - the fitting coefficients had terms of order -8e138 which fpe in float */ 00449 long int magic_found, 00450 magic_expect = 70513; 00451 /* special case, new data file from Oak Ridge project - 00452 * call init routine and return - data is always read in when large H2 is included - 00453 * but data are only used (for now, mid 2005) when command 00454 * atom H2 He OLD (Meudon) NEW (Stancil) and OFF given */ 00455 /*H2_He_coll_init receives the name of the file that contains the fitting coefficients 00456 * of all transitions and read into 3d vectors. It outputs 'test.out' to test the arrays*/ 00457 chFilename = "H2_coll_He_ORNL.dat"; 00458 if( (magic_found = H2_He_coll_init( chFilename )) != magic_expect ) 00459 { 00460 fprintf(ioQQQ,"The H2 - He collision data file H2_coll_He_ORNL.dat does not have the correct magic number.\n"); 00461 fprintf(ioQQQ,"I found %li but expected %li\n", magic_found , magic_expect ); 00462 00463 cdEXIT(EXIT_FAILURE); 00464 } 00465 return; 00466 } 00467 else 00468 { 00469 magic_expect = 20429; 00470 /* the Le Bourlot et al data */ 00471 chFilename = "H2_coll_He_LeBourlot.dat"; 00472 } 00473 } 00474 else 00475 { 00476 /* files with only one version */ 00477 magic_expect = 20429; 00478 chFilename = cdDATAFILE[nColl]; 00479 } 00480 00481 ioDATA = open_data( chFilename, "r" ); 00482 00483 /* read the first line and check that magic number is ok */ 00484 if( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) == NULL ) 00485 { 00486 fprintf( ioQQQ, " H2_CollidRateRead could not read first line of %s\n", chFilename ); 00487 cdEXIT(EXIT_FAILURE); 00488 } 00489 00490 /* level 1 magic number */ 00491 n1 = atoi( chLine ); 00492 00493 /* magic number 00494 * the following is the set of numbers that appear at the start of level1.dat 01 08 10 */ 00495 if( n1 != magic_expect ) 00496 { 00497 fprintf( ioQQQ, 00498 " H2_CollidRateRead: the version of %s is not the current version.\n", chFilename ); 00499 fprintf( ioQQQ, 00500 " I expected to find the number %li and got %li instead.\n" , 00501 magic_expect , n1 ); 00502 fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine ); 00503 cdEXIT(EXIT_FAILURE); 00504 } 00505 00506 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL ) 00507 { 00508 if( chLine[0]=='#' ) 00509 continue; 00510 double a[3]; 00511 sscanf(chLine,"%li\t%li\t%li\t%li\t%le\t%le\t%le", 00512 &iVibHi ,&iRotHi , &iVibLo , &iRotLo , &a[0],&a[1],&a[2] ); 00513 /*fprintf(ioQQQ,"DEBUG %li %li %li %li \n", 00514 iVibHi , iRotHi , iVibLo , iRotLo);*/ 00515 ASSERT( iRotHi >= 0 && iVibHi >= 0 && iRotLo >= 0 && iVibLo >=0 ); 00516 00517 /* check that we actually included the levels in the model representation 00518 * depending on the potential surface, the collision date set 00519 * may not agree with our adopted model - skip those */ 00520 if( iVibHi > VIB_COLLID || 00521 iVibLo > VIB_COLLID || 00522 iRotHi > h2.nRot_hi[0][iVibHi] || 00523 iRotLo > h2.nRot_hi[0][iVibLo]) 00524 { 00525 iVibHi = -1; 00526 continue; 00527 } 00528 00529 /* some collision rates have the same upper and lower indices - skip them */ 00530 if( !( (iVibHi == iVibLo) && (iRotHi == iRotLo )) ) 00531 { 00532 /* this is downward transition - make sure that the energy difference is positive */ 00533 ASSERT( (energy_wn[0][iVibHi][iRotHi] - energy_wn[0][iVibLo][iRotLo] ) > 0. ); 00534 for( i=0; i<3; ++i ) 00535 { 00536 CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][i] = (realnum)a[i]; 00537 } 00538 00539 /* this prints all levels with rates 00540 fprintf(ioQQQ,"no\t%li\t%li\t%li\t%li\t%.2e\t%.2e\t%.2e\n", 00541 iVibHi,iRotHi,iVibLo,iRotLo,a[0],a[1],a[2]);*/ 00542 } 00543 } 00544 fclose( ioDATA ); 00545 00546 return; 00547 } 00548 /* end of H2_CollidRateRead */ 00549 00550 /* This code was written by Terry Yun, 2005 */ 00551 /* H2_He_coll_init receives the name of the file that contains the fitting 00552 * coefficients of all transitions and read into 3d vectors. 00553 * It returns magic number*/ 00554 long int H2_He_coll_init(const char FILE_NAME_IN[]) 00555 { 00556 00557 int i, j; 00558 long int magic; 00559 double par[N_H2_HE_FIT_PAR]; 00560 char line[INPUT_LINE_LENGTH]; 00561 int h2_i, h2_f, he_i, he_f;/* target, projectile initial and final indices */ 00562 char quality, space; 00563 /*'space' variable is for reading spaces between characters */ 00564 /* scanf can skip spaces(or tap) when it reads numbers, but not for characters. */ 00565 double error; 00566 00567 FILE *ifp; 00568 00569 DEBUG_ENTRY( "H2_He_coll_init()" ); 00570 00571 /* create space for two arrays - H2_He_coll_fit_par will contain the 00572 * fit parameters in form [init state][final state][param] 00573 * lgDefn_H2He_coll is flag saying whether data exists */ 00574 H2_He_coll_fit_par = (realnum ***)MALLOC(sizeof(realnum **)*(unsigned)nLevels_per_elec[0] ); 00575 lgDefn_H2He_coll = (bool**)MALLOC(sizeof(bool *)*(unsigned)nLevels_per_elec[0] ); 00576 for( i=1; i<nLevels_per_elec[0]; ++i ) 00577 { 00578 lgDefn_H2He_coll[i] = (bool*)MALLOC(sizeof(bool)*(unsigned)i ); 00579 H2_He_coll_fit_par[i] = (realnum **)MALLOC(sizeof(realnum *)*(unsigned)i ); 00580 for( j=0; j<i; ++j) 00581 H2_He_coll_fit_par[i][j] = (realnum *)MALLOC(sizeof(realnum)*(unsigned)N_H2_HE_FIT_PAR ); 00582 } 00583 00584 /* set a flag saying whether data exits: initially everything is '0' */ 00585 for( i=1; i<nLevels_per_elec[0]; i++ ) 00586 { 00587 for( j=0; j<i; j++ ) 00588 { 00589 lgDefn_H2He_coll[i][j] = false; 00590 } 00591 } 00592 00593 /*open the input file and put into 3d arrays*/ 00594 ifp = open_data( FILE_NAME_IN, "r" ); 00595 00596 /*read the magic number*/ 00597 if( read_whole_line(line, (int)sizeof(line), ifp)==NULL ) 00598 { 00599 printf("DISASTER H2_He_coll_init can't read first line of %s\n", FILE_NAME_IN); 00600 cdEXIT(EXIT_FAILURE); 00601 } 00602 sscanf(line, "%li", &magic); 00603 00604 /* read the file until the end of line */ 00605 while( read_whole_line(line, (int)sizeof(line), ifp) != NULL ) 00606 { 00607 /* skip any line that starts with '#' */ 00608 if( line[0]!='#' ) 00609 { 00610 sscanf(line, "%i%i%i%i%c%c%c%c%lf%lf%lf%lf%lf%lf%lf%lf%lf", &h2_i, 00611 &h2_f, &he_i, &he_f,&space, &space, &space, &quality, 00612 &error, &par[0], &par[1], &par[2], &par[3], &par[4], 00613 &par[5], &par[6], &par[7]); 00614 00615 /* >>chng 07 may 28, extensive testing showed rate > 1e38 for 00616 * temperatures where fitting equation hit a pole. Exchange 00617 * with Phillip Stancil: 00618 "You are correct. The problem is the parameter d which is 00619 negative in all the problem cases you found. A pole occurs 00620 when d*T/1000 + 1 =0. 00621 00622 I found that there are a total of 11 transitions with d<1. 00623 A quick solution is to take abs(d). I checked one case and the 00624 resulting function went smoothly through the pole. On either 00625 side of the pole the two functions looked the same on a log-log plot. 00626 00627 LeBourlot used a similar function, but had d=1. So, they never 00628 got a pole. I guess I tried to make the function a little too 00629 flexible. 00630 00631 There is a similar term with g*T/1000 + 1. There are no fits 00632 with a negative g, but you might take abs(g) to be on the space 00633 side for the future." 00634 */ 00635 if( par[3] < 0. ) 00636 { 00637 /*fprintf(ioQQQ,"DEBUG negative [3]=%.2e\n", par[3]);*/ 00638 par[3] = -par[3]; 00639 } 00640 if( par[6] < 0. ) 00641 { 00642 /*fprintf(ioQQQ,"DEBUG negative [6]=%.2e\n", par[6]);*/ 00643 par[6] = -par[6]; 00644 } 00645 /* input file is on physics or Fortran scale so lowest energy 00646 * index is 1. Store on C scale */ 00647 --h2_f; 00648 --h2_i; 00649 /* set true when the indices are defined */ 00650 lgDefn_H2He_coll[h2_i][h2_f] = true; 00651 for( i=0; i<N_H2_HE_FIT_PAR; i++ ) 00652 /* assigning the parameters to 3d array */ 00653 H2_He_coll_fit_par[h2_i][h2_f][i] = (realnum)par[i]; 00654 } 00655 } 00656 fclose(ifp); 00657 return magic; 00658 } 00659 00660 /* This code was written by Terry Yun, 2005 */ 00661 /* H2_He_coll Interpolate the rate coefficients 00662 * It receives initial and final transition and temperature 00663 * The range of the temperature is between 2K - 1e8K */ 00664 double H2_He_coll(int init, int final, double temp) 00665 { 00666 00667 double k, t3Plus1, b2, c2, f2, t3; 00668 00669 DEBUG_ENTRY( "H2_He_coll()" ); 00670 00671 /* Invalid entries returns '-1':the initial indices are smaller than the final indices */ 00672 if( temp<2 || temp > 1e8 ) 00673 { 00674 k = -1; 00675 } 00676 else if( init <= final ) 00677 { 00678 k = -1; 00679 } 00680 /* Invalid returns '-1': the indices are greater than 302 or smaller than 0 */ 00681 else if( init < 0 || init >302 || final < 0 || final > 302 ) 00682 { 00683 k = -1; 00684 } 00685 /* Undefined indices returns '0' */ 00686 else if( !lgDefn_H2He_coll[init][final] ) 00687 { 00688 k = -1; 00689 } 00690 /* defined indices */ 00691 else if( lgDefn_H2He_coll[init][final] ) 00692 { 00693 /* the fitting equation we used: 00694 * k_(vj,v'j') = 10^(a + b/(d*T/10^3+1) + c/t^2) 00695 * + 10^(e + f/(g*T/10^3+1)**h) 00696 * T in K, k in cm3/s. */ 00697 double sum1 , sum2; 00698 00699 /* this kludge is to bypass errors in fitting formula - there 00700 * is a pole possible at around 1e6 K and rates can diverge 00701 * and high temperatures */ 00703 temp = MIN2( 1e4 , temp ); 00704 00705 t3 = temp/1e3; 00706 /* t3Plus1 = T/1000 + 1*/ 00707 t3Plus1 = t3+1; 00708 b2 = H2_He_coll_fit_par[init][final][1]/(H2_He_coll_fit_par[init][final][3]*t3+1); 00709 c2 = H2_He_coll_fit_par[init][final][2]/(t3Plus1*t3Plus1); 00710 { 00711 enum {DEBUG_LOC=false}; 00712 if( DEBUG_LOC ) 00713 { 00714 fprintf(ioQQQ,"bug H2 He coll\t%i %i %.3e %.3e %.3e \n", 00715 init,final, 00716 H2_He_coll_fit_par[init][final][5], 00717 H2_He_coll_fit_par[init][final][6]*t3+1, 00718 H2_He_coll_fit_par[init][final][7] 00719 /*pow(H2_He_coll_fit_par[init][final][6]*t3+1,H2_He_coll_fit_par[init][final][7])*/ 00720 ); 00721 } 00722 } 00723 /* this is log of f2 - see whether it is within bounds */ 00724 sum1 = H2_He_coll_fit_par[init][final][7] * log10( H2_He_coll_fit_par[init][final][6]*t3+1. ); 00725 /* this protects against overflow */ 00726 if( fabs(sum1)< 38. ) 00727 { 00728 /* >>chng 06 jun 07, Mitchell Martin, from following to equivalent below. This was basically a 00729 * bug in the pgcc compiler that caused h2_pdr_leiden_f1 to throw an fpe 00730 * the changed code is equivalent */ 00731 /*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]);*/ 00732 f2 = H2_He_coll_fit_par[init][final][5]/pow(10. , sum1 ); 00733 } 00734 else 00735 f2= 0.; 00736 { 00737 enum {DEBUG_LOC=false }; 00738 if( DEBUG_LOC ) 00739 { 00740 fprintf(ioQQQ,"bug H2 He coll\t%i %i %.3e %.3e %.3e %.3e %.3e sum %.3e %.3e \n", 00741 init,final, 00742 H2_He_coll_fit_par[init][final][0], 00743 b2, 00744 c2, 00745 H2_He_coll_fit_par[init][final][4], 00746 f2 , 00747 H2_He_coll_fit_par[init][final][0]+b2+c2, 00748 H2_He_coll_fit_par[init][final][4]+f2); 00749 } 00750 } 00751 00752 sum1 = H2_He_coll_fit_par[init][final][0]+b2+c2; 00753 sum2 = H2_He_coll_fit_par[init][final][4]+f2; 00754 k = 0.; 00755 if( fabs(sum1) < 38. ) 00756 k += pow(10., sum1 ); 00757 if( fabs(sum2) < 38. ) 00758 k += pow(10., sum2 ); 00759 00760 if( k>1e-6 ) 00761 { 00762 /* rate of 1e-6 is largest possible rate according to Phillip 00763 * Stancil email 07 may 27 */ 00764 fprintf(ioQQQ,"PROBLEM H2-He rate coefficient hit cap, upper index of %i" 00765 " lower index of %i rate was %.2e resetting to 1e-6, Te=%e\n", 00766 init , final , k , phycon.te ); 00767 k = 1e-6; 00768 } 00769 { 00770 enum {DEBUG_LOC=false}; 00771 if( DEBUG_LOC ) 00772 { 00773 fprintf(ioQQQ,"DEBUG H2 He rate %.3e \n", 00774 k); 00775 } 00776 } 00777 } 00778 /*unknown invalid entries returns '99'*/ 00779 else 00780 { 00781 k = 99; 00782 } 00783 return k; 00784 }