00001
00002
00003 #include "cddefines.h"
00004 #include "physconst.h"
00005 #include "h2.h"
00006 #include "h2_priv.h"
00007 #include "hmi.h"
00008 #include "rfield.h"
00009 #include "thirdparty.h"
00010 #include "ipoint.h"
00011 #include "mole.h"
00012
00013 STATIC void GetDissociationRateCoeff( diss_tran& tran );
00014
00015 void diatomics::Read_Mol_Diss_cross_sections()
00016 {
00017 DEBUG_ENTRY( "diatomics::Read_Mol_Diss_cross_sections()" );
00018
00019 FILE *ioDATA;
00020
00021 char chPath[FILENAME_PATH_LENGTH_2],
00022 chDirectory[FILENAME_PATH_LENGTH_2],
00023 chLine[FILENAME_PATH_LENGTH_2];
00024
00025
00026 const int ipNUM_FILES = 1;
00027
00028 char chFileNames[ipNUM_FILES][17] =
00029 {"cont_diss.dat"} ;
00030
00031
00032 ASSERT( lgEnabled );
00033
00034
00035
00036 Cont_Dissoc_Rate.reserve( n_elec_states );
00037 for( int iElecLo=0; iElecLo<n_elec_states; ++iElecLo )
00038 {
00039 Cont_Dissoc_Rate.reserve( iElecLo, nVib_hi[iElecLo]+1 );
00040 for( int iVibLo=0; iVibLo<=nVib_hi[iElecLo]; ++iVibLo )
00041 {
00042 Cont_Dissoc_Rate.reserve( iElecLo, iVibLo, nRot_hi[iElecLo][iVibLo]+1 );
00043 }
00044 }
00045 Cont_Dissoc_Rate.alloc();
00046
00047
00048 # ifdef _WIN32
00049 strcpy( chDirectory, "h2\\" );
00050 # else
00051 strcpy( chDirectory, "h2/" );
00052 # endif
00053
00054
00055 strcpy( chPath, chDirectory );
00056 strcat( chPath, chFileNames[0] );
00057 ioDATA = open_data( chPath, "r" );
00058
00059 Diss_Trans.clear();
00060
00061
00062 while( read_whole_line( chLine , (int)sizeof(chLine) , ioDATA ) != NULL )
00063 {
00064 static bool skipData = false;
00065 long ini=0, inf, iv, ij;
00066
00067
00068
00069 if(sizeof(chLine) >= 2 && chLine[0] == '#' && chLine[1] == '!' && chLine[4] == 'n' && chLine[5] =='e' )
00070 {
00071 sscanf(chLine,"#! nei=%li, nef=%li, vi= %li, ji= %li",&ini, &inf, &iv, &ij);
00072
00073 if( ini > n_elec_states )
00074 skipData = true;
00075 else if( iv > nVib_hi[ini] )
00076 skipData = true;
00077 else if( ij > nRot_hi[ini][iv] )
00078 skipData = true;
00079 else
00080 {
00081 skipData = false;
00082 diss_level Leveli = {ini, iv, ij};
00083 diss_level Levelf = {inf, -1, -1};
00084 diss_tran thisTran( Leveli, Levelf );
00085 Diss_Trans.push_back( thisTran );
00086 }
00087 }
00088
00089 if( skipData )
00090 continue;
00091
00092
00093
00094 if(chLine[0] != '#')
00095 {
00096 realnum energy, xsection;
00097 const realnum AngSquared = 1e-16f;
00098 sscanf(chLine,"%f,%f",&energy, &xsection);
00099
00100
00101 Diss_Trans.back().energies.push_back( energy*WAVNRYD );
00102
00103
00104 Diss_Trans.back().xsections.push_back( xsection*AngSquared );
00105 }
00106 }
00107
00108 fclose(ioDATA);
00109 return;
00110 }
00111
00112 double diatomics::MolDissocOpacity( const diss_tran& tran, const double& Mol_Ene )
00113 {
00114 DEBUG_ENTRY( "diatomics::MolDissocOpacity()" );
00115
00116 double cross_section = MolDissocCrossSection( tran, Mol_Ene ) *
00117 states[ ipEnergySort[ tran.initial.n ][ tran.initial.v ][ tran.initial.j ] ].Pop();
00118 return cross_section;
00119 }
00120
00121 double MolDissocCrossSection( const diss_tran& tran, const double& Mol_Ene )
00122 {
00123 DEBUG_ENTRY( "MolDissocCrossSection()" );
00124
00125 double photodiss_cs = 0.;
00126
00127
00128 if (Mol_Ene < tran.energies[0])
00129 photodiss_cs = 0.;
00130
00131
00132 else if(Mol_Ene > tran.energies.back())
00133 photodiss_cs = tran.xsections.back()/sqrt(powi(Mol_Ene/tran.energies.back(),7));
00134
00135
00136 else
00137 {
00138 ASSERT( Mol_Ene > tran.energies[0] && Mol_Ene < tran.energies.back() );
00139 photodiss_cs = linint(&tran.energies[0],
00140 &tran.xsections[0],
00141 tran.xsections.size(),
00142 Mol_Ene);
00143 }
00144
00145 return photodiss_cs;
00146 }
00147
00148 void diatomics::Mol_Photo_Diss_Rates( void )
00149 {
00150 DEBUG_ENTRY( "diatomics::Mol_Photo_Diss_Rates()" );
00151
00152
00153
00154
00155 ASSERT( lgEnabled && mole_global.lgStancil );
00156
00157
00158 Cont_Dissoc_Rate.zero();
00159 Cont_Dissoc_Rate_H2g = 0.;
00160 Cont_Dissoc_Rate_H2s = 0.;
00161
00162 for_each( Diss_Trans.begin(), Diss_Trans.end(), GetDissociationRateCoeff );
00163 for( vector< diss_tran >::const_iterator dt = Diss_Trans.begin(); dt != Diss_Trans.end(); ++dt )
00164 {
00165 double rate = (*this).GetDissociationRate( *dt );
00166
00167 Cont_Dissoc_Rate[dt->initial.n][dt->initial.v][dt->initial.j] += dt->rate_coeff;
00168 if( states[ ipEnergySort[dt->initial.n][dt->initial.v][dt->initial.j] ].energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
00169 {
00170
00171
00172 Cont_Dissoc_Rate_H2s += rate;
00173 }
00174 else
00175 {
00176
00177
00178 Cont_Dissoc_Rate_H2g += rate;
00179 }
00180 }
00181
00182
00183
00184
00185 Cont_Dissoc_Rate_H2g = Cont_Dissoc_Rate_H2g / SDIV(H2_den_g);
00186 Cont_Dissoc_Rate_H2s = Cont_Dissoc_Rate_H2s / SDIV(H2_den_s);
00187
00188 return;
00189 }
00190
00191 STATIC void GetDissociationRateCoeff( diss_tran& tran )
00192 {
00193 DEBUG_ENTRY( "GetDissociationRateCoeff()" );
00194
00195 long index_min = ipoint( tran.energies[0] );
00196 long index_max = ipoint( tran.energies.back() );
00197 index_max = MIN2( index_max, rfield.nflux-1 );
00198
00199 tran.rate_coeff = 0.;
00200 for(long i = index_min; i <= index_max; ++i)
00201 {
00202
00203 double photodiss_cs = MolDissocCrossSection( tran, rfield.anu[i] );
00204
00205
00206 tran.rate_coeff += photodiss_cs*( rfield.flux[0][i] + rfield.ConInterOut[i]+
00207 rfield.outlin[0][i]+ rfield.outlin_noplot[i] );
00208 }
00209 }
00210
00211 double diatomics::GetDissociationRate( const diss_tran& tran )
00212 {
00213 DEBUG_ENTRY( "diatomics::GetDissociationRate()" );
00214
00215 long iElecLo = tran.initial.n;
00216 long iVibLo = tran.initial.v;
00217 long iRotLo = tran.initial.j;
00218
00219
00220 return states[ ipEnergySort[iElecLo][iVibLo][iRotLo] ].Pop() * tran.rate_coeff;
00221 }
00222
00223 double diatomics::Cont_Diss_Heat_Rate( void )
00224 {
00225 DEBUG_ENTRY( "diatomics::Cont_Diss_Heat_Rate()" );
00226
00227
00228
00229
00230
00231
00232 if( !mole_global.lgStancil || !lgEnabled )
00233 return 0.;
00234
00235 Mol_Photo_Diss_Rates();
00236
00237 double Cont_Dissoc_Heat_Rate = 0.0;
00238 for( vector< diss_tran >::const_iterator dt = Diss_Trans.begin(); dt != Diss_Trans.end(); ++dt )
00239 Cont_Dissoc_Heat_Rate += (*this).GetHeatRate( *dt );
00240
00241 return Cont_Dissoc_Heat_Rate;
00242 }
00243
00244 double diatomics::GetHeatRate( const diss_tran& tran )
00245 {
00246 DEBUG_ENTRY( "diatomics::GetHeatRate()" );
00247
00248 long index_min = ipoint( tran.energies[0] );
00249 long index_max = ipoint( tran.energies.back() );
00250 index_max = MIN2( index_max, rfield.nflux-1 );
00251
00252 long iElecLo = tran.initial.n;
00253 long iVibLo = tran.initial.v;
00254 long iRotLo = tran.initial.j;
00255 double rate = 0.;
00256
00257 for( long i = index_min; i<= index_max; ++i )
00258 {
00259
00260 double E_H2 = states[ ipEnergySort[iElecLo][iVibLo][iRotLo] ].energy().Ryd();
00261
00262
00263 double E_Reactants = rfield.anu[i] + E_H2;
00264
00265
00266
00267 double E_Products = 0.75;
00268
00269
00270 double density = states[ ipEnergySort[iElecLo][iVibLo][iRotLo] ].Pop();
00271
00272
00273
00274
00275 double Rate_Coeff = tran.rate_coeff;
00276
00277
00278
00279
00280
00281 rate += EN1RYD*(E_Reactants - E_Products) * Rate_Coeff * density;
00282 }
00283
00284 return rate;
00285 }
00286