parse_blackbody.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 */
/*ParseBlackbody parse parameters off black body command */
#include "cddefines.h"
#include "physconst.h"
#include "optimize.h"
#include "input.h" 
#include "rfield.h"
#include "radius.h"
#include "parse.h"
#include "parser.h"

void ParseBlackbody(
  /* input command line, already changed to caps */
  Parser &p)
{
        bool 
                lgIntensitySet=false;
        double a, 
          dil, 
          rlogl;

        char chParamType[20] = "";
        int nParam = 0;

        DEBUG_ENTRY( "ParseBlackbody()" );

        set_NaN( rlogl );

        /* type is blackbody */
        strcpy( rfield.chSpType[rfield.nShape], "BLACK" );
        strcpy( rfield.chSpNorm[p.m_nqh], "LUMI" );

        /* these two are not used for this continuum shape */
        rfield.cutoff[rfield.nShape][0] = 0.;
        rfield.cutoff[rfield.nShape][1] = 0.;

        /* get the blackbody temperature */
        rfield.slope[rfield.nShape] = p.FFmtRead();
        if( p.lgEOL() )
                p.NoNumb("blackbody temperature");

        /* this is the temperature - make sure its linear in the end
         * there are two keys, LINEAR and LOG, that could be here,
         * else choose which is here by which side of 10 */
        if( (rfield.slope[rfield.nShape] <= 10. && !p.nMatch("LINE")) || 
                p.nMatch(" LOG") )
        {
                /* log option */
                rfield.slope[rfield.nShape] = pow(10.,rfield.slope[rfield.nShape]);
        }

        /* check that temp is not too low - could happen if log misused */
        if( rfield.slope[rfield.nShape] < 1e4 )
        {
                fprintf( ioQQQ, " Is T(star)=%10.2e correct???\n", 
                  rfield.slope[rfield.nShape] );
        }

        /* now check that temp not too low - would peak below low
         * energy limit of the code
         * factor is temperature of 1 Ryd, egamry is high-energy limit of code */
        if( rfield.slope[rfield.nShape]/TE1RYD < rfield.emm )
        {
                fprintf( ioQQQ, " This temperature is very low - the blackbody will have significant flux low the low energy limit of the code, presently %10.2e Ryd.\n", 
                  rfield.emm );
                fprintf( ioQQQ, " Was this intended?\n" );
        }

        /* now check that temp not too high - would extend beyond high
         * energy limit of the code
         * factor is temperature of 1 Ryd, egamry is high-energy limit of code */
        if( rfield.slope[rfield.nShape]/TE1RYD*2. > rfield.egamry )
        {
                fprintf( ioQQQ, " This temperature is very high - the blackbody will have significant flux above the high-energy limit of the code,%10.2e Ryd.\n", 
                  rfield.egamry );
                fprintf( ioQQQ, " Was this intended?\n" );
        }

        /* also possible to input log(total luminosity)=real log(l) */
        a = p.FFmtRead();

        /* there was not a second number on the line; check if LTE or STE */
        if( p.nMatch(" LTE") || p.nMatch("LTE ") ||
            p.nMatch(" STE") || p.nMatch("STE ") )
        {
                /* set energy density to the STE - strict thermodynamic equilibrium - value */
                strcpy( chParamType , "STE" );
                nParam = 1;

                if( !p.lgEOL() )
                {
                        fprintf(ioQQQ,"PROBLEM the luminosity was specified on "
                                "the BLACKBODY K STE command.\n");
                        fprintf(ioQQQ,"Do not specify the luminosity since STE does this.\n");
                        fprintf( ioQQQ, " Sorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                /* use blackbody relations to get intensity from temperature */
                rlogl = log10(4.*STEFAN_BOLTZ) + 4.*log10(rfield.slope[rfield.nShape]);

                /* set radius to very large value if not already set */
                if( !radius.lgRadiusKnown )
                {
                        radius.Radius = pow(10.,radius.rdfalt);
                }

                strcpy( rfield.chRSpec[p.m_nqh], "SQCM" );
                lgIntensitySet = true;
        }

        /* a second number was entered, what was it? */
        else if( p.nMatch("LUMI") )
        {
                strcpy( chParamType , "LUMINOSITY" );
                nParam = 2;
                rlogl = a;
                strcpy( rfield.chRSpec[p.m_nqh], "4 PI" );
                if( p.lgEOL() )
                        p.NoNumb("luminosity" );
                lgIntensitySet = true;
        }

        else if( p.nMatch("RADI") )
        {
                strcpy( chParamType , "RADIUS" );
                nParam = 2;
                /* radius was entered, convert to total luminosity */
                rlogl = -3.147238 + 2.*a + 4.*log10(rfield.slope[rfield.nShape]);
                strcpy( rfield.chRSpec[p.m_nqh], "4 PI" );
                if( p.lgEOL() )
                        p.NoNumb("radius" );
                lgIntensitySet = true;
        }

        else if( p.nMatch("DENS") )
        {
                strcpy( chParamType , "ENERGY DENSITY" );
                nParam = 2;
                /* number was temperature to deduce energy density
                 * number is linear if greater than 10, or if LINEAR appears on line
                 * want number to be log of temperature at end of this */
                if( !p.nMatch(" LOG") && (p.nMatch("LINE") || a > 10.) )
                {
                        a = log10(a);
                }
                rlogl = log10(4.*STEFAN_BOLTZ) + 4.*a;
                /* set radius to very large value if not already set */
                if( !radius.lgRadiusKnown )
                {
                        radius.Radius = pow(10.,radius.rdfalt);
                }
                strcpy( rfield.chRSpec[p.m_nqh], "SQCM" );
                if( p.lgEOL() )
                        p.NoNumb("energy density");
                lgIntensitySet = true;
        }

        else if( p.nMatch("DILU") )
        {
                strcpy( chParamType , "DILUTION FACTOR" );
                nParam = 2;
                /* number is dilution factor, if negative then its log */
                if( a <= 0. )
                        dil = a;
                else
                        dil = log10(a);

                if( dil > 0. )
                        fprintf( ioQQQ, "PROBLEM Is the dilution factor > 1 on this "
                        "blackbody command physical?\n" );

                /* intensity from black body relations and temperature */
                rlogl = log10(4.*STEFAN_BOLTZ) + 4.*log10(rfield.slope[rfield.nShape]);

                /* add on dilution factor */
                rlogl += dil;

                /* set radius to very large value if not already set */
                if( !radius.lgRadiusKnown )
                {
                        radius.Radius = pow(10.,radius.rdfalt);
                }
                strcpy( rfield.chRSpec[p.m_nqh], "SQCM" );
                if( p.lgEOL() )
                        p.NoNumb("dilution factor" );
                lgIntensitySet = true;
        }

        else if( p.nMatch("DISK") )
        {
                if( p.lgEOL() )
                        p.NoNumb("disk Te" );

                strcpy( chParamType , "DISK" );
                nParam = 2;
                
                rfield.cutoff[rfield.nShape][0] = a;
                /* this is the temperature - make sure its linear in the end
                 * there are two keys, LINEAR and LOG, that could be here,
                 * else choose which is here by which side of 10 */
                if( (rfield.cutoff[rfield.nShape][0] <= 10. && !p.nMatch("LINE")) ||
                        p.nMatch(" LOG") )
                {
                        /* log option */
                        rfield.cutoff[rfield.nShape][0] = pow(10.,rfield.cutoff[rfield.nShape][0]);
                }
                a = log10( rfield.cutoff[rfield.nShape][0] );

                strcpy( rfield.chSpType[rfield.nShape], "DISKB" );
                lgIntensitySet = false;
        }

        if( lgIntensitySet )
        {
                /* a luminosity option was specified
                 * check that stack of shape and luminosity specifications
                 * is parallel, stop if not - this happens is background comes
                 * BETWEEN another set of shape and luminosity commands */
                if( rfield.nShape != p.m_nqh )
                {
                        fprintf( ioQQQ, " This command has come between a previous ordered pair of continuum shape and luminosity commands.\n Reorder the commands to complete each continuum specification before starting another.\n" );
                        fprintf( ioQQQ, " Sorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                rfield.range[p.m_nqh][0] = rfield.emm;
                rfield.range[p.m_nqh][1] = rfield.egamry;
                rfield.totpow[p.m_nqh] = rlogl;
                ++p.m_nqh;
        }
        /* vary option */
        if( optimize.lgVarOn )
        {
                /* this test no option on blackbody command */
                if( strcmp(chParamType,"") == 0 )
                {
                        /* no luminosity options on vary */
                        optimize.nvarxt[optimize.nparm] = 1;
                        strcpy( optimize.chVarFmt[optimize.nparm], "BLACKbody= %f LOG" );
                }
                else
                {
                        char chHold[100];
                        /* there was an option - honor it */
                        if( nParam==1 )
                        {
                                optimize.nvarxt[optimize.nparm] = 1;
                                strcpy( chHold , "BLACKbody= %f LOG ");
                                strcat( chHold , chParamType );
                        }
                        else if( nParam==2 )
                        {
                                optimize.nvarxt[optimize.nparm] = 2;
                                optimize.vparm[1][optimize.nparm] = (realnum)a;
                                strcpy( chHold , "BLACKbody= %f LOG %f ");
                                strcat( chHold , chParamType );
                        }
                        else
                                TotalInsanity();
                        strcpy( optimize.chVarFmt[optimize.nparm], chHold );
                }

                /* pointer to where to write */
                optimize.nvfpnt[optimize.nparm] = input.nRead;
                /* log of temp stored here  */
                optimize.vparm[0][optimize.nparm] = (realnum)log10(rfield.slope[rfield.nShape]);
                /* the increment in the first steps away from the original value */
                optimize.vincr[optimize.nparm] = 0.5f;
                ++optimize.nparm;
        }

        /* increment SED indices */
        ++rfield.nShape;
        if( rfield.nShape >= LIMSPC )
        {
                fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
                cdEXIT(EXIT_FAILURE);
        }

        return;
}