/home66/gary/public_html/cloudy/c13_branch/source/mole_h2_create.cpp

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 */
/*H2_Create create variables for the H2 molecule, called by ContCreatePointers after continuum
 * mesh has been set up */
#include "cddefines.h" 
#include "collision.h" 
#include "physconst.h" 
#include "taulines.h"
#include "lines_service.h"
#include "opacity.h" 
#include "hmi.h" 
#include "ipoint.h"
#include "grainvar.h"
#include "h2.h"
#include "h2_priv.h"
#include "mole.h"
#include "dense.h"

/* this integer is added to rotation quantum number J for the test of whether
 * a particular J state is ortho or para - the state is ortho if J+below is odd,
 * and para if J+below is even */
int H2_nRot_add_ortho_para[N_ELEC] = {0 , 1 , 1 , 0, 1, 1 , 0};

/* if this is set true then code will print energies and stop */
enum {DEBUG_ENER=false};

/* this is equation 8 of Takahashi 2001, clearer definition is given in
 * equation 5 and following discussion of
 * >>refer      H2      formation       Takahashi, J., & Uehara, H., 2001, ApJ, 561, 843-857
 * 0.27eV, convert into wavenumbers */
static double XVIB[H2_TOP] = { 0.70 , 0.60 , 0.20 };
static double Xdust[H2_TOP] = { 0.04 , 0.10 , 0.40 };

/* this is energy difference between bottom of potential well and 0,0
 * the Takahashi energy scale is from the bottom,
 * 2201.9 wavenumbers  */
static const double energy_off = 0.273*FREQ_1EV/SPEEDLIGHT;

STATIC double EH2_eval(  int ipH2, double DissocEnergy, double energy_wn )
{
        double EH2_here;
        double Evm = DissocEnergy * XVIB[ipH2] + energy_off;

        double Ev = (energy_wn+energy_off);
        /* equation 9 of Takahashi 2001 which is only an approximation
         * equation 1, 2 of 
         * Takahashi, Junko, & Uehara, Hideya, 2001, ApJ, 561, 843-857,
         * this is heat deposited on grain by H2 formation in this state */
        double Edust = DissocEnergy * Xdust[ipH2] *
                ( 1. - ( (Ev - Evm) / (DissocEnergy+energy_off-Evm)) *
                ( (1.-Xdust[ipH2])/2.) );
        ASSERT( Edust >= 0. );

        /* energy is total binding energy less energy lost on grain surface 
         * and energy offset */
        EH2_here = DissocEnergy +energy_off - Edust;
        ASSERT( EH2_here >= 0.);

        return EH2_here;
}

/*H2_vib_dist evaluates the vibration distribution for H2 formed on grains */
STATIC double H2_vib_dist( int ipH2 , double EH2, double DissocEnergy, double energy_wn)
{
        double G1[H2_TOP] = { 0.3 , 0.4 , 0.9 };
        double G2[H2_TOP] = { 0.6 , 0.6 , 0.4 };
        double Evm = DissocEnergy * XVIB[ipH2] + energy_off;
        double Fv;
        if( (energy_wn+energy_off) <= Evm )
        {
                /* equation 4 of Takahashi 2001 */
                Fv = sexp( POW2( (energy_wn+energy_off - Evm)/(G1[ipH2]* Evm ) ) );
        }
        else
        {
                /* equation 5 of Takahashi 2001 */
                Fv = sexp( POW2( (energy_wn+energy_off - Evm)/(G2[ipH2]*(EH2 - Evm ) ) ) );
        }
        return Fv;
}

STATIC bool lgRadiative( const TransitionList::iterator &tr )
{
        //quantumState_diatoms *Hi = static_cast<quantumState_diatoms*>( tr.Hi );       
        qList::iterator Lo = (*tr).Lo() ;       
        //if( Lo->n==0 && lgH2_radiative[ ipEnergySort[(*Hi).n][(*Hi).v][(*Hi).J] ][ ipEnergySort[Lo->n][Lo->v][Lo->J] ] )
        if( (*Lo).n()==0 && (*tr).Emis().Aul() > 1.01e-30 )
                return true;
        else
                return false;
}

STATIC bool compareEmis( const TransitionList::iterator& tr1, const TransitionList::iterator &tr2 )
{
        if( lgRadiative( tr1 ) && !lgRadiative( tr2 ) )
                return true;
        else 
                return false; 
}

bool compareEnergies( qStateProxy st1, qStateProxy st2 )
{
        if( st1.energy().Ryd() <= st2.energy().Ryd() )
                return true;
        else 
                return false; 
}

/* create variables for the H2 molecule, called by
 * ContCreatePointers after continuum mesh has been set up */
