/home66/gary/public_html/cloudy/c13_branch/source/species2.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 */
#include "cddefines.h"
#include "atmdat.h"
#include "phycon.h"
#include "taulines.h"
#include "mole.h"
#include "mole_priv.h"
#include "atoms.h"
#include "string.h"
#include "thirdparty.h"
#include "dense.h"
#include "conv.h"
#include "h2.h"
#include "physconst.h"
#include "secondaries.h"
#include "thermal.h"
#include "cooling.h"
#include "lines_service.h"
#include "atmdat.h"

STATIC double LeidenCollRate(long, long, const TransitionProxy& ,double);
STATIC double StoutCollRate(long ipSpecies, long ipCollider, const TransitionProxy&, double ftemp);
STATIC double ChiantiCollRate(long ipSpecies, long ipCollider, const TransitionProxy&, double ftemp);

static const bool DEBUGSTATE = false;

static double *g, *ex, *pops, *depart, *source, *sink;
static double **AulEscp, **col_str, **AulDest, **AulPump, **CollRate;

/*Solving for the level populations*/

void dBase_solve(void )
{
        realnum abund;
        DEBUG_ENTRY( "dBase_solve()" );
        static bool lgFirstPass = true;
        static long maxNumLevels = 1;
        double totalHeating = 0.;

        if( nSpecies==0 )
                return;

        if( lgFirstPass )
        {
                for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
                        maxNumLevels = MAX2( maxNumLevels, dBaseSpecies[ipSpecies].numLevels_max );

                g = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
                ex = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
                pops = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
                depart = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
                source = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));
                sink = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));

                AulEscp = (double **)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double *)); 
                col_str = (double **)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double *)); 
                AulDest = (double **)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double *)); 
                AulPump = (double **)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double *));  
                CollRate = (double **)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double *)); 

                for( long j=0; j< maxNumLevels; j++ )
                {
                        AulEscp[j] = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double)); 
                        col_str[j] = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double)); 
                        AulDest[j] = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double)); 
                        AulPump[j] = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double));  
                        CollRate[j] = (double *)MALLOC( (unsigned long)(maxNumLevels)*sizeof(double)); 
                }

                lgFirstPass = false;
        }

        // zero all of these values
        memset( g, 0, (unsigned long)(maxNumLevels)*sizeof(double) );
        memset( ex, 0, (unsigned long)(maxNumLevels)*sizeof(double) );
        memset( pops, 0, (unsigned long)(maxNumLevels)*sizeof(double) );
        memset( depart, 0, (unsigned long)(maxNumLevels)*sizeof(double) );
        memset( source, 0, (unsigned long)(maxNumLevels)*sizeof(double) );
        memset( sink, 0, (unsigned long)(maxNumLevels)*sizeof(double) );
        for( long j=0; j< maxNumLevels; j++ )
        {
                memset( AulEscp[j], 0, (unsigned long)(maxNumLevels)*sizeof(double) );
                memset( col_str[j], 0, (unsigned long)(maxNumLevels)*sizeof(double) );
                memset( AulDest[j], 0, (unsigned long)(maxNumLevels)*sizeof(double) );
                memset( AulPump[j], 0, (unsigned long)(maxNumLevels)*sizeof(double) );
                memset( CollRate[j], 0, (unsigned long)(maxNumLevels)*sizeof(double) );
        }


        for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
        {
                const char *spName = dBaseSpecies[ipSpecies].chLabel;
                double cooltl, coolder;
                int nNegPop;
                bool lgZeroPop, lgDeBug = false;

                dBaseSpecies[ipSpecies].CoolTotal = 0.;

#if 0
                //limit for now to small number of levels
                dBaseSpecies[ipSpecies].numLevels_local = MIN2( dBaseSpecies[ipSpecies].numLevels_local, 10 );
#endif

                /* first find current density (cm-3) of species */
                if( dBaseSpecies[ipSpecies].lgMolecular )
                {
                        molezone *SpeciesCurrent;
                        if( (SpeciesCurrent = findspecieslocal(dBaseSpecies[ipSpecies].chLabel)) == null_molezone )
                        {
                                /* did not find the species - print warning for now */
                                if( !conv.nTotalIoniz )
                                        fprintf(ioQQQ," PROBLEM dBase_solve did not find molecular species %li\n",ipSpecies);
                        }
                        abund = (realnum)SpeciesCurrent->den;
                }
                else
                {
                        /* an atom or ion */
                        ASSERT( dBaseStates[ipSpecies][0].nelem()<=LIMELM && dBaseStates[ipSpecies][0].IonStg()<=dBaseStates[ipSpecies][0].nelem()+1 );
                        abund = dense.xIonDense[ dBaseStates[ipSpecies][0].nelem()-1 ][ dBaseStates[ipSpecies][0].IonStg()-1 ];
                }

                abund *= dBaseSpecies[ipSpecies].fracType * dBaseSpecies[ipSpecies].fracIsotopologue;

                // initialization at start of each iteration
                if( conv.nTotalIoniz == 0)
                        dBaseSpecies[ipSpecies].lgActive = true;

                bool lgMakeInActive = (abund <= 1e-20 * dense.xNucleiTotal);
                if( lgMakeInActive && dBaseSpecies[ipSpecies].lgActive )
                {
                        // zero out populations and intensities, if previously not set
                        dBaseStates[ipSpecies][0].Pop() = 0.;
                        for(long ipHi = 1; ipHi < dBaseSpecies[ipSpecies].numLevels_max; ipHi++ )
                        {       
                                dBaseStates[ipSpecies][ipHi].Pop() = 0.;
# ifdef USE_NLTE7
                                dBaseStates[ipSpecies][ipHi].DestCollBB() = 0.;
                                dBaseStates[ipSpecies][ipHi].DestPhotoBB() = 0.;
# endif
                        }
                        for (TransitionList::iterator tr=dBaseTrans[ipSpecies].begin(); 
                                  tr != dBaseTrans[ipSpecies].end(); ++tr)
                        {
                                (*tr).Emis().phots() = 0.;
                                (*tr).Emis().xIntensity() = 0.;
                                (*tr).Coll().col_str() = 0.;
                                (*tr).Coll().cool() = 0.;
                                (*tr).Coll().heat() = 0.;
                        }
                        dBaseSpecies[ipSpecies].lgActive = false;
                }

                if( !lgMakeInActive )
                        dBaseSpecies[ipSpecies].lgActive = true;

                if( !dBaseSpecies[ipSpecies].lgActive )
                        continue;

                // we always hit search phase first, reset number of levels
                if( conv.lgSearch )
                        dBaseSpecies[ipSpecies].numLevels_local = dBaseSpecies[ipSpecies].numLevels_max;
                for( long ipLo = 0; ipLo < dBaseSpecies[ipSpecies].numLevels_local; ipLo++ )
                {
                        /* statistical weights & Excitation Energies*/
                        g[ipLo] = dBaseStates[ipSpecies][ipLo].g() ;
                        // parts of the code assert that ground is at zero energy - this is
                        // not true for the stored molecular data - so rescale to zero
                        ex[ipLo] = dBaseStates[ipSpecies][ipLo].energy().WN() - 
                                          dBaseStates[ipSpecies][0].energy().WN();
                        /* zero some rates */   
                        source[ipLo] = 0.;
                        sink[ipLo] = 0.;
                }

                // non-zero was due to roundoff errors on 32-bit
                if( ex[0] <= dBaseStates[ipSpecies][0].energy().WN()* 10. *DBL_EPSILON )
                        ex[0] = 0.;
                else
                        TotalInsanity();

                for( long ipHi= 0; ipHi<dBaseSpecies[ipSpecies].numLevels_local; ipHi++)
                {
                        for( long ipLo= 0; ipLo<dBaseSpecies[ipSpecies].numLevels_local; ipLo++)
                        {
                                if (ipHi > ipLo)
                                {
                                        AulEscp[ipHi][ipLo] = SMALLFLOAT;
                                        AulDest[ipHi][ipLo] = SMALLFLOAT;
                                        AulPump[ipLo][ipHi] = SMALLFLOAT;
                                }
                                else
                                {
                                        AulEscp[ipHi][ipLo] = 0.;
                                        AulDest[ipHi][ipLo] = 0.;
                                        AulPump[ipLo][ipHi] = 0.;
                                }
                        }
                }

                for(TransitionList::iterator tr = dBaseTrans[ipSpecies].begin();
                         tr != dBaseTrans[ipSpecies].end(); ++tr)
                {
                        int ipHi = (*tr).ipHi();
                        if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local || (*tr).ipCont() <= 0)
                                continue;
                        int ipLo = (*tr).ipLo();
                        AulEscp[ipHi][ipLo] = (*tr).Emis().Aul()*
                                ((*tr).Emis().Pesc() + (*tr).Emis().Pelec_esc());
                        AulDest[ipHi][ipLo] = (*tr).Emis().Aul()*(*tr).Emis().Pdest();
                        AulPump[ipLo][ipHi] = (*tr).Emis().pump();
                }

                /*Setting all the collision strengths and collision rate to zero*/
                for( long ipHi= 0; ipHi<dBaseSpecies[ipSpecies].numLevels_local; ipHi++)
                {
                        for( long ipLo= 0; ipLo<dBaseSpecies[ipSpecies].numLevels_local; ipLo++)
                        {
                                col_str[ipHi][ipLo] = 0.;
                                CollRate[ipHi][ipLo] = 0.;
                        }
                }

                /*Setting all the collision strengths and collision rate to zero*/
                for(TransitionList::iterator tr = dBaseTrans[ipSpecies].begin();
                         tr != dBaseTrans[ipSpecies].end(); ++tr)
                {
                        int ipHi = (*tr).ipHi();
                        if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local)
                                continue;
                        CollisionZero( (*tr).Coll() );
                        for( long k=0; k<ipNCOLLIDER; ++k )
                                (*tr).Coll().rate_coef_ul_set()[k] = 0.f;
                }

                /* update the collision rates */
                /* molecule */
                if( dBaseSpecies[ipSpecies].lgLAMDA )
                {
                        for(TransitionList::iterator tr = dBaseTrans[ipSpecies].begin();
                                 tr != dBaseTrans[ipSpecies].end(); ++tr)
                        {
                                int ipHi = (*tr).ipHi();
                                if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local)
                                        continue;
                                for( long intCollNo=0; intCollNo<ipNCOLLIDER; intCollNo++)
                                {
                                        /*using the collision rate coefficients directly*/
                                        (*tr).Coll().rate_coef_ul_set()[intCollNo] =
                                                (realnum)LeidenCollRate(ipSpecies, intCollNo, *tr, phycon.te);
                                }
                        }
                }
                /* Chianti */
                else if( dBaseTrans[ipSpecies].chLabel() == "Chianti" )
                {
                        for(TransitionList::iterator tr = dBaseTrans[ipSpecies].begin();
                                 tr != dBaseTrans[ipSpecies].end(); ++tr)
                        {
                                if ((*tr).ipHi() >= dBaseSpecies[ipSpecies].numLevels_local)
                                        continue;
                                for( long intCollNo=0; intCollNo<ipNCOLLIDER; intCollNo++)
                                {
                                        (*tr).Coll().rate_coef_ul_set()[intCollNo] =
                                                (realnum)ChiantiCollRate(ipSpecies, intCollNo, *tr, phycon.te);
                                }
                        }
                }
                /* Stout */
                else if( dBaseTrans[ipSpecies].chLabel() == "Stout" )
                {
                        for(TransitionList::iterator tr = dBaseTrans[ipSpecies].begin();
                                 tr != dBaseTrans[ipSpecies].end(); ++tr)
                        {
                                if ((*tr).ipHi() >= dBaseSpecies[ipSpecies].numLevels_local)
                                        continue;
                                for( long intCollNo=0; intCollNo<ipNCOLLIDER; intCollNo++)
                                {
                                        (*tr).Coll().rate_coef_ul_set()[intCollNo] =
                                                        (realnum)StoutCollRate(ipSpecies, intCollNo, *tr, phycon.te);
                                }
                        }
                }
                else
                        TotalInsanity();
                
                /* guess some missing data */
                for(TransitionList::iterator tr = dBaseTrans[ipSpecies].begin();
                         tr != dBaseTrans[ipSpecies].end(); ++tr)
                {
                        int ipHi = (*tr).ipHi();
                        if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local)
                                continue;
                        int ipLo = (*tr).ipLo();
                        const CollisionProxy &coll_temp = (*tr).Coll();

                                /* make educated guesses for some missing data */
                        if( dBaseSpecies[ipSpecies].lgMolecular )
                        {
                                ASSERT( dBaseSpecies[ipSpecies].lgLAMDA );
                                /*The collision rate coefficients for helium should not be present and that for molecular hydrogen should be present*/
                                if( AtmolCollRateCoeff[ipSpecies][ipATOM_HE].temps.size() == 0 &&
                                        AtmolCollRateCoeff[ipSpecies][ipH2].temps.size() != 0 )
                                {
                                        coll_temp.rate_coef_ul_set()[ipATOM_HE] = 0.7f * coll_temp.rate_coef_ul()[ipH2];
                                }

                                /* Put in something for hydrogen collisions if not in database */
                                if( AtmolCollRateCoeff[ipSpecies][ipATOM_H].temps.size() == 0 )
                                {
                                        if( AtmolCollRateCoeff[ipSpecies][ipATOM_HE].temps.size() != 0 ) //He0
                                        {
                                                coll_temp.rate_coef_ul_set()[ipATOM_H] = 2.0f * coll_temp.rate_coef_ul()[ipATOM_HE];
                                        }
                                        else if( AtmolCollRateCoeff[ipSpecies][ipH2_ORTHO].temps.size() != 0 ) //ortho-H2
                                        {
                                                coll_temp.rate_coef_ul_set()[ipATOM_H] = 1.4f * coll_temp.rate_coef_ul()[ipH2_ORTHO];
                                        }
                                        else if( AtmolCollRateCoeff[ipSpecies][ipH2_PARA].temps.size() != 0 ) //para-H2
                                        {
                                                coll_temp.rate_coef_ul_set()[ipATOM_H] = 1.4f * coll_temp.rate_coef_ul()[ipH2_PARA];
                                        }
                                        else if( AtmolCollRateCoeff[ipSpecies][ipH2].temps.size() != 0 ) // total H2
                                        {
                                                coll_temp.rate_coef_ul_set()[ipATOM_H] = 1.4f * coll_temp.rate_coef_ul()[ipH2];
                                        }
                                        else
                                                coll_temp.rate_coef_ul_set()[ipATOM_H] = 1e-13f * (realnum)g[ipLo];
                                }

                                /* Put in something for proton collisions if not in database */
                                if( AtmolCollRateCoeff[ipSpecies][ipPROTON].temps.size() == 0 )
                                {
                                        if( AtmolCollRateCoeff[ipSpecies][ipHE_PLUS].temps.size() != 0 ) //He+
                                        {
                                                coll_temp.rate_coef_ul_set()[ipPROTON] = 2.0f * coll_temp.rate_coef_ul()[ipHE_PLUS];
                                        }
                                        else
                                                coll_temp.rate_coef_ul_set()[ipPROTON] = 1e-13f * (realnum)g[ipLo];

                                }
                                
#if     0
                                /* if nothing else has been done, just put a small rate coefficient in */
                                for( long intCollNo=0; intCollNo<ipNCOLLIDER; intCollNo++)
                                {
                                        if( coll_temp.rate_coef_ul()[intCollNo] == 0. )
                                                coll_temp.rate_coef_ul_set()[intCollNo] = 1e-13;
                                }
#endif
                        }
                        else
                        {
                                /* test for transitions without collision data */
                                if( (*tr).Coll().rate_coef_ul_set()[ipELECTRON] == 0. )
                                {
                                        /* Specific transitions of Fe 3,4,and 5 have data that are not in the Chianti/Kurucz files*/
                                        if( strcmp(dBaseSpecies[ipSpecies].chLabel,"Fe 3") == 0 && ipHi < 14 )
                                        {
                                                coll_temp.col_str() = Fe3_cs(ipLo,ipHi);
                                        }
                                        else if( strcmp(dBaseSpecies[ipSpecies].chLabel,"Fe 4") == 0 && ipHi < 12 )
                                        {
                                                coll_temp.col_str() = Fe4_cs(ipLo,ipHi);
                                        }
                                        else if( strcmp(dBaseSpecies[ipSpecies].chLabel,"Fe 5") == 0 && ipHi < 14 )
                                        {
                                                coll_temp.col_str() = Fe5_cs(ipLo,ipHi);
                                        }
                                        else
                                        {
                                                /* All other transitions should use gbar */
                                                MakeCS(*tr);
                                                /* Sometimes gf is zero and then gbar = 0
                                                 * If this happens use a default value */
                                                ASSERT( (*tr).Emis().Aul() != 0. );
                                                if( coll_temp.col_str() == 0. )
                                                {
                                                        coll_temp.col_str() = atmdat.collstrDefault;
                                                }
                                        }

                                        coll_temp.rate_coef_ul_set()[ipELECTRON] = (COLL_CONST*coll_temp.col_str())/
                                                        (g[ipHi]*phycon.sqrte);
                                }
                        }
                }

                /*Updating the CollRate*/
                for(TransitionList::iterator tr = dBaseTrans[ipSpecies].begin();
                         tr != dBaseTrans[ipSpecies].end(); ++tr)
                {
                        int ipHi = (*tr).ipHi();
                        if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local)
                                continue;
                        int ipLo = (*tr).ipLo();
                        CollRate[ipHi][ipLo] = (*tr).Coll().ColUL( colliders );
                        
                        CollRate[ipLo][ipHi] = CollRate[ipHi][ipLo] * g[ipHi] / g[ipLo] *
                                dsexp( (*tr).EnergyK() / phycon.te );
                        
                        /* now add in excitations resulting from cosmic ray secondaries */
                        // \todo 2 add branch to do forbidden transitions       
                        // this g-bar only works for permitted lines
                        if( (*tr).ipCont() > 0 )
                        {
                                /* get secondaries for all permitted lines by scaling LyA 
                                 * excitation by ratio of cross section (oscillator strength/energy) 
                                 * Born approximation or plane-wave approximation */
                                (*tr).Coll().rate_lu_nontherm_set() = secondaries.x12tot *
                                        ((*tr).Emis().gf()/(*tr).EnergyWN()) /
                                        (iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,0).Emis().gf()/
                                        iso_sp[ipH_LIKE][ipHYDROGEN].trans(ipH2p,0).EnergyWN());
                                        
                                CollRate[ipLo][ipHi] += (*tr).Coll().rate_lu_nontherm();
                                CollRate[ipHi][ipLo] += (*tr).Coll().rate_lu_nontherm() * (*(*tr).Lo()).g() / (*(*tr).Hi()).g();
                        }
                }

                multi_arr<double,2> Cool(dBaseSpecies[ipSpecies].numLevels_local, dBaseSpecies[ipSpecies].numLevels_local);

                /* solve the n-level atom */
                atom_levelN(
                        /* dBaseSpecies[ipSpecies].numLevels_local is the number of levels to compute*/ 
                        dBaseSpecies[ipSpecies].numLevels_local, 
                        /* ABUND is total abundance of species, used for nth equation
                         * if balance equations are homogeneous */
                        abund, 
                        /* g(dBaseSpecies[ipSpecies].numLevels_local) is statistical weight of levels */
                        g, 
                        /* EX(dBaseSpecies[ipSpecies].numLevels_local) is excitation potential of levels, deg K or wavenumbers
                         * 0 for lowest level, all are energy rel to ground NOT d(ENER) */
                        ex, 
                        /* this is 'K' for ex[] as Kelvin deg, is 'w' for wavenumbers */
                        'w',
                        /* populations [cm-3] of each level as deduced here */
                        pops, 
                        /* departure coefficient, derived below */
                        depart,
                        /* net transition rate, A * esc prob, s-1 */
                        &AulEscp, 
                        /* col str from up to low */
                        &col_str, 
                        /* AulDest[ihi][ilo] is destruction rate, trans from ihi to ilo, A * dest prob,
                         * asserts confirm that [ihi][ilo] is zero */
                        &AulDest, 
                        /* AulPump[ilo][ihi] is pumping rate from lower to upper level (s^-1), (hi,lo) must be zero  */
                        &AulPump, 
                        /* collision rates (s^-1), evaluated here and returned for cooling by calling function,
                         * unless following flag is true.  If true then calling function has already filled
                         * in these rates.  CollRate[ipSpecies][j] is rate from ipSpecies to j */
                        &CollRate,
                        /* this is an additional creation rate from continuum, normally zero, units cm-3 s-1 */
                        source,
                        /* this is an additional destruction rate to continuum, normally zero, units s-1 */
                        sink,
                        // flag saying whether CollRate already done (true), or we need to do it here (false),
                        true,
                        /* total cooling and its derivative, set here but nothing done with it*/
                        &cooltl, 
                        &coolder, 
                        /* string used to identify calling program in case of error */
                        spName, 
                        /* nNegPop flag indicating what we have done
                         * positive if negative populations occurred
                         * zero if normal calculation done
                         * negative if too cold (for some atoms other routine will be called in this case) */
                        &nNegPop,
                        /* true if populations are zero, either due to zero abundance of very low temperature */
                        &lgZeroPop,
                        /* option to print debug information */
                        lgDeBug,
                        /* option to do the molecule in LTE */
                        dBaseSpecies[ipSpecies].lgLTE,
                        /* cooling per line */
                        &Cool);

                if( nNegPop > 0 )
                {
                        /* negative populations occurred */
                        fprintf(ioQQQ," PROBLEM in dBase_solve, atom_levelN returned negative population .\n");
                        cdEXIT( EXIT_FAILURE );
                }

                // highest levels may have no population
                while( (pops[dBaseSpecies[ipSpecies].numLevels_local-1]<=0 ) &&
                        (dBaseSpecies[ipSpecies].numLevels_local > 1) )
                                --dBaseSpecies[ipSpecies].numLevels_local;

                for( long j=0;j< dBaseSpecies[ipSpecies].numLevels_local; j++ )
                        dBaseStates[ipSpecies][j].Pop() = MAX2(pops[j],SMALLFLOAT);

                for( long j=dBaseSpecies[ipSpecies].numLevels_local;
                        j< dBaseSpecies[ipSpecies].numLevels_max; j++ )
                                dBaseStates[ipSpecies][j].Pop() = 0.;

