h2_priv.h

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2022 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */

#ifndef H2_PRIV_H_
#define H2_PRIV_H_

#include "transition.h"
#include "container_classes.h"

typedef void (*linefunc)(const TransitionProxy& t, 
                                                                 bool lgShield_this_zone,
                                                                 realnum pestrk,
                                                                 realnum DopplerWidth);
class Parser;
class molecule;
class CollRateCoeffArray;

const int N_X_COLLIDER = 5;
const int chN_X_COLLIDER = 10;
// colliders are (see h2.cpp diatoms_init)
// 0 H
// 1 He
// 2 H2 ortho
// 3 H2 para
// 4 H+

const int nTE_HMINUS = 7;
                
const int N_ELEC = 7;

const realnum H2_logte_hminus[nTE_HMINUS] = {1.,1.47712,2.,2.47712,3.,3.47712,4.};

struct diss_level 
{
        long n, v, j;
};

class diss_tran 
{
public:
        explicit diss_tran( diss_level a, diss_level b )
        {
                initial = a;
                final = b;
                energies.clear();
                xsections.clear();
                rate_coeff = 0.;
        };
        diss_level initial, final;
        vector<double> energies;        
        vector<double> xsections;
        double rate_coeff;      
};

struct t_coll_source
{
        t_coll_source()
        {
                magic = 0;
                filename = "";
        };
        long magic;
        string filename;
};

class diatomics
{
public:
        double Abund() const
        {
                return *dense_total;
        }
        void GetIndices( long& ipHi, long& ipLo, const char* chLine, long& i ) const;
        void CalcPhotoionizationRate(void);
        double (*photoion_opacity_fun)( double energy );
        long OpacityCreate( double *stack );
        double GetHeatRate( const diss_tran& tran );
        double GetDissociationRate( const diss_tran& tran );
        double MolDissocOpacity( const diss_tran& tran, const double& Mol_Ene );
        double Cont_Diss_Heat_Rate( void );
        void Mol_Photo_Diss_Rates( void );
        void Read_Mol_Diss_cross_sections(void);
        void SolveExcitedElectronicLevels(void);
        void SolveSomeGroundElectronicLevels(void);
        double GetExcitedElecDensity(void);
        realnum GetXColden( long iVib, long iRot );

        long int getLine( long iElecHi, long iVibHi, long iRotHi, long iElecLo, long iVibLo, long iRotLo, double *relint, double *absint );

        /* compute rate coefficient for a single quenching collision */
        realnum H2_CollidRateEvalOne( long iVibHi, long iRotHi, long iVibLo, long iRotLo, long ipHi, long ipLo, long nColl, double temp_K );

        void H2_Calc_Average_Rates( void );

        void H2_X_sink_and_source( void );
        /*H2_X_coll_rate_evaluate find collisional rates within X - 
         * this is one time upon entry into H2_LevelPops */
        void H2_X_coll_rate_evaluate( void );

        /*H2_Level_low_matrix evaluate lower populations within X */
        /* total abundance within matrix */
        void H2_Level_low_matrix(realnum abundance );

        /*H2_Read_Cosmicray_distribution read distribution function for H2 population following cosmic ray collisional excitation
        void H2_Read_Cosmicray_distribution(void); */

        // read energies for all electronic levels
        void H2_ReadEnergies();
        void H2_ReadEnergies( long int nelec, vector<int>& n, vector<int>& v, vector<int>&J, vector<double>& eWN );

        void H2_ReadDissprob( long int nelec );

        void H2_CollidRateEvalAll( void );

        void H2_CollidRateRead( long int nColl );

        void H2_ReadTransprob( long int nelec, TransitionList &trans );

        void H2_Read_hminus_distribution(void);

        void mole_H2_form( void );

        void mole_H2_LTE( void );

        void H2_Solomon_rate( void );

        double gs_rate( void );

        void H2_zero_pops_too_low( void );

        void init(void);

        void H2_ContPoint( void );

        double H2_DR(void);

        double H2_Accel(void);

        void H2_RT_OTS( void );

        double H2_RadPress(void);

        void H2_LinesAdd(void);

        void H2_Reset( void );

        double H2_InterEnergy(void);

        void H2_Colden( const char *chLabel );

        void H2_Cooling(void);

        double LTE_Cooling_per_H2();

        void H2_Punch_line_data(
                FILE* ioPUN ,
                bool lgDoAll );

        void H2_PunchLineStuff( FILE * io , realnum xLimit  , long index);

        void H2_RT_diffuse(void);

        void H2_RTMake( linefunc line_one );

        void H2_RT_tau_inc(void);

        void H2_Prt_Zone(void);

        // print departure coefficients for all X levels
        void H2_PrtDepartCoef(void);

        void H2_LineZero( void );

        void H2_RT_tau_reset( void );

