/home66/gary/public_html/cloudy/c08_branch/source/parse_element.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 */
/*ParseElement parse options on element command */
#include "cddefines.h"
#include "optimize.h"
#include "abund.h"
#include "dense.h"
#include "input.h"
#include "iso.h"
#include "hmi.h"
#include "mole.h"
#include "elementnames.h"
#include "parse.h"

void ParseElement(char *chCard )
{
        /* this will be used to remember how many levels are active in any element we turn off,
         * so that we can retain state if turned back on  */
        static bool lgFirst = true;
        static long levels[NISO][LIMELM];
        char chCap[INPUT_LINE_LENGTH];
        bool lgEOL, 
          lgEnd, 
          lgHIT;

        bool lgForceLog=false, lgForceLinear=false;

        long int i, 
          nelem, 
          j, 
          k;
        double param;

        DEBUG_ENTRY( "ParseElement()" );

        /* zero out this array if first call */
        if( lgFirst )
        {
                lgFirst = false;
                for( i=0; i<NISO; ++i )
                {
                        for( nelem=0; nelem<LIMELM; ++nelem )
                        {
                                levels[i][nelem] = 0;
                        }
                }
        }
        /* say that abundances have been changed */
        abund.lgAbnSolar = false;

        /* read in list of elements for abundances command */
        if( nMatch("READ",chCard) )
        {
                abund.npSolar = 0;
                for( i=1; i < LIMELM; i++ )
                {
                        input_readarray(chCard,&lgEnd);

                        if( lgEnd )
                        {
                                fprintf( ioQQQ, " Hit EOF while reading element list; use END to end list.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }

                        strcpy( chCap, chCard );
                        caps(chCap);

                        /* this would be a line starting with END to say end on list */
                        if( strncmp(chCap,"END",3) == 0 )
                        {
                                return;
                        }

                        j = 1;
                        lgHIT = false;
                        while( j < LIMELM && !lgHIT )
                        {
                                j += 1;

                                if( strncmp( chCap,elementnames.chElementNameShort[j-1] , 4) == 0 )
                                {
                                        abund.npSolar += 1;
                                        abund.ipSolar[abund.npSolar-1] = j;
                                        lgHIT = true;
                                }
                        }

                        if( !lgHIT )
                        {
                                fprintf( ioQQQ, "%80.80s\n", chCard );
                                fprintf( ioQQQ, " Sorry, but I did not recognize element name on this line.\n" );
                                fprintf( ioQQQ, " Here is the list of names I recognize.\n" );
                                fprintf( ioQQQ, " " );

                                for( k=2; k <= LIMELM; k++ )
                                {
                                        fprintf( ioQQQ, "%4.4s\n", elementnames.chElementNameShort[k-1] );
                                }

                                cdEXIT(EXIT_FAILURE);
                        }
                }

                /* fell through, make sure one more line, the end line */
                input_readarray(chCard,&lgEnd);
                strcpy( chCap, chCard );
                caps(chCap);

                if( strncmp(chCap,"END",3) == 0 )
                {
                        return;
                }

                else
                {
                        fprintf( ioQQQ, " Too many elements were entered.\n" );
                        fprintf( ioQQQ, " I only know about%3d elements.\n", 
                          LIMELM );
                        fprintf( ioQQQ, " Sorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        /* find which element - will be used in remainder of routine to 
         * adjust aspects of this element */
        nelem = GetElem(chCard );

        bool lgElementSet = false;
        if( nelem >= ipHYDROGEN )
        {
                lgElementSet = true;
                /* nelem is now the atomic number of the element, and must be correct */
                ASSERT( nelem>=0 && nelem < LIMELM );
        }

        /* look for log or linear keywords */
        if( nMatch(" LOG",chCard) )
                lgForceLog = true;
        else if( nMatch("LINE",chCard) )
                lgForceLinear = true;

        i = 4;
        param = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

        if( nMatch("SCAL",chCard) )
        {
                /* enter abundance as scale factor, relative to what is in now */
                if( lgEOL )
                {
                        fprintf( ioQQQ, " There must be a number on this line.\n" );
                        fprintf( ioQQQ, " Sorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                else
                {

                        if( !lgElementSet )
                        {
                                fprintf( ioQQQ, 
                                        " ParseElement did not find an element on the following line:\n" );
                                fprintf( ioQQQ, 
                                        "%80.80s\n", chCard );
                                cdEXIT(EXIT_FAILURE);
                        }
                        /* interpret as log unless forced linear */
                        if( (lgForceLog || param <= 0.) && !lgForceLinear )
                        {
                                /* option for log of scale factor */
                                param = pow(10.,param);
                        }
                        abund.ScaleElement[nelem] = (realnum)param;
                }
        }

        else if( nMatch("ABUN",chCard) )
        {
                /* log of actual abundance */
                if( lgEOL )
                {
                        fprintf( ioQQQ, " There must be a number on this line.\n" );
                        fprintf( ioQQQ, " Sorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                else
                {

                        if( !lgElementSet )
                        {
                                fprintf( ioQQQ, 
                                        " ParseElement did not find an element on the following line:\n" );
                                fprintf( ioQQQ, 
                                        "%80.80s\n", chCard );
                                cdEXIT(EXIT_FAILURE);
                        }
                        if( lgForceLinear )
                        {
                                abund.solar[nelem] = (realnum)param;
                        }
                        else
                        {
                                abund.solar[nelem] = (realnum)pow(10.,param);
                        }

                        if( abund.solar[nelem]> 1. )
                        {
                                fprintf( ioQQQ, 
                                        " Please check the abundance of this element.  It seems high to me.\n" );
                        }
                }
        }

        else if( nMatch(" OFF",chCard) )
        {
                /* keyword LIMIT sets lower limit to abundance of element
                 * that will be included */
                /* option to turn off this element, set abundance to zero */
                if( nMatch( "LIMI",chCard) )
                {
                        if( lgEOL )
                        {
                                fprintf( ioQQQ, " There must be an abundances on the ELEMENT OFF LIMIT command.\n" );
                                fprintf( ioQQQ, " Sorry.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                        dense.AbundanceLimit = (realnum)pow(10., param);
                }
                else
                {
                        if( !lgElementSet )
                        {
                                fprintf( ioQQQ, 
                                        " ParseElement did not find an element on the following line:\n" );
                                fprintf( ioQQQ, 
                                        "%80.80s\n", chCard );
                                cdEXIT(EXIT_FAILURE);
                        }
                        /* specify element name */
                        dense.lgElmtOn[nelem] = false;
                        /* another flag that element is off */
                        dense.IonHigh[nelem] = -1;
                        dense.lgSetIoniz[nelem] = false;
                        /* set no levels for all elements that are turned off, except for helium itself, which always exists */
                        if( nelem > ipHELIUM )
                        {
                                levels[ipHYDROGEN][nelem] = iso.numLevels_max[ipH_LIKE][nelem];
                                levels[ipHELIUM][nelem] = iso.numLevels_max[ipHE_LIKE][nelem];
                                iso.numLevels_max[ipH_LIKE][nelem] = 0;
                                iso.numLevels_max[ipHE_LIKE][nelem] = 0;
                        }

                        if( nelem == ipHYDROGEN )
                        {
                                fprintf( ioQQQ, " It is not possible to turn hydrogen off.\n" );
                                fprintf( ioQQQ, " Sorry.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                }
        }

        /* specify an ionization distribution */
        else if( nMatch("IONI",chCard) )
        {
                bool lgLogSet = false;
                int ion;
                int ihi , low;


                if( !lgElementSet )
                {
                        fprintf( ioQQQ, 
                                " ParseElement did not find an element on the following line:\n" );
                        fprintf( ioQQQ, 
                                "%80.80s\n", chCard );
                        cdEXIT(EXIT_FAILURE);
                }
                if( !dense.lgElmtOn[nelem] )
                {
                        fprintf(ioQQQ,"Sorry, you cannot set the ionization of %s since it has been turned off.\n" ,
                                elementnames.chElementName[nelem] );
                        cdEXIT(EXIT_FAILURE);
                }

                i = 4;
                /* option to specify ionization distribution for this element */
                dense.lgSetIoniz[nelem] = true;
                ion = 0;
                while( ion<nelem+2 )
                {
                        /* the ionization fractions that are set when above set true,
                         * gas phase abundance is this times total abundance
                         * Ionization fraction for [nelem][ion] */
                        dense.SetIoniz[nelem][ion] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

                        if( lgEOL ) 
                                break;

                        /* all are log if any are negative */
                        if( dense.SetIoniz[nelem][ion] < 0. )
                                lgLogSet = true;
                        ++ion;
                }

                /* zero out ones not specified */
                for( i=ion; i<nelem+2; ++i )
                        dense.SetIoniz[nelem][i] = 0.;

                /* convert rest to linear if any were negative */
                if( lgLogSet )
                {
                        for( i=0; i<ion; ++i )
                                dense.SetIoniz[nelem][i] = (realnum)pow((realnum)10.f, dense.SetIoniz[nelem][i]);
                }

                /* now check that no zero abundances exist between lowest and highest non-zero values */
                low = 0;
                while( dense.SetIoniz[nelem][low]==0 && low < nelem+1 )
                        ++low;
                ihi = nelem+1;
                while( dense.SetIoniz[nelem][ihi]==0 && ihi > low )
                        --ihi;

                if( ihi==low && dense.SetIoniz[nelem][ihi]==0 )
                {
                        fprintf(ioQQQ," element ionization command has all zero ionization fractions.  This is not possible.\n Sorry\n");
                        cdEXIT(EXIT_FAILURE);
                }
                for(i=low; i<=ihi; ++i )
                {
                        if( dense.SetIoniz[nelem][i]==0 )
                        {
                                fprintf(ioQQQ," element abundance command has zero abundance between positive values.  This is not possible.\n Sorry\n");
                                cdEXIT(EXIT_FAILURE);
                        }
                }
                /* >>chng 05 mar 13, add this
                 * finally turn off any part of the chemical networks that might use this element,
                 * since chemistry must come into equilibrium with the ionization, and vice vese -
                 * ion distribution cannot come into equilibrium with chemistry when ions are set,
                 * so much turn off chemistry */
                if( nelem==ipHYDROGEN && !hmi.lgNoH2Mole )
                {
                        fprintf(ioQQQ," Hydrogen ionization has been set, H chemistry is disabled.\n");
                        /* turn off only H2 */
                        hmi.lgNoH2Mole = true;
                }
                /* loop over all atoms in co chem net - these occur as the last NUM_ELEMENTS 
                 * elements in the full set of mole.num_heavy_molec */
                /* >>chng 06 mar 19, do not print if network already turned off */
                if( mole.lgElem_in_chemistry[nelem] && !co.lgNoCOMole )
                {
                        fprintf(ioQQQ," The ionization of an element in the CO chemistry network has been set, CO chemistry is disabled.\n");
                        /* turn off CO network */
                        co.lgNoCOMole = true;
                }
        }

        /* option to turn element on - some ini files turn off elements,
         * this can turn them back on */
        else if( nMatch(" ON ",chCard) )
        {

                if( !lgElementSet )
                {
                        fprintf( ioQQQ, 
                                " ParseElement did not find an element on the following line:\n" );
                        fprintf( ioQQQ, 
                                "%80.80s\n", chCard );
                        cdEXIT(EXIT_FAILURE);
                }
                /* option to turn off this element, set abundance to zero */
                dense.lgElmtOn[nelem] = true;
                /* reset levels to default if they were ever turned off with element off command */
                if( levels[ipHYDROGEN][nelem] )
                {
                        iso.numLevels_max[ipH_LIKE][nelem] = levels[ipHYDROGEN][nelem];
                        iso.numLevels_max[ipHE_LIKE][nelem] = levels[ipHELIUM][nelem];
                }
        }

        else if( nMatch("TABL",chCard) )
        {

                if( !lgElementSet )
                {
                        fprintf( ioQQQ, 
                                " ParseElement did not find an element on the following line:\n" );
                        fprintf( ioQQQ, 
                                "%80.80s\n", chCard );
                        cdEXIT(EXIT_FAILURE);
                }
                /* read in table of position-dependent abundances for a particular element. */
                abund.lgAbunTabl[nelem] = true;

                /* general flag saying this option turned on */
                abund.lgAbTaON = true;

                /* does the table give depth or radius?  keyword DEPTH
                 * says it is depth, default is radius */
                if( nMatch("DEPT",chCard) )
                        abund.lgAbTaDepth[nelem] = true;
                else
                        abund.lgAbTaDepth[nelem] = false;

                /* make sure not trying to change hydrogen */
                if( nelem == ipHYDROGEN )
                {
                        fprintf( ioQQQ, " cannot change abundance of hydrogen.\n" );
                        fprintf( ioQQQ, " Sorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                /* read pair giving radius and abundance */
                input_readarray(chCard,&lgEnd);
                i = 1;
                abund.AbTabRad[0][nelem] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                abund.AbTabFac[0][nelem] = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

                if( lgEOL )
                {
                        fprintf( ioQQQ, " no pairs entered - cannot interpolate\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                /* number of points in the table */
                abund.nAbunTabl = 2;

                lgEnd = false;
                /* LIMTAB is limit to number of points we can store and is
                 * set to 500 in abundances */
                while( !lgEnd && abund.nAbunTabl < LIMTABD )
                {
                        input_readarray(chCard,&lgEnd);
                        if( !lgEnd )
                        {
                                /* convert first 4 char to caps, into chCap */
                                cap4(chCap , chCard);
                                if( strncmp(chCap,"END",3) == 0 )
                                        lgEnd = true;
                        }

                        /* lgEnd may have been set within above if, if end line encountered*/
                        if( !lgEnd )
                        {
                                i = 1;
                                abund.AbTabRad[abund.nAbunTabl-1][nelem] = 
                                        (realnum)FFmtRead(chCard ,&i,INPUT_LINE_LENGTH,&lgEOL);

                                abund.AbTabFac[abund.nAbunTabl-1][nelem] = 
                                        (realnum)FFmtRead(chCard ,&i,INPUT_LINE_LENGTH,&lgEOL);
                                abund.nAbunTabl += 1;
                        }
                }

                abund.nAbunTabl -= 1;

                /* now check that radii are in increasing order */
                for( i=1; i < abund.nAbunTabl; i++ )
                {
                        /* the radius values are assumed to be strictly increasing */
                        if( abund.AbTabRad[i][nelem] <= abund.AbTabRad[i-1][nelem] )
                        {
                                fprintf( ioQQQ, "ParseElement: TABLE ELEMENT TABLE radii "
                                        "must be in increasing order\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                }
        }

        else
        {
                fprintf( ioQQQ, "ParseElement: ELEMENT command - there must be "
                        "a keyword on this line.\n" );
                fprintf( ioQQQ, " The keys I know about are TABLE, SCALE, _OFF, "
                        "_ON_, IONIZATION, and ABUNDANCE.\n" );
                fprintf( ioQQQ, " Sorry.\n" );
                cdEXIT(EXIT_FAILURE);
        }

        /* vary option */
        if( optimize.lgVarOn )
        {
                optimize.nvarxt[optimize.nparm] = 1;
                /* pointer to where to write */
                optimize.nvfpnt[optimize.nparm] = input.nRead;

                if( nMatch("SCAL",chCard) )
                {
                        /* vary scale factor */
                        sprintf( optimize.chVarFmt[optimize.nparm], "ELEMENT %4.4s SCALE %%f LOG", 
                          elementnames.chElementNameShort[nelem] );

                        /* param is linear scale factor */
                        optimize.vparm[0][optimize.nparm] = (realnum)log10(param);
                        optimize.vincr[optimize.nparm] = 0.2f;
                }

                else if( nMatch("ABUN",chCard) )
                {
                        /* vary absolute abundance */
                        sprintf( optimize.chVarFmt[optimize.nparm], "ELEMENT %4.4s ABUND %%f LOG ", 
                          elementnames.chElementNameShort[nelem] );

                        /* param is log of abundance by number relative to hydrogen */
                        optimize.vparm[0][optimize.nparm] = (realnum)param;
                        optimize.vincr[optimize.nparm] = 0.2f;
                }
                ++optimize.nparm;
        }
        return;
}

---