00001
00002
00003
00004 #include "cddefines.h"
00005 #include "grainvar.h"
00006 #include "phycon.h"
00007 #include "mole.h"
00008 #include "hmi.h"
00009 #include "dense.h"
00010 #include "h2.h"
00011 #include "h2_priv.h"
00012 #include "deuterium.h"
00013
00014
00015 void diatomics::mole_H2_form( void )
00016 {
00017 DEBUG_ENTRY( "mole_H2_form()" );
00018
00019
00020
00021
00022 for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
00023 {
00024 for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
00025 {
00026
00027
00028
00029
00030 H2_X_formation[iVib][iRot] = 0.;
00031 H2_X_Hmin_back[iVib][iRot] = 0.;
00032 }
00033 }
00034
00035
00036
00037
00038 hmi.H2_forms_grains = 0.;
00039 hmi.H2star_forms_grains = 0.;
00040
00041 double sum_check = 0.;
00042
00043
00044 for( size_t nd=0; nd < gv.bin.size(); ++nd )
00045 {
00046 int ipH2 = (int)gv.which_H2distr[gv.bin[nd]->matType];
00047 for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
00048 {
00049 for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
00050 {
00051
00052 double one =
00053
00054 H2_X_grain_formation_distribution[ipH2][iVib][iRot] *
00055
00056 gv.bin[nd]->rate_h2_form_grains_used;
00057
00058 sum_check += one;
00059
00060
00061
00062
00063 H2_X_formation[iVib][iRot] += (realnum)one*dense.xIonDense[ipHYDROGEN][0];
00064
00065
00066
00067
00068
00069 if( states[ ipEnergySort[0][iVib][iRot] ].energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
00070 {
00071 hmi.H2star_forms_grains += one;
00072 }
00073 else
00074 {
00075
00076 hmi.H2_forms_grains += one;
00077 }
00078 }
00079 }
00080 }
00081
00082 ASSERT( fp_equal_tol( sum_check, gv.rate_h2_form_grains_used_total, 1e-5 * (sum_check + SMALLFLOAT) ) );
00083
00084
00085
00086 hmi.H2star_forms_hminus = 0.;
00087 hmi.H2_forms_hminus = 0.;
00088 double frac_lo , frac_hi;
00089 long ipT;
00090
00091
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
00107 for( ipT=0; ipT<nTE_HMINUS-1; ++ipT )
00108 {
00109 if( H2_logte_hminus[ipT+1]>phycon.alogte )
00110 break;
00111 }
00112 frac_hi = (phycon.alogte-H2_logte_hminus[ipT])/(H2_logte_hminus[ipT+1]-H2_logte_hminus[ipT]);
00113 frac_lo = 1.-frac_hi;
00114 }
00115
00116
00117
00118
00119
00120 double rate = (mole.findrk("H,H-=>H2,e-")+mole.findrk("H,H-=>H2*,e-")) * dense.xIonDense[ipHYDROGEN][0] * findspecieslocal("H-")->den;
00121
00122 double rateback = (mole.findrk("H2,e-=>H,H-")+mole.findrk("H2*,e-=>H,H-"))*dense.eden;
00123 double oldrate = 0.;
00124
00125
00126 for( long iVib=0; iVib<=nVib_hi[0]; ++iVib )
00127 {
00128 for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
00129 {
00130
00131
00132 double rate_interp =
00133 frac_lo*H2_X_hminus_formation_distribution[ipT][iVib][iRot] +
00134 frac_hi*H2_X_hminus_formation_distribution[ipT+1][iVib][iRot];
00135
00136
00137
00138 double one = rate * rate_interp;
00139
00140
00141 H2_X_Hmin_back[iVib][iRot] = (realnum)(rateback * rate_interp);
00142
00143
00144 H2_X_formation[iVib][iRot] += (realnum)one;
00145
00146 oldrate += rate_interp;
00147
00148
00149 if( states[ ipEnergySort[0][iVib][iRot] ].energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
00150 {
00151 hmi.H2star_forms_hminus += one;
00152 }
00153 else
00154 {
00155
00156 hmi.H2_forms_hminus += one;
00157 }
00158 }
00159 }
00160
00161 ASSERT( fabs(1.-oldrate)<1e-4 );
00162
00163
00164
00165
00166
00167
00168 rate = mole.findrk("H,H2+=>H+,H2*") * dense.xIonDense[ipHYDROGEN][0] * findspecieslocal("H2+")->den;
00169
00170 H2_X_formation[4][0] += (realnum)rate;
00171
00172 fixit();
00173 ASSERT( this==&h2 || this==&hd );
00174 if( this != &h2 )
00175 {
00176 for( long iVib=0; iVib<=nVib_hi[0]; ++iVib )
00177 {
00178 for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
00179 {
00180
00181
00182 H2_X_formation[iVib][iRot] *= deut.gas_phase/dense.gas_phase[ipHYDROGEN];
00183 }
00184 }
00185 }
00186
00187 return;
00188 }