        void H2_LevelPops( bool &lgPopsConverged, double &old_value, double &new_value );

        void H2_PunchDo( FILE* io , char chJOB[] , const char chTime[] , long int ipPun );

        void H2_ParseSave( Parser &p, ostringstream& chHeader );

        double H2_itrzn( void );

        void H2_Prt_column_density( FILE *ioMEAN );

        void set_numLevelsMatrix( long numLevels );

        void H2_ReadDissocEnergies( void );

        bool lgREAD_DATA;

        double photoionize_rate;
        double photo_heat_soft;
        double photo_heat_hard;
        double photodissoc_BigH2_H2s;
        double photodissoc_BigH2_H2g;

        /* total spontaneous dissociation rate [s-1],
         * summed over excited electronic states, weighted by pops */ 
        double spon_diss_tot;

        double Solomon_dissoc_rate_g;
        double Solomon_dissoc_rate_s;
        
        double Solomon_elec_decay_g;
        double Solomon_elec_decay_s;

        double rate_grain_op_conserve;
        double rate_grain_J1_to_J0;

        double Cont_Dissoc_Rate_H2s;
        double Cont_Dissoc_Rate_H2g;
        multi_arr<double,3> Cont_Dissoc_Rate;
        
        double rel_pop_LTE_g;
        double rel_pop_LTE_s;

        double average_energy_g;
        double average_energy_s;

        double HeatDiss;
        double HeatDexc;
        double HeatDexc_old;
        double HeatDexc_deriv;
        double HeatChangeOld, HeatChange;

        double Average_A;
        double Average_collH2_deexcit;
        double Average_collH_deexcit;
        double Average_collH2_excit;
        double Average_collH_excit;
        double Average_collH_dissoc_g;
        double Average_collH_dissoc_s;
        double Average_collH2_dissoc_g;
        double Average_collH2_dissoc_s;

        bool lgEvaluated;

        long ip_photo_opac_thresh;
        long ip_photo_opac_offset;

        t_coll_source coll_source[N_X_COLLIDER];

        double ortho_density,
                para_density;

        // single precision versions of the above       
        realnum ortho_density_f,
                para_density_f;

        double ortho_colden ,
                para_colden;

        /* old and older ortho - para ratios, used to determine whether soln is converged */
        double ortho_para_old, ortho_para_older, ortho_para_current;

        // these remember the largest and smallest factors needed to
        // renormalize the H2 chemistry 
        double renorm_max ,
                renorm_min;

        // this will say how many times the large H2 molecule has been called in this zone -
        // if not called (due to low H2 abundance) then not need to update its line arrays 
        long int nCall_this_zone;

        // flag saying whether to bother with the large molecule at all,
        // default is false, set true with atom h2 on command 
        bool lgEnabled;

        // this is the number of electronic levels to include in the output - default is 1,
        // only X.  changed with PRINT LINES H2 ELECTRONIC and option on PUNCH H2 LINES commands 
        int nElecLevelOutput;

        /* true to use 2007 set of H2 - H collision rate, false use 1999 */
        bool lgH2_H_coll_07;

        bool lgColl_gbar;

        bool lgColl_deexec_Calc;

        bool lgColl_dissoc_coll;

        // include collision rates that come from real calculations,
        // off with atom h2 collisions off command 
        bool lgH2_grain_deexcitation;

        bool lgLTE;

        bool lgH2_ortho_para_coll_on;

        // which set of He - H2 collisions to use? default is ORNL, other
        // is Le BOURlet 
        bool lgH2_He_ORNL;

        // flag saying whether (true) or not to use ORNL H2 - H2 collisions
        bool lgH2_ORH2_ORNL;
        bool lgH2_PAH2_ORNL;

        bool lgH2_NOISE;
        bool lgH2_NOISECOSMIC;
                        
        long int loop_h2_oscil;
        long int nzoneEval;
        
        double xMeanNoise , xSTDNoise;

        double H2_to_H_limit;

        realnum mass_amu;

        int nTRACE;
        
        int n_trace_final , 
                n_trace_iterations , 
                n_trace_full,
                n_trace_matrix;

        // the number of electronic quantum states to include.
        // To do both Lyman and Werner bands want nelec = 3 
        long int n_elec_states;
                
        /* this is fraction of population that is within levels done with matrix */
        double frac_matrix;

        /* used to recall the temperature used for last set of Boltzmann factors */
        double TeUsedBoltz;
        double TeUsedColl;
        
        explicit diatomics( const string& a, const double& e_star, const double* const abund, double (*fun)(double) ) ;

        molecule *sp;   
        molecule *sp_star;
        qList states;   
        TransitionList trans;   
        TransitionList::iterator rad_end;
        vector< diss_tran > Diss_Trans;

private:
        string label;
        string shortlabel;
        string path;

