/home66/gary/public_html/cloudy/c08_branch/source/nemala2.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
* others.  For conditions of distribution and use see copyright notice in license.txt */
#include "cddefines.h"
#include "phycon.h"
#include "taulines.h"
#include "mole.h"
#include "atoms.h"
#include "string.h"
#include "thirdparty.h"
#include "dense.h"
#include "conv.h"
#include "h2.h"
#include "physconst.h"

void states_popfill( void);
STATIC double LeidenCollRate(long, long, long, long,double);
STATIC double Chianti_Upsilon(long, long, long, long,double);

#define DEBUGSTATE false


/*Solving for the level populations*/

void atmol_popsolve(void )
{
        realnum abund;
        DEBUG_ENTRY( "atmol_popsolve()" );

        for( long i=0; i<nSpecies; i++ )
        {
                const char *spName = Species[i].chptrSpName;
                double *g, *ex, *pops, *depart;
                double **AulEscp, **col_str, **AulDest, **AulPump, **CollRate,fcolldensity[9];
                double *source, *sink;
                double cooltl, coolder;
                double fupsilon;
                long ipLo,ipHi,intCollNo;
                int nNegPop;
                bool lgZeroPop, lgDeBug = false;
                long nlev = Species[i].numLevels_max;

                /* first find current density (cm-3) of species */
                if( lgSpeciesMolecule[i] )
                {
                        struct molecule *SpeciesCurrent;
                        if( (SpeciesCurrent = findspecies(Species[i].chptrSpName)) == &null_mole )
                        {
                                /* did not find the species - print warning for now */
                                fprintf(ioQQQ," PROBLEM atmol_popsolve did not find molecular species %li\n",i);
                        }
                        abund = SpeciesCurrent->hevmol;
                }
                else
                {
                        /* an atom or ion */
                        ASSERT( Species[i].intAtNo<LIMELM && Species[i].intIonStage<=Species[i].intAtNo);
                        abund = dense.xIonDense[Species[i].intAtNo][Species[i].intIonStage];
                }

                if(abund<=SMALLFLOAT)
                {
                        /* zero out populations and intensities */
                        /* NB - can't use states_popfill here, because it zeros ALL species */
                        atmolStates[i][0].Pop = 0.;
                        for(long ipHi = 1; ipHi < nlev; ipHi++ )
                        {       
                                atmolStates[i][ipHi].Pop = 0.;
                                for(long ipLo = 0; ipLo < ipHi; ipLo++ )
                                {       
                                        if( atmolTrans[i][ipHi][ipLo].ipCont > 0 )
                                        {
                                                atmolTrans[i][ipHi][ipLo].Emis->phots = 0.;
                                                atmolTrans[i][ipHi][ipLo].Emis->xIntensity = 0.;
                                        }
                                }
                        }
                        continue;
                }

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

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

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

                for( ipLo = 0; ipLo < nlev; ipLo++ )
                {
                        /*Filling in statistical weights & Excitation Energies*/
                        g[ipLo] = atmolStates[i][ipLo].g ;
                        ex[ipLo] = atmolStates[i][ipLo].energy;
                        /* zero some rates */   
                        source[ipLo] = 0.;
                        sink[ipLo] = 0.;
                        AulEscp[ipLo][ipLo] = 0.;
                        AulDest[ipLo][ipLo] = 0.;
                        AulPump[ipLo][ipLo] = 0.;
                }

                for( ipHi=1; ipHi<nlev; ipHi++ )
                {
                        for( ipLo=0; ipLo<ipHi; ipLo++ )
                        {
                                if( atmolTrans[i][ipHi][ipLo].ipCont > 0 )
                                {       
                                        AulEscp[ipHi][ipLo] = atmolTrans[i][ipHi][ipLo].Emis->Aul*
                                                (atmolTrans[i][ipHi][ipLo].Emis->Pesc +
                                                atmolTrans[i][ipHi][ipLo].Emis->Pelec_esc);
                                        AulDest[ipHi][ipLo] = atmolTrans[i][ipHi][ipLo].Emis->Aul*
                                                atmolTrans[i][ipHi][ipLo].Emis->Pdest;
                                        AulPump[ipLo][ipHi] = atmolTrans[i][ipHi][ipLo].Emis->pump;
                                }
                                else
                                {
                                        AulEscp[ipHi][ipLo] = SMALLFLOAT;
                                        AulDest[ipHi][ipLo] = SMALLFLOAT;
                                        AulPump[ipLo][ipHi] = SMALLFLOAT;
                                }

                                AulEscp[ipLo][ipHi] = 0.;
                                AulDest[ipLo][ipHi] = 0.;
                                AulPump[ipHi][ipLo] = 0.;
                        }
                }

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

                /*Updating the CollRatesArray*/
                /*You need to do this for a species indexed by i*/
                for( intCollNo=0; intCollNo<NUM_COLLIDERS; intCollNo++)
                {
                        for( ipHi=1; ipHi<nlev; ipHi++)
                        {
                                for( ipLo=0; ipLo<ipHi; ipLo++)
                                {
                                        /* molecule */
                                        if( lgSpeciesMolecule[i] )
                                        {
                                                if(CollRatesArray[i][intCollNo]!=NULL)
                                                {
                                                        /*using the collision rate coefficients directly*/
                                                        CollRatesArray[i][intCollNo][ipHi][ipLo] = 
                                                                LeidenCollRate(i, intCollNo, ipHi, ipLo, phycon.te);
                                                }
                                        }
                                        /* atom or ion */
                                        else
                                        {
                                                if(CollRatesArray[i][intCollNo]!=NULL)
                                                {
                                                        fupsilon = Chianti_Upsilon(i, intCollNo, ipHi, ipLo, phycon.te);
                                                        /*converting the collision strength to a collision rate coefficient*/
                                                        /*This formula converting collision strength to collision rate coefficent works fine for the electrons*/
                                                        /*For any other collider the mass would be different*/
                                                        CollRatesArray[i][intCollNo][ipHi][ipLo] = ((8.6291e-6)*fupsilon)/(g[ipHi]*phycon.sqrte);
                                                }
                                        }

                                        /* now get the corresponding excitation rate */
                                        if(CollRatesArray[i][intCollNo]!=NULL)
                                        {
                                                CollRatesArray[i][intCollNo][ipLo][ipHi] = 
                                                        CollRatesArray[i][intCollNo][ipHi][ipLo] * g[ipHi] / g[ipLo] *
                                                        sexp( atmolTrans[i][ipHi][ipLo].EnergyK / phycon.te );
                                        }
                                }
                        }
                }
                /*Get the densities seperately*/
                /*Electron*/
                fcolldensity[0] = dense.eden;
                /*Proton*/
                fcolldensity[1] = dense.xIonDense[1][1];
                /*Atomic Hydrogen*/
                fcolldensity[2] = dense.xIonDense[1][0];
                /*Helium*/
                fcolldensity[3] = dense.xIonDense[2][0];
                /*Helium+*/
                fcolldensity[4] = dense.xIonDense[2][1];
                /*Helium ++*/
                fcolldensity[5] = dense.xIonDense[2][2];
                /*Molecular Hydrogen*/
                fcolldensity[8] = findspecies("H2")->hevmol;
                /*Ortho(6) and Para(7) Mol Hydrogen*/
                if(h2.lgH2ON)
                {
                        fcolldensity[6] = h2.ortho_density;
                        fcolldensity[7] = h2.para_density;
                }
                else
                {               
                        fcolldensity[6] = (fcolldensity[8])*(0.75);
                        fcolldensity[7] = (fcolldensity[8])*(0.25);
                }
                
                /*Updating the CollRate*/
                for( ipHi=0; ipHi<nlev; ipHi++ )
                {       
                        for( ipLo=0; ipLo<nlev; ipLo++ )
                        {
                                if( ipHi == ipLo )
                                {
                                        CollRate[ipHi][ipLo] = 0.;
                                        continue;
                                }

                                for( intCollNo=0; intCollNo<NUM_COLLIDERS; intCollNo++ )
                                {
                                        if(CollRatesArray[i][intCollNo]!=NULL)
                                        {
                                                CollRate[ipHi][ipLo] += CollRatesArray[i][intCollNo][ipHi][ipLo]*fcolldensity[intCollNo];
                                        }
                                }
                                /*Correction for Helium*/
                                /*To include the effects of collision with Helium,the user must multiply the density by 1.14*/
                                /*pg 5,Schoier et al 2004*/
                                if(lgSpeciesMolecule[i])
                                {
                                        /*The collision rate coefficients for helium should not be present and that for molecular hydrogen should be present*/
                                        if(CollRatesArray[i][3]==NULL && CollRatesArray[i][8]!=NULL )
                                        {
                                                CollRate[ipHi][ipLo] += 0.7 *CollRatesArray[i][8][ipHi][ipLo]*fcolldensity[3];
                                        }
                                }
                        }
                }
                /*Level Lowering as seen in mole_co_atom.cpp*/
                /*Excitation energies need to be converted to temperature*/
                /*ex[nUsed]/phycon.te cannot be much larger than 25, or matrix inversion will fail*/
                while ( ((ex[nlev-1]*T1CM) > phycon.te*25.) && (nlev > 1) )
                {
                        --nlev;
                }
                Species[i].numLevels_local = nlev;

                /* build some arrays */
                atom_levelN(
                        /* nlev is the number of levels to compute*/ 
                        nlev, 
                        /* ABUND is total abundance of species, used for nth equation
                         * if balance equations are homogeneous */
                        abund, 
                        /* G(nlev) is stat weight of levels */
                        g, 
                        /* EX(nlev) 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[i][j] is rate from i 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),
                         * this is stored in data)[ihi][ilo] as either downward rate or collis strength*/
                        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 );

                if( nNegPop == 0 )
                {
                        ASSERT( lgZeroPop == false );

                        for( long j=0; j< nlev; j++ )
                        {
                                atmolStates[i][j].Pop = pops[j];
                        }
                }
                else if( nNegPop > 0 )
                {
                        /* something should be done here. */
                        /* try another solver? */
                        /* probably don't want to quit */
                        fprintf(ioQQQ," PROBLEM, atom_levelN returned negative population .\n");
                        cdEXIT( EXIT_FAILURE ); 
                }
                else
                {
                        /*In search phase we dont lower the levels*/
                        if(conv.lgSearch!= true)
                        {
                                for(long j=Species[i].numLevels_max-1; j>0; j--)
                                {
                                        if((atmolStates[i][j].Pop/abund)<POPTHRES)
                                        {
                                                nlev--;
                                        }
                                        else
                                                break;
                                }
                                
                                ASSERT( nlev >= 1 );
                                /*Setting the number of energy levels to new limit*/
                                Species[i].numLevels_local = nlev;

                                if( nlev == 1 )
                                {
                                        atmolStates[i][0].Pop = abund;
                                        for( long j=1; j< Species[i].numLevels_max; j++ )
                                        {
                                                atmolStates[i][j].Pop = 0.;
                                        }
                                }
                                else
                                {
                                        /*Solving for level populations*/
                                        atom_levelN(nlev,abund,g,ex,'w',pops,depart,&AulEscp,&col_str,
                                                                        &AulDest, &AulPump, &CollRate, source, sink,true,&cooltl, 
                                                                        &coolder, spName, &nNegPop,     &lgZeroPop, lgDeBug );
                                        
                                        for( long j=0; j< nlev; j++ )
                                        {
                                                atmolStates[i][j].Pop = pops[j];
                                        }
                                        for( long j=nlev; j< Species[i].numLevels_max; j++ )
                                        {
                                                atmolStates[i][j].Pop = 0.;
                                        }
                                }
                        }
                }

                /*Atmol  line*/
                for(long ipHi = 1; ipHi < nlev; ipHi++ )
                {       
                        for(long ipLo = 0; ipLo < ipHi; ipLo++ )
                        {       
                                if( atmolTrans[i][ipHi][ipLo].ipCont > 0 )
                                {
                                        atmolTrans[i][ipHi][ipLo].Emis->phots =
                                                atmolTrans[i][ipHi][ipLo].Emis->Aul * 
                                                atmolTrans[i][ipHi][ipLo].Hi->Pop *
                                                (atmolTrans[i][ipHi][ipLo].Emis->Pesc + 
                                                atmolTrans[i][ipHi][ipLo].Emis->Pelec_esc);

                                        atmolTrans[i][ipHi][ipLo].Emis->xIntensity = 
                                                atmolTrans[i][ipHi][ipLo].Emis->phots * 
                                                atmolTrans[i][ipHi][ipLo].EnergyErg;

                                        /* population of lower level rel to ion, corrected for stim em */

                                        atmolTrans[i][ipHi][ipLo].Emis->PopOpc =
                                                atmolTrans[i][ipHi][ipLo].Lo->Pop - atmolTrans[i][ipHi][ipLo].Hi->Pop*
                                                atmolTrans[i][ipHi][ipLo].Lo->g/atmolTrans[i][ipHi][ipLo].Hi->g;
                                        fixit(); /* need to do something with these! */
                                        atmolTrans[i][ipHi][ipLo].Emis->ColOvTot = 0.;
                                        atmolTrans[i][ipHi][ipLo].Coll.cool = 0.;
                                        atmolTrans[i][ipHi][ipLo].Coll.heat = 0.;
                                }
                        }
                }

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

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

                /* free vectors */
                free( g );
                free( ex );
                free( pops );
                free( depart );
                free( source );
                free( sink );

                /* free arrays */
                for( long j=0; j< nlev; j++ )
                {
                        free( AulEscp[j] );
                        free( col_str[j] ); 
                        free( AulDest[j] ); 
                        free( AulPump[j] );  
                        free( CollRate[j] ); 
                }
                free( AulEscp );
                free( col_str ); 
                free( AulDest ); 
                free( AulPump );  
                free( CollRate );
        }
        return;
}

