/home66/gary/public_html/cloudy/c08_branch/source/parse_blackbody.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 */
/*ParseBlackbody parse parameters off black body command */
#include "cddefines.h"
#include "physconst.h"
#include "rfield.h"
#include "radius.h"
#include "parse.h"

void ParseBlackbody(char *chCard,       /* input command line, already changed to caps */
  /* counter for which continuum source this is */
  long int *nqh,
  /* optional area that might be set here */
  realnum *ar1)
{
        bool lgEOL;
        long int i;
        double a, 
          dil, 
          rlogl;

        DEBUG_ENTRY( "ParseBlackbody()" );

        /* type is blackbody */
        strcpy( rfield.chSpType[rfield.nspec], "BLACK" );

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

        /* get the temperature */
        i = 5;
        rfield.slope[rfield.nspec] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);

        /* there must be at least one number, the temperature */
        if( lgEOL )
        {
                fprintf( ioQQQ, " There must be at least 1 number on a BLACKBODY command line.  Sorry.\n" );
                cdEXIT(EXIT_FAILURE);
        }

        /* 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.nspec] <= 10. && !nMatch("LINE",chCard)) || 
                nMatch(" LOG",chCard) )
        {
                /* log option */
                rfield.slope[rfield.nspec] = pow(10.,rfield.slope[rfield.nspec]);
        }

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

        /* 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.nspec]/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.nspec]/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" );
        }

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

        /* also possible to input log(total luminosity)=real log(l) */
        a = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
        if( lgEOL )
        {
                /* there was not a second number on the line; check if LTE or STE */
                if( nMatch(" LTE",chCard) || nMatch("LTE ",chCard) ||
                    nMatch(" STE",chCard) || nMatch("STE ",chCard) )
                {
                        /* set energy density to the STE - strict thermodynamic equilibrium - value */
                        strcpy( rfield.chSpNorm[*nqh], "LUMI" );

                        /* need nspec-1 since was incremented above */
                        a = log10(rfield.slope[rfield.nspec-1]);
                        rlogl = log10(4.*STEFAN_BOLTZ) + 4.*a;

                        /* set radius to very large value if not already set */
                        /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
                        if( !radius.lgRadiusKnown )
                        {
                                *ar1 = (realnum)radius.rdfalt;
                                radius.Radius = pow(10.,radius.rdfalt);
                        }
                        strcpy( rfield.chRSpec[*nqh], "SQCM" );
                        rfield.range[*nqh][0] = rfield.emm;
                        rfield.range[*nqh][1] = rfield.egamry;
                        rfield.totpow[*nqh] = rlogl;
                        *nqh = MIN2(*nqh+1,LIMSPC);

                        /* 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.nspec != *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);
                        }
                }

                return;
        }

        strcpy( rfield.chSpNorm[*nqh], "LUMI" );

        /* a second number was entered, what was it? */
        if( nMatch("LUMI",chCard) )
        {
                rlogl = a;
                strcpy( rfield.chRSpec[*nqh], "4 PI" );
        }

        else if( nMatch("RADI",chCard) )
        {
                /* radius was entered, convert to total lumin */
                rlogl = -3.147238 + 2.*a + 4.*log10(rfield.slope[rfield.nspec-1]);
                strcpy( rfield.chRSpec[*nqh], "4 PI" );
        }

        else if( nMatch("DENS",chCard) )
        {
                /* 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( nMatch("LINE",chCard) || a > 10. )
                {
                        a = log10(a);
                }
                rlogl = log10(4.*STEFAN_BOLTZ) + 4.*a;
                /* set radius to very large value if not already set */
                /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
                if( !radius.lgRadiusKnown )
                {
                        *ar1 = (realnum)radius.rdfalt;
                        radius.Radius = pow(10.,radius.rdfalt);
                }
                strcpy( rfield.chRSpec[*nqh], "SQCM" );
        }

        else if( nMatch("DILU",chCard) )
        {
                /* number is dilution factor, if negative then its log */
                if( a < 0. )
                {
                        dil = a;
                }
                else
                {
                        dil = log10(a);
                }
                if( dil > 0. )
                {
                        fprintf( ioQQQ, " Is a dilution factor greater than one physical?\n" );
                }

                /* this is LTE energy density */
                a = log10(rfield.slope[rfield.nspec-1]);
                rlogl = log10(4.*STEFAN_BOLTZ) + 4.*a;

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

                /* set radius to very large value if not already set */
                /* >>chng 01 jul 24, from Radius == 0 to this, as per PvH comments */
                if( !radius.lgRadiusKnown )
                {
                        *ar1 = (realnum)radius.rdfalt;
                        radius.Radius = pow(10.,radius.rdfalt);
                }
                strcpy( rfield.chRSpec[*nqh], "SQCM" );

        }
        else
        {
                fprintf( ioQQQ, " I do not understand what the second number was- keywords are LUMINOSITY, RADIUS, DILUTION, and DENSITY.\n" );
                cdEXIT(EXIT_FAILURE);
        }

        rfield.range[*nqh][0] = rfield.emm;
        rfield.range[*nqh][1] = rfield.egamry;
        rfield.totpow[*nqh] = rlogl;
        /* the second number will be some sort of luminosity */
        *nqh = MIN2(*nqh+1,LIMSPC);

        /* 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.nspec != *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);
        }

        return;
}

---