# ifdef USE_NLTE7
                // do sums of totals out (destruction) and into (creation) level
                for( int lvl = 0; lvl < dBaseSpecies[ipSpecies].numLevels_local; lvl++)
                {
                        for( int lvl2 = 0; lvl2 < dBaseSpecies[ipSpecies].numLevels_local; lvl2++)
                        {
                                if( lvl != lvl2 )
                                {
                                        dBaseStates[ipSpecies][lvl].DestCollBB() += CollRate[lvl][lvl2];
                                        dBaseStates[ipSpecies][lvl].CreatCollBB() += CollRate[lvl2][lvl]*pops[lvl2];
                                        if( lvl2 < lvl)
                                                dBaseStates[ipSpecies][lvl].DestPhotoBB() += AulDest[lvl][lvl2];
                                }
                        }
                }
# endif

                /*Atmol  line*/
                for(TransitionList::iterator tr = dBaseTrans[ipSpecies].begin();
                         tr != dBaseTrans[ipSpecies].end(); ++tr)
                {
                        int ipHi = (*tr).ipHi();
                        if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local)
                                continue;
                        int ipLo = (*tr).ipLo();
                        (*tr).Coll().cool() = max(Cool[ipHi][ipLo],0.);
                        (*tr).Coll().heat() = max(-Cool[ipHi][ipLo],0.);
                        
                        if ( (*tr).ipCont() > 0 )
                        {
                        
                                (*tr).Emis().phots() =  (*tr).Emis().Aul() * (*(*tr).Hi()).Pop() *
                                                ((*tr).Emis().Pesc() + (*tr).Emis().Pelec_esc());

                                (*tr).Emis().xIntensity() = (*tr).Emis().phots() * (*tr).EnergyErg();

                                        /* population of lower level rel to ion, corrected for stim em */
                                (*tr).Emis().PopOpc() = (*(*tr).Lo()).Pop() - (*(*tr).Hi()).Pop()*
                                        (*(*tr).Lo()).g()/(*(*tr).Hi()).g();

                                /* it's possible for this sum to be zero.  Set ratio to zero in that case. */
                                if( CollRate[ipLo][ipHi]+AulPump[ipLo][ipHi] > 0. )
                                {
                                        (*tr).Emis().ColOvTot() = CollRate[ipLo][ipHi]/
                                                (CollRate[ipLo][ipHi]+AulPump[ipLo][ipHi]);
                                }
                                else
                                        (*tr).Emis().ColOvTot() = 0.;
                                
                                // this is only used for save line printout.  Maybe colliders may be involved, but
                                // this simple approximation of a "collision strength" should be good enough for
                                // the purposes of that printout.
                                (*tr).Coll().col_str() = (realnum)( (*tr).Coll().rate_coef_ul()[ipELECTRON] *
                                        (g[ipHi]*phycon.sqrte)/COLL_CONST);
                        }
                        else
                        {
                                (*tr).Coll().col_str() = 0.;
                        }
                }

                for(TransitionList::iterator tr = dBaseTrans[ipSpecies].begin();
                         tr != dBaseTrans[ipSpecies].end(); ++tr)
                {
                        int ipHi = (*tr).ipHi();
                        if (ipHi < dBaseSpecies[ipSpecies].numLevels_local)
                                continue;
                        EmLineZero( (*tr).Emis() );             
                }

                dBaseSpecies[ipSpecies].CoolTotal = cooltl;
                CoolAdd( dBaseSpecies[ipSpecies].chLabel, 0., max(0.,dBaseSpecies[ipSpecies].CoolTotal) );
                if( dBaseSpecies[ipSpecies].lgMolecular )
                        thermal.elementcool[LIMELM] += max(0.,dBaseSpecies[ipSpecies].CoolTotal);
                else
                        thermal.elementcool[dBaseStates[ipSpecies][0].nelem()-1] += max(0.,dBaseSpecies[ipSpecies].CoolTotal);
                totalHeating += max(0., -dBaseSpecies[ipSpecies].CoolTotal);
                thermal.dCooldT += coolder;

                /* option to print departure coefficients */
                {
                        enum {DEBUG_LOC=false};

                        if( DEBUG_LOC )
                        {
                                fprintf( ioQQQ, " Departure coefficients for species %li\n", ipSpecies );
                                for( long j=0; j< dBaseSpecies[ipSpecies].numLevels_local; j++ )
                                {
                                        fprintf( ioQQQ, " level %li \t Depar Coef %e\n", j, depart[j] );
                                }
                        }
                }
        }

        // total heating for all dBase dBaseSpecies     
        thermal.heating[0][27] = totalHeating;

        return;
}