void diatomics::init(void)
{
        /* this is flag set above - when true h2 code is not executed - this is way to
         * avoid this code when it is not working */
        /* only malloc vectors one time per core load */
        if( lgREAD_DATA || !lgEnabled )
                return;

        DEBUG_ENTRY( "H2_Create()" );

        /* print string if H2 debugging is enabled */
        if( nTRACE )
                fprintf(ioQQQ," H2_Create called in DEBUG mode.\n");
                
        // find the species in the chemistry network    
        sp = findspecies( label.c_str() );
        ASSERT( sp != null_mole );

        // find the same species with the excitation marker     
        string strSpStar = shortlabel + "*";
        sp_star = findspecies( strSpStar.c_str() );
        ASSERT( sp != null_mole );
                
        mass_amu = sp->mole_mass / ATOMIC_MASS_UNIT;
        ASSERT( mass_amu > 0. );

        H2_ReadDissocEnergies();

        /* This does more than just read energies... it actually defines the states */
        H2_ReadEnergies();

        // now sort states into energy order
        fixit(); // this doesn't work!
        //sort( states.begin(), states.end(), compareEnergies );
        
        ASSERT( n_elec_states > 0 );
        /* create a vector of sorted energies */
        ipEnergySort.reserve(n_elec_states);
        for( long iElecHi=0; iElecHi<n_elec_states; ++iElecHi )
        {
                ipEnergySort.reserve(iElecHi,nVib_hi[iElecHi]+1);
                for( long iVibHi = 0; iVibHi <= nVib_hi[iElecHi]; ++iVibHi )
                {
                        ipEnergySort.reserve(iElecHi,iVibHi,nRot_hi[iElecHi][iVibHi]+1);
                }
        }
        ipEnergySort.alloc();

        /* create arrays for energy sorted referencing of e, v, J */
        ipElec_H2_energy_sort.resize( states.size() );
        ipVib_H2_energy_sort.resize( states.size() );
        ipRot_H2_energy_sort.resize( states.size() );
        for( unsigned nEner = 0; nEner < states.size(); ++nEner )
        {
                long iElec = states[nEner].n();
                long iVib = states[nEner].v();
                long iRot = states[nEner].J();
                ipElec_H2_energy_sort[nEner] = iElec;
                ipVib_H2_energy_sort[nEner] = iVib;
                ipRot_H2_energy_sort[nEner] = iRot;
                ipEnergySort[iElec][iVib][iRot] = nEner;
                /* following will print quantum indices and energies */
                /*fprintf(ioQQQ,"%li\t%li\t%li\t%.3e\n", iElec, iVib, iRot, 
                        states[nEner].energy().WN() );*/
        }
        
        H2_stat.alloc( ipEnergySort.clone() );
        H2_lgOrtho.alloc( ipEnergySort.clone() );

        /* set all statistical weights - ours is total statistical weight - 
         * including nuclear spin */
        for( qList::iterator st = states.begin(); st != states.end(); ++st )
        {
                /* unlike atoms, for H2 nuclear spin is taken into account - so the
                 * statistical weight of even and odd J states differ by factor of 3 - see page 166, sec par
                 * >>>refer     H2      H2_stat wght    Shull, J.M., & Beckwith, S., 1982, ARAA, 20, 163-188 */
                if( is_odd( (*st).J() + H2_nRot_add_ortho_para[(*st).n()]) )
                {
                        /* ortho */
                        H2_lgOrtho[(*st).n()][(*st).v()][(*st).J()] = true;
                        H2_stat[(*st).n()][(*st).v()][(*st).J()] = 3.f*(2.f*(*st).J()+1.f);
                }
                else
                {
                        /* para */
                        H2_lgOrtho[(*st).n()][(*st).v()][(*st).J()] = false;
                        H2_stat[(*st).n()][(*st).v()][(*st).J()] = (2.f*(*st).J()+1.f);
                }
                (*st).g() = H2_stat[(*st).n()][(*st).v()][(*st).J()];
        }

        if( nTRACE >= n_trace_full ) 
        {
                for( long iElec=0; iElec<n_elec_states; ++iElec)
                {
                        /* print the number of levels within iElec */
                        fprintf(ioQQQ,"\t(%li %li)", iElec ,  nLevels_per_elec[iElec] );
                }
                fprintf(ioQQQ,
                        " H2_Create: there are %li electronic levels, in each level there are",
                        n_elec_states);
                fprintf(ioQQQ,
                                  " for a total of %li levels.\n", (long int) states.size() );
        }

        /* now find number of levels in H2g */
        for( long nEner=0; nEner<nLevels_per_elec[0]; ++nEner )
        {
                if( states[nEner].energy().WN() > ENERGY_H2_STAR && hmi.lgLeiden_Keep_ipMH2s )
                        break;
                nEner_H2_ground = nEner;
        }
        /* need to increment it so that this is the number of levels, not the index
         * of the highest level */
        ++nEner_H2_ground;

        /* this is the number of levels to do with the matrix - set with the
         * atom h2 matrix command, keyword ALL means to do all of X in the matrix
         * but number of levels within X was not known when the command was parsed,
         * so this was set to -1 to defer setting to all until now */
        if( nXLevelsMatrix<0 )
        {
                nXLevelsMatrix = nLevels_per_elec[0];
        }
        else if( nXLevelsMatrix > nLevels_per_elec[0] )
        {
                fprintf( ioQQQ, 
                        " 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",
                        nXLevelsMatrix ,
                        nLevels_per_elec[0]);
                cdEXIT(EXIT_FAILURE);
        }

        /* at this stage the full electronic, vibration, and rotation energies have been defined,
         * this is an option to print the energies */
        {
                /* set following true to get printout, false to not print energies */
                if( DEBUG_ENER )
                {
                        /* print title for quantum numbers and energies */
                        /*fprintf(ioQQQ,"elec\tvib\trot\tenergy\n");*/
                        for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
                        {
                                fprintf(ioQQQ,"%li\t%li\t%li\t%.5e\n", (*st).n(), (*st).v(), (*st).J(), (*st).energy().WN() );
                        }
                        /* this will exit the program after printing the level energies */
                        cdEXIT(EXIT_SUCCESS);
                }
        }

        /* this will prevent data from being read twice */
        lgREAD_DATA = true;

        /* create space for the electronic levels */
        H2_populations_LTE.reserve(n_elec_states);
        pops_per_vib.reserve(n_elec_states);
        H2_dissprob.reserve(n_elec_states);

        for( long iElec = 0; iElec<n_elec_states; ++iElec )
        {

                if( nTRACE >= n_trace_full )
                        fprintf(ioQQQ,"elec %li highest vib= %li\n", iElec , nVib_hi[iElec] );

                ASSERT( nVib_hi[iElec] > 0 );

                /* nVib_hi is now the highest vibrational level before dissociation,
                 * now allocate space to hold the number of rotation levels */
                H2_populations_LTE.reserve(iElec,nVib_hi[iElec]+1);
                pops_per_vib.reserve(iElec,nVib_hi[iElec]+1);

                /* now loop over all vibrational levels, and find out how many rotation levels there are */
                /* ground is special since use tabulated data - there are 14 vib states,
                 * ivib=14 is highest */
                for( long iVib = 0; iVib <= nVib_hi[iElec]; ++iVib )
                {
                        /* lastly create the space for the rotation quantum number */
                        H2_populations_LTE.reserve(iElec,iVib,nRot_hi[iElec][iVib]+1);
                }
        }

        H2_populations_LTE.alloc();
        H2_old_populations.alloc( H2_populations_LTE.clone() );
        H2_Boltzmann.alloc( H2_populations_LTE.clone() );
        H2_rad_rate_out.alloc( H2_populations_LTE.clone() );
        H2_dissprob.alloc( H2_populations_LTE.clone() );
        H2_disske.alloc( H2_populations_LTE.clone() );
        // this has a different geometry than the above
        pops_per_vib.alloc();

        H2_dissprob.zero();
        H2_disske.zero();

        /* set this one time, will never be set again, but might be printed */
        H2_rad_rate_out.zero();

        /* these do not have electronic levels - all within X */
        H2_ipPhoto.reserve(nVib_hi[0]+1);

        /* space for the vibration levels */
        for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
        {
                /* space for the rotation quantum number */
                H2_ipPhoto.reserve(iVib,nRot_hi[0][iVib]+1);
        }

        H2_ipPhoto.alloc();
        H2_col_rate_in.alloc( H2_ipPhoto.clone() );
        H2_col_rate_out.alloc( H2_ipPhoto.clone() );
        H2_rad_rate_in.alloc( H2_ipPhoto.clone() );
        H2_coll_dissoc_rate_coef.alloc( H2_ipPhoto.clone() );
        H2_coll_dissoc_rate_coef_H2.alloc( H2_ipPhoto.clone() );
        H2_X_colden.alloc( H2_ipPhoto.clone() );
        H2_X_rate_from_elec_excited.alloc( H2_ipPhoto.clone() );
        H2_X_rate_to_elec_excited.alloc( H2_ipPhoto.clone() );
        H2_X_colden_LTE.alloc( H2_ipPhoto.clone() );
        H2_X_formation.alloc( H2_ipPhoto.clone() );
        H2_X_Hmin_back.alloc( H2_ipPhoto.clone() );

        for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
        {
                for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
                {
                        /* >>chng 04 jun 14, set these to bad numbers */
                        H2_rad_rate_in[iVib][iRot] = -BIGFLOAT;
                        H2_coll_dissoc_rate_coef[iVib][iRot] = -BIGFLOAT;
                        H2_coll_dissoc_rate_coef_H2[iVib][iRot] = -BIGFLOAT;
                }
        }
        /* zero out the matrices */
        H2_X_colden.zero();
        H2_X_colden_LTE.zero();
        H2_X_formation.zero();
        H2_X_Hmin_back.zero();
        /* rates [cm-3 s-1] from elec excited states into X only vib and rot */
        H2_X_rate_from_elec_excited.zero();
        /* rates [s-1] to elec excited states from X only vib and rot */
        H2_X_rate_to_elec_excited.zero();

        /* distribution function for populations following formation from H minus H- */
        H2_X_hminus_formation_distribution.reserve(nTE_HMINUS);
        for( long i=0; i<nTE_HMINUS; ++i )
        {
                H2_X_hminus_formation_distribution.reserve(i,nVib_hi[0]+1);
                /* space for the vibration levels */
                for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
                {
                        H2_X_hminus_formation_distribution.reserve(i,iVib,nRot_hi[0][iVib]+1);
                }
        }
        H2_X_hminus_formation_distribution.alloc();
        H2_X_hminus_formation_distribution.zero();
        H2_Read_hminus_distribution();

        /* >>chng 05 jun 20, do not use this, which is highly processed - use ab initio
         * rates of excitation to electronic levels instead */
        /* read in cosmic ray distribution information
        H2_Read_Cosmicray_distribution(); */

        /* grain formation matrix */
        H2_X_grain_formation_distribution.reserve(H2_TOP);
        for( long ipH2=0; ipH2<(int)H2_TOP; ++ipH2 )
        {
                H2_X_grain_formation_distribution.reserve(ipH2,nVib_hi[0]+1);

                /* space for the vibration levels */
                for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
                {
                        H2_X_grain_formation_distribution.reserve(ipH2,iVib,nRot_hi[0][iVib]+1);
                }
        }
        H2_X_grain_formation_distribution.alloc();
        H2_X_grain_formation_distribution.zero();

        for( long iElec=0; iElec<n_elec_states; ++iElec )
        {
                /* get dissociation probabilities and energies - ground state is stable */
                if( iElec > 0 )
                        H2_ReadDissprob(iElec);
        }

        /* >>02 oct 18, add photodissociation, H2 + hnu => 2H + KE */
        /* we now have ro-vib energies, now set up threshold array offsets
         * for photodissociation */
        for( qList::const_iterator st = states.begin(); st != states.end(); ++st )
        {
                long iElec = (*st).n();
                if( iElec > 0 ) continue;
                long iVib = (*st).v();
                long iRot = (*st).J();
                /* this is energy needed to get up to n=3 electronic continuum 
                 * H2 cannot dissociate following absorption of a continuum photon into the
                 * continuum above X (which would require little energy) 
                 * because that process violates momentum conservation 
                 * these would be the triplet states - permitted are into singlets
                 * the effective full wavelength range of this process is from Lya to 
                 * Lyman limit in shielded regions 
                 * tests show limits are between 850A and 1220A - so Lya is included */
                /*>>KEYWORD     Allison & Dalgarno; continuum dissociation; */
                double thresh = (H2_DissocEnergies[1] - (*st).energy().WN())*WAVNRYD;
                /*fprintf(ioQQQ,"DEBUG\t%.2f\t%f\n", RYDLAM/thresh , thresh);*/
                /* in theory we should be able to assert that thesh just barely reaches
                 * lya, but actual numbers reach down to 0.749 ryd */
                ASSERT( thresh > 0.74 );
                H2_ipPhoto[iVib][iRot] = ipoint(thresh);
                fixit(); // this needs to be generalized
        }

        CollRateCoeff.reserve( nLevels_per_elec[0] );
        for( long j = 1; j < nLevels_per_elec[0]; ++j )
        {
                CollRateCoeff.reserve( j, j );
                for( long k = 0; k < j; ++k )
                {
                        CollRateCoeff.reserve( j, k, N_X_COLLIDER );
                }
        }
        CollRateCoeff.alloc();
        CollRateCoeff.zero();

        fixit(); // Does RateCoefTable only need to be N_X_COLLIDER long?
        RateCoefTable.resize(ipNCOLLIDER);
        /* now read in the various sets of collision data */
        for( long nColl=0; nColl<N_X_COLLIDER; ++nColl )
        {
                /* ground state has tabulated data */
                H2_CollidRateRead(nColl);
        }

        CollRateErrFac.alloc( CollRateCoeff.clone() );
        CollRateErrFac.zero();

        /* option to add gaussian random mole */
        if( lgH2_NOISE )
        {
                for( long ipHi = 1; ipHi < nLevels_per_elec[0]; ++ipHi )
                {
                        for( long ipLo = 0; ipLo < ipHi; ++ipLo )
                        {
                                for( long nColl=0; nColl<N_X_COLLIDER; ++nColl )
                                {
                                        /* this returns the log of the random noise */
                                        realnum r = (realnum)RandGauss( xMeanNoise , xSTDNoise );
                                        CollRateErrFac[ipHi][ipLo][nColl] = pow((realnum)10.f,r);
                                }
                        }
                }
        }

        /* this will be total collision rate from an upper to a lower level within X */
        H2_X_source.resize( nLevels_per_elec[0] );
        H2_X_sink.resize( nLevels_per_elec[0] );

        H2_X_coll_rate.reserve(nLevels_per_elec[0]);
        /* now expand out to include all lower levels as lower state */
        for( long i=1; i<nLevels_per_elec[0]; ++i )
        {
                H2_X_coll_rate.reserve(i,i);
        }
        H2_X_coll_rate.alloc();

        ipTransitionSort.reserve( states.size() );
        for( unsigned nEner = 1; nEner < states.size(); ++nEner )
                ipTransitionSort.reserve( nEner, nEner );
        ipTransitionSort.alloc();
        ipTransitionSort.zero();
        lgH2_radiative.alloc( ipTransitionSort.clone() );
        lgH2_radiative.zero();
        
        trans.resize( (states.size() * (states.size()-1))/2 );
        AllTransitions.push_back(trans);
        qList initStates(1);
        TransitionList initlist("H2InitList",&initStates);
        vector<TransitionList::iterator> initptrs;
        initlist.resize(trans.size());
        initlist.states() = &states;
        initptrs.resize(trans.size());

        {
                long lineIndex = 0;
                TransitionList::iterator tr = initlist.begin();
                for( unsigned ipHi=1; ipHi< states.size(); ++ipHi )
                {
                        for( unsigned ipLo=0; ipLo<ipHi; ++ipLo )
                        {
                                (*tr).Junk();
                                (*tr).setHi(ipHi);
                                (*tr).setLo(ipLo);
                                (*tr).Zero();
                                initptrs[lineIndex] = tr;
                                ipTransitionSort[ipHi][ipLo] = lineIndex;
                                lineIndex++;
                                ++tr;
                        }
                }
        }

        /* create the main array of lines */
        H2_SaveLine.reserve(n_elec_states);
        for( long iElecHi=0; iElecHi<n_elec_states; ++iElecHi )
        {
                H2_SaveLine.reserve(iElecHi,nVib_hi[iElecHi]+1);
                for( long iVibHi=0; iVibHi<=nVib_hi[iElecHi]; ++iVibHi )
                {
                        H2_SaveLine.reserve(iElecHi,iVibHi,nRot_hi[iElecHi][iVibHi]+1);
                        for( long iRotHi=Jlowest[iElecHi]; iRotHi<=nRot_hi[iElecHi][iVibHi]; ++iRotHi )
                        {
                                /* now the lower levels */
                                /* NB - X is the only lower level considered here, since we are only 
                                 * concerned with excited electronic levels as a photodissociation process
                                 * code exists to relax this assumption - simply change following to iElecHi */
                                long int lim_elec_lo = 0;
                                H2_SaveLine.reserve(iElecHi,iVibHi,iRotHi,lim_elec_lo+1);
                                for( long iElecLo=0; iElecLo<=lim_elec_lo; ++iElecLo )
                                {
                                        H2_SaveLine.reserve(iElecHi,iVibHi,iRotHi,iElecLo,nVib_hi[iElecLo]+1);
                                        for( long iVibLo=0; iVibLo<=nVib_hi[iElecLo]; ++iVibLo )
                                        {
                                                H2_SaveLine.reserve(iElecHi,iVibHi,iRotHi,iElecLo,iVibLo,nRot_hi[iElecLo][iVibLo]+1);
                                        }
                                }
                        }
                }
        }

        H2_SaveLine.alloc();
        
        /* zero out array used to save emission line intensities */
        H2_SaveLine.zero();

        /* space for the energy vector is now malloced, must read trans probs from table */
        for( long iElec=0; iElec<n_elec_states; ++iElec )
        {
                /* ground state has tabulated data */
                H2_ReadTransprob(iElec,initlist);
        }

        // sort sys.trans so that radiative lines are at the beginning
        stable_sort( initptrs.begin(), initptrs.end(), compareEmis );
        rad_end = trans.begin();
        // rad_end will be used for the end of the range (non-inclusive) for operations on radiative lines
        {
                TransitionList::iterator tr = trans.begin();
                vector<TransitionList::iterator>::iterator ptr = initptrs.begin();
                for (size_t i=0; i < initptrs.size(); ++i, ++tr, ++ptr)
                {
                        (*tr).copy(*(*ptr));
                        if (lgRadiative(tr))
                        {
                                rad_end = tr;
                        }
                }
        }
        ++rad_end;
        ASSERT( rad_end != trans.end() ); 
        // after above sorting, ipTransitionSort is now invalid.  Fix here
        for( unsigned i = 0; i < trans.size(); ++i )
        {
                qList::iterator Hi = trans[i].Hi();
                qList::iterator Lo = trans[i].Lo();
                ipTransitionSort[ ipEnergySort[(*Hi).n()][(*Hi).v()][(*Hi).J()] ][ ipEnergySort[(*Lo).n()][(*Lo).v()][(*Lo).J()] ] = i;
                trans[i].ipHi() = ipEnergySort[(*Hi).n()][(*Hi).v()][(*Hi).J()];
                trans[i].ipLo() = ipEnergySort[(*Lo).n()][(*Lo).v()][(*Lo).J()];
        } 

        findspecieslocal("H2")->levels = &states;
        findspecieslocal("H2")->lines = &trans;

        // after above sorting, ipTransitionSort is now invalid.  Fix here
        for( unsigned i = 0; i < trans.size(); ++i )
        {
                qList::iterator Hi = trans[i].Hi();
                qList::iterator Lo = trans[i].Lo();
                ipTransitionSort[ ipEnergySort[(*Hi).n()][(*Hi).v()][(*Hi).J()] ][ ipEnergySort[(*Lo).n()][(*Lo).v()][(*Lo).J()] ] = i;
                //trans[i].ipHi() = ipEnergySort[(*Hi).n()][(*Hi).v()][(*Hi).J()];
                //trans[i].ipLo() = ipEnergySort[(*Lo).n()][(*Lo).v()][(*Lo).J()];
        }

        // this loop is over all transitions
        for( TransitionList::iterator tr = trans.begin(); tr != trans.end(); ++tr )
        {
                (*tr).EnergyWN() = (realnum)((*(*tr).Hi()).energy().WN() - (*(*tr).Lo()).energy().WN());
                /*wavelength of transition in Angstroms */
                if( (*tr).EnergyWN() > SMALLFLOAT)
                        (*tr).WLAng() = (realnum)(1.e8f/(*tr).EnergyWN() / RefIndex( (*tr).EnergyWN() ) );

                (*tr).Coll().col_str() = 0.;
        }
        
        // this loop is over only radiative transitions.  Notice the end iterator.
        for( TransitionList::iterator tr = trans.begin(); tr != rad_end; ++tr )
        {
                /* line redistribution function - will use complete redistribution */
                /* >>chng 04 mar 26, should include damping wings, especially for electronic
                 * transitions, had used doppler core only */
                (*tr).resetEmis();
                (*tr).Emis().iRedisFun() = ipCRDW;
                /* line optical depths in direction towards source of ionizing radiation */
                (*tr).Emis().TauIn() = opac.taumin;
                (*tr).Emis().TauCon() = opac.taumin;
                /* outward optical depth */
                (*tr).Emis().TauTot() = 1e20f;

                (*tr).Emis().dampXvel() = (realnum)( (*tr).Emis().Aul()/(*tr).EnergyWN()/PI4);
                (*tr).Emis().gf() = (realnum)(GetGF( (*tr).Emis().Aul(),(*tr).EnergyWN(), (*(*tr).Hi()).g() ) );

                /* derive the absorption coefficient, call to function is gf, wl (A), g_low */
                (*tr).Emis().opacity() = (realnum)( abscf( (*tr).Emis().gf(), (*tr).EnergyWN(), (*(*tr).Lo()).g()) );
                
                qList::iterator Hi = (*tr).Hi() ;       

                if( (*Hi).n() == 0 )
                {
                        /* the ground electronic state, most excitations are not direct pumping 
                         * (rather indirect, which does not count for ColOvTot) */
                        (*tr).Emis().ColOvTot() = 1.;
                }
                else
                {
                        /* these are excited electronic states, mostly pumped, except for supras */
                        (*tr).Emis().ColOvTot() = 0.;
                }

                fixit(); // why not include this for excitations within X as well?
                if( (*Hi).n() > 0 )
                {
                        /* cosmic ray and non-thermal suprathermal excitation 
                         * to singlet state of H2 (B,C,B',D)
                         * cross section is equation 5 of 
                         *>>refer       H2      cs      Liu, W. & Dalgarno, A. 1994, ApJ, 428, 769
                         * relative to H I Lya cross section 
                         * this is used to derive H2 electronic excitations
                         * from the H I Lya rate
                         * the following is dimensionless scale factor for excitation 
                         * relative to H I Lya
                         */
                        (*tr).Coll().col_str() = (realnum)(
                                pow3( (*tr).WLAng()*1e-8 ) *
                                ( (*(*tr).Hi()).g()/(*(*tr).Lo()).g()) *
                                (*tr).Emis().Aul() *
                                log(3.)*HPLANCK/(160.f*pow3(PI)*0.5*1e-8*EN1EV)/6.e-17);
                        ASSERT( (*tr).Coll().col_str()>0.);
                }
        }

        /*  Read continuum photodissociation cross section files */
        if( mole_global.lgStancil )
                Read_Mol_Diss_cross_sections();

        /* define branching ratios for deposition of H2 formed on grain surfaces,
         * set true to use Takahashi distribution, false to use Draine & Bertoldi */

        /* loop over all types of grain surfaces */
        /* >>chng 02 oct 08, resolved grain types */
        /* number of different grain types H2_TOP is set in grainvar.h,
         * types are ice, silicate, graphite */
        for( long ipH2=0; ipH2<(int)H2_TOP; ++ipH2 )
        {
                realnum sum = 0., sumj = 0., sumv = 0., sumo = 0., sump = 0.;
                
                /* first is Draine distribution */
                if( hmi.chGrainFormPump == 'D' )
                {
                        long iElec = 0;
                        /* H2 formation temperature, for equation 19, page 271, of
                        * >>refer       H2      formation distribution  Draine, B.T., & Bertoldi, F., 1996, ApJ, 468, 269-289
                        */
                        double T_H2_FORM = 50000.;
                        for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
                        {
                                for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
                                {
                                        /* no distinction between grain surface composition */
                                        H2_X_grain_formation_distribution[ipH2][iVib][iRot] = 
                                                /* first term is nuclear H2_stat weight */
                                                (1.f+2.f*H2_lgOrtho[iElec][iVib][iRot]) * (1.f+iVib) *
                                                (realnum)sexp( states[ ipEnergySort[iElec][iVib][iRot] ].energy().K()/T_H2_FORM );
                                        sum += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
                                        sumj += iRot * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
                                        sumv += iVib * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
                                        if( H2_lgOrtho[iElec][iVib][iRot] )
                                        {
                                                sumo += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
                                        }
                                        else
                                        {
                                                /* >>chng 02 nov 14, [0][iVib][iRot] -> [ipH2][iVib][iRot], PvH */
                                                sump += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
                                        }
                                }
                        }
                }
                else if( hmi.chGrainFormPump == 'T' )
                {
                        /* Takahashi 2001 distribution */
                        double Xrot[H2_TOP] = { 0.14 , 0.15 , 0.15 };
                        double Xtrans[H2_TOP] = { 0.12 , 0.15 , 0.25 };
                        /* first normalize the vibration distribution function */
                        double sumvib = 0.;
                        double EH2;
                        long iElec = 0;

                        for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
                        {
                                double vibdist;
                                EH2 = EH2_eval( ipH2, H2_DissocEnergies[0], states[ ipEnergySort[0][iVib][0] ].energy().WN() );
                                vibdist = H2_vib_dist( ipH2 , EH2, H2_DissocEnergies[0], states[ ipEnergySort[0][iVib][0] ].energy().WN() );
                                sumvib += vibdist;
                        }
                        /* this branch, use distribution function from
                        * >>refer       grain   physics Takahashi, Junko, 2001, ApJ, 561, 254-263 */
                        for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
                        {
                                double Ev = states[ ipEnergySort[iElec][iVib][0] ].energy().WN()+energy_off;
                                double Fv;
                                /* equation 10 of Takahashi 2001, extra term is energy offset between bottom of potential
                                 * the 0,0 level */
                                double Erot;
                                /*fprintf(ioQQQ," Evvvv\t%i\t%li\t%.3e\n", ipH2 ,iVib , Ev*WAVNRYD*EVRYD);*/

                                EH2 = EH2_eval( ipH2, H2_DissocEnergies[0], states[ ipEnergySort[0][iVib][0] ].energy().WN() );

                                /* equation 3 of Taktahashi & Uehara */
                                Erot = (EH2 - Ev) * Xrot[ipH2] / (Xrot[ipH2] + Xtrans[ipH2]);

                                /* email exchange with Junko Takahashi - 
                                Thank you for your E-mail.
                                I did not intend to generate negative Erot.
                                I cut off the populations if their energy levels are negative, and made the total
                                population be unity by using normalization factors (see, e.g., Eq. 12).

                                I hope that my answer is of help to you and your work is going well.
                                With best wishes,
                                Junko

                                >Thanks for the reply.  By cutting off the population, should we set the
                                >population to zero when Erot becomes negative, or should we set Erot to
                                >a small positive number? 

                                I just set the population to zero when Erot becomes negative.
                                Our model is still a rough one for the vibration-rotation distribution function
                                of H2 newly formed on dust, because we have not yet had any exact
                                experimental or theoretical data about it.
                                With best wishes,
                                Junko

                                 */

                                if( Erot > 0. )
                                {
                                        /* the vibrational distribution */
                                        Fv = H2_vib_dist( ipH2 , EH2, H2_DissocEnergies[0], states[ ipEnergySort[0][iVib][0] ].energy().WN() ) / sumvib;
                                        /*fprintf(ioQQQ," vibbb\t%li\t%.3e\n", iVib , Fv );*/

                                        for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
                                        {
                                                double deltaE = states[ ipEnergySort[iElec][iVib][iRot] ].energy().WN() -
                                                        states[ ipEnergySort[iElec][iVib][0] ].energy().WN();
                                                /* equation 6 of Takahashi 2001 */
                                                double gaussian = sexp( POW2( (deltaE - Erot) / (0.5 * Erot) ) );
                                                /* equation 7 of Takahashi 2001 */
                                                double thermal_dist = sexp( deltaE / Erot );

                                                /* take the mean of the two */
                                                double aver = ( gaussian + thermal_dist ) / 2.;
                                                /*fprintf(ioQQQ,"rottt\t%i\t%li\t%li\t%.3e\t%.3e\t%.3e\t%.3e\n",
                                                        ipH2,iVib,iRot,
                                                        deltaE*WAVNRYD*EVRYD,
                                                        gaussian, thermal_dist , aver );*/

                                                /* thermal_dist does become > 1 since Erot can become negative */
                                                ASSERT( gaussian <= 1. /*&& thermal_dist <= 10.*/ );

                                                H2_X_grain_formation_distribution[ipH2][iVib][iRot] = (realnum)(
                                                        /* first term is nuclear H2_stat weight */
                                                        (1.f+2.f*H2_lgOrtho[iElec][iVib][iRot]) * Fv * (2.*iRot+1.) * aver );

                                                sum += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
                                                sumj += iRot * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
                                                sumv += iVib * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
                                                if( H2_lgOrtho[iElec][iVib][iRot] )
                                                {
                                                        sumo += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
                                                }
                                                else
                                                {
                                                        sump += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
                                                }

                                        }
                                }
                                else
                                {
                                        /* this branch Erot is non-positive, so no distribution */
                                        for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
                                        {
                                                H2_X_grain_formation_distribution[ipH2][iVib][iRot] = 0.;
                                        }
                                }
                        }
                }
                else if( hmi.chGrainFormPump == 't' )
                { 
                        /* thermal distribution at 1.5 eV, as suggested by Amiel & Jaques */
                        /* thermal distribution, upper right column of page 239 of
                         *>>refer       H2      formation       Le Bourlot, J, 1991, A&A, 242, 235 
                         * set with command
                         * set h2 grain formation pumping thermal */
                        double T_H2_FORM = 17329.;
                        long iElec = 0;
                        for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
                        {
                                for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
                                {
                                        /* no distinction between grain surface composition */
                                        H2_X_grain_formation_distribution[ipH2][iVib][iRot] = 
                                                /* first term is nuclear H2_stat weight */
                                                H2_stat[0][iVib][iRot] *
                                                (realnum)sexp( states[ ipEnergySort[0][iVib][iRot] ].energy().K()/T_H2_FORM );
                                        sum += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
                                        sumj += iRot * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
                                        sumv += iVib * H2_X_grain_formation_distribution[ipH2][iVib][iRot];
                                        if( H2_lgOrtho[iElec][iVib][iRot] )
                                        {
                                                sumo += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
                                        }
                                        else
                                        {
                                                /* >>chng 02 nov 14, [0][iVib][iRot] -> [ipH2][iVib][iRot], PvH */
                                                sump += H2_X_grain_formation_distribution[ipH2][iVib][iRot];
                                        }
                                }
                        }
                }
                else
                        TotalInsanity();

                if( nTRACE >= n_trace_full )
                        fprintf(ioQQQ, "H2 form grains mean J= %.3f mean v = %.3f ortho/para= %.3f\n", 
                                sumj/sum , sumv/sum , sumo/sump );

                //iElec = 0;
                /* now rescale so that integral is unity */
                for( long iVib = 0; iVib <= nVib_hi[0]; ++iVib )
                {
                        for( long iRot=Jlowest[0]; iRot<=nRot_hi[0][iVib]; ++iRot )
                        {
                                H2_X_grain_formation_distribution[ipH2][iVib][iRot] /= sum;
                                /* print the distribution function */
                                /*if( states[ ipEnergySort[iElec][iVib][iRot] ].energy().WN() < 5200. )
                                fprintf(ioQQQ,"disttt\t%i\t%li\t%li\t%li\t%.4e\t%.4e\t%.4e\t%.4e\n",
                                        ipH2, iVib , iRot, (long)H2_stat[0][iVib][iRot] ,
                                        states[ ipEnergySort[iElec][iVib][iRot] ].energy().WN(), 
                                        states[ ipEnergySort[iElec][iVib][iRot] ].energy().K(), 
                                        H2_X_grain_formation_distribution[ipH2][iVib][iRot],
                                        H2_X_grain_formation_distribution[ipH2][iVib][iRot]/H2_stat[0][iVib][iRot]
                                        );*/
                        }
                }
        }

        return;
}

---