/home66/gary/public_html/cloudy/c08_branch/source/parse_constant.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 */
/*ParseConstant parse parameters from the 'constant ...' command */
#include "cddefines.h"
#include "physconst.h"
#include "optimize.h"
#include "thermal.h"
#include "dense.h"
#include "pressure.h"
#include "phycon.h"
#include "input.h"
#include "parse.h"

void ParseConstant(char *chCard )
{
        bool lgEOL;
        long int i;

        DEBUG_ENTRY( "ParseConstant()" );

        if( nMatch("GRAI",chCard) && nMatch("TEMP",chCard) )
        {
                /* constant grain temperature command */
                i = 5;
                thermal.ConstGrainTemp = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

                /* if linear option is not on the line, convert to exponent if <= 10 */
                if( !nMatch("LINE",chCard) )
                {
                        if( thermal.ConstGrainTemp <= 10. )
                                thermal.ConstGrainTemp = (realnum)pow((realnum)10.f,thermal.ConstGrainTemp);
                }

                if( lgEOL )
                {
                        NoNumb(chCard);
                }
        }

        else if( nMatch("TEMP",chCard) )
        {
                /* a constant temperature model */
                thermal.lgTemperatureConstant = true;
                thermal.lgTemperatureConstantCommandParsed = true;

                /* this is an option to specify the temperature in different units
                 * keV, eV now supported */
                realnum convert_to_Kelvin = 1;
                if( nMatch(" EV ",chCard) )
                {
                        convert_to_Kelvin = (realnum)EVDEGK;
                }
                else if( nMatch(" KEV",chCard) )
                {
                        convert_to_Kelvin = (realnum)(EVDEGK * 1000.);
                }

                i = 5;
                /* this is the "force" temperature.  same var used in force temp
                 * command, but lgTSetOn is not set so then allowed to vary 
                 * so constant temperature requires both lgTSetOn true and ConstTemp > 0 */
                thermal.ConstTemp = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                if( lgEOL )
                        NoNumb(chCard);

                /* if linear option is not on the line and T<=10, assume number is log */
                if( !nMatch("LINE",chCard) )
                {
                        if( thermal.ConstTemp <= 10. )
                                thermal.ConstTemp = (realnum)pow((realnum)10.f,thermal.ConstTemp);
                }
                /* do units conversion here */
                thermal.ConstTemp *= convert_to_Kelvin;

                /* check that temperature is not less than 3K */
                if( thermal.ConstTemp < 3. )
                {
                        fprintf( ioQQQ, " PROBLEM TE<3K, reset to 3K.\n" );
                        thermal.ConstTemp = 3.;
                }

                /* bail if temperature is too high */
                /* >>chng 06 feb 14, can cannot do a temperature of 1e10 - gaunt factors overrun
                 * energy array.  temp must be less than 1e10 
                 * >>chng 06 feb 27, change 1e10 to phycon.TEMP_LIMIT_HIGH - which is the same */
                if( thermal.ConstTemp > phycon.TEMP_LIMIT_HIGH )
                {
                        fprintf( ioQQQ, " The constant temperature (%.3e K) is > %.3e K, Cloudy cannot handle this temperature.\n",
                                thermal.ConstTemp,
                                phycon.TEMP_LIMIT_HIGH);
                        fprintf( ioQQQ, " The upper temperature limit is %.3e K.\n" ,
                                phycon.TEMP_LIMIT_HIGH);
                        fprintf( ioQQQ, " Sorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                /* set the real electron temperature to the forced value */
                TempChange(thermal.ConstTemp , false);

                /* vary option */
                if( optimize.lgVarOn )
                {
                        /*  no luminosity options on vary */
                        optimize.nvarxt[optimize.nparm] = 1;
                        strcpy( optimize.chVarFmt[optimize.nparm], "CONStant TEMP %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(phycon.te);
                        optimize.vincr[optimize.nparm] = 0.1f;
                        ++optimize.nparm;
                }
        }

        else if( nMatch("DENS",chCard) )
        {
                /* constant density */
                strcpy( dense.chDenseLaw, "CDEN" );
                /* turn off radiation pressure */
                pressure.lgPres_radiation_ON = false;
                pressure.lgPres_magnetic_ON = false;
                pressure.lgPres_ram_ON = false;
        }

        else if( nMatch("PRES",chCard) )
        {
                /* constant pressure  */
                strcpy( dense.chDenseLaw, "CPRE" );

                /* >>chng 06 jun 20, add reset option, to reset density to keep 
                 * initial pressure itself constant from iteration to iteration, 
                 * rather than initial density */
                if( nMatch("RESE",chCard) )
                {
                        /* this says not to keep initial density constant, 
                         * reset it to keep pressure const */
                        dense.lgDenseInitConstant = false;
                }
                else
                {
                        /* this is default, says keep initial density constant, 
                         * so pressure from iter to iter not really const */
                        dense.lgDenseInitConstant = true;
                }

                if( nMatch(" GAS",chCard) )
                {
                        /*  constant gas pressure (no radiation)
                         *  turn off radiation pressure */
                        pressure.lgPres_radiation_ON = false;

                        /*  turn off incident continuum */
                        pressure.lgContRadPresOn = false;

                        /* turn off magnetic and ram pressure */
                        pressure.lgPres_magnetic_ON = false;
                        pressure.lgPres_ram_ON = false;

                        /*  optional second number is power law index */
                        i = 4;
                        pressure.PresPowerlaw = (realnum)FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                }

                else
                {
                        /*  constant total pressure, gas+rad+incident continuum
                         *  turn on radiation pressure */
                        pressure.lgPres_radiation_ON = true;
                        pressure.lgPres_magnetic_ON = true;
                        pressure.lgPres_ram_ON = true;

                        /*  option to turn off continuum pressure */
                        if( nMatch("O CO",chCard) )
                        {
                                pressure.lgContRadPresOn = false;
                        }
                        else
                        {
                                pressure.lgContRadPresOn = true;
                        }

                        /*  option to not abort when too much radiation pressure, no abort */
                        if( nMatch("O AB",chCard) )
                        {
                                pressure.lgRadPresAbortOK = false;
                        }
                        else
                        {
                                pressure.lgRadPresAbortOK = true;
                        }
                        /*  there is no optional power law option for constant total pressure */
                        pressure.PresPowerlaw = 0.;

                        /* option to set pressure */
                        if( nMatch(" SET",chCard) )
                        {
                                /* number on line is log of nT - option to specify initial pressure */
                                pressure.lgPressureInitialSpecified = true;
                                /* this is log of nT product - if not present then set zero */
                                i = 5;
                                pressure.PressureInitialSpecified = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
                                if( lgEOL )
                                        NoNumb( chCard );
                                else
                                        /* pressure in nkT units */
                                        pressure.PressureInitialSpecified = pow(10.,pressure.PressureInitialSpecified ) *
                                        BOLTZMANN;
                        }
                        else
                                pressure.lgPressureInitialSpecified = false;
                }

                /* vary option */
                if( optimize.lgVarOn &&pressure.lgPressureInitialSpecified )
                {
                        /*  no options on vary */
                        optimize.nvarxt[optimize.nparm] = 1;
                        strcpy( optimize.chVarFmt[optimize.nparm], "CONStant PRESsure SET %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(pressure.PressureInitialSpecified);
                        optimize.vincr[optimize.nparm] = 0.1f;
                        ++optimize.nparm;
                }
        }

        else
        {
                /* no keys were recognized */
                fprintf( ioQQQ, " The keyword should be TEMPerature, DENSity, GAS or PRESsure, sorry.\n" );
                cdEXIT(EXIT_FAILURE);
        }
        return;
}

---