/*This function fills the population of the states to a dummy value of 0*/
void states_popfill( void)
{
        DEBUG_ENTRY( "states_popfill()" );

        for( long i=0; i<nSpecies; i++)
        {
                for( long j=0; j<Species[i].numLevels_max; j++)
                {
                        atmolStates[i][j].Pop = 0.;
                }
        }
        return;
}

/*Leiden*/
double LeidenCollRate(long ipSpecies, long ipCollider, long ipHi, long ipLo, double ftemp)
{
        double ret_collrate;
        int inttemps;
        DEBUG_ENTRY("LeidenCollRate()");
        inttemps = AtmolCollRateCoeff[ipSpecies][ipCollider].ntemps;

        if( AtmolCollRateCoeff[ipSpecies][ipCollider].temps == NULL )
        {
                return 0.;
        }

        if(ftemp < AtmolCollRateCoeff[ipSpecies][ipCollider].temps[0])
        {
                ret_collrate = AtmolCollRateCoeff[ipSpecies][ipCollider].collrates[ipHi][ipLo][0];
        }
        else if(ftemp > AtmolCollRateCoeff[ipSpecies][ipCollider].temps[inttemps-1])
        {
                ret_collrate = AtmolCollRateCoeff[ipSpecies][ipCollider].collrates[ipHi][ipLo][inttemps-1];
        }
        else
        {       
                ret_collrate = linint(AtmolCollRateCoeff[ipSpecies][ipCollider].temps,
                        AtmolCollRateCoeff[ipSpecies][ipCollider].collrates[ipHi][ipLo],
                        AtmolCollRateCoeff[ipSpecies][ipCollider].ntemps,
                        ftemp);
        }

        if(DEBUGSTATE)
                printf("\nThe collision rate at temperature %f is %e",ftemp,ret_collrate);

        ASSERT( !isnan( ret_collrate ) );
        return(ret_collrate);

}

