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_Create create variables for the H2 molecule, called by ContCreatePointers after continuum 00004 * mesh has been set up */ 00005 #include "cddefines.h" 00006 #include "physconst.h" 00007 #include "mole.h" 00008 #include "taulines.h" 00009 #include "lines_service.h" 00010 #include "opacity.h" 00011 #include "hmi.h" 00012 #include "ipoint.h" 00013 #include "grainvar.h" 00014 #include "h2.h" 00015 #include "h2_priv.h" 00016 00017 /* if this is set true then code will print energies and stop */ 00018 /*@-redef@*/ 00019 enum {DEBUG_ENER=false}; 00020 /*@+redef@*/ 00021 00022 /* this is equation 8 of Takahashi 2001, clearer definition is given in 00023 * equation 5 and following discussion of 00024 * >>refer H2 formation Takahashi, J., & Uehara, H., 2001, ApJ, 561, 843-857 00025 * 0.27eV, convert into wavenumbers */ 00026 static double XVIB[H2_TOP] = { 0.70 , 0.60 , 0.20 }; 00027 static double Xdust[H2_TOP] = { 0.04 , 0.10 , 0.40 }; 00028 00029 /* this is energy difference between bottom of potential well and 0,0 00030 * the Takahashi energy scale is from the bottom, 00031 * 2201.9 wavenumbers */ 00032 static const double energy_off = 0.273*FREQ_1EV/SPEEDLIGHT; 00033 00034 STATIC double EH2_eval( long int iVib , int ipH2 ) 00035 { 00036 double EH2_here; 00037 double Evm = H2_DissocEnergies[0]* XVIB[ipH2] + energy_off; 00038 00039 double Ev = (energy_wn[0][iVib][0]+energy_off); 00040 /* equation 9 of Takahashi 2001 which is only an approximation 00041 double EH2 = H2_DissocEnergies[0] * (1. - Xdust[ipH2] ); */ 00042 /* equation 1, 2 of 00043 * Takahashi, Junko, & Uehara, Hideya, 2001, ApJ, 561, 843-857, 00044 * this is heat deposited on grain by H2 formation in this state */ 00045 double Edust = H2_DissocEnergies[0] * Xdust[ipH2] * 00046 ( 1. - ( (Ev - Evm) / (H2_DissocEnergies[0]+energy_off-Evm)) * 00047 ( (1.-Xdust[ipH2])/2.) ); 00048 ASSERT( Edust >= 0. ); 00049 00050 /* energy is total binding energy less energy lost on grain surface 00051 * and energy offset */ 00052 EH2_here = H2_DissocEnergies[0]+energy_off - Edust; 00053 ASSERT( EH2_here >= 0.); 00054 00055 return EH2_here; 00056 } 00057 00058 /*H2_vib_dist evaluates the vibration distribution for H2 formed on grains */ 00059 STATIC double H2_vib_dist( long int iVib , int ipH2 , double EH2) 00060 { 00061 double G1[H2_TOP] = { 0.3 , 0.4 , 0.9 }; 00062 double G2[H2_TOP] = { 0.6 , 0.6 , 0.4 }; 00063 double Evm = H2_DissocEnergies[0]* XVIB[ipH2] + energy_off; 00064 double Fv; 00065 if( (energy_wn[0][iVib][0]+energy_off) <= Evm ) 00066 { 00067 /* equation 4 of Takahashi 2001 */ 00068 Fv = sexp( POW2( (energy_wn[0][iVib][0]+energy_off - Evm)/(G1[ipH2]* Evm ) ) ); 00069 } 00070 else 00071 { 00072 /* equation 5 of Takahashi 2001 */ 00073 Fv = sexp( POW2( (energy_wn[0][iVib][0]+energy_off - Evm)/(G2[ipH2]*(EH2 - Evm ) ) ) ); 00074 } 00075 return Fv; 00076 } 00077 00078 00079 /*H2_Create create variables for the H2 molecule, called by 00080 * ContCreatePointers after continuum mesh has been set up */ 00081 void H2_Create(void) 00082 { 00083 long int i , iElecHi , iElecLo; 00084 long int iVibHi , iVibLo; 00085 long int iRotHi , iRotLo; 00086 long int iElec, iVib , iRot; 00087 long int nColl, 00088 nlines; 00089 int ier; 00090 int nEner; 00091 /* >>chng 03 nov 26, from enum H2_type to int */ 00092 int ipH2; 00093 realnum sum , sumj , sumv , sumo , sump; 00094 00095 /* this is flag set above - when true h2 code is not executed - this is way to 00096 * avoid this code when it is not working */ 00097 /* only malloc vectors one time per core load */ 00098 if( lgH2_READ_DATA || !h2.lgH2ON ) 00099 return; 00100 00101 DEBUG_ENTRY( "H2_Create()" ); 00102 00103 /* print string if H2 debugging is enabled */ 00104 if( mole.nH2_TRACE ) 00105 fprintf(ioQQQ," H2_Create called in DEBUG mode.\n"); 00106 00107 /* this was option to print all electronic states in the main printout - but 00108 * number of electronic states was not known at initialization so set to -1, 00109 * will set properly now */ 00110 if( h2.nElecLevelOutput < 1 ) 00111 h2.nElecLevelOutput = mole.n_h2_elec_states; 00112 00113 /* this var is in h2.h and prevents h2 from being change once committed here */ 00114 lgH2_READ_DATA = true; 00115 00116 /* create special vector that saves collision rates within ground */ 00117 /* this will contain a vector for collisions within the X ground electronic state, 00118 * CollRateFit[vib_up][rot_up][vib_lo][rot_lo][coll_type][3] */ 00119 /* N_X_COLLIDER is number of different species that collide */ 00120 iElecHi = 0; 00121 /* the current data set is limited to vib hi <= 3 */ 00122 /* will define collision rates for all possible transitions within X */ 00123 CollRateFit.reserve(h2.nVib_hi[iElecHi]+1); 00124 H2_CollRate.reserve(h2.nVib_hi[iElecHi]+1); 00125 for( iVibHi = 0; iVibHi <= h2.nVib_hi[iElecHi]; ++iVibHi ) 00126 { 00127 CollRateFit.reserve(iVibHi,h2.nRot_hi[iElecHi][iVibHi]+1); 00128 H2_CollRate.reserve(iVibHi,h2.nRot_hi[iElecHi][iVibHi]+1); 00129 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 00130 { 00131 CollRateFit.reserve(iVibHi,iRotHi,h2.nVib_hi[iElecHi]+1); 00132 H2_CollRate.reserve(iVibHi,iRotHi,h2.nVib_hi[iElecHi]+1); 00133 for( iVibLo=0; iVibLo<(h2.nVib_hi[iElecHi]+1); ++iVibLo ) 00134 { 00135 CollRateFit.reserve(iVibHi,iRotHi,iVibLo,h2.nRot_hi[iElecHi][iVibLo]+1); 00136 H2_CollRate.reserve(iVibHi,iRotHi,iVibLo,h2.nRot_hi[iElecHi][iVibLo]+1); 00137 for( iRotLo=0; iRotLo<=h2.nRot_hi[iElecHi][iVibLo]; ++iRotLo ) 00138 { 00139 H2_CollRate.reserve(iVibHi,iRotHi,iVibLo,iRotLo,N_X_COLLIDER); 00140 CollRateFit.reserve(iVibHi,iRotHi,iVibLo,iRotLo,N_X_COLLIDER); 00141 for( nColl=0; nColl<N_X_COLLIDER; ++nColl ) 00142 { 00143 /* the last one - the three coefficients */ 00144 CollRateFit.reserve(iVibHi,iRotHi,iVibLo,iRotLo,nColl,3); 00145 } 00146 } 00147 } 00148 } 00149 } 00150 00151 CollRateFit.alloc(); 00152 H2_CollRate.alloc(); 00153 00154 /* zero out the collisional rates since only a minority of them are known*/ 00155 CollRateFit.zero(); 00156 H2_CollRate.zero(); 00157 00158 /* create space for the electronic levels */ 00159 H2_populations.reserve(mole.n_h2_elec_states); 00160 pops_per_vib.reserve(mole.n_h2_elec_states); 00161 H2_dissprob.reserve(mole.n_h2_elec_states); 00162 00163 for( iElec = 0; iElec<mole.n_h2_elec_states; ++iElec ) 00164 { 00165 00166 if( mole.nH2_TRACE >= mole.nH2_trace_full) 00167 fprintf(ioQQQ,"elec %li highest vib= %li\n", iElec , h2.nVib_hi[iElec] ); 00168 00169 ASSERT( h2.nVib_hi[iElec] > 0 ); 00170 00171 /* h2.nVib_hi is now the highest vibrational level before dissociation, 00172 * now allocate space to hold the number of rotation levels */ 00173 H2_populations.reserve(iElec,h2.nVib_hi[iElec]+1); 00174 pops_per_vib.reserve(iElec,h2.nVib_hi[iElec]+1); 00175 if( iElec > 0 ) 00176 H2_dissprob.reserve(iElec,h2.nVib_hi[iElec]+1); 00177 00178 /* now loop over all vibrational levels, and find out how many rotation levels there are */ 00179 /* ground is special since use tabulated data - there are 14 vib states, 00180 * ivib=14 is highest */ 00181 for( iVib = 0; iVib <= h2.nVib_hi[iElec]; ++iVib ) 00182 { 00183 /* lastly create the space for the rotation quantum number */ 00184 H2_populations.reserve(iElec,iVib,h2.nRot_hi[iElec][iVib]+1); 00185 if( iElec > 0 ) 00186 H2_dissprob.reserve(iElec,iVib,h2.nRot_hi[iElec][iVib]+1); 00187 } 00188 } 00189 00190 H2_populations.alloc(); 00191 H2_populations_LTE.alloc( H2_populations.clone() ); 00192 H2_old_populations.alloc( H2_populations.clone() ); 00193 H2_Boltzmann.alloc( H2_populations.clone() ); 00194 H2_stat.alloc( H2_populations.clone() ); 00195 energy_wn.alloc( H2_populations.clone() ); 00196 H2_rad_rate_out.alloc( H2_populations.clone() ); 00197 H2_lgOrtho.alloc( H2_populations.clone() ); 00198 00199 pops_per_vib.alloc(); 00200 00201 H2_dissprob.alloc(); 00202 H2_disske.alloc( H2_dissprob.clone() ); 00203 00204 /* set this one time, will never be set again, but might be printed */ 00205 H2_rad_rate_out.zero(); 00206 00207 /* these do not have electronic levels - all within X */ 00208 H2_ipPhoto.reserve(h2.nVib_hi[0]+1); 00209 00210 /* space for the vibration levels */ 00211 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib ) 00212 { 00213 /* space for the rotation quantum number */ 00214 H2_ipPhoto.reserve(iVib,h2.nRot_hi[0][iVib]+1); 00215 } 00216 00217 H2_ipPhoto.alloc(); 00218 H2_col_rate_in.alloc( H2_ipPhoto.clone() ); 00219 H2_col_rate_out.alloc( H2_ipPhoto.clone() ); 00220 H2_rad_rate_in.alloc( H2_ipPhoto.clone() ); 00221 H2_coll_dissoc_rate_coef.alloc( H2_ipPhoto.clone() ); 00222 H2_coll_dissoc_rate_coef_H2.alloc( H2_ipPhoto.clone() ); 00223 H2_X_colden.alloc( H2_ipPhoto.clone() ); 00224 H2_X_rate_from_elec_excited.alloc( H2_ipPhoto.clone() ); 00225 H2_X_rate_to_elec_excited.alloc( H2_ipPhoto.clone() ); 00226 H2_X_colden_LTE.alloc( H2_ipPhoto.clone() ); 00227 H2_X_formation.alloc( H2_ipPhoto.clone() ); 00228 H2_X_Hmin_back.alloc( H2_ipPhoto.clone() ); 00229 00230 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib ) 00231 { 00232 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot ) 00233 { 00234 /* >>chng 04 jun 14, set these to bad numbers */ 00235 H2_rad_rate_in[iVib][iRot] = -BIGFLOAT; 00236 H2_coll_dissoc_rate_coef[iVib][iRot] = -BIGFLOAT; 00237 H2_coll_dissoc_rate_coef_H2[iVib][iRot] = -BIGFLOAT; 00238 } 00239 } 00240 /* zero out the matrices */ 00241 H2_X_colden.zero(); 00242 H2_X_colden_LTE.zero(); 00243 H2_X_formation.zero(); 00244 H2_X_Hmin_back.zero(); 00245 /* rates [cm-3 s-1] from elec excited states into X only vib and rot */ 00246 H2_X_rate_from_elec_excited.zero(); 00247 /* rates [s-1] to elec excited states from X only vib and rot */ 00248 H2_X_rate_to_elec_excited.zero(); 00249 00250 /* distribution function for populations following formation from H minus H- */ 00251 H2_X_hminus_formation_distribution.reserve(nTE_HMINUS); 00252 for( i=0; i<nTE_HMINUS; ++i ) 00253 { 00254 H2_X_hminus_formation_distribution.reserve(i,h2.nVib_hi[0]+1); 00255 /* space for the vibration levels */ 00256 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib ) 00257 { 00258 H2_X_hminus_formation_distribution.reserve(i,iVib,h2.nRot_hi[0][iVib]+1); 00259 } 00260 } 00261 H2_X_hminus_formation_distribution.alloc(); 00262 H2_X_hminus_formation_distribution.zero(); 00263 H2_Read_hminus_distribution(); 00264 00265 /* >>chng 05 jun 20, do not use this, which is highly processed - use ab initio 00266 * rates of excitation to electronic levels instead */ 00267 /* read in cosmic ray distribution information 00268 H2_Read_Cosmicray_distribution(); */ 00269 00270 /* grain formation matrix */ 00271 H2_X_grain_formation_distribution.reserve(H2_TOP); 00272 for( ipH2=0; ipH2<(int)H2_TOP; ++ipH2 ) 00273 { 00274 H2_X_grain_formation_distribution.reserve(ipH2,h2.nVib_hi[0]+1); 00275 00276 /* space for the vibration levels */ 00277 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib ) 00278 { 00279 H2_X_grain_formation_distribution.reserve(ipH2,iVib,h2.nRot_hi[0][iVib]+1); 00280 } 00281 } 00282 H2_X_grain_formation_distribution.alloc(); 00283 H2_X_grain_formation_distribution.zero(); 00284 00285 /* space for the energy vector is now malloced, must fill it in, 00286 * defines array energy_wn[nelec][iVib][iRot] */ 00287 for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec ) 00288 { 00289 /* get energies out of files into array energy_wn[nelec][iVib][iRot] */ 00290 H2_ReadEnergies(iElec); 00291 00292 /* get dissociation probabilities and energies - ground state is stable */ 00293 if( iElec > 0 ) 00294 H2_ReadDissprob(iElec); 00295 } 00296 00297 /* >>02 oct 18, add photodissociation, H2 + hnu => 2H + KE */ 00298 /* we now have ro-vib energies, now set up threshold array offsets 00299 * for photodissociation */ 00300 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib ) 00301 { 00302 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot ) 00303 { 00304 /* this is energy needed to get up to n=3 electronic continuum 00305 * H2 cannot dissociate following absorption of a continuum photon into the 00306 * continuum above X, which would require little energy, and would be given 00307 * by H2_DissocEnergies[0], because that process violates momentum conservation 00308 * these would be the triplet states - permitted are into singlets 00309 * the effective full wavelength range of this process is from Lya to 00310 * Lyman limit in shielded regions 00311 * tests show limits are between 850A and 1220A - so Lya is included */ 00313 /*>>KEYWORD Allison & Dalgarno; continuum dissociation; */ 00314 double thresh = (H2_DissocEnergies[1] - energy_wn[0][iVib][iRot])*WAVNRYD; 00315 /*fprintf(ioQQQ,"DEBUG\t%.2f\t%f\n", RYDLAM/thresh , thresh);*/ 00316 /* in theory we should be able to assert that thesh just barely reaches 00317 * lya, but actual numbers reach down to 0.749 ryd */ 00318 ASSERT( thresh > 0.74 ); 00319 H2_ipPhoto[iVib][iRot] = ipoint(thresh); 00320 } 00321 } 00322 00323 nH2_energies = 0; 00324 for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec) 00325 { 00326 /* the number of levels within the molecule */ 00327 nH2_energies += nLevels_per_elec[iElec]; 00328 } 00329 00330 if( mole.nH2_TRACE >= mole.nH2_trace_full ) 00331 { 00332 for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec) 00333 { 00334 /* print the number of levels within iElec */ 00335 fprintf(ioQQQ,"\t(%li %li)", iElec , nLevels_per_elec[iElec] ); 00336 } 00337 fprintf(ioQQQ, 00338 " H2_Create: there are %li electronic levels, in each level there are", 00339 mole.n_h2_elec_states); 00340 fprintf(ioQQQ, 00341 " for a total of %li levels.\n", nH2_energies ); 00342 } 00343 00344 /* now read in the various sets of collision data */ 00345 for( nColl=0; nColl<N_X_COLLIDER; ++nColl ) 00346 { 00347 /* ground state has tabulated data */ 00348 H2_CollidRateRead(nColl); 00349 } 00350 /* option to add gaussian random mole */ 00351 if( mole.lgH2_NOISE ) 00352 { 00353 iElecHi = 0; 00354 /* loop over all transitions */ 00355 for( iVibHi = 0; iVibHi <= VIB_COLLID; ++iVibHi ) 00356 { 00357 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 00358 { 00359 for( iVibLo=0; iVibLo<(VIB_COLLID+1); ++iVibLo ) 00360 { 00361 for( iRotLo=0; iRotLo<=h2.nRot_hi[iElecHi][iVibLo]; ++iRotLo ) 00362 { 00363 /* first set of expressions are series that adds to log of rate, 00364 * so we will add the gaussian noise */ 00365 /* >>chng 05 dec 13, GS, last two fits are different, 00366 * loop had been to N_X_COLLIDER-1 and so included Stancil data 00367 * with noise became negative */ 00368 for( nColl=0; nColl<N_X_COLLIDER-2; ++nColl ) 00369 { 00370 /* the gaussian random number, many possible collision rates 00371 * have no data, and CollRateFit[][][][][][0] is zero - do not 00372 * scramble these, only scramble the non-zero rates */ 00373 if( CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][0] != 0. ) 00374 { 00375 /* this returns the log of the random noise */ 00376 realnum r = (realnum)RandGauss( mole.xMeanNoise , mole.xSTDNoise ); 00377 /* check that coefficient 0 is the one we want to hit with the mole */ 00378 /* these are used at line 2990 below, */ 00379 CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][nColl][0] += r; 00380 } 00381 } 00382 /* >>chng 04 feb 19, break out last one which is linear and must be treated separately */ 00383 /* for late one is linear so use pow */ 00384 if( CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][N_X_COLLIDER-2][0] != 0. ) 00385 { 00386 /* this returns the log of the random noise */ 00387 realnum r = (realnum)RandGauss( mole.xMeanNoise , mole.xSTDNoise ); 00388 /* check that coefficient 0 is the one we want to hit with the mole */ 00389 /* these are used at line 2990 below, */ 00390 CollRateFit[iVibHi][iRotHi][iVibLo][iRotLo][N_X_COLLIDER-2][0] *= pow((realnum)10.f,r); 00391 } 00392 } 00393 } 00394 } 00395 } 00396 } 00397 00398 00399 /* create arrays for energy sorted referencing of e, v, J */ 00400 H2_energies = (realnum*)MALLOC(sizeof(realnum)*(unsigned)nH2_energies ); 00401 H2_ipX_ener_sort = (long int*)MALLOC(sizeof(long int)*(unsigned)nH2_energies ); 00402 ipElec_H2_energy_sort = (long int*)MALLOC(sizeof(long int)*(unsigned)nH2_energies ); 00403 ipVib_H2_energy_sort = (long int*)MALLOC(sizeof(long int)*(unsigned)nH2_energies ); 00404 ipRot_H2_energy_sort = (long int*)MALLOC(sizeof(long int)*(unsigned)nH2_energies ); 00405 00406 /* this will be total collision rate from an upper to a lower level within X */ 00407 H2_X_source = (realnum*)MALLOC(sizeof(realnum)*(unsigned)nLevels_per_elec[0] ); 00408 H2_X_sink = (realnum*)MALLOC(sizeof(realnum)*(unsigned)nLevels_per_elec[0] ); 00409 00410 H2_X_coll_rate.reserve(nLevels_per_elec[0]); 00411 /* now expand out to include all lower levels as lower state */ 00412 for( i=1; i<nLevels_per_elec[0]; ++i ) 00413 { 00414 H2_X_coll_rate.reserve(i,i); 00415 } 00416 H2_X_coll_rate.alloc(); 00417 00418 /* create a vector of sorted energies for X */ 00419 ipEnergySort.reserve(mole.n_h2_elec_states); 00420 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 00421 { 00422 ipEnergySort.reserve(iElecHi,h2.nVib_hi[iElecHi]+1); 00423 for( iVibHi = 0; iVibHi <= h2.nVib_hi[iElecHi]; ++iVibHi ) 00424 { 00425 ipEnergySort.reserve(iElecHi,iVibHi,h2.nRot_hi[iElecHi][iVibHi]+1); 00426 } 00427 } 00428 ipEnergySort.alloc(); 00429 00430 nEner = 0; 00431 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 00432 { 00433 /* get set of energies */ 00434 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 00435 { 00436 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 00437 { 00438 H2_energies[nEner] = (realnum)energy_wn[iElecHi][iVibHi][iRotHi]; 00439 ipElec_H2_energy_sort[nEner] = iElecHi; 00440 ipVib_H2_energy_sort[nEner] = iVibHi; 00441 ipRot_H2_energy_sort[nEner] = iRotHi; 00442 ipEnergySort[iElecHi][iVibHi][iRotHi] = -1; 00443 ++nEner; 00444 } 00445 } 00446 } 00447 00448 ASSERT( nH2_energies == nEner ); 00449 00450 /* sort the energy levels so that we can do top-down trickle of states */ 00451 /*spsort netlib routine to sort array returning sorted indices */ 00452 spsort( 00453 /* input array to be sorted */ 00454 H2_energies, 00455 /* number of values in the molecule */ 00456 nH2_energies, 00457 /* permutation output array */ 00458 H2_ipX_ener_sort, 00459 /* flag saying what to do - 1 sorts into increasing order, not changing 00460 * the original vector, -1 sorts into decreasing order. 2, -2 change vector */ 00461 1, 00462 /* error condition, should be 0 */ 00463 &ier); 00464 00465 /* now loop over the energies confirming the order */ 00466 for( nEner=0; nEner<nH2_energies; ++nEner ) 00467 { 00468 if( nEner+1 < nLevels_per_elec[0] ) 00469 ASSERT( H2_energies[H2_ipX_ener_sort[nEner]] < 00470 H2_energies[H2_ipX_ener_sort[nEner+1]] ); 00471 /* following will print quantum indices and energies */ 00472 /*fprintf(ioQQQ,"%li\t%li\t%.3e\n", 00473 ipVib_H2_energy_sort[H2_ipX_ener_sort[nEner]], 00474 ipRot_H2_energy_sort[H2_ipX_ener_sort[nEner]], 00475 H2_energies[H2_ipX_ener_sort[nEner]]);*/ 00476 i = H2_ipX_ener_sort[nEner]; 00477 iElec = ipElec_H2_energy_sort[i]; 00478 iRot = ipRot_H2_energy_sort[i]; 00479 iVib = ipVib_H2_energy_sort[i]; 00480 /* this allows v,J to map into energy sorted array */ 00481 ipEnergySort[iElec][iVib][iRot] = nEner; 00482 } 00483 00484 /* now find number of levels in H2g */ 00485 for( nEner=0; nEner<nLevels_per_elec[0]; ++nEner ) 00486 { 00487 i = H2_ipX_ener_sort[nEner]; 00488 iRot = ipRot_H2_energy_sort[i]; 00489 iVib = ipVib_H2_energy_sort[i]; 00490 if( energy_wn[0][iVib][iRot] > ENERGY_H2_STAR ) 00491 break; 00492 nEner_H2_ground = nEner; 00493 } 00494 /* need to increment it so that this is the number of levels, not the index 00495 * of the highest level */ 00496 ++nEner_H2_ground; 00497 00498 /* this is the number of levels to do with the matrix - set with the 00499 * atom h2 matrix command, keyword ALL means to do all of X in the matrix 00500 * but number of levels within X was not known when the command was parsed, 00501 * so this was set to -1 to defer setting to all until now */ 00502 if( nXLevelsMatrix<0 ) 00503 { 00504 nXLevelsMatrix = nLevels_per_elec[0]; 00505 } 00506 else if( nXLevelsMatrix > nLevels_per_elec[0] ) 00507 { 00508 fprintf( ioQQQ, 00509 " The total number of levels used in the matrix solver was set to %li but there are only %li levels in X.\n Sorry.\n", 00510 nXLevelsMatrix , 00511 nLevels_per_elec[0]); 00512 cdEXIT(EXIT_FAILURE); 00513 } 00514 00515 /* create the main array of lines */ 00516 H2Lines.reserve(mole.n_h2_elec_states); 00517 00518 nlines = 0; 00519 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 00520 { 00521 H2Lines.reserve(iElecHi,h2.nVib_hi[iElecHi]+1); 00522 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 00523 { 00524 H2Lines.reserve(iElecHi,iVibHi,h2.nRot_hi[iElecHi][iVibHi]+1); 00525 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 00526 { 00527 /* now the lower levels */ 00528 /* NB - X is the only lower level considered here, since we are only 00529 * concerned with excited electronic levels as a photodissociation process 00530 * code exists to relax this assumption - simply change following to iElecHi */ 00531 long int lim_elec_lo = 0; 00532 H2Lines.reserve(iElecHi,iVibHi,iRotHi,1); 00533 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo ) 00534 { 00535 /* want to include all vib states in lower level if 00536 * different elec level, but only lower vib levels if 00537 * same elec level */ 00538 long int nv = h2.nVib_hi[iElecLo]; 00539 /* within X, no transitions v_hi < v_lo transitions exist */ 00540 if( iElecLo==iElecHi ) 00541 nv = iVibHi; 00542 H2Lines.reserve(iElecHi,iVibHi,iRotHi,iElecLo,nv+1); 00543 for( iVibLo=0; iVibLo<=nv; ++iVibLo ) 00544 { 00545 long nr = h2.nRot_hi[iElecLo][iVibLo]; 00546 if( iElecLo==iElecHi && iVibHi==iVibLo ) 00547 /* max because cannot malloc 0 bytes */ 00548 nr = MAX2(1,iRotHi-1); 00549 H2Lines.reserve(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,nr+1); 00550 nlines += nr+1; 00551 } 00552 } 00553 } 00554 } 00555 } 00556 00557 H2Lines.alloc(); 00558 H2_SaveLine.alloc( H2Lines.clone() ); 00559 lgH2_line_exists.alloc( H2Lines.clone() ); 00560 00561 /* set flag saying that space exists */ 00562 H2_SaveLine.zero(); 00563 lgH2_line_exists.zero(); 00564 00565 if( mole.nH2_TRACE >= mole.nH2_trace_full ) 00566 fprintf(ioQQQ," There are a total of %li lines in the entire H2 molecule.\n", nlines ); 00567 00568 /* junk the transitions */ 00569 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 00570 { 00571 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 00572 { 00573 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 00574 { 00575 /* NB - X is the only lower level considered here, since we are only 00576 * concerned with excited electronic levels as a photodissociation process 00577 * code exists to relax this assumption - simply change following to iElecHi */ 00578 long int lim_elec_lo = 0; 00579 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo ) 00580 { 00581 /* want to include all vib states in lower level if different elec level, 00582 * but only lower vib levels if same elec level */ 00583 long int nv = h2.nVib_hi[iElecLo]; 00584 if( iElecLo==iElecHi ) 00585 nv = iVibHi; 00586 for( iVibLo=0; iVibLo<=nv; ++iVibLo ) 00587 { 00588 long nr = h2.nRot_hi[iElecLo][iVibLo]; 00589 if( iElecLo==iElecHi && iVibHi==iVibLo ) 00590 nr = iRotHi-1; 00591 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 00592 { 00593 TransitionJunk( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] ); 00594 } 00595 } 00596 } 00597 } 00598 } 00599 } 00600 00601 /* now set up state pointers and zero out the transitions */ 00602 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 00603 { 00604 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 00605 { 00606 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 00607 { 00608 H2Lines[iElecHi][iVibHi][iRotHi][0][0][0].Hi = AddState2Stack(); 00609 00610 /* NB - X is the only lower level considered here, since we are only 00611 * concerned with excited electronic levels as a photodissociation process 00612 * code exists to relax this assumption - simply change following to iElecHi */ 00613 long int lim_elec_lo = 0; 00614 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo ) 00615 { 00616 /* want to include all vib states in lower level if different elec level, 00617 * but only lower vib levels if same elec level */ 00618 long int nv = h2.nVib_hi[iElecLo]; 00619 if( iElecLo==iElecHi ) 00620 nv = iVibHi; 00621 for( iVibLo=0; iVibLo<=nv; ++iVibLo ) 00622 { 00623 long nr = h2.nRot_hi[iElecLo][iVibLo]; 00624 if( iElecLo==iElecHi && iVibHi==iVibLo ) 00625 nr = iRotHi-1; 00626 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 00627 { 00628 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi = 00629 H2Lines[iElecHi][iVibHi][iRotHi][0][0][0].Hi; 00630 /* lower level is higher level of a previous transition. */ 00631 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo = 00632 H2Lines[iElecLo][iVibLo][iRotLo][0][0][0].Hi; 00633 00634 /* set initial values for each line structure */ 00635 TransitionZero( &H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] ); 00636 } 00637 } 00638 } 00639 } 00640 } 00641 } 00642 00643 /* space for the energy vector is now malloced, must read trans probs from table */ 00644 for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec ) 00645 { 00646 /* ground state has tabulated data */ 00647 H2_ReadTransprob(iElec); 00648 } 00649 00650 /* set all statistical weights - ours is total statistical weight - 00651 * including nuclear spin */ 00652 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 00653 { 00654 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 00655 { 00656 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 00657 { 00658 /* unlike atoms, for H2 nuclear spin is taken into account - so the 00659 * statistical weight of even and odd J states differ by factor of 3 - see page 166, sec par 00660 * >>>refer H2 H2_stat wght Shull, J.M., & Beckwith, S., 1982, ARAA, 20, 163-188 */ 00661 if( is_odd(iRotHi+H2_nRot_add_ortho_para[iElecHi]) ) 00662 { 00663 /* ortho */ 00664 H2_lgOrtho[iElecHi][iVibHi][iRotHi] = true; 00665 H2_stat[iElecHi][iVibHi][iRotHi] = 3.f*(2.f*iRotHi+1.f); 00666 } 00667 else 00668 { 00669 /* para */ 00670 H2_lgOrtho[iElecHi][iVibHi][iRotHi] = false; 00671 H2_stat[iElecHi][iVibHi][iRotHi] = (2.f*iRotHi+1.f); 00672 } 00673 } 00674 } 00675 } 00676 00677 /* set up transition parameters */ 00678 for( iElecHi=0; iElecHi<mole.n_h2_elec_states; ++iElecHi ) 00679 { 00680 for( iVibHi=0; iVibHi<=h2.nVib_hi[iElecHi]; ++iVibHi ) 00681 { 00682 for( iRotHi=h2.Jlowest[iElecHi]; iRotHi<=h2.nRot_hi[iElecHi][iVibHi]; ++iRotHi ) 00683 { 00684 /* now the lower levels */ 00685 /* NB - X is the only lower level considered here, since we are only 00686 * concerned with excited electronic levels as a photodissociation process 00687 * code exists to relax this assumption - simply change following to iElecHi */ 00688 long int lim_elec_lo = 0; 00689 for( iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo ) 00690 { 00691 /* want to include all vib states in lower level if different elec level, 00692 * but only lower vib levels if same elec level */ 00693 long int nv = h2.nVib_hi[iElecLo]; 00694 if( iElecLo==iElecHi ) 00695 nv = iVibHi; 00696 for( iVibLo=0; iVibLo<=nv; ++iVibLo ) 00697 { 00698 long nr = h2.nRot_hi[iElecLo][iVibLo]; 00699 00700 if( iElecLo==iElecHi && iVibHi==iVibLo ) 00701 nr = iRotHi-1; 00702 for( iRotLo=h2.Jlowest[iElecLo]; iRotLo<=nr; ++iRotLo ) 00703 { 00704 /* this is special flag for H2 - these are used in velset (in tfidle.c) to 00705 * set doppler velocities for species */ 00706 /* NB this must be kept parallel with nelem and ionstag in H2Lines transition struc, 00707 * since that struc expects to find the abundances here - abund set in hmole.c */ 00708 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->nelem = LIMELM+3; 00709 /* this does not mean anything for a molecule */ 00710 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->IonStg = 1; 00711 00712 /* statistical weights of lower and upper levels */ 00713 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->g = H2_stat[iElecLo][iVibLo][iRotLo]; 00714 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g = H2_stat[iElecHi][iVibHi][iRotHi]; 00715 00716 /* energy of the transition in wavenumbers */ 00717 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN = 00718 (realnum)(energy_wn[iElecHi][iVibHi][iRotHi] - energy_wn[iElecLo][iVibLo][iRotLo]); 00719 00720 /*wavelength of transition in Angstroms */ 00721 if( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN > SMALLFLOAT) 00722 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].WLAng = 00723 (realnum)(1.e8f/H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN / 00724 RefIndex( H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN)); 00725 00726 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyK = 00727 (realnum)(T1CM)*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN; 00728 00729 /* energy of photon in ergs */ 00730 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyErg = 00731 (realnum)(ERG1CM)*H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN; 00732 00733 /* only do this if radiative transition exists */ 00734 if( lgH2_line_exists[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo] ) 00735 { 00736 /* line redistribution function - will use complete redistribution */ 00737 /* >>chng 04 mar 26, should include damping wings, especially for electronic 00738 * transitions, had used doppler core only */ 00739 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->iRedisFun = ipCRDW; 00740 00741 /* line optical depths in direction towards source of ionizing radiation */ 00742 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->TauIn = opac.taumin; 00743 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->TauCon = opac.taumin; 00744 /* outward optical depth */ 00745 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->TauTot = 1e20f; 00746 00747 00748 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->dampXvel = 00749 (realnum)(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul/ 00750 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN/PI4); 00751 00752 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->gf = 00753 (realnum)(GetGF(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul, 00754 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN, 00755 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g ) ); 00756 00757 /* derive the absorption coefficient, call to function is gf, wl (A), g_low */ 00758 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->opacity = (realnum)( 00759 abscf(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->gf, 00760 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].EnergyWN, 00761 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->g)); 00762 00763 if( iElecHi > 0 ) 00764 { 00765 /* cosmic ray and non-thermal suprathermal excitation 00766 * to singlet state of H2 (B,C,B',D) from 00767 *>>refer H2 cr exc Dalgarno, Yan, & Liu 1999 ApJs; 00768 * cross section is scaled from transition probability 00769 * and equation 5 of 00770 *>>refer H2 cr excit Dalgarno, A., Yan, Min, & Liu, Weihong 1999, ApJS, 125, 237 00771 * in terms of hydrogen ionization cross-section 00772 *>>refer H2 cr exc Shemansky et al. 1985 00773 * this is used in mole_h2.cpp at place where Solomon 00774 * rate is added 00775 * cosmic ray excitation of triplets is 00776 * done 00777 */ 00778 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Coll.cs = (realnum)( 00779 pow(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].WLAng*1e-8,3)* 00780 (H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Hi->g/ 00781 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Lo->g)* 00782 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Emis->Aul* 00783 log(3.)*HPLANCK/(160.f*pow(PI,3)*0.5*1e-8*EN1EV)); 00784 ASSERT(H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Coll.cs>0.); 00785 /*fprintf(ioQQQ,"DEBUG %3li %3li %3li %3li %3li %3li %.2e\n", 00786 iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,iRotLo, 00787 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].cs 00788 );*/ 00789 } 00790 else 00791 { 00792 /* excitation within X - not treated this way */ 00793 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Coll.cs = 0.; 00794 } 00795 } 00796 else 00797 { 00798 /* Aul is zero but cosmic ray collisions are not 00799 * zero in this case - this is the Aul = 0 branch */ 00800 H2Lines[iElecHi][iVibHi][iRotHi][iElecLo][iVibLo][iRotLo].Coll.cs = 0.; 00801 } 00802 } 00803 } 00804 } 00805 } 00806 } 00807 } 00808 00809 /* define branching ratios for deposition of H2 formed on grain surfaces, 00810 * set true to use Takahashi distribution, false to use Draine & Bertoldi */ 00811 00812 /* loop over all types of grain surfaces */ 00813 /* >>chng 02 oct 08, resolved grain types */ 00814 /* number of different grain types H2_TOP is set in grainvar.h, 00815 * types are ice, silicate, graphite */ 00816 for( ipH2=0; ipH2<(int)H2_TOP; ++ipH2 ) 00817 { 00818 sum = 0.; 00819 sumj = 0.; 00820 sumv = 0.; 00821 sumo = 0.; 00822 sump = 0.; 00823 iElec = 0; 00824 /* first is Draine distribution */ 00825 if( hmi.chGrainFormPump == 'D' ) 00826 { 00827 /* H2 formation temperature, for equation 19, page 271, of 00828 * >>refer H2 formation distribution Draine, B.T., & Bertoldi, F., 1996, ApJ, 468, 269-289 00829 */ 00830 # define T_H2_FORM 50000. 00831 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib ) 00832 { 00833 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot ) 00834 { 00835 /* no distinction between grain surface composition */ 00836 H2_X_grain_formation_distribution[ipH2][iVib][iRot] = 00837 /* first term is nuclear H2_stat weight */ 00838 (1.f+2.f*H2_lgOrtho[iElec][iVib][iRot]) * (1.f+iVib) * 00839 (realnum)sexp( energy_wn[iElec][iVib][iRot]*T1CM/T_H2_FORM ); 00840 sum += H2_X_grain_formation_distribution[ipH2][iVib][iRot]; 00841 sumj += iRot * H2_X_grain_formation_distribution[ipH2][iVib][iRot]; 00842 sumv += iVib * H2_X_grain_formation_distribution[ipH2][iVib][iRot]; 00843 if( H2_lgOrtho[iElec][iVib][iRot] ) 00844 { 00845 sumo += H2_X_grain_formation_distribution[ipH2][iVib][iRot]; 00846 } 00847 else 00848 { 00849 /* >>chng 02 nov 14, [0][iVib][iRot] -> [ipH2][iVib][iRot], PvH */ 00850 sump += H2_X_grain_formation_distribution[ipH2][iVib][iRot]; 00851 } 00852 } 00853 } 00854 } 00855 else if( hmi.chGrainFormPump == 'T' ) 00856 { 00857 /* Takahashi 2001 distribution */ 00858 double Xrot[H2_TOP] = { 0.14 , 0.15 , 0.15 }; 00859 double Xtrans[H2_TOP] = { 0.12 , 0.15 , 0.25 }; 00860 /* first normalize the vibration distribution function */ 00861 double sumvib = 0.; 00862 double EH2; 00863 00864 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib ) 00865 { 00866 double vibdist; 00867 EH2 = EH2_eval( iVib , ipH2 ); 00868 vibdist = H2_vib_dist( iVib , ipH2 , EH2); 00869 sumvib += vibdist; 00870 } 00871 /* this branch, use distribution function from 00872 * >>refer grain physics Takahashi, Junko, 2001, ApJ, 561, 254-263 */ 00873 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib ) 00874 { 00875 double Ev = (energy_wn[iElec][iVib][0]+energy_off); 00876 double Fv; 00877 /* equation 10 of Takahashi 2001, extra term is energy offset between bottom of potential 00878 * the 0,0 level */ 00879 double Erot; 00880 /*fprintf(ioQQQ," Evvvv\t%i\t%li\t%.3e\n", ipH2 ,iVib , Ev*WAVNRYD*EVRYD);*/ 00881 00882 EH2 = EH2_eval( iVib , ipH2 ); 00883 00884 /* equation 3 of Taktahashi & Uehara */ 00885 Erot = (EH2 - Ev) * Xrot[ipH2] / (Xrot[ipH2] + Xtrans[ipH2]); 00886 00887 /* email exchange with Junko Takahashi - 00888 Thank you for your E-mail. 00889 I did not intend to generate negative Erot. 00890 I cut off the populations if their energy levels are negative, and made the total 00891 population be unity by using normalization factors (see, e.g., Eq. 12). 00892 00893 I hope that my answer is of help to you and your work is going well. 00894 With best wishes, 00895 Junko 00896 00897 >Thanks for the reply. By cutting off the population, should we set the 00898 >population to zero when Erot becomes negative, or should we set Erot to 00899 >a small positive number? 00900 00901 I just set the population to zero when Erot becomes negative. 00902 Our model is still a rough one for the vibration-rotation distribution function 00903 of H2 newly formed on dust, because we have not yet had any exact 00904 experimental or theoretical data about it. 00905 With best wishes, 00906 Junko 00907 00908 */ 00909 00910 if( Erot > 0. ) 00911 { 00912 /* the vibrational distribution */ 00913 Fv = H2_vib_dist( iVib , ipH2 , EH2) / sumvib; 00914 /*fprintf(ioQQQ," vibbb\t%li\t%.3e\n", iVib , Fv );*/ 00915 00916 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot ) 00917 { 00918 /* equation 6 of Takahashi 2001 */ 00919 double gaussian = 00920 sexp( POW2( (energy_wn[iElec][iVib][iRot] - energy_wn[iElec][iVib][0] - Erot) / 00921 (0.5 * Erot) ) ); 00922 /* equation 7 of Takahashi 2001 */ 00923 double thermal_dist = 00924 sexp( (energy_wn[iElec][iVib][iRot] - energy_wn[iElec][iVib][0]) / 00925 Erot ); 00926 00927 /* take the mean of the two */ 00928 double aver = ( gaussian + thermal_dist ) / 2.; 00929 /*fprintf(ioQQQ,"rottt\t%i\t%li\t%li\t%.3e\t%.3e\t%.3e\t%.3e\n", 00930 ipH2,iVib,iRot, 00931 (energy_wn[iElec][iVib][iRot]+energy_off)*WAVNRYD*EVRYD, 00932 gaussian, thermal_dist , aver );*/ 00933 00934 /* thermal_dist does become > 1 since Erot can become negative */ 00935 ASSERT( gaussian <= 1. /*&& thermal_dist <= 10.*/ ); 00936 00937 H2_X_grain_formation_distribution[ipH2][iVib][iRot] = (realnum)( 00938 /* first term is nuclear H2_stat weight */ 00939 (1.f+2.f*H2_lgOrtho[iElec][iVib][iRot]) * Fv * (2.*iRot+1.) * aver ); 00940 00941 sum += H2_X_grain_formation_distribution[ipH2][iVib][iRot]; 00942 sumj += iRot * H2_X_grain_formation_distribution[ipH2][iVib][iRot]; 00943 sumv += iVib * H2_X_grain_formation_distribution[ipH2][iVib][iRot]; 00944 if( H2_lgOrtho[iElec][iVib][iRot] ) 00945 { 00946 sumo += H2_X_grain_formation_distribution[ipH2][iVib][iRot]; 00947 } 00948 else 00949 { 00950 sump += H2_X_grain_formation_distribution[ipH2][iVib][iRot]; 00951 } 00952 00953 } 00954 } 00955 else 00956 { 00957 /* this branch Erot is non-positive, so no distribution */ 00958 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot ) 00959 { 00960 H2_X_grain_formation_distribution[ipH2][iVib][iRot] = 0.; 00961 } 00962 } 00963 } 00964 } 00965 # undef T_H2_FORM 00966 else if( hmi.chGrainFormPump == 't' ) 00967 { 00968 /* thermal distribution at 1.5 eV, as suggested by Amiel & Jaques */ 00969 /* thermal distribution, upper right column of page 239 of 00970 *>>refer H2 formation Le Bourlot, J, 1991, A&A, 242, 235 00971 * set with command 00972 * set h2 grain formation pumping thermal */ 00973 # define T_H2_FORM 17329. 00974 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib ) 00975 { 00976 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot ) 00977 { 00978 /* no distinction between grain surface composition */ 00979 H2_X_grain_formation_distribution[ipH2][iVib][iRot] = 00980 /* first term is nuclear H2_stat weight */ 00981 H2_stat[0][iVib][iRot] * 00982 (realnum)sexp( energy_wn[0][iVib][iRot]*T1CM/T_H2_FORM ); 00983 sum += H2_X_grain_formation_distribution[ipH2][iVib][iRot]; 00984 sumj += iRot * H2_X_grain_formation_distribution[ipH2][iVib][iRot]; 00985 sumv += iVib * H2_X_grain_formation_distribution[ipH2][iVib][iRot]; 00986 if( H2_lgOrtho[iElec][iVib][iRot] ) 00987 { 00988 sumo += H2_X_grain_formation_distribution[ipH2][iVib][iRot]; 00989 } 00990 else 00991 { 00992 /* >>chng 02 nov 14, [0][iVib][iRot] -> [ipH2][iVib][iRot], PvH */ 00993 sump += H2_X_grain_formation_distribution[ipH2][iVib][iRot]; 00994 } 00995 } 00996 } 00997 } 00998 else 00999 TotalInsanity(); 01000 01001 if( mole.nH2_TRACE >= mole.nH2_trace_full ) 01002 fprintf(ioQQQ, "H2 form grains mean J= %.3f mean v = %.3f ortho/para= %.3f\n", 01003 sumj/sum , sumv/sum , sumo/sump ); 01004 01005 iElec = 0; 01006 /* now rescale so that integral is unity */ 01007 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib ) 01008 { 01009 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot ) 01010 { 01011 H2_X_grain_formation_distribution[ipH2][iVib][iRot] /= sum; 01012 /* print the distribution function */ 01013 /*if( energy_wn[iElec][iVib][iRot] < 5200. ) 01014 fprintf(ioQQQ,"disttt\t%i\t%li\t%li\t%li\t%.4e\t%.4e\t%.4e\t%.4e\n", 01015 ipH2, iVib , iRot, (long)H2_stat[0][iVib][iRot] , 01016 energy_wn[iElec][iVib][iRot] , 01017 energy_wn[iElec][iVib][iRot]*T1CM , 01018 H2_X_grain_formation_distribution[ipH2][iVib][iRot], 01019 H2_X_grain_formation_distribution[ipH2][iVib][iRot]/H2_stat[0][iVib][iRot] 01020 );*/ 01021 } 01022 } 01023 } 01024 01025 /* at this stage the full electronic, vibration, and rotation energies have been defined, 01026 * this is an option to print the energies */ 01027 { 01028 /* set following true to get printout, false to not print energies */ 01029 if( DEBUG_ENER ) 01030 { 01031 /* print title for quantum numbers and energies */ 01032 /*fprintf(ioQQQ,"elec\tvib\trot\tenergy\n");*/ 01033 for( iElec=0; iElec<mole.n_h2_elec_states; ++iElec ) 01034 { 01035 /* now must specify the number of rotation levels within the vib levels */ 01036 for( iVib=0; iVib<=h2.nVib_hi[iElec]; ++iVib ) 01037 { 01038 for( iRot=0; iRot<=h2.nRot_hi[iElec][iVib]; ++iRot ) 01039 { 01040 fprintf(ioQQQ,"%li\t%li\t%li\t%.5e\n", 01041 iElec, iVib, iRot , 01042 energy_wn[iElec][iVib][iRot]); 01043 } 01044 } 01045 } 01046 /* this will exit the program after printing the level energies */ 01047 cdEXIT(EXIT_SUCCESS); 01048 } 01049 } 01050 01051 return; 01052 }