h2_priv.h

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2013 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 "parser.h"
#include "mole.h"

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);

        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( void );

        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_Prt_line_tau(void);

        void H2_ParseSave( Parser &p ,
                                           char *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) ) : trans(a, &states), ENERGY_H2_STAR (e_star), dense_total(abund)
        {
                fixit(); //should probably force path lower and label upper case.
                path = a;
                label = a;
                {
                        unsigned j;
                        for( j = 0; j < a.size(); ++j )
                                label[j] = toupper( label[j] );
                        shortlabel = label;
                        for( j = a.size(); j < 4; ++j )
                                label += ' ';
                        label += '\0';
                }
                /* option to turn off or on gbar collisions of the collision rate,
                * default is to have it on */
                /* turn h2.lgColl_gbar on/off with atom h2 gbar on off */
                lgColl_gbar = true;
                /* include collision rates that come from real calculations,
                * off with atom h2 collisions off command */
                lgColl_deexec_Calc = true;
                lgColl_dissoc_coll = true;
                lgEnabled = false;
                lgEvaluated = false;
                /* option to turn off H2 - grain collision & deexcitation,
                * atom h2 grain collision on/off */
                lgH2_grain_deexcitation = false;
                /* option to scramble collision data */
                lgH2_NOISE = false;
                lgH2_NOISECOSMIC = false; 
                /* option to turn off ortho-para collisions, command ATOM H2 COLLISIONS ORTHO PARA OFF */
                lgH2_ortho_para_coll_on = true;
                /* which set of H2 - He collision data to use?  default is ORNL data set
                 * also Le Bourlot Meudon set available,
                 * set to ORNL with command
                 * atom H2 He collisions ORNL */
                lgH2_He_ORNL = true;
                /* use 1999 or 2007 H2 - H collision set atomic data mole H2 h */
                lgH2_H_coll_07 = true;
                /*>>chng 08 feb 27, GS, ORNL or Meudon H2 - H2 collision data 
                 * >>chng 09 may 11, make it the default */
                lgH2_ORH2_ORNL = true;
                lgH2_PAH2_ORNL = true;
                /* flag to force using LTE level populations */
                lgLTE = false;
                lgREAD_DATA = false;
                loop_h2_oscil = -1;
                HeatDiss = 0.;
                HeatDexc = 0.;
                HeatDexc_old = 0.;
                HeatDexc_deriv = 0.;
                HeatChangeOld = 0.;
                HeatChange = 0.;
                photo_heat_soft = 0.;
                photo_heat_hard = 0.;
                photodissoc_BigH2_H2s = 0.;
                photodissoc_BigH2_H2g = 0.;
                photoion_opacity_fun = fun;
                Solomon_dissoc_rate_s = 0.;
                Solomon_dissoc_rate_g = 0.;
                spon_diss_tot = 0.;
                rate_grain_op_conserve = 0.;
                rate_grain_J1_to_J0 = 0.;
                Average_A = 0.;
                Average_collH2_deexcit = 0.;
                Average_collH_deexcit = 0.;
                Average_collH2_excit = 0.;
                Average_collH_excit = 0.;
                Average_collH_dissoc_g = 0.;
                Average_collH_dissoc_s = 0.;
                Average_collH2_dissoc_g = 0.;
                Average_collH2_dissoc_s = 0.;

                /* these remember the largest and smallest factors needed to
                 * renormalize the H2 chemistry */
                renorm_max = 1.;
                renorm_min = 1.;
                /* ortho and para column densities */
                ortho_colden = 0.;
                para_colden = 0.;
                ortho_para_old = 0.;
                ortho_para_older = 0.;
                ortho_para_current = 0.;
                TeUsedBoltz = -1.;
                TeUsedColl = -1.;

                /* counters used by H2_itrzn to find number of calls of h2 per zone */
                nH2_pops  = 0;
                nH2_zone = 0;
                /* these are used to set trace levels of output by options on atom h2 trace command 
                 * first is minimum level of trace, keyword is FINAL */
                n_trace_final = 1;
                /* intermediate level of trace, info per iteration, key ITERATION */
                n_trace_iterations = 2;
                /* full trace, keyword is FULL */
                n_trace_full = 3;
                /* print matrices used in solving X */
                n_trace_matrix = 4;
                /* holds options for trace set with atom h2 command */
                nTRACE = false;
                /* this is number of electronic levels to include in the print and save output 
                * changed with the PRINT LINE H2 ELECTRONIC and PUNCH H2 commands 
                * by default only include X */
                nElecLevelOutput = 1;
                /* the number of electronic quantum states to include.
                 * To do both Lyman and Werner bands want nelec = 3,
                 * default is to do all bands included */
                n_elec_states = N_ELEC;
                nCall_this_zone = 0;
                /* the number of levels used in the matrix solution
                 * of the level populations - set with atom h2 matrix nlevel,
                 * >>chng 04 oct 05, make default 30 levels 
                 * >>chng 04 dec 23, make default 70 levels */
                nXLevelsMatrix = 70;
                ndimMalloced = 0;

                levelAsEval = -1;
                lgFirst = true;
                nzone_eval = -1;
                nzoneEval = -1;
                iteration_evaluated = -1;

                /* this is used to establish zone number for evaluation of number of levels in matrix */
                nzone_nlevel_set = -1;
                nzoneAsEval = -1;
                /* this is used to establish zone number for evaluation of number of levels in matrix */
                nzone_nlevel_set = 0;

                /* this is the smallest ratio of H2/H where we will bother with the large H2 molecule
                 * this value was chosen so that large mole used at very start of TH85 standard PDR,
                 * NB - this appears in headinfo and must be updated there if changed here */
                /* >>chng 03 jun 02, from 1e-6 to 1e-8 - in orion veil large H2 turned on half way
                 * across, and Solomon process was very fast since all lines optically thin.  correct
                 * result has some shielding, so reset to lower value so that H2 comes on sooner. */
                H2_to_H_limit = 1e-8;
                iterationAsEval = -1; 

                strcpy( chH2ColliderLabels[0] , "H0" );
                strcpy( chH2ColliderLabels[1] , "He" );
                strcpy( chH2ColliderLabels[2] , "H2 o" );
                strcpy( chH2ColliderLabels[3] , "H2 p" );
                strcpy( chH2ColliderLabels[4] , "H+" );
        };

        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_Boltzmann;
        multi_arr<double,3> H2_populations_LTE;
        multi_arr<realnum,3> H2_stat;
        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,
                **col_str,
                **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;

};

/* 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_ */