parse_constant.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2013 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"
#include "parser.h"
 
void ParseConstant(Parser &p )
{
        DEBUG_ENTRY( "ParseConstant()" );
 
        if( p.nMatch("GRAI") && p.nMatch("TEMP") )
        {
                /* constant grain temperature command */
                thermal.ConstGrainTemp = (realnum)p.FFmtRead();
 
                /* if linear option is not on the line, convert to exponent if <= 10 */
                if( !p.nMatch("LINE") )
                {
                        if( thermal.ConstGrainTemp <= 10. )
                                thermal.ConstGrainTemp = (realnum)pow((realnum)10.f,thermal.ConstGrainTemp);
                }
 
                if( p.lgEOL() )
                {
                        p.NoNumb("grain temperature");
                }
        }
 
        else if( p.nMatch("TEMP") )
        {
                /* 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( p.nMatch(" EV ") )
                {
                        convert_to_Kelvin = (realnum)EVDEGK;
                }
                else if( p.nMatch(" KEV") )
                {
                        convert_to_Kelvin = (realnum)(EVDEGK * 1000.);
                }
 
                /* 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)p.FFmtRead();
                if( p.lgEOL() )
                        p.NoNumb("temperature");
 
                /* if linear option is not on the line and T<=10, assume number is log */
                if( p.nMatch(" LOG") || (thermal.ConstTemp <= 10. && !p.nMatch("LINE")) )
                {
                        if( thermal.ConstTemp > log10(BIGFLOAT) )
                        {
                                fprintf(ioQQQ," PROBLEM temperature entered as a log but is too large "\
                                        "for this processor.  I am interpreting it as the linear temperature.\n");
                        }
                        else
                                thermal.ConstTemp = (realnum)pow((realnum)10.f,thermal.ConstTemp);
                }
                /* do units conversion here */
                thermal.ConstTemp *= convert_to_Kelvin;
 
                /* check temperature bounds */
                if( thermal.ConstTemp < phycon.TEMP_LIMIT_LOW )
                {
                        thermal.ConstTemp = (realnum)(1.0001*phycon.TEMP_LIMIT_LOW);
                        fprintf( ioQQQ, " PROBLEM Te too low, reset to %g K.\n",
                                 thermal.ConstTemp );
                }
                if( thermal.ConstTemp > phycon.TEMP_LIMIT_HIGH )
                {
                        thermal.ConstTemp = (realnum)(0.9999*phycon.TEMP_LIMIT_HIGH);
                        fprintf( ioQQQ, " PROBLEM Te too high, reset to %g K.\n",
                                 thermal.ConstTemp );
                }
 
                /* 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;
                        // the keyword LOG is not used above, but is checked elsewhere
                        strcpy( optimize.chVarFmt[optimize.nparm], "CONSTANT TEMP %f LOG" );
 
                        /*  pointer to where to write */
                        optimize.nvfpnt[optimize.nparm] = input.nRead;
 
                        /*  log of temp will be pointer */
                        optimize.vparm[0][optimize.nparm] = (realnum)log10(thermal.ConstTemp);
                        optimize.vincr[optimize.nparm] = 0.1f;
                        optimize.varang[optimize.nparm][0] = (realnum)log10(1.00001*phycon.TEMP_LIMIT_LOW);
                        optimize.varang[optimize.nparm][1] = (realnum)log10(0.99999*phycon.TEMP_LIMIT_HIGH);
                        ++optimize.nparm;
                }
        }
 
        else if( p.nMatch("DENS") )
        {
                /* 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( p.nMatch("PRES") )
        {
                /* 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( p.nMatch("RESE") )
                {
                        /* 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( p.nMatch("TIME") )
                {
                        // pressure varies as a function of time
                        /* this says not to keep initial density constant,
                         * reset it to keep pressure const */
                        dense.lgDenseInitConstant = false;
                        dense.lgPressureVaryTime = true;
 
                        //  required number is timescale for time variation
                        dense.PressureVaryTimeTimescale = (realnum)p.FFmtRead();
                        if( dense.PressureVaryTimeTimescale<SMALLFLOAT )
                        {
                                // need two numbers
                                fprintf(ioQQQ," PROBLEM the constant pressure time command requires"
                                        " a positive timescale.\n");
                                cdEXIT(EXIT_FAILURE);
                        }
                        //  optional number is index for time variation
                        dense.PressureVaryTimeIndex = (realnum)p.FFmtRead();
                        // pressure will be initial pressure - pressure.PresTotlInit - multiplied by
                        // (time / time scale ) ^ index
                        if( p.lgEOL() )
                        {
                                // need two numbers
                                fprintf(ioQQQ," PROBLEM the constant pressure time command requires"
                                        " two numbers, the timescale for the variation and an index.\n");
                                cdEXIT(EXIT_FAILURE);
                        }
                }
 
                if( p.nMatch(" GAS") )
                {
                        /*  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 number is power law index */
                        pressure.PresPowerlaw = (realnum)p.FFmtRead();
                }
 
                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( p.nMatch("NO CO") )
                        {
                                pressure.lgContRadPresOn = false;
                        }
                        else
                        {
                                pressure.lgContRadPresOn = true;
                        }
 
                        /*  option to not abort when too much radiation pressure, no abort */
                        if( p.nMatch("NO AB") )
                        {
                                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( p.nMatch(" SET") )
                        {
                                /* 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 */
                                pressure.PressureInitialSpecified = p.FFmtRead();
                                if( p.lgEOL() )
                                        p.NoNumb("initial pressure" );
                                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;
                        // the keyword LOG is not used above, but is checked elsewhere
                        strcpy( optimize.chVarFmt[optimize.nparm], "CONSTANT PRESSURE SET %f LOG" );
                        if( !dense.lgDenseInitConstant )
                                strcat( optimize.chVarFmt[optimize.nparm], " RESET" );
                        if( !pressure.lgContRadPresOn )
                                strcat( optimize.chVarFmt[optimize.nparm], " NO CONTINUUM" );
                        if( !pressure.lgRadPresAbortOK )
                                strcat( optimize.chVarFmt[optimize.nparm], " NO ABORT" );
 
                        /*  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/BOLTZMANN);
                        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;
}