/home66/gary/public_html/cloudy/c08_branch/source/parse_abundances.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 */
/*ParseAbundances parse and read in composition as set by abundances command */
#include "cddefines.h"
#include "grains.h"
#include "abund.h"
#include "phycon.h"
#include "called.h"
#include "elementnames.h"
#include "input.h"
#include "parse.h"

void ParseAbundances(char *chCard,
                        /* following set true by grains command,
                        * so this will not set more grains if grains already on. */
                     bool lgDSet)
{
        bool lgEOF, 
          lgEOL, 
          lgLog;
        long int i, 
          j, 
          k;
        double absav[LIMELM], 
          chk;
        GrainPar gp;

        DEBUG_ENTRY( "ParseAbundances()" );

        j = 5;
        absav[0] = FFmtRead(chCard,&j,INPUT_LINE_LENGTH,&lgEOL);
        /* check whether it scanned off end of line on first try
         * at least one number on the line if it passes this test
         * no numbers means a keyword */
        if( !lgEOL )
        {
                absav[1] = FFmtRead(chCard,&j,INPUT_LINE_LENGTH,&lgEOL);
                if( lgEOL )
                {
                        /* this branch, we found one, but not two, numbers -
                         * must be one of a few special cases */
                        if( nMatch(" ALL",chCard) )
                        {
                                /* special option, all abundances will be this number */
                                if( absav[0] <= 0. )
                                {
                                        absav[0] = pow(10.,absav[0]);
                                }
                                for( i=1; i < LIMELM; i++ )
                                {
                                        abund.solar[i] = (realnum)absav[0];
                                }

                        }
                        else if( nMatch("OLD ",chCard) && nMatch("SOLA",chCard) )
                        {
                                i = (int)absav[0];
                                /* old solar - better be the number "84" on the line */
                                if( i!=84 )
                                {
                                        fprintf( ioQQQ, 
                                                " The only old abundance set I have is for version 84 - %3ld was requested.  Sorry.\n", 
                                                i );
                                        cdEXIT(EXIT_FAILURE);
                                }
                                for( i=1; i < LIMELM; i++ )
                                {
                                        /* these are the old abundances */
                                        abund.SolarSave[i] = abund.OldSolar84[i];
                                        abund.solar[i] = abund.OldSolar84[i];
                                }
                        }

                        else if( nMatch("STAR",chCard) )
                        {
                                /* Fred Hamonn's star burst galaxy mixture */
                                abund_starburst(chCard);
                        }
                        else
                        {
                                fprintf( ioQQQ, 
                                        " I did not recognize a sub-keyword - options are ALL and OLD SOLAR 84.  Sorry.\n");
                                cdEXIT(EXIT_FAILURE);
                        }

                        /* normal return */
                        return;
                }

                /* we get here if there is a second number - read in all abundances */
                for( i=2; i < abund.npSolar; i++ )
                {
                        absav[i] = FFmtRead(chCard,&j,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                /*   read CONTINUE line if J scanned off end before reading all abund */
                                input_readarray(chCard,&lgEOF);
                                if( lgEOF )
                                {
                                        fprintf( ioQQQ, " There MUST be%3ld abundances entered, there were only%3ld.  Sorry.\n", 
                                                abund.npSolar, i );
                                        cdEXIT(EXIT_FAILURE);
                                }

                                /* option to ignore all lines starting with #, *, or % */
                                while( lgInputComment(chCard) )
                                {
                                        if( lgEOF )
                                        {
                                                fprintf( ioQQQ, " There MUST be%3ld abundances entered, there were only%3ld.  Sorry.\n", 
                                                        abund.npSolar, i );
                                                cdEXIT(EXIT_FAILURE);
                                        }
                                        /* read in another line */
                                        input_readarray(chCard,&lgEOF);
                                }

                                /* we have the line image, print it */
                                if( called.lgTalk )
                                {
                                        fprintf( ioQQQ, "                       * ");
                                        k=0;
                                        while( chCard[k]!='\0' )
                                        {
                                                fprintf(ioQQQ,"%c",chCard[k]);
                                                ++k;
                                        }
                                        while( k<80 )
                                        {
                                                fprintf(ioQQQ,"%c",' ');
                                                ++k;
                                        }
                                        fprintf( ioQQQ,"*\n"); 
                                }

                                /* now convert to caps */
                                caps(chCard);
                                if( strncmp(chCard,"CONT",4) != 0 )
                                {
                                        fprintf( ioQQQ, " There MUST be%3ld abundances entered, there were only%3ld.  Sorry.\n", 
                                                abund.npSolar, i );
                                        cdEXIT(EXIT_FAILURE);
                                }
                                else
                                {
                                        j = 1;
                                        absav[i] = FFmtRead(chCard,&j,INPUT_LINE_LENGTH,&lgEOL);
                                        if( lgEOF )
                                        {
                                                fprintf( ioQQQ, " There MUST be%3ld abundances entered, there were only%3ld.  Sorry.\n", 
                                                        abund.npSolar, i);
                                                cdEXIT(EXIT_FAILURE);
                                        }
                                }
                        }
                }

                /* fell through to here after reading in N abundances for N elements
                 * check that there are no more abundances on the line - that would
                 * be an error - a typo */
                chk = FFmtRead(chCard,&j,INPUT_LINE_LENGTH,&lgEOL);
                if( !lgEOL || (chk!=0.) )
                {
                        /* got another number, not lgEOL, so too many numbers */
                        fprintf( ioQQQ, " There were more than %ld abundances entered\n", 
                                abund.npSolar );
                        fprintf( ioQQQ, " Could there have been a typo somewhere?\n" );
                }

                /* are numbers scale factors, or log of abund rel to H?? */
                lgLog = false;
                for( i=0; i < abund.npSolar; i++ )
                {
                        if( absav[i] < 0. )
                                lgLog = true;
                }

                if( lgLog )
                {
                        /* entered as log of number rel to hydrogen */
                        for( i=0; i < abund.npSolar; i++ )
                        {
                                abund.solar[abund.ipSolar[i]-1] = (realnum)pow(10.,absav[i]);
                        }
                }
                else
                {
                        /* scale factors relative to solar */
                        for( i=0; i < abund.npSolar; i++ )
                        {
                                abund.solar[abund.ipSolar[i]-1] *= (realnum)absav[i];
                        }
                }

                /* check whether the abundances are reasonable */
                if( abund.solar[1] > 0.2 )
                {
                        fprintf( ioQQQ, " Is an abundance of%10.3e for %2.2s reasonable?\n", 
                                abund.solar[1], elementnames.chElementSym[1] );
                }

                for( i=2; i < LIMELM; i++ )
                {
                        if( abund.solar[i] > 0.1 )
                        {
                                fprintf( ioQQQ, " Is an abundance of%10.3e for %2.2s reasonable?\n", 
                                        abund.solar[i], elementnames.chElementSym[i] );
                        }
                }
                return;
        }

        /* this branch if no numbers on the line - option to use one of several
         * special stored abundance sets */
        if( nMatch(" AGB",chCard) || nMatch("AGB ",chCard) || nMatch("PLAN",chCard) )
        {
                /* AGB/planetary nebula abundances */
                /* only turn on grains if "no grains" is not present */
                if( !nMatch("O GR",chCard) )
                {
                        /* turn on grains if not already done */
                        if( !lgDSet )
                        {
                                /* return bins allocated by previous abundances ... commands */
                                ReturnGrainBins();
                                /* now allocate new grain bins */
                                gp.dep = 1.;
                                gp.lgAbunVsDepth = false;
                                gp.lgRequestQHeating = true;
                                gp.lgForbidQHeating = false;
                                gp.lgGreyGrain = false;

                                /* NO QHEAT option to turn off quantum heating for grains */
                                if( nMatch("O QH",chCard) ) 
                                {
                                        gp.lgForbidQHeating = true;
                                        gp.lgRequestQHeating = false;
                                        phycon.lgPhysOK = false;
                                }

                                /* actually set up the grains */
                                mie_read_opc("graphite_ism_10.opc",gp);
                                mie_read_opc("silicate_ism_10.opc",gp);
                        }
                }

                for( i=0; i < LIMELM; i++ )
                {
                        abund.solar[i] = abund.apn[i];
                }
        }

        else if( nMatch("HII ",chCard) || nMatch("H II",chCard) || nMatch("ORIO",chCard) )
        {
                /* H II region abundances */
                /* only turn on grains if "no grains" is not present */
                if( !nMatch("O GR",chCard) )
                {
                        /* option to turn on grains */
                        if( !lgDSet )
                        {
                                /* return bins allocated by previous abundances ... commands */
                                ReturnGrainBins();
                                /* now allocate new grain bins */
                                gp.dep = 1.;
                                gp.lgAbunVsDepth = false;
                                gp.lgRequestQHeating = true;
                                gp.lgForbidQHeating = false;
                                gp.lgGreyGrain = false;

                                /* NO QHEAT option to turn off quantum heating for grains */
                                if( nMatch("O QH",chCard) ) 
                                {
                                        gp.lgForbidQHeating = true;
                                        gp.lgRequestQHeating = false;
                                        phycon.lgPhysOK = false;
                                }
                                /* This scales the Orion grain abundance so that the observed 
                                 * dust to gas ratio that Cloudy predicts is in agreement with
                                 * that observed in the Veil,
                                 *>>refer       grain   Abel, N., Brogan, C., Ferland, G., O'Dell, C.R., 
                                 *>>refercon    Shaw, G., Troland, T., 2004, ApJ, submitted */
                                gp.dep *= 0.85;

                                mie_read_opc("graphite_orion_10.opc",gp);
                                mie_read_opc("silicate_orion_10.opc",gp);
                        }
                }

                for( i=0; i < LIMELM; i++ )
                {
                        abund.solar[i] = abund.ahii[i];
                }
        }

        else if( nMatch("ISM ",chCard) || nMatch(" ISM",chCard) )
        {
                /* ISM abundances from Cowie and Songaila Ann Rev '86 */
                /* only turn on grains if "no grains" is not present */
                if( !nMatch("O GR",chCard) )
                {
                        if( !lgDSet )
                        {
                                /* return bins allocated by previous abundances ... commands */
                                ReturnGrainBins();
                                /* now allocate new grain bins */
                                gp.dep = 1.;
                                gp.lgAbunVsDepth = false;
                                gp.lgRequestQHeating = true;
                                gp.lgForbidQHeating = false;
                                gp.lgGreyGrain = false;

                                /* NO QHEAT option to turn off quantum heating */
                                if( nMatch("O QH",chCard) ) 
                                {
                                        gp.lgForbidQHeating = true;
                                        gp.lgRequestQHeating = false;
                                        phycon.lgPhysOK = false;
                                }

                                /* actually set up the grains */
                                mie_read_opc("graphite_ism_10.opc",gp);
                                mie_read_opc("silicate_ism_10.opc",gp);
                        }
                }

                for( i=0; i < LIMELM; i++ )
                {
                        abund.solar[i] = abund.aism[i];
                }
        }

        else if( nMatch("NOVA",chCard) )
        {
                /* Nova Cyg abundances */
                for( i=0; i < LIMELM; i++ )
                {
                        abund.solar[i] = abund.anova[i];
                }
        }

        else if( nMatch("PRIM",chCard) )
        {
                /* roughly primordial abundances: He/H=.072 */
                for( i=0; i < LIMELM; i++ )
                {
                        abund.solar[i] = abund.aprim[i];
                }

                /* now turn off the heavy elements */
                for( i=4; i < LIMELM; i++ )
                {
                        /* turn off heavy elements - do this way to make sure,
                         * that H-like and He-like level limits are handled properly */
                        char chDUMMY[LIMELM];
                        sprintf(chDUMMY,"element %s off ", elementnames.chElementName[i] );
                        /* now convert to caps */
                        caps(chDUMMY);
                        ParseElement( chDUMMY );
                }
        }

        else if( nMatch("CAME",chCard) )
        {
                /* mix from Cameron 1982, "Essays on Nuclear Astrophysics" */
                for( i=0; i < LIMELM; i++ )
                {
                        abund.solar[i] = abund.camern[i];
                }
        }

        else
        {
                fprintf( ioQQQ, 
                        " ABUND must have PLAN, H II, CAMERON, NOVA, ALL, STARBURST, OLD SOLAR 84 or PRIMORDIAL.  Sorry.\n" );
                cdEXIT(EXIT_FAILURE);
        }
        return;
}

---