49 double Edust = DissocEnergy *
Xdust[
ipH2] *
50 ( 1. - ( (Ev - Evm) / (DissocEnergy+
energy_off-Evm)) *
65 double G1[
H2_TOP] = { 0.3 , 0.4 , 0.9 };
66 double G2[
H2_TOP] = { 0.6 , 0.6 , 0.4 };
87 if( (*Lo).n()==0 && (*tr).Emis().Aul() > 1.01e-30 )
121 bool lgDebugPrt =
false;
125 fprintf(
ioQQQ,
" H2_Create called in DEBUG mode.\n");
154 for(
long iVibHi = 0; iVibHi <=
nVib_hi[iElecHi]; ++iVibHi )
165 for(
unsigned nEner = 0; nEner <
states.size(); ++nEner )
167 long iElec =
states[nEner].n();
168 long iVib =
states[nEner].v();
169 long iRot =
states[nEner].J();
192 H2_lgOrtho[(*st).n()][(*st).v()][(*st).J()] =
true;
193 H2_stat[(*st).n()][(*st).v()][(*st).J()] = 3.f*(2.f*(*st).J()+1.f);
198 H2_lgOrtho[(*st).n()][(*st).v()][(*st).J()] =
false;
199 H2_stat[(*st).n()][(*st).v()][(*st).J()] = (2.f*(*st).J()+1.f);
201 (*st).g() =
H2_stat[(*st).n()][(*st).v()][(*st).J()];
212 " H2_Create: there are %li electronic levels, in each level there are",
215 " for a total of %li levels.\n", (
long int)
states.size() );
240 " 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",
256 fprintf(
ioQQQ,
"%li\t%li\t%li\t%.5e\n", (*st).n(), (*st).v(), (*st).J(), (*st).energy().WN() );
275 fprintf(
ioQQQ,
"elec %li highest vib= %li\n", iElec ,
nVib_hi[iElec] );
287 for(
long iVib = 0; iVib <=
nVib_hi[iElec]; ++iVib )
313 for(
long iVib = 0; iVib <=
nVib_hi[0]; ++iVib )
332 for(
long iVib = 0; iVib <=
nVib_hi[0]; ++iVib )
358 for(
long iVib = 0; iVib <=
nVib_hi[0]; ++iVib )
379 for(
long iVib = 0; iVib <=
nVib_hi[0]; ++iVib )
399 long iElec = (*st).n();
400 if( iElec > 0 )
continue;
401 long iVib = (*st).v();
402 long iRot = (*st).J();
426 for(
long k = 0; k < j; ++k )
451 for(
long ipLo = 0; ipLo < ipHi; ++ipLo )
476 for(
unsigned nEner = 1; nEner <
states.size(); ++nEner )
487 vector<TransitionList::iterator> initptrs;
490 initptrs.resize(
trans.size());
495 for(
unsigned ipHi=1; ipHi<
states.size(); ++ipHi )
497 for(
unsigned ipLo=0; ipLo<ipHi; ++ipLo )
503 initptrs[lineIndex] = tr;
516 for(
long iVibHi=0; iVibHi<=
nVib_hi[iElecHi]; ++iVibHi )
519 for(
long iRotHi=
Jlowest[iElecHi]; iRotHi<=
nRot_hi[iElecHi][iVibHi]; ++iRotHi )
525 long int lim_elec_lo = 0;
526 H2_SaveLine.reserve(iElecHi,iVibHi,iRotHi,lim_elec_lo+1);
527 for(
long iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
530 for(
long iVibLo=0; iVibLo<=
nVib_hi[iElecLo]; ++iVibLo )
552 stable_sort( initptrs.begin(), initptrs.end(),
compareEmis );
557 vector<TransitionList::iterator>::iterator ptr = initptrs.begin();
558 for (
size_t i=0; i < initptrs.size(); ++i, ++tr, ++ptr)
570 for(
unsigned i = 0; i <
trans.size(); ++i )
584 for(
unsigned i = 0; i <
trans.size(); ++i )
596 (*tr).EnergyWN() = (
realnum)((*(*tr).Hi()).energy().WN() - (*(*tr).Lo()).energy().WN());
599 (*tr).WLAng() = (
realnum)(1.e8f/(*tr).EnergyWN() /
RefIndex( (*tr).EnergyWN() ) );
601 (*tr).Coll().col_str() = 0.;
611 (*tr).Emis().iRedisFun() =
ipCRDW;
613 (*tr).Emis().TauIn() =
opac.taumin;
614 (*tr).Emis().TauCon() =
opac.taumin;
616 (*tr).Emis().TauTot() = 1e20f;
618 (*tr).Emis().dampXvel() = (
realnum)( (*tr).Emis().Aul()/(*tr).EnergyWN()/
PI4);
619 (*tr).Emis().gf() = (
realnum)(
GetGF( (*tr).Emis().Aul(),(*tr).EnergyWN(), (*(*tr).Hi()).g() ) );
622 (*tr).Emis().opacity() = (
realnum)(
abscf( (*tr).Emis().gf(), (*tr).EnergyWN(), (*(*tr).Lo()).g()) );
630 (*tr).Emis().ColOvTot() = 1.;
636 (*tr).Emis().ColOvTot() = 0.;
652 (*tr).Coll().col_str() = (
realnum)(
653 pow3( (*tr).WLAng()*1e-8 ) *
654 ( (*(*tr).Hi()).g()/(*(*tr).Lo()).g()) *
657 ASSERT( (*tr).Coll().col_str()>0.);
674 realnum sum = 0., sumj = 0., sumv = 0., sumo = 0., sump = 0.;
677 if(
hmi.chGrainFormPump ==
'D' )
683 double T_H2_FORM = 50000.;
684 for(
long iVib = 0; iVib <=
nVib_hi[0]; ++iVib )
691 (1.f+2.f*
H2_lgOrtho[iElec][iVib][iRot]) * (1.f+iVib) *
708 else if(
hmi.chGrainFormPump ==
'T' )
711 double Xrot[
H2_TOP] = { 0.14 , 0.15 , 0.15 };
712 double Xtrans[
H2_TOP] = { 0.12 , 0.15 , 0.25 };
718 for(
long iVib = 0; iVib <=
nVib_hi[0]; ++iVib )
727 for(
long iVib = 0; iVib <=
nVib_hi[0]; ++iVib )
739 Erot = (EH2 - Ev) * Xrot[
ipH2] / (Xrot[
ipH2] + Xtrans[
ipH2]);
775 double gaussian =
sexp(
POW2( (deltaE - Erot) / (0.5 * Erot) ) );
777 double thermal_dist =
sexp( deltaE / Erot );
780 double aver = ( gaussian + thermal_dist ) / 2.;
791 (1.f+2.f*
H2_lgOrtho[iElec][iVib][iRot]) * Fv * (2.*iRot+1.) * aver );
817 else if(
hmi.chGrainFormPump ==
't' )
824 double T_H2_FORM = 17329.;
826 for(
long iVib = 0; iVib <=
nVib_hi[0]; ++iVib )
854 fprintf(
ioQQQ,
"H2 form grains mean J= %.3f mean v = %.3f ortho/para= %.3f\n",
855 sumj/sum , sumv/sum , sumo/sump );
859 for(
long iVib = 0; iVib <=
nVib_hi[0]; ++iVib )
sys_float sexp(sys_float x)
double RandGauss(double xMean, double s)
NORETURN void TotalInsanity(void)
#define DEBUG_ENTRY(funcname)
TransitionProxy::iterator iterator
void resize(size_t newsize)
multi_arr< realnum, 2 > H2_X_coll_rate
multi_arr< double, 2 > H2_col_rate_out
multi_arr< bool, 3 > H2_lgOrtho
multi_arr< realnum, 3 > H2_X_hminus_formation_distribution
multi_arr< realnum, 3 > H2_dissprob
multi_arr< realnum, 2 > H2_X_formation
valarray< long > ipElec_H2_energy_sort
void H2_ReadDissprob(long int nelec)
multi_arr< realnum, 2 > H2_coll_dissoc_rate_coef
void H2_Read_hminus_distribution(void)
multi_arr< long int, 2 > ipTransitionSort
valarray< long > ipVib_H2_energy_sort
multi_arr< realnum, 3 > CollRateErrFac
multi_arr< realnum, 3 > H2_stat
double H2_DissocEnergies[N_ELEC]
multi_arr< double, 2 > H2_col_rate_in
multi_arr< realnum, 6 > H2_SaveLine
vector< CollRateCoeffArray > RateCoefTable
multi_arr< double, 3 > H2_populations_LTE
valarray< long > ipRot_H2_energy_sort
multi_arr< realnum, 3 > CollRateCoeff
multi_arr< double, 3 > H2_Boltzmann
const double ENERGY_H2_STAR
void Read_Mol_Diss_cross_sections(void)
multi_arr< double, 2 > H2_X_rate_from_elec_excited
multi_arr< int, 2 > H2_ipPhoto
void H2_ReadDissocEnergies(void)
void H2_ReadTransprob(long int nelec, TransitionList &trans)
multi_arr< bool, 2 > lgH2_radiative
multi_arr< double, 2 > pops_per_vib
valarray< realnum > H2_X_source
multi_arr< realnum, 2 > H2_coll_dissoc_rate_coef_H2
multi_arr< realnum, 2 > H2_X_colden_LTE
multi_arr< realnum, 2 > H2_X_colden
multi_arr< realnum, 2 > H2_X_Hmin_back
multi_arr< realnum, 3 > H2_X_grain_formation_distribution
valarray< long > nRot_hi[N_ELEC]
multi_arr< double, 3 > H2_old_populations
multi_arr< double, 3 > H2_rad_rate_out
TransitionList::iterator rad_end
multi_arr< realnum, 3 > H2_disske
multi_arr< double, 2 > H2_rad_rate_in
multi_arr< double, 2 > H2_X_rate_to_elec_excited
long int nLevels_per_elec[N_ELEC]
void H2_CollidRateRead(long int nColl)
valarray< realnum > H2_X_sink
multi_arr< long int, 3 > ipEnergySort
ProxyIterator< qStateConstProxy, qStateConstProxy > const_iterator
ProxyIterator< qStateProxy, qStateConstProxy > iterator
long ipoint(double energy_ryd)
UNUSED const realnum BIGFLOAT
double RefIndex(double EnergyWN)
double abscf(double gf, double enercm, double gl)
double GetGF(double trans_prob, double enercm, double gup)
t_mole_global mole_global
molecule * findspecies(const char buf[])
static const double energy_off
STATIC double EH2_eval(int ipH2, double DissocEnergy, double energy_wn)
static double XVIB[H2_TOP]
int H2_nRot_add_ortho_para[N_ELEC]
STATIC bool compareEmis(const TransitionList::iterator &tr1, const TransitionList::iterator &tr2)
STATIC double H2_vib_dist(int ipH2, double EH2, double DissocEnergy, double energy_wn)
bool compareEnergies(qStateProxy st1, qStateProxy st2)
static double Xdust[H2_TOP]
STATIC bool lgRadiative(const TransitionList::iterator &tr)
UNUSED const double HPLANCK
UNUSED const double SPEEDLIGHT
UNUSED const double FREQ_1EV
UNUSED const double EN1EV
UNUSED const double ATOMIC_MASS_UNIT
UNUSED const double WAVNRYD
vector< TransitionList > AllTransitions