        /* >> chng 05 jul 15, TE, H2g = sum (v=0, J=0,1) */
        /* >>chng 05 jul 29, to 0.5 eV, this goes up to J=8 for v=0 */
        /* >>chng 05 aug 03, slight upward change in energy to include the J=8 level,
         * also give energy in waveumbers for simplicity (save h2 levels give energy in ryd) */
        /*#define       ENERGY_H2_STAR  (0.5/EVRYD/WAVNRYD)*/
        /* energy of v=0, J=8 is 4051.73, J=9 is 5001.97
         * v=1, J=0 is 4161.14 */
public:
        const double ENERGY_H2_STAR;

private:
        // pointer to the density of the species        
        const double* const dense_total;        
        
        char chH2ColliderLabels[N_X_COLLIDER][chN_X_COLLIDER];

        /* these vars are private for H2 but uses same style as all other header files -
         * the extern is extern in all except cddefines */

        long int nEner_H2_ground;

        multi_arr<double,2> pops_per_vib;

        double H2_renorm_chemistry;

        multi_arr<realnum,2> H2_X_coll_rate;

        //int H2_nRot_add_ortho_para[N_ELEC];
        double H2_DissocEnergies[N_ELEC];
        long int nVib_hi[N_ELEC];
        valarray<long> nRot_hi[N_ELEC];
        long int Jlowest[N_ELEC];
        long int nLevels_per_elec[N_ELEC];
        double pops_per_elec[N_ELEC];
        multi_arr<realnum,3> CollRateCoeff;
        multi_arr<realnum,3> CollRateErrFac; 
        vector<CollRateCoeffArray> RateCoefTable;

        // quantities dealing with chemistry 
#if 1
#endif

        // these are quantities for each state
#if 1

        multi_arr<realnum,3> H2_dissprob;
        multi_arr<realnum,3> H2_disske;
        multi_arr<double,3> H2_rad_rate_out;

        multi_arr<double,3> H2_old_populations;
        multi_arr<double,3> H2_populations_LTE;
        multi_arr<bool,3> H2_lgOrtho;
#endif

        long int nzoneAsEval , iterationAsEval;

        // these are quantities for states with X
#if 1
        multi_arr<int,2> H2_ipPhoto;
        multi_arr<double,2> H2_col_rate_in;
        multi_arr<double,2> H2_col_rate_out;
        multi_arr<double,2> H2_rad_rate_in;
        multi_arr<realnum,2> H2_X_formation;
        multi_arr<realnum,2> H2_X_Hmin_back;
        multi_arr<realnum,2> H2_X_colden;
        multi_arr<realnum,2> H2_X_colden_LTE;
        multi_arr<double,2> H2_X_rate_from_elec_excited;
        multi_arr<double,2> H2_X_rate_to_elec_excited;
        multi_arr<realnum,2> H2_coll_dissoc_rate_coef;

        multi_arr<realnum,2> H2_coll_dissoc_rate_coef_H2;
#endif

        valarray<realnum> H2_X_source;
        valarray<realnum> H2_X_sink;

        multi_arr<realnum,3> H2_X_grain_formation_distribution;

        double H2_den_s , H2_den_g;

        multi_arr<realnum,3> H2_X_hminus_formation_distribution;

        valarray<long> ipVib_H2_energy_sort;
        valarray<long> ipElec_H2_energy_sort;
        valarray<long> ipRot_H2_energy_sort;
        multi_arr<long int,3> ipEnergySort;
        multi_arr<long int,2> ipTransitionSort;

        long int nXLevelsMatrix;
        long int ndimMalloced;
        double **AulEscp,
                **AulDest, 
                **AulPump,
                **CollRate_levn;
        vector<double> pops, create, destroy, depart, stat_levn, excit;
        
        long int levelAsEval;
        bool lgFirst;
        long int nzone_eval;
        long int iteration_evaluated;

        multi_arr<realnum,6> H2_SaveLine;

        multi_arr<bool,2> lgH2_radiative;

        long int nH2_pops;
        long int nH2_zone;

        long int nzone_nlevel_set;

        long int nCall_this_iteration;

        /* LTE cooling per molecule */
        vector<double> LTE_Temp, LTE_cool;

        bool lgPrtMatrix;

public:
        /* Read LTE cooling per molecule */
        void H2_Read_LTE_cooling_per_H2();

        /* interpolate over LTE cooling array */
        double interpolate_LTE_Cooling( double Temp );
};

/* compute H2 continuum dissociation cross sections */
double MolDissocCrossSection( const diss_tran& tran, const double& Mol_Ene );

double Yan_H2_CS( double energy_ryd /* photon energy in ryd */);

/* compute H2 continuum dissoication opacities */
//double MolDissocOpacity( const diss_tran& tran, const double& Mol_Ene );

#endif /* H2_PRIV_H_ */