/*Chianti*/
double Chianti_Upsilon(long ipSpecies, long ipCollider, long ipHi, long ipLo, double ftemp)
{
        double fdeltae,fscalingparam,fkte,fxt,fsups,fups;
        double *xs,*spl,*spl2;
        int i,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();

        /*Creating spline points array*/
        if(intsplinepts == 5)
        {
                xs = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double));
                spl = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double));
                spl2 = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double));
                for(intxs=0;intxs<5;intxs++)
                {
                        xs[intxs] = 0.25*intxs;
                        spl[intxs] = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].collspline[intxs];
                        if(DEBUGSTATE)
                        {
                                printf("The xs and spl values are %f and %f \n",xs[intxs],spl[intxs]);
                                getchar();
                        }
                }
        }
        else if(intsplinepts == 9)
        {
                xs = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double));
                spl = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double));
                spl2 = (double *)MALLOC((unsigned long)intsplinepts*sizeof(double));
                for( intxs=0; intxs<9; intxs++ )
                {
                        xs[intxs] = 0.125*intxs;
                        spl[intxs] = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].collspline[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*/
        for( i=0; i<intsplinepts; i++)
        {
                spl2[i] = AtmolCollSplines[ipSpecies][ipHi][ipLo][ipCollider].SplineSecDer[i];
        }

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

        /*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();
        }

        free( xs );
        free( spl );
        free( spl2 );

        ASSERT(fups>=0);
        return(fups);
}

---