00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "cddefines.h"
00012 #include "phycon.h"
00013 #include "hmi.h"
00014 #include "taulines.h"
00015 #include "h2.h"
00016 #include "h2_priv.h"
00017 #include "rfield.h"
00018 #include "thermal.h"
00019 #include "gammas.h"
00020 #include "ionbal.h"
00021
00022
00023
00024 void diatomics::H2_Solomon_rate( void )
00025 {
00026 DEBUG_ENTRY( "H2_Solomon_rate()" );
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 Solomon_dissoc_rate_g = 0.;
00037 Solomon_dissoc_rate_s = 0.;
00038
00039
00040 Solomon_elec_decay_g = 0.;
00041 Solomon_elec_decay_s = 0.;
00042
00043
00044
00045
00046
00047 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
00048 {
00049 qList::iterator Hi = (*tr).Hi() ;
00050 if( (*Hi).n() < 1 ) continue;
00051 double factor = (double)H2_dissprob[(*Hi).n()][(*Hi).v()][(*Hi).J()]/H2_rad_rate_out[(*Hi).n()][(*Hi).v()][(*Hi).J()];
00052
00053
00054
00055 double rate_up_cont = (*(*tr).Lo()).Pop() * (*tr).Emis().pump() * factor;
00056
00057
00058 double elec_decay = (*(*tr).Hi()).Pop() * (*tr).Emis().Aul() *
00059 ((*tr).Emis().Pesc() + (*tr).Emis().Pdest() + (*tr).Emis().Pelec_esc());
00060 if( (*(*tr).Lo()).energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
00061 {
00062
00063
00064 Solomon_dissoc_rate_s += rate_up_cont;
00065
00066 Solomon_elec_decay_s += elec_decay;
00067 }
00068 else
00069 {
00070
00071
00072 Solomon_dissoc_rate_g += rate_up_cont;
00073
00074 Solomon_elec_decay_g += elec_decay;
00075 }
00076 }
00077
00078 double H2_sum_excit_elec_den = GetExcitedElecDensity();
00079
00080
00081
00082
00083 if( *dense_total > SMALLFLOAT )
00084 {
00085 Solomon_elec_decay_g /= SDIV( H2_sum_excit_elec_den );
00086 Solomon_elec_decay_s /= SDIV( H2_sum_excit_elec_den );
00087
00088
00089 Solomon_dissoc_rate_s /= SDIV(H2_den_s);
00090
00091
00092 Solomon_dissoc_rate_g /= SDIV(H2_den_g);
00093
00094 }
00095 else
00096 {
00097 Solomon_dissoc_rate_s = 0.;
00098 Solomon_dissoc_rate_g = 0.;
00099
00100 }
00101
00102
00103
00104
00105
00106
00107 return;
00108 }
00109
00110
00111 double diatomics::gs_rate( void )
00112 {
00113 DEBUG_ENTRY( "diatomics::gs_rate()" );
00114
00115
00116 double ground_to_star_rate = 0.;
00117
00118
00119 for( long ipLoX=0; ipLoX < nEner_H2_ground; ++ipLoX )
00120 {
00121
00122
00123 for( long iElecHi=1; iElecHi<n_elec_states; ++iElecHi )
00124 {
00125 for( long iVibHi=0; iVibHi<=nVib_hi[iElecHi]; ++iVibHi )
00126 {
00127 for( long iRotHi=Jlowest[iElecHi]; iRotHi<=nRot_hi[iElecHi][iVibHi]; ++iRotHi )
00128 {
00129 long ipHi = ipEnergySort[iElecHi][iVibHi][iRotHi];
00130 if( lgH2_radiative[ipHi][ipLoX] )
00131 {
00132
00133
00134
00135 double rate_up_cont =
00136 states[ipLoX].Pop() *
00137 trans[ ipTransitionSort[ipHi][ipLoX] ].Emis().pump();
00138 double decay_star = H2_rad_rate_out[iElecHi][iVibHi][iRotHi] - H2_dissprob[iElecHi][iVibHi][iRotHi];
00139
00140
00141
00142 for( long ipOther=0; ipOther < nEner_H2_ground; ++ipOther )
00143 {
00144 if( lgH2_radiative[ipHi][ipOther] )
00145 {
00146 EmissionProxy em = trans[ ipTransitionSort[ipHi][ipOther] ].Emis();
00147 decay_star -= em.Aul() * ( em.Pesc() + em.Pdest() + em.Pelec_esc() );
00148 }
00149 }
00150
00151
00152 decay_star = MAX2(0., decay_star)/SDIV(H2_rad_rate_out[iElecHi][iVibHi][iRotHi]);
00153 ground_to_star_rate += rate_up_cont*decay_star;
00154
00155 }
00156 }
00157 }
00158 }
00159 }
00160
00161
00162 ground_to_star_rate /= SDIV( H2_den_g );
00163 return ground_to_star_rate;
00164 }
00165
00166 long diatomics::OpacityCreate( double *stack )
00167 {
00168 DEBUG_ENTRY( "diatomics::OpacityCreate()" );
00169
00170 ASSERT( photoion_opacity_fun != NULL );
00171
00172 for( long i=ip_photo_opac_thresh-1; i < rfield.nupper; i++ )
00173 {
00174
00175
00176 stack[ i - ip_photo_opac_thresh + ip_photo_opac_offset ] =
00177 photoion_opacity_fun( rfield.AnuOrg[i] );
00178
00179 }
00180
00181
00182 return rfield.nupper - ip_photo_opac_thresh + 1;
00183 }
00184
00185
00186
00187 void diatomics::H2_zero_pops_too_low( void )
00188 {
00189 DEBUG_ENTRY( "H2_zero_pops_too_low()" );
00190
00191
00192 fill_n( pops_per_elec, N_ELEC, 0. );
00193 pops_per_vib.zero();
00194 for( qList::iterator st = states.begin(); st != states.end(); ++st )
00195 {
00196 double pop = H2_populations_LTE[ (*st).n() ][ (*st).v() ][ (*st).J() ] * (*dense_total);
00197 H2_old_populations[ (*st).n() ][ (*st).v() ][ (*st).J() ] = pop;
00198 (*st).Pop() = pop;
00199 }
00200
00201
00202 for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
00203 {
00204
00205 (*tr).Emis().PopOpc() = (*(*tr).Lo()).Pop() - (*(*tr).Hi()).Pop() * (*(*tr).Lo()).g() / (*(*tr).Hi()).g();
00206
00207
00208 (*tr).Coll().cool() = 0.;
00209 (*tr).Coll().heat() = 0.;
00210
00211
00212 (*tr).Emis().xIntensity() = 0.;
00213 (*tr).Emis().phots() = 0.;
00214 (*tr).Emis().ots() = 0.;
00215 }
00216
00217 photodissoc_BigH2_H2s = 0.;
00218 photodissoc_BigH2_H2g = 0.;
00219 HeatDiss = 0.;
00220 HeatDexc = 0.;
00221 HeatDexc_deriv = 0.;
00222 Solomon_dissoc_rate_g = 0.;
00223 Solomon_dissoc_rate_s = 0.;
00224 return;
00225 }
00226
00227
00228 void diatomics::mole_H2_LTE( void )
00229 {
00230 DEBUG_ENTRY( "mole_H2_LTE()" );
00231
00232
00233 if( ! fp_equal( phycon.te, TeUsedBoltz ) )
00234 {
00235 double part_fun = 0.;
00236 TeUsedBoltz = phycon.te;
00237
00238 for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
00239 {
00240 long iElec = (*st).n();
00241 long iVib = (*st).v();
00242 long iRot = (*st).J();
00243 H2_Boltzmann[iElec][iVib][iRot] =
00244
00245
00246 sexp( (*st).energy().K() / phycon.te );
00247
00248 part_fun += H2_Boltzmann[iElec][iVib][iRot] * (*st).g();
00249 ASSERT( part_fun > 0 );
00250 }
00251
00252 for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
00253 {
00254 long iElec = (*st).n();
00255 long iVib = (*st).v();
00256 long iRot = (*st).J();
00257
00258
00259 H2_populations_LTE[iElec][iVib][iRot] =
00260 H2_Boltzmann[iElec][iVib][iRot] *
00261 (*st).g() / part_fun;
00262 }
00263 if( nTRACE >= n_trace_full )
00264 fprintf(ioQQQ,
00265 "mole_H2_LTE set H2_Boltzmann factors, T=%.2f, partition function is %.2f\n",
00266 phycon.te,
00267 part_fun);
00268 }
00269
00270 return;
00271 }
00272
00273
00274 void diatomics::H2_Reset( void )
00275 {
00276
00277 DEBUG_ENTRY( "H2_Reset()" );
00278
00279 if(nTRACE)
00280 fprintf(ioQQQ,
00281 "\n***************H2_Reset called, resetting nCall_this_iteration, zone %.2f iteration %li\n",
00282 fnzone,
00283 iteration );
00284
00285
00286
00287 nCall_this_iteration = 0;
00288
00289
00290
00291 renorm_max = 1.;
00292 renorm_min = 1.;
00293
00294
00295 nH2_pops = 0;
00296 nH2_zone = 0;
00297
00298 nzone_nlevel_set = 0;
00299
00300 nzoneAsEval = -1;
00301 iterationAsEval = -1;
00302
00303 TeUsedColl = -1;
00304 TeUsedBoltz = -1;
00305
00306 lgEvaluated = false;
00307
00308
00309 H2_SaveLine.zero();
00310
00311
00312
00313
00314 if( nElecLevelOutput < 1 )
00315 nElecLevelOutput = n_elec_states;
00316
00317 return;
00318 }
00319
00320
00321
00322
00323
00324
00325 double Yan_H2_CS( double energy_ryd )
00326 {
00327
00328 double energy_keV;
00329 double cross_section;
00330 double x;
00331 double xsqrt , x15 , x2;
00332 double energy = energy_ryd * EVRYD;
00333
00334 DEBUG_ENTRY( "Yan_H2_CS()" );
00335
00336
00337 x = energy / 15.4;
00338 energy_keV = energy/1000.0;
00339 xsqrt = sqrt(x);
00340 x15 = x * xsqrt;
00341 x2 = x*x;
00342
00343 if( energy < 15.4 )
00344 {
00345 cross_section = 0.;
00346 }
00347
00348 else if(energy >= 15.4 && energy < 18.)
00349 {
00350 cross_section = 1e7 * (1 - 197.448/xsqrt + 438.823/x - 260.481/x15 + 17.915/x2);
00351
00352 cross_section = MAX2( 0. , cross_section );
00353 }
00354
00355 else if(energy >= 18. && energy <= 30.)
00356 {
00357 cross_section = (-145.528 +351.394*xsqrt - 274.294*x + 74.320*x15)/pow(energy_keV,3.5);
00358 }
00359
00360 else if(energy > 30. && energy <= 85.)
00361 {
00362 cross_section = (65.304 - 91.762*xsqrt + 51.778*x - 9.364*x15)/pow(energy_keV,3.5);
00363 }
00364
00365
00366 else
00367 {
00368 cross_section = 45.57*(1 - 2.003/xsqrt - 4.806/x + 50.577/x15 - 171.044/x2 + 231.608/(xsqrt*x2) - 81.885/(x*x2))/pow(energy_keV,3.5);
00369 }
00370
00371 return( cross_section * 1e-24 );
00372 }
00373
00374 void diatomics::CalcPhotoionizationRate(void)
00375 {
00376 DEBUG_ENTRY( "diatomics::CalcPhotoionizationRate()" );
00377
00378
00379 if( ( nzone_eval!=nzone || iteration_evaluated!=iteration ) || !nzone )
00380 {
00381 t_phoHeat photoHeat;
00382
00383 photoionize_rate =
00384 GammaK( ip_photo_opac_thresh,
00385 rfield.nupper,
00386 ip_photo_opac_offset,1.,
00387 &photoHeat)*
00388 ionbal.lgPhotoIoniz_On +
00389
00390
00391
00392 2.*ionbal.CompRecoilIonRate[ipHYDROGEN][0];
00393
00394
00395
00396 photo_heat_soft = photoHeat.HeatLowEnr * ionbal.lgPhotoIoniz_On;
00397 photo_heat_hard = photoHeat.HeatHiEnr * ionbal.lgPhotoIoniz_On;
00398 nzone_eval = nzone;
00399 iteration_evaluated = iteration;
00400 }
00401
00402 return;
00403 }
00404