00001 /* This file is part of Cloudy and is copyright (C)1978-2010 by Gary J. Ferland and 00002 * others. For conditions of distribution and use see copyright notice in license.txt */ 00003 /*mole_H2_form find state specific rates grains and H- form H2 */ 00004 #include "cddefines.h" 00005 #include "grainvar.h" 00006 #include "phycon.h" 00007 #include "hmi.h" 00008 #include "dense.h" 00009 #include "h2.h" 00010 #include "h2_priv.h" 00011 00012 /*mole_H2_form find state specific rates grains and H- form H2 */ 00013 void mole_H2_form( void ) 00014 { 00015 00016 long int iVib, iRot; 00017 long int ipT; 00018 double rate, oldrate , 00019 frac_lo , frac_hi; 00020 00021 DEBUG_ENTRY( "mole_H2_form()" ); 00022 00023 /* rate of entry into X from H- and formation on grain surfaces 00024 * will one of several distribution functions derived elsewhere 00025 * first zero out formation rates and rates others collide into particular level */ 00026 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib ) 00027 { 00028 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot ) 00029 { 00030 /* this will be the rate formation (s-1) of H2 due to 00031 * both formation on grain surfaces and the H minus route, 00032 * also H2+ + H => H2 + H+ into one vJ level */ 00033 /* units cm-3 s-1 */ 00034 H2_X_formation[iVib][iRot] = 0.; 00035 H2_X_Hmin_back[iVib][iRot] = 0.; 00036 } 00037 } 00038 00039 /* loop over all grain types, finding total formation rate into each ro-vib level, 00040 * also keeps trace of total formation into H2 ground and star, as defined by Tielens & Hollenbach, 00041 * these are used in the H molecular network */ 00042 hmi.H2_forms_grains = 0.; 00043 hmi.H2star_forms_grains = 0.; 00044 /* >>chng 02 oct 08, resolved grain types */ 00045 for( size_t nd=0; nd < gv.bin.size(); ++nd ) 00046 { 00047 int ipH2 = (int)gv.which_H2distr[gv.bin[nd]->matType]; 00048 for( iVib = 0; iVib <= h2.nVib_hi[0]; ++iVib ) 00049 { 00050 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot ) 00051 { 00052 /* >>chng 02 nov 14, changed indexing into H2_X_grain_formation_distribution and gv.bin, PvH */ 00053 realnum one = 00054 /* H2_X_grain_formation_distribution is normalized to a summed total of unity */ 00055 H2_X_grain_formation_distribution[ipH2][iVib][iRot] * 00056 /* units of following are s-1 */ 00057 (realnum)gv.bin[nd]->rate_h2_form_grains_used; 00058 /* final units are s-1 */ 00059 /* units cm-3 s-1 */ 00060 /* >>chng 04 may 05, added atomic hydrogen density, units cm-3 s-1 */ 00061 H2_X_formation[iVib][iRot] += one*dense.xIonDense[ipHYDROGEN][0]; 00062 00063 /* save rates for formation into "H2" and "H2*" in the chemical 00064 * network - it resolves the H2 into two species, as in 00065 * Hollenbach / Tielens work - these rates will be used in the 00066 * chemistry solver to get H2 and H2* densities */ 00067 if( energy_wn[0][iVib][iRot] < ENERGY_H2_STAR ) 00068 { 00069 /* unit s-1, h2 means h2g*/ 00070 hmi.H2_forms_grains += one; 00071 } 00072 else 00073 { 00074 hmi.H2star_forms_grains += one; 00075 } 00076 } 00077 } 00078 } 00079 00080 /* formation of H2 in excited states from H- H minus */ 00081 /* >>chng 02 oct 17, temp dependent fit to rate, updated reference, 00082 * about 40% larger than before */ 00083 /* rate in has units cm-3 s-1 */ 00084 rate = hmi.Hmolec[ipMHm] * dense.xIonDense[ipHYDROGEN][0] * hmi.assoc_detach; 00085 /*rate = hmi.hminus*1.35e-9f;*/ 00086 /* convert to dimensionless factors that add to unity */ 00087 /* >>chng 02 oct 17, use proper distribution function */ 00088 hmi.H2star_forms_hminus = 0.; 00089 hmi.H2_forms_hminus = 0.; 00090 oldrate = 0.; 00091 /* which temperature point to use? */ 00092 if( phycon.alogte<=1. ) 00093 { 00094 ipT = 0; 00095 frac_lo = 1.; 00096 frac_hi = 0.; 00097 } 00098 else if( phycon.alogte>=4. ) 00099 { 00100 ipT = nTE_HMINUS-2; 00101 frac_lo = 0.; 00102 frac_hi = 1.; 00103 } 00104 else 00105 { 00106 /* find the temp */ 00107 for( ipT=0; ipT<nTE_HMINUS-1; ++ipT ) 00108 { 00109 if( H2_te_hminus[ipT+1]>phycon.alogte ) 00110 break; 00111 } 00112 frac_hi = (phycon.alogte-H2_te_hminus[ipT])/(H2_te_hminus[ipT+1]-H2_te_hminus[ipT]); 00113 frac_lo = 1.-frac_hi; 00114 } 00115 00116 /* we now know how to interpolate, now fill in H- formation sites */ 00117 for( iVib=0; iVib<=h2.nVib_hi[0]; ++iVib ) 00118 { 00119 for( iRot=h2.Jlowest[0]; iRot<=h2.nRot_hi[0][iVib]; ++iRot ) 00120 { 00121 /* the temperature-interpolated distribution function, adds up to unity, 00122 * dimensionless */ 00123 double rate_interp = 00124 frac_lo*H2_X_hminus_formation_distribution[ipT][iVib][iRot] + 00125 frac_hi*H2_X_hminus_formation_distribution[ipT+1][iVib][iRot]; 00126 00127 /* above rate was set, had dimensions cm-3 s-1 00128 * rate is product of parent densities and total formation rate */ 00129 realnum one = (realnum)(rate * rate_interp); 00130 00131 /* save this rate [cm3 s-1] for back reaction in levels solver for v,J */ 00132 H2_X_Hmin_back[iVib][iRot] = (realnum)(rate_interp * hmi.assoc_detach); 00133 00134 /* units cm-3 s-1 */ 00135 H2_X_formation[iVib][iRot] += one; 00136 00137 oldrate += rate_interp; 00138 00139 /* save rates to pass back into molecule network */ 00140 if( energy_wn[0][iVib][iRot] < ENERGY_H2_STAR ) 00141 { 00142 /* unit cm-3s-1, h2 means h2g*/ 00143 hmi.H2_forms_hminus += one; 00144 } 00145 else 00146 { 00147 hmi.H2star_forms_hminus += one; 00148 } 00149 } 00150 } 00151 /* confirm that shape function is normalized correctly */ 00152 ASSERT( fabs(1.-oldrate)<1e-4 ); 00153 00154 /* >>chng 03 feb 10, add this population process */ 00155 /* H2+ + H => H2 + H+, 00156 * >>refer H2 population Krstic, P.S., preprint 00157 * all goes into v=4 but no J information, assume into J = 0 */ 00158 /* >>chng 04 may 05, add density at end */ 00159 rate = hmi.bh2h2p * hmi.Hmolec[ipMH2p] * dense.xIonDense[ipHYDROGEN][0]; 00160 iVib = 4; 00161 iRot = 0; 00162 /* units cm-3 s-1 */ 00163 H2_X_formation[iVib][iRot] += (realnum)rate; 00164 00165 return; 00166 }