/home66/gary/public_html/cloudy/c08_branch/source/parse_powerlawcontinuum.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 */
/*ParsePowerlawContinuum parse the power law continuum command */
#include "cddefines.h"
#include "rfield.h"
#include "optimize.h"
#include "input.h"
#include "parse.h"
#include "physconst.h"

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

        DEBUG_ENTRY( "ParsePowerlawContinuum()" );

        /* power law with cutoff and X-ray continuum */
        strcpy( rfield.chSpType[rfield.nspec], "POWER" );

        /* first parameter is slope of continuum, probably should be negative */
        i = 5;
        rfield.slope[rfield.nspec] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,&lgEOL);
        if( lgEOL )
        {
                fprintf( ioQQQ, " There should have been a number on this line.   Sorry.\n" );
                cdEXIT(EXIT_FAILURE);
        }

        if( rfield.slope[rfield.nspec] >= 0. )
        {
                fprintf( ioQQQ, " Is the slope of this power law correct?\n" );
        }

        /* second optional parameter is high energy cut off */
        rfield.cutoff[rfield.nspec][0] = FFmtRead(chCard,&i,INPUT_LINE_LENGTH,
          &lgEOL);

        /* no cutoff if eof hit */
        if( lgEOL )
        {
                /* no extra parameters at all, so put in extreme cutoffs */
                rfield.cutoff[rfield.nspec][0] = 1e4;
                rfield.cutoff[rfield.nspec][1] = 1e-4;
        }
        else
        {
                /* first cutoff was present, check for second */
                rfield.cutoff[rfield.nspec][1] = FFmtRead(chCard,&i,
                  INPUT_LINE_LENGTH,&lgEOL);
                if( lgEOL )
                        rfield.cutoff[rfield.nspec][1] = 1e-4;
        }

        /* check that energies were entered in the correct order */
        if( rfield.cutoff[rfield.nspec][1] > rfield.cutoff[rfield.nspec][0] )
        {
                fprintf( ioQQQ, " The optional cutoff energies do not appear to be in the correct order.  Check Hazy.\n" );
                cdEXIT(EXIT_FAILURE);
        }

        /* check for keyword KELVIN to interpret cutoff energies as degrees */
        if( nMatch("KELV",chCard) )
        {
                /* temps are log if first le 10 */
                if( rfield.cutoff[rfield.nspec][0] <= 10. )
                        rfield.cutoff[rfield.nspec][0] = pow(10.,rfield.cutoff[rfield.nspec][0]);
                if( rfield.cutoff[rfield.nspec][1] <= 10. )
                        rfield.cutoff[rfield.nspec][1] = pow(10.,rfield.cutoff[rfield.nspec][1]);
                rfield.cutoff[rfield.nspec][0] /= TE1RYD;
                rfield.cutoff[rfield.nspec][1] /= TE1RYD;
        }

        if( rfield.cutoff[rfield.nspec][0] < 0. || 
                rfield.cutoff[rfield.nspec][1] < 0. )
        {
                fprintf( ioQQQ, " A negative cutoff energy is not physical.  Sorry.\n" );
                cdEXIT(EXIT_FAILURE);
        }

        if( rfield.cutoff[rfield.nspec][1] == 0. && 
                rfield.slope[rfield.nspec] <= -1. )
        {
                fprintf( ioQQQ, " A power-law with this slope, and no low energy cutoff, may have an unphysically large\n brightness temperature in the radio.\n" );
        }

        /* vary option */
        if( optimize.lgVarOn )
        {
                /* pointer to where to write */
                optimize.nvfpnt[optimize.nparm] = input.nRead;
                if( nMatch("ARYB",chCard) )
                {
                        /* this test is for key "varyb", meaning to vary second parameter
                         * the cutoff temperature
                         * this is the number of parameters to feed onto the input line */
                        optimize.nvarxt[optimize.nparm] = 2;
                        /* vary ?? */
                        sprintf( optimize.chVarFmt[optimize.nparm], 
                                "POWER LAW %f KELVIN%%f %%f", 
                          rfield.slope[rfield.nspec] );
                        /* param is linear scale factor */
                        optimize.vparm[0][optimize.nparm] = (realnum)log10(rfield.cutoff[rfield.nspec][0]*
                          TE1RYD);
                        optimize.vparm[1][optimize.nparm] = (realnum)log10(rfield.cutoff[rfield.nspec][1]*
                          TE1RYD);
                        optimize.vincr[optimize.nparm] = 0.2f;
                }
                else if( nMatch("ARYC",chCard) )
                {
                        /* the keyword was "varyc"
                         * this is the number of parameters to feed onto the input line */
                        optimize.nvarxt[optimize.nparm] = 1;

                        /* vary ?? */
                        sprintf( optimize.chVarFmt[optimize.nparm], "POWER LAW%f %f KELVIN %%f )", 
                          rfield.slope[rfield.nspec], 
                          log10(rfield.cutoff[rfield.nspec][0]* TE1RYD) );

                        optimize.vparm[0][optimize.nparm] = (realnum)log10(rfield.cutoff[rfield.nspec][1]*
                          TE1RYD);
                        optimize.vincr[optimize.nparm] = 0.2f;
                }
                else
                {
                        /* vary the first parameter only, but still are two more
                         * this is the number of parameters to feed onto the input line */
                        optimize.nvarxt[optimize.nparm] = 3;
                        strcpy( optimize.chVarFmt[optimize.nparm], 
                                "POWER LAW KELVIN%f %f %f" );
                        /* param is slope of the power law */
                        optimize.vparm[0][optimize.nparm] = (realnum)rfield.slope[rfield.nspec];
                        optimize.vparm[1][optimize.nparm] = (realnum)log10(rfield.cutoff[rfield.nspec][0]*
                          TE1RYD);
                        optimize.vparm[2][optimize.nparm] = (realnum)log10(rfield.cutoff[rfield.nspec][1]*
                          TE1RYD);
                        optimize.vincr[optimize.nparm] = 0.2f;
                }
                ++optimize.nparm;
        }

        /*>>chng 06 nov 10, BUGFIX, nspec was incremented before previous branch
         * and so crashed with log10 0 since nspec was beyond set values
         * caused fpe domain function error with log 0
         * fpe caught by Pavel Abolmasov */
        ++rfield.nspec;
        if( rfield.nspec >= LIMSPC )
        {
                fprintf( ioQQQ, " Too many continua entered; increase LIMSPC\n" );
                cdEXIT(EXIT_FAILURE);
        }

        return;
}

---