/*Leiden*/
STATIC double LeidenCollRate(long ipSpecies, long ipCollider, const TransitionProxy& tr, double ftemp)
{
        DEBUG_ENTRY("LeidenCollRate()");
        double ret_collrate = InterpCollRate( AtmolCollRateCoeff[ipSpecies][ipCollider], tr.ipHi(), tr.ipLo(), ftemp);
        return ret_collrate;
}

/*STOUT*/
STATIC double StoutCollRate(long ipSpecies, long ipCollider, const TransitionProxy& tr, double ftemp)
{
        DEBUG_ENTRY("StoutCollRate()");

        double rate = 0.;
        // deexcitation rate or collision strength?
        bool lgIsRate = StoutCollData[ipSpecies][tr.ipHi()][tr.ipLo()][ipCollider].lgIsRate;
        int n = StoutCollData[ipSpecies][tr.ipHi()][tr.ipLo()][ipCollider].ntemps;
        if( n < 2)
                return 0.;

        double *x = (double*)MALLOC(n*sizeof(double));
        double *y = (double*)MALLOC(n*sizeof(double));

        double fupsilon = 0.;
        for(int j = 0; j < n; j ++)
        {
                x[j] = StoutCollData[ipSpecies][tr.ipHi()][tr.ipLo()][ipCollider].temps[j];
                y[j] = StoutCollData[ipSpecies][tr.ipHi()][tr.ipLo()][ipCollider].collstrs[j];
                ASSERT( x[j] > 0. && y[j] > 0.);
        }
        //If the temperature is above or below the temperature range, use the CS from the closest temperature.
        //Otherwise, do the linear interpolation.
        if( ftemp < x[0] )
        {
                fupsilon = y[0];
        }
        else if( ftemp > x[n-1] )
        {
                fupsilon = y[n-1];
        }
        else
        {
                fupsilon = linint(x,y,n,ftemp);
        }

        free(x);
        free(y);
        ASSERT(fupsilon > 0);

        /* We can deal with derexcitation rates and collision strengths currently */
        if( lgIsRate )
        {
                rate = fupsilon;
        }
        else
        {
                /* convert the collision strength to a collision rate coefficient */
                /* This formula converting collision strength to collision rate coefficient works fine for the electrons*/
                /* For any other collider the mass would be different*/
                rate = (COLL_CONST*fupsilon)/((*tr.Hi()).g()*sqrt(ftemp));
        }

        return rate;
}

/*CHIANTI*/
STATIC double ChiantiCollRate(long ipSpecies, long ipCollider, const TransitionProxy& tr, double ftemp)
{
        DEBUG_ENTRY("ChiantiCollRate()");

        double rate = 0.;
        double fupsilon = CHIANTI_Upsilon(ipSpecies, ipCollider, tr.ipHi(), tr.ipLo(), ftemp);

        /* NB NB - if proton colliders, the upsilons returned here are actually already rate coefficients. */
        /* these are designated by a collider index and a transition type */
        if( ipCollider == ipPROTON && AtmolCollSplines[ipSpecies][tr.ipHi()][tr.ipLo()][ipCollider].intTranType == 6 )
        {
                rate = fupsilon;
        }
        else if( ipCollider == ipELECTRON )
        {
                /* convert the collision strength to a collision rate coefficient */
                /*This formula converting collision strength to collision rate coefficient works fine for the electrons*/
                /*For any other collider the mass would be different*/
                rate = (COLL_CONST*fupsilon)/((*tr.Hi()).g()*sqrt(ftemp));
        }
        else
                rate = 0.;

        return rate;
}

double CHIANTI_Upsilon(long ipSpecies, long ipCollider, long ipHi, long ipLo, double ftemp)
{
        double fdeltae,fscalingparam,fkte,fxt,fsups,fups;
        int intxs,inttype,intsplinepts;

        DEBUG_ENTRY("CHIANTI_Upsilon()");

        if( AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].collspline == NULL )
        {
                return 0.;
        }

        intsplinepts = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].nSplinePts;
        inttype = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].intTranType;
        fdeltae = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].EnergyDiff;
        fscalingparam = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].ScalingParam;

        fkte = ftemp/fdeltae/1.57888e5;

        /*Way the temperature is scaled*/
        /*Burgess&Tully 1992:Paper gives only types 1 to 4*/
        /*Found that the expressions were the same for 5 & 6 from the associated routine DESCALE_ALL*/
        /*What about 7,8&9?*/
        if( inttype ==1 || inttype==4 )
        {
                fxt = 1-(log(fscalingparam)/(log(fkte+fscalingparam)));
        }
        else if(inttype  == 2 || inttype == 3||inttype == 5 || inttype == 6)
        {
                fxt = fkte/(fkte+fscalingparam);
        }
        else
                TotalInsanity();

        double xs[9],*spl,*spl2;
        /*Creating spline points array*/
        if(intsplinepts == 5)
        {
                spl = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].collspline;
                for(intxs=0;intxs<5;intxs++)
                {
                        xs[intxs] = 0.25*intxs;
                        if(DEBUGSTATE)
                        {
                                printf("The xs and spl values are %f and %f \n",xs[intxs],spl[intxs]);
                                getchar();
                        }
                }
        }
        else if(intsplinepts == 9)
        {
                spl = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].collspline;   
                for( intxs=0; intxs<9; intxs++ )
                {
                        xs[intxs] = 0.125*intxs;
                        if(DEBUGSTATE)
                        {
                                printf("The xs and spl values are %f and %f \n",xs[intxs],spl[intxs]);
                                getchar();
                        }
                }
        }
        else
        {
                TotalInsanity();
        }

        /*Finding the second derivative*/
        spl2 = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].SplineSecDer;

        if(DEBUGSTATE)
        {
                printf("\n");
                for(intxs=0;intxs<intsplinepts;intxs++)
                {
                        printf("The %d value of 2nd derivative is %f \n",intxs+1,spl2[intxs]);
                }
        }

        /*Extracting out the value*/
        //splint(xs,spl,spl2,intsplinepts,fxt,&fsups);

        fsups = linint( xs, spl, intsplinepts, fxt);

        /*Finding upsilon*/
        if(inttype == 1)
        {
                fups = fsups*log(fkte+exp(1.0));
        }
        else if(inttype == 2)
        {
                fups = fsups;
        }
        else if(inttype == 3)
        {
                fups = fsups/(fkte+1.0) ;
        }
        else if(inttype == 4)
        {
                fups = fsups*log(fkte+fscalingparam) ;
        }
        else if(inttype == 5)
        {
                fups = fsups/fkte ;
        }
        else if(inttype == 6)
        {
                fups = pow(10.0,fsups) ;
        }
        else
        {
                TotalInsanity();
        }

        if( fups < 0. ) 
        {
                fprintf( ioQQQ," WARNING: Negative upsilon in species %s, collider %li, indices %4li %4li, Te = %e.\n",
                                dBaseSpecies[ipSpecies].chLabel, ipCollider, ipHi, ipLo, ftemp );
                fups = 0.;
        }
        ASSERT(fups>=0);
        return(fups);
}

---