/home66/gary/public_html/cloudy/c08_branch/source/parse_stop.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 */
/*ParseStop parse the stop command */
#include "cddefines.h"
#include "physconst.h"
#include "optimize.h"
#include "phycon.h"
#include "rfield.h"
#include "radius.h"
#include "geometry.h"
#include "iterations.h"
#include "stopcalc.h"
#include "input.h"
#include "parse.h"

void ParseStop(char *chCard)
{
        bool lgEOL;

        long int i, 
          j;

        double a, 
          effcol, 
          tread;

        DEBUG_ENTRY( "ParseStop()" );

        /* first scan off a generic number, all subcommands use at least one */
        i = 5;
        a = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

        /* >>chng 04 nov 09, introduce off option */
        if( lgEOL && !nMatch(" OFF",chCard) )
        {
                NoNumb(chCard);
        }

        if( nMatch("TEMP",chCard) )
        {
                /* >>chng 04 nov 09, introduce off option to disable this stopping criterion */
                if( lgEOL && nMatch(" OFF",chCard) )
                {
                        /* this is special case for ending temperature - do not use -
                         * but will still stop if Te falls below TeLowest, the lowest
                         * possible temperature */
                        StopCalc.tend = -1.f;
                }
                else
                {
                        /* lowest electron temperature allowed before stopping
                         * assumed to be the log of the temperature if <10
                         * optional keyword LINEAR forces linear */
                        if( a <= 10. && !nMatch("LINE",chCard) )
                        {
                                tread = pow(10.,a);
                        }
                        else
                        {
                                tread = a;
                        }

                        /* tread is linear temperature*/
                        if( tread < phycon.TEMP_LIMIT_LOW )
                        {
                                fprintf( ioQQQ, 
                                        " Temperatures below %.2e K not allowed. Reset to lowest value."
                                        "   I am doing this myself.\n" ,
                                        phycon.TEMP_LIMIT_LOW );
                                /* set slightly off extreme limit for safety */
                                tread = phycon.TEMP_LIMIT_LOW*1.01;
                        }
                        else if( tread > phycon.TEMP_LIMIT_HIGH )
                        {
                                fprintf( ioQQQ, 
                                        " Temperatures is above %.2e K not allowed. Reset to highest value."
                                        "   I am doing this myself.\n" ,
                                        phycon.TEMP_LIMIT_HIGH);
                                /* set slightly off extreme limit for safety */
                                tread = phycon.TEMP_LIMIT_HIGH*0.99;
                        }

                        if( nMatch("EXCE",chCard) )
                        {
                                /* option for this to be highest allowed temperature,
                                 * stop temperate exceeds */
                                StopCalc.T2High = (realnum)tread;
                        }
                        else
                        {
                                /* this is ending temperature - we stop if kinetic temperature
                                 * falls below this */
                                StopCalc.tend = (realnum)tread;
                        }
                }
        }

        /* stop at 21cm line center optical depth  */
        else if( nMatch("OPTI",chCard) && nMatch("21CM",chCard) )
        {
                /* default is for number to be log of optical depth */
                bool lgLOG = true;
                if( nMatch("LINE",chCard) )
                {
                        /* force it to be linear not log */
                        lgLOG = false;
                }
                i = 4;
                j = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                if( j!=21 )
                {
                        fprintf( ioQQQ, " First number on STOP 21CM OPTICAL DEPTH command must be 21\n" );
                        cdEXIT(EXIT_FAILURE);
                }
                /* now get the next number, which is the optical depth */
                a = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

                /* tau must be a log, and second number is energy where tau specified */
                if( lgLOG )
                {
                        StopCalc.tauend = (realnum)pow(10.,a);
                }
                else
                {
                        StopCalc.tauend = (realnum)a;
                }
                /* this flag says that 21cm line optical depth is the stop quantity */
                StopCalc.lgStop21cm = true;
        }
        /* stop optical depth at some energy */
        else if( nMatch("OPTI",chCard) )
        {
                /* default is for number to be log of optical depth */
                bool lgLOG = true;
                if( nMatch("LINE",chCard) )
                {
                        /* force it to be linear not log */
                        lgLOG = false;
                }

                if( a > 37. )
                {
                        fprintf( ioQQQ, " optical depth too big\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                /* tau must be a log, and second number is energy where tau specified */
                if( lgLOG )
                {
                        StopCalc.tauend = (realnum)pow(10.,a);
                }
                else
                {
                        StopCalc.tauend = (realnum)a;
                }

                /* energy where tau specified */
                StopCalc.taunu = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

                if( lgEOL )
                {
                        if( nMatch("LYMA",chCard) )
                        {
                                /* largest Lyman limit optical depth */
                                StopCalc.taunu = 1.;
                        }
                        else if( nMatch("BALM",chCard) )
                        {
                                /* stop at this Balmer continuum optical depth */
                                StopCalc.taunu = 0.25;
                        }
                        else
                        {
                                fprintf( ioQQQ, " There must be a second number, the energy in Ryd.  Sorry.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                }

                else
                {
                        /* if second number is negative then log of energy in rydbergs */
                        if( StopCalc.taunu < 0. )
                        {
                                StopCalc.taunu = (realnum)pow((realnum)10.f,StopCalc.taunu);
                        }

                        /* check that energy is within bounds of code */
                        if( StopCalc.taunu < rfield.emm || StopCalc.taunu > rfield.egamry )
                        {
                                fprintf( ioQQQ, " The energy must be in the range %10.2e to %10.2e.  It was %10.2e. Sorry.\n", 
                                  rfield.emm, rfield.egamry, StopCalc.taunu );
                                cdEXIT(EXIT_FAILURE);
                        }
                }
        }

        /* stop optical depth at extinction in V filter */
        else if( nMatch(" AV ",chCard) )
        {
                /* default is for number to be A_V, log if negative */
                if( a<=0. )
                {
                        a = pow(10.,a);
                }
                /* A_V can be for either point or extended source, difference is (1-g) multiplied by scat opacity
                 * if keyword point occurs then for point source, otherwise for extended source */
                if( nMatch("EXTE" , chCard ) )
                {
                        StopCalc.AV_extended = (realnum)a;
                }
                else
                {
                        /* default is point, as measured in ism work */
                        StopCalc.AV_point = (realnum)a;
                }
        }

        /* stop when a fraction of molecules frozen out on grain surfaces is reached */
        else if( nMatch("MOLE",chCard) && nMatch("DEPL",chCard) )
        {
                if( a <= 0. )
                {
                        StopCalc.StopDepleteFrac = (realnum)pow(10.,a);
                }
                else
                {
                        StopCalc.StopDepleteFrac = (realnum)a;
                }
        }

        /* stop when absolute value of flow velocity falls below this value */
        else if( nMatch("VELO",chCard) )
        {
                /* entered in km/s but stored as cm/s */
                StopCalc.StopVelocity = (realnum)(a*1e5);
        }

        /* stop at a given computed mass */
        else if( nMatch("MASS",chCard) )
        {
                /* number of log of mass in gm if inner radius is specified,
                 * mass per unit area, gm cm-2 if not
                 * leave it as a log since dare not deal with linear mass */
                StopCalc.xMass = (realnum)a;
                /* NB 0 is sentinel for not set, if a is zero we must reset it */
                if( StopCalc.xMass == 0 )
                        StopCalc.xMass = SMALLFLOAT;
        }

        /* stop thickness command, also stop depth, this must come after stop
         * optical depth, since do not want to trigger on depth in optical depth */
        else if( nMatch("THIC",chCard) || nMatch("DEPT",chCard) )
        {
                /* decide whether linear or log,
                 * this branch if key linear appears */
                if( nMatch("LINE",chCard) )
                {
                        radius.router[0] = a;
                }
                else
                {
                        /* this is default branch, no key, so log */
                        if( a > 37. )
                        {
                                fprintf( ioQQQ, " thickness too large\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                        radius.router[0] = pow(10.,a);
                }

                /* optional parsec unit */
                if( nMatch("PARS",chCard) )
                {
                        radius.router[0] *= PARSEC;
                }

                /* can stop at different thickness on each iteration */
                for( j=1; j < iterations.iter_malloc; j++ )
                {
                        a = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( lgEOL )
                        {
                                radius.router[j] = radius.router[j-1];
                        }
                        else
                        {
                                if( nMatch("LINE",chCard) )
                                {
                                        radius.router[j] = a;
                                }
                                else
                                {
                                        if( a > 37. )
                                        {
                                                fprintf( ioQQQ, " thickness too large\n" );
                                                cdEXIT(EXIT_FAILURE);
                                        }
                                        radius.router[j] = pow(10.,a);
                                }
                                if( nMatch("PARS",chCard) )
                                {
                                        radius.router[j] *= PARSEC;
                                }
                        }
                }

                /* vary option */
                if( optimize.lgVarOn )
                {
                        optimize.nvarxt[optimize.nparm] = 1;
                        strcpy( optimize.chVarFmt[optimize.nparm], "STOP THICKNESS %f" );

                        /* pointer to where to write */
                        optimize.nvfpnt[optimize.nparm] = input.nRead;

                        /* log of temp will be pointer */
                        optimize.vparm[0][optimize.nparm] = (realnum)log10(radius.router[0]);
                        optimize.vincr[optimize.nparm] = 0.5;

                        ++optimize.nparm;
                }
        }

        /* stop at a particular zone, for each iteration */
        else if( nMatch("ZONE",chCard) )
        {
                /* stop after computing this zone */
                /* >>chng 03 jun 06, do not let fall below 1, stop zone 0 has same effect
                 * as stop zone 1, bug caught by Joop Schaye */
                geometry.nend[0] = (long)MAX2(1.,a);
                geometry.lgZoneSet = true;

                /* this tells code that we intend to stop at this zone, so caution not generated*/
                geometry.lgEndDflt = false;

                for( j=1; j < iterations.iter_malloc; j++ )
                {
                        geometry.nend[j] = (long)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        /* if eol on this iteration, set to previous.  In most cases
                         * all will be equal to the first */
                        if( lgEOL )
                        {
                                geometry.nend[j] = geometry.nend[j-1];
                        }
                        else
                        {
                                /* >>chng 03 jun 06, do not let fall below 1, stop zone 0 has same effect
                                 * as stop zone 1, bug caught by Joop Schaye */
                                geometry.nend[j] = MAX2( 1 , geometry.nend[j] );
                        }
                }
        }

        /* stop at this electron fraction, relative to hydrogen */
        else if( nMatch("EFRA",chCard) )
        {
                if( a <= 0. )
                {
                        StopCalc.StopElecFrac = (realnum)pow(10.,a);
                }
                else
                {
                        StopCalc.StopElecFrac = (realnum)a;
                }
        }

        /* stop at a hydrogen molecular fraction, relative to total hydrogen,
         * this is 2H_2 / H_total*/
        else if( nMatch("MFRA",chCard) )
        {
                if( a <= 0. )
                {
                        StopCalc.StopH2MoleFrac = (realnum)pow(10.,a);
                }
                else
                {
                        StopCalc.StopH2MoleFrac = (realnum)a;
                }
        }

        /* stop at a ionized hydrogen fraction, relative to total hydrogen,
         * this is H+ / H_total */
        else if( nMatch("PFRA",chCard) )
        {
                if( a <= 0. )
                {
                        StopCalc.StopHPlusFrac = (realnum)pow(10.,a);
                }
                else
                {
                        StopCalc.StopHPlusFrac = (realnum)a;
                }
        }

        /* stop at a particular column density */
        else if( nMatch("COLU",chCard) )
        {
                /* check for linear option, if present take log since a being
                 * log column density is default */
                if( nMatch( "LINE" , chCard ) )
                        a = log10(a);

                /*  stop at an effective column density */
                if( nMatch("EFFE",chCard) )
                {
                        /* actually stop at certain optical depth at 1keV */
                        effcol = pow(10.,a);
                        StopCalc.tauend = (realnum)(effcol*2.14e-22);
                        StopCalc.taunu = 73.5;
                        /* vary option */
                        if( optimize.lgVarOn )
                        {
                                optimize.nvarxt[optimize.nparm] = 1;
                                strcpy( optimize.chVarFmt[optimize.nparm], "STOP EFFECTIVE COLUMN DENSITY %f" );
                                /* pointer to where to write */
                                optimize.nvfpnt[optimize.nparm] = input.nRead;
                                /* log of temp will be pointer */
                                optimize.vparm[0][optimize.nparm] = (realnum)log10(effcol);
                                optimize.vincr[optimize.nparm] = 0.5;
                                ++optimize.nparm;
                        }
                }

                else if( nMatch("IONI",chCard) )
                {
                        /*  this is ionized column */
                        if( a > 37. )
                        {
                                fprintf( ioQQQ, " column too big\n" );
                                cdEXIT(EXIT_FAILURE);
                        }

                        StopCalc.colpls = (realnum)pow(10.,a);

                        /*  vary option */
                        if( optimize.lgVarOn )
                        {
                                optimize.nvarxt[optimize.nparm] = 1;
                                strcpy( optimize.chVarFmt[optimize.nparm], "STOP IONIZED COLUMN DENSITY %f" );
                                /*  pointer to where to write */
                                optimize.nvfpnt[optimize.nparm] = input.nRead;
                                /*  log of temp will be pointer */
                                optimize.vparm[0][optimize.nparm] = (realnum)log10(StopCalc.colpls);
                                optimize.vincr[optimize.nparm] = 0.5;
                                ++optimize.nparm;
                        }
                }

                /*  stop at a neutral column */
                else if( nMatch("NEUT",chCard) )
                {
                        StopCalc.colnut = (realnum)pow(10.,a);

                        /*  vary option */
                        if( optimize.lgVarOn )
                        {
                                optimize.nvarxt[optimize.nparm] = 1;
                                strcpy( optimize.chVarFmt[optimize.nparm], "STOP NEUTRAL COLUMN DENSITY %f");
                                /*  pointer to where to write */
                                optimize.nvfpnt[optimize.nparm] = input.nRead;
                                /*  log of temp will be pointer */
                                optimize.vparm[0][optimize.nparm] = (realnum)log10(StopCalc.colnut);
                                optimize.vincr[optimize.nparm] = 0.5;
                                ++optimize.nparm;
                        }
                }

                /* >>chng 03 apr 15, add this option 
                 *  stop at a molecular hydrogen column density, input parameter
                 * is log of the H2 column density */
                else if( nMatch(" H2 ",chCard) )
                {
                        /* this command has a 2 in the H2 label - must not parse the two by
                         * accident.  Get the first number off the line image, and confirm that
                         * it is a 2 */
                        i = 5;
                        j = (long int)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        if( j != 2 )
                        {
                                fprintf( ioQQQ, " Something is wrong with the order of the numbers on this line.\n" );
                                fprintf( ioQQQ, " The first number I encounter should be the 2 in H2.\n Sorry.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }
                        a = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                        StopCalc.col_h2 = (realnum)pow(10.,a);

                        /*  vary option */
                        if( optimize.lgVarOn )
                        {
                                optimize.nvarxt[optimize.nparm] = 1;
                                strcpy( optimize.chVarFmt[optimize.nparm], "STOP H2 COLUMN DENSITY %f");
                                /*  pointer to where to write */
                                optimize.nvfpnt[optimize.nparm] = input.nRead;
                                /*  log of temp will be pointer */
                                optimize.vparm[0][optimize.nparm] = (realnum)log10(StopCalc.col_h2);
                                optimize.vincr[optimize.nparm] = 0.5;
                                ++optimize.nparm;
                        }
                }

                else if( nMatch("ATOM",chCard) )
                {
                        StopCalc.col_h2_nut = (realnum)pow(10.,a);
                        /*  vary option */
                        if( optimize.lgVarOn )
                        {
                                optimize.nvarxt[optimize.nparm] = 1;
                                strcpy( optimize.chVarFmt[optimize.nparm], "STOP ATOMIC COLUMN DENSITY %f");
                                /*  pointer to where to write */
                                optimize.nvfpnt[optimize.nparm] = input.nRead;
                                /*  log of temp will be pointer */
                                optimize.vparm[0][optimize.nparm] = (realnum)log10(StopCalc.col_h2_nut);
                                optimize.vincr[optimize.nparm] = 0.5;
                                ++optimize.nparm;
                        }
                }

                else if( nMatch("H/TS",chCard) )
                {
                        /* >> 05 jan 09, add stop integrated n(H0) / Tspin */
                        StopCalc.col_H0_ov_Tspin = (realnum)pow(10.,a);
                        /*  vary option */
                        if( optimize.lgVarOn )
                        {
                                optimize.nvarxt[optimize.nparm] = 1;
                                strcpy( optimize.chVarFmt[optimize.nparm], "STOP H/TSPIN COLUMN DENSITY %f");
                                /*  pointer to where to write */
                                optimize.nvfpnt[optimize.nparm] = input.nRead;
                                /*  log of temp will be pointer */
                                optimize.vparm[0][optimize.nparm] = (realnum)log10(StopCalc.col_H0_ov_Tspin);
                                optimize.vincr[optimize.nparm] = 0.5;
                                ++optimize.nparm;
                        }
                }

                else if( nMatch(" CO ",chCard) )
                {
                        /* chng << 03 Oct. 27--Nick Abel, add this option */
                        /* stop at a carbon monoxide column density */
                        StopCalc.col_monoxco = (realnum)pow(10.,a);
                        /*  vary option */
                        if( optimize.lgVarOn )
                        {
                                optimize.nvarxt[optimize.nparm] = 1;
                                strcpy( optimize.chVarFmt[optimize.nparm], "STOP CO COLUMN DENSITY %f");
                                /*  pointer to where to write */
                                optimize.nvfpnt[optimize.nparm] = input.nRead;
                                /*  log of temp will be pointer */
                                optimize.vparm[0][optimize.nparm] = (realnum)log10(StopCalc.col_monoxco);
                                optimize.vincr[optimize.nparm] = 0.5;
                                ++optimize.nparm;
                        }
                }

                /* fall through default is total hydrogen column density */
                else
                {
                        /*  both HII and HI */
                        if( a > 37. )
                        {
                                fprintf( ioQQQ, " column too big\n" );
                                cdEXIT(EXIT_FAILURE);
                        }

                        StopCalc.HColStop = (realnum)pow(10.,a);

                        /*  vary option */
                        if( optimize.lgVarOn )
                        {
                                optimize.nvarxt[optimize.nparm] = 1;
                                strcpy( optimize.chVarFmt[optimize.nparm], "STOP COLUMN DENSITY %f" );
                                /*  pointer to where to write */
                                optimize.nvfpnt[optimize.nparm] = input.nRead;
                                /*  log of temp will be pointer */
                                optimize.vparm[0][optimize.nparm] = (realnum)log10(StopCalc.HColStop);
                                optimize.vincr[optimize.nparm] = 0.5;
                                ++optimize.nparm;
                        }
                }
        }

        /* stop when electron density falls below this value, linear or log */
        else if( nMatch("EDEN",chCard) )
        {
                /* stop if electron density falls below this value
                 * LINEAR option */
                if( nMatch("LINE",chCard) )
                {
                        StopCalc.StopElecDensity = (realnum)a;
                }
                else
                {
                        StopCalc.StopElecDensity = (realnum)pow(10.,a);
                }
        }

        /* stop at a particular line ratio - this must come last since many commands
         * have linear option - don't want to trigger on that */
        else if( nMatch("LINE",chCard) )
        {
                char chLabel[5];
                /* first line wavelength, then intensity relative to Hbeta for stop
                 * if third number is entered, it is wl of line in denominator */

                /* get label for the line - must do this first so we clear the label string before
                 * trying to read the wavelength */
                GetQuote( chLabel , chCard , true );

                /* copy first four char of label into caps, and null terminate*/
                strncpy( StopCalc.chStopLabel1[StopCalc.nstpl], chLabel , 4 );
                StopCalc.chStopLabel1[StopCalc.nstpl][4] = 0;

                i = 5;
                /* get line wavelength */
                StopCalc.StopLineWl1[StopCalc.nstpl] = 
                        (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH, &lgEOL);

                /* line was entered, look for possible micron or cm label */
                if( input.chCARDCAPS[i-1] == 'M' )
                {
                        /* microns */
                        StopCalc.StopLineWl1[StopCalc.nstpl] *= 1e4;
                }
                else if( input.chCARDCAPS[i-1] == 'C' )
                {
                        /* microns */
                        StopCalc.StopLineWl1[StopCalc.nstpl] *= 1e8;
                }

                /* get relative intensity */
                StopCalc.stpint[StopCalc.nstpl] = 
                        (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH, &lgEOL);
                if( lgEOL )
                {
                        fprintf( ioQQQ, " There MUST be a relative intensity  entered "
                                "for first line in STOP LINE command.  Sorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                /* check for second line - use Hbeta is not specified */
                j = i;
                a = FFmtRead(chCard,&j,INPUT_LINE_LENGTH, &lgEOL);

                if( lgEOL )
                {
                        /* hit the eol, normalization line will be H beta */
                        strncpy( StopCalc.chStopLabel2[StopCalc.nstpl], "TOTL" , 4 );
                        StopCalc.chStopLabel2[StopCalc.nstpl][4] = 0;
                        StopCalc.StopLineWl2[StopCalc.nstpl] = 4861.f;
                }
                else
                {
                        /* get label for the line - must do this first so we clear the label string before
                         * trying to read the wavelength */
                        GetQuote( chLabel , chCard , true );
                        /* copy first four char of label into caps, and null terminate*/
                        strncpy( StopCalc.chStopLabel2[StopCalc.nstpl], chLabel , 4 );
                        StopCalc.chStopLabel2[StopCalc.nstpl][4] = 0;

                        /* wavelength of second line, may be absent and so zero -
                         * we will use Hbeta if not specified */
                        StopCalc.StopLineWl2[StopCalc.nstpl] = 
                                (realnum)FFmtRead(chCard,&i, INPUT_LINE_LENGTH,&lgEOL);

                        /* line was entered, look for possible micron or cm label */
                        if( input.chCARDCAPS[i-1] == 'M' )
                        {
                                /* microns */
                                StopCalc.StopLineWl2[StopCalc.nstpl] *= 1e4;
                        }
                        else if( input.chCARDCAPS[i-1] == 'C' )
                        {
                                /* microns */
                                StopCalc.StopLineWl2[StopCalc.nstpl] *= 1e8;
                        }
                }
                /* increment number of stop lines commands entered */
                StopCalc.nstpl = MIN2(StopCalc.nstpl+1,MXSTPL-1);
        }

        else if( nMatch("NTOT" , chCard ) && nMatch("ALIO" , chCard ) )
        {
                /* this counter is incremented at end of conv_base - total number of
                 * times conv base was called 
                 * this is a debugging aid - code aborts */
                StopCalc.nTotalIonizStop = (long)a;
        }

        /* oops! no keyword that we could find */
        else
        {
                fprintf( ioQQQ, " I did not recognize a keyword on this STOP line, line image follows;\n" );
                fprintf( ioQQQ, " \"%s\"\n Sorry.\n" , chCard);
                cdEXIT(EXIT_FAILURE);
        }
        return;
}

---