root/branches/newmole/source/cont_ffun.cpp

/* 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 */
/*ffun evaluate total flux for sum of all continuum sources */
/*ffun1 derive flux at a specific energy, for one continuum */
/*ReadTable called by TABLE READ to read in continuum from PUNCH TRANSMITTED CONTINUUM */
#include "cddefines.h"
#include "physconst.h"
#include "rfield.h"
#include "ipoint.h"
#include "opacity.h"
#include "continuum.h"

/*ReadTable called by TABLE READ to read in continuum from PUNCH TRANSMITTED CONTINUUM */
STATIC void ReadTable(const string& fnam);

/* evaluate sum of all individual continua at one energy, return is
* continuum intensity */
double ffun(
                        /* the energy in Rydbergs where the continuum will be evaluated */
                        double anu )
{
        double frac_beam_time;
        /* fraction of beamed continuum that is constant */
        double frac_beam_const;
        /* fraction of continuum that is isotropic */
        double frac_isotropic;
        double a;

        DEBUG_ENTRY( "ffun()" );

        /* call real function */
        a = ffun( anu , &frac_beam_time , &frac_beam_const , &frac_isotropic );
        return a;
}

/* evaluate sum of all individual continua at one energy, return is
 * continuum intensity */
double ffun(
        /* the energy in Rydbergs where the continuum will be evaluated */
        double anu , 
        /* fraction of beamed continuum that is varies with time */
        double *frac_beam_time,
        /* fraction of beamed continuum that is constant */
        double *frac_beam_const,
        /* fraction of continuum that is isotropic */
        double *frac_isotropic )
{
        double ffun_v;
        static bool lgWarn = false;
        double flx_beam_time , flx_beam_const , flx_isotropic;

        DEBUG_ENTRY( "ffun()" );

        /* This routine, ffun, returns the sum of photons per unit time, area, energy,
         * for all continua in the calculation.  ffun1 is called and returns 
         * a single continuum.  We loop over all nspec continuum sources - 
         * ipspec points to each */
        ffun_v = 0.;
        flx_beam_time = 0.;
        flx_beam_const = 0.;
        flx_isotropic = 0.;
        for( rfield.ipspec=0; rfield.ipspec < rfield.nspec; rfield.ipspec++ )
        {
                double one = ffun1(anu)*rfield.spfac[rfield.ipspec];
                ffun_v += one;

                /* find fraction of total that is constant vs variable and 
                 * isotropic vs beamed */
                if( rfield.lgBeamed[rfield.ipspec] )
                {
                        if( rfield.lgTimeVary[rfield.ipspec] )
                                flx_beam_time += one;
                        else
                                flx_beam_const += one;
                }
                else
                        flx_isotropic += one;
        }

        /* at this point rfield.flux is the sum of the following three continua
         * now keep track of the three different types */
        if( ffun_v < SMALLFLOAT )
        {
                *frac_beam_time = 0.;
                *frac_beam_const = 1.;
                *frac_isotropic = 0.;
        }
        else
        {
                /* fraction of beamed continuum that varies with time */
                *frac_beam_time = flx_beam_time / ffun_v;
                /* part of beamed continuum that is constant */
                *frac_beam_const = flx_beam_const / ffun_v;
                /* the constant isotropic continuum */
                *frac_isotropic = flx_isotropic / ffun_v;
        }
        ASSERT( *frac_beam_time >=0. && *frac_beam_time<=1.+3.*DBL_EPSILON );
        ASSERT( *frac_beam_const >=0.&& *frac_beam_const<=1.+3.*DBL_EPSILON );
        ASSERT( *frac_isotropic >=0. && *frac_isotropic<=1.+3.*DBL_EPSILON );
        ASSERT( fabs( 1.-*frac_beam_time-*frac_beam_const-*frac_isotropic)<
                10.*DBL_EPSILON);

        if( ffun_v > BIGFLOAT && !lgWarn )
        {
                lgWarn = true;
                fprintf( ioQQQ, " FFUN:  The net continuum is very intense.\n" );
                fprintf( ioQQQ, " I will try to press on, but may have problems.\n" );
        }
        return( ffun_v );
}

/*ffun1 derive flux at a specific energy, for one continuum */
double ffun1(double xnu)
{
        char chKey[6];
        long int i;
        double fac, 
          ffun1_v, 
          y;
        static bool lgWarn = false;

        DEBUG_ENTRY( "ffun1()" );


        /* confirm that pointer is within range */
        ASSERT( rfield.ipspec >= 0);
        ASSERT( rfield.ipspec < rfield.nspec );

        /* FFUN1 returns photons per unit time, area, energy, for one continuum
         * ipspec is the pointer to the continuum source, in the order
         * entered on the command lines */

        /*begin sanity check */
        ASSERT( xnu >= rfield.emm*0.99 );
        ASSERT( xnu <= rfield.egamry*1.01 );
        /*end sanity check */

        strcpy( chKey, rfield.chSpType[rfield.ipspec] );

        if( strcmp(chKey,"AGN  ") == 0 )
        {
                /* power law with cutoff at both ends
                 * nomenclature all screwed up - slope is cutoff energy in Ryd,
                 * cutoff[1][i] is ratio of two continua from alpha ox
                 * cutoff[2][i] is slope of bb temp */
                ffun1_v = pow(xnu,-1. + rfield.cutoff[rfield.ipspec][1])*
                  sexp(xnu/rfield.slope[rfield.ipspec])*sexp(0.01/
                  xnu);
                /* only add on x-ray for energies above 0.1 Ryd */
                if( xnu > 0.1 )
                {
                        if( xnu < 7350. )
                        {
                                /* cutoff is actually normalization constant
                                 * below 100keV continuum is nu-1 */
                                ffun1_v += pow(xnu,rfield.cutoff[rfield.ipspec][2] - 
                                  1.)*rfield.cutoff[rfield.ipspec][0]*sexp(1./
                                  xnu);
                        }
                        else
                        {
                                ffun1_v += pow(7350.,rfield.cutoff[rfield.ipspec][2] - 
                                  1.)*rfield.cutoff[rfield.ipspec][0]/
                                  POW3(xnu/7350.);
                        }
                }

        }
        else if( strcmp(chKey,"POWER") == 0 )
        {
                /* power law with cutoff at both ends */
                ffun1_v = pow(xnu,-1. + rfield.slope[rfield.ipspec])*
                  sexp(xnu/rfield.cutoff[rfield.ipspec][0]+rfield.cutoff[rfield.ipspec][1]/
                  xnu);

        }
        else if( strcmp(chKey,"BLACK") == 0 )
        {
                const double db_log = log(DBL_MAX);

                /* black body */
                fac = TE1RYD*xnu/rfield.slope[rfield.ipspec];
                /* >>>chng 00 apr 13 from 80 to log(dbl_max) */
                if( fac > db_log )
                {
                        ffun1_v = 0.;
                }
                else if( fac > 1.e-5 )
                {
                        ffun1_v = xnu*xnu/(exp(fac) - 1.);
                }
                else
                {
                        ffun1_v = xnu*xnu/(fac*(1. + fac/2.));
                }

        }
        else if( strcmp(chKey,"INTER") == 0 )
        {
                /* interpolate on tabulated input spectrum, factor of 1.0001 to 
                 * make sure that requested energy is within bounds of array */
                if( xnu >= rfield.tNuRyd[rfield.ipspec][0]*1.000001 )
                {
                        /* loop starts at second array element, [1], since want to 
                         * find next continuum energy greater than desired point */
                        i = 1;
                        /* up to NCELL tabulated points may be read in.  Very fine
                         * continuum mesh such as that output by stellar atmospheres 
                         * can have very large number of points */
                        while( i< NCELL-1 && rfield.tNuRyd[rfield.ipspec][i]>0. )
                        {
                                if( xnu < rfield.tNuRyd[rfield.ipspec][i] )
                                {
                                        /* the energy xnu is between points rfield.tNuRyd[rfield.ipspec][i-1]
                                         * and rfield.tNuRyd[rfield.ipspec][i] - do linear 
                                         * interpolation in log log space */
                                        y = rfield.tFluxLog[rfield.ipspec][i-1] + 
                                          rfield.tslop[rfield.ipspec][i-1]*
                                          log10(xnu/rfield.tNuRyd[rfield.ipspec][i-1]);

                                        /* return value is photon density, div by energy */
                                        ffun1_v = pow(10.,y);

                                        /* this checks that overshoots did not occur - interpolated
                                         * value must be between lowest and highest point */
#                                       ifndef NDEBUG
                                        double ys1 = MIN2( rfield.tFluxLog[rfield.ipspec][i-1],rfield.tFluxLog[rfield.ipspec][i]);
                                        double ys2 = MAX2( rfield.tFluxLog[rfield.ipspec][i-1],rfield.tFluxLog[rfield.ipspec][i]);
                                        ys1 = pow( 10. , ys1 );
                                        ys2 = pow( 10. , ys2 );
                                        ASSERT( ffun1_v >= ys1/(1.+100.*FLT_EPSILON) );
                                        ASSERT( ffun1_v <= ys2*(1.+100.*FLT_EPSILON) );
#                                       endif
                                        /* return value is photon density, div by energy */
                                        return( ffun1_v/xnu );
                                }
                                ++i;
                        }
                        /* energy above highest in table */
                        ffun1_v = 0.;
                }
                else
                {
                        /* energy below lowest on table */
                        ffun1_v = 0.;
                }
        }

        else if( strcmp(chKey,"BREMS") == 0 )
        {
                /* brems continuum, rough gaunt factor */
                fac = TE1RYD*xnu/rfield.slope[rfield.ipspec];
                ffun1_v = sexp(fac)/pow(xnu,1.2);

        }
        else if( strcmp(chKey,"LASER") == 0 )
        {
                const double BIG = 1.e10;
                const double SMALL = 1.e-25;
                /* a laser, mostly one frequency */
                /* >>chng 01 jul 01, was hard-wired 0.05 rel frac, change to optional
                 * second parameter, with default of 0.05 */
                /*if( xnu > 0.95*rfield.slope[rfield.ipspec] && xnu < 
                  1.05*rfield.slope[rfield.ipspec] )*/
                if( xnu > (1.-rfield.cutoff[rfield.ipspec][0])*rfield.slope[rfield.ipspec] && 
                        xnu < (1.+rfield.cutoff[rfield.ipspec][0])*rfield.slope[rfield.ipspec] )
                {
                        ffun1_v = BIG;
                }
                else
                {
                        ffun1_v = SMALL;
                }

        }
        else if( strcmp(chKey,"READ ") == 0 )
        {
                /* implement the table read command
                 * if this is first time called then we need to read in file */
                if( rfield.tslop[rfield.ipspec][0] == 0. )
                {
                        /* >>chng 01 nov 01, array index for below was bogus, only worked for a single continuum */
                        ReadTable(rfield.ioTableRead[rfield.ipspec]);
                        rfield.tslop[rfield.ipspec][0] = 1.;
                }
                /* use array of values read in on TABLE READ command */
                if( xnu >= rfield.egamry )
                {
                        ffun1_v = 0.;
                }
                else
                {
                        i = ipoint(xnu);
                        if( i > rfield.nupper || i < 1 )
                        {
                                ffun1_v = 0.;
                        }
                        else
                        {
                                ffun1_v = rfield.ConTabRead[i-1]/rfield.anu[i-1]/rfield.anu[i-1];
                        }
                }
        }

        /* >>chng 06 jul 10, retired TABLE STARBURST command, PvH */
        /* >>chng 05 nov 30, retired TABLE TLUSTY command, PvH */

        else if( strcmp(chKey,"VOLK ") == 0 )
        {
                /* use array of values read in from Kevin Volk's rebinning of
                 * large atlas grids */
                if( xnu >= rfield.egamry )
                {
                        ffun1_v = 0.;
                }
                else
                {
                        i = ipoint(xnu);
                        if( i > rfield.nupper )
                        {
                                fprintf( ioQQQ, " ffun1: Too many points - increase ncell\n" );
                                fprintf( ioQQQ, " cell needed=%4ld ncell=%4ld\n", 
                                  i, rfield.nupper );
                                cdEXIT(EXIT_FAILURE);
                        }
                        if( i > rfield.nupper || i < 1 )
                        {
                                ffun1_v = 0.;
                        }
                        else
                        {
                                /* bug fixed Jul 9 93: FFUN1 = TSLOP(IPSPEC,I) / ANU(I) / ANU(I)
                                 *   i has value 939 */
                                ffun1_v = rfield.tslop[rfield.ipspec][i-1]/ rfield.anu[i-1];
                        }
                }
        }
        else
        {
                fprintf( ioQQQ, " ffun1: I do not understand continuum label \"%s\" for continuum %li.\n", 
                  chKey , rfield.ipspec);
                cdEXIT(EXIT_FAILURE);
        }

        if( ffun1_v > 1e35 && !lgWarn )
        {
                lgWarn = true;
                fprintf( ioQQQ, " FFUN1:  Continuum %ld is very intense.\n", 
                  rfield.ipspec );
                fprintf( ioQQQ, " I will try to press on, but may have problems.\n" );
        }
        return( ffun1_v );
}

/*ReadTable called by TABLE READ to read in continuum from PUNCH TRANSMITTED CONTINUUM */
STATIC void ReadTable(const string& fnam)
{
        char chLine[INPUT_LINE_LENGTH];
        long int i, 
          nPoints;
        double Differ, 
          EnerLast;
        FILE *io;

        DEBUG_ENTRY( "ReadTable()" );

        io = open_data( fnam.c_str(), "r", AS_LOCAL_ONLY );

        /* read in first line of header */
        if( read_whole_line( chLine , (int)sizeof(chLine) , io )==NULL )
        {
                fprintf( ioQQQ, " error 1 reading input continuum file.\n" );
                cdEXIT(EXIT_FAILURE);
        }

        /* now read in the file of numbers */
        i = 0;
        /* keep reading until we hit eol or run out of room in the array */
        while( (read_whole_line(chLine, (int)sizeof(chLine),io)!=NULL) && (i<rfield.nupper) )
        {
                double help[2];
                sscanf( chLine, "%lf%lf ", &help[0], &help[1] );
                opac.tmn[i] = (realnum)help[0];
                rfield.ConTabRead[i] = (realnum)help[1];
                ++i;
        }
        /* put pointer at last good value */
        nPoints = i - 1;

        /* check that sane number of values entered */
        if( nPoints < 1 )
        {
                fprintf( ioQQQ, " ReadTable, file for TABLE READ has too few points, =%5ld\n", 
                  nPoints );
                cdEXIT(EXIT_FAILURE);
        }

        /* check on units of energy scale, convert to Rydbergs if necessary
         * tmn is scale read in from punch file, anu is correct scale
         * EnerLast is energy of last point in Rydbergs */
        EnerLast = opac.tmn[nPoints];
        if( fabs(opac.tmn[0]-rfield.anu[0])/rfield.anu[0] > 1e-3 )
        {
                /* first energy does not appear to have been in Rydbergs */
                if( opac.tmn[0] <= 0. )
                {
                        fprintf( ioQQQ, 
                                " DISASTER ReadTable:  the first energy in table read file is not positive.  Something is wrong.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
                else if( fabs(opac.tmn[0]-RYDLAM/rfield.anu[0]*1e-4)/opac.tmn[0] < 
                  1e-3 )
                {
                        /* wavelength in microns */
                        EnerLast = RYDLAM/opac.tmn[nPoints]/1e4;
                }
                else if( fabs(opac.tmn[0]-RYDLAM/rfield.anu[0])/opac.tmn[0] < 
                  1e-3 )
                {
                        /* wavelength in Angstroms */
                        EnerLast = RYDLAM/opac.tmn[nPoints];
                }
                else if( fabs(opac.tmn[0]-rfield.anu[0]*EVRYD*1e-3)/opac.tmn[0] < 
                  1e-3 )
                {
                        /* wavelength in keV */
                        EnerLast = opac.tmn[nPoints]/EVRYD/1e-3;
                }
                else if( fabs(opac.tmn[0]-rfield.anu[0]*EVRYD)/opac.tmn[0] < 
                  1e-3 )
                {
                        /* wavelength in eV */
                        EnerLast = opac.tmn[nPoints]/EVRYD;
                }
                else
                {
                        fprintf( ioQQQ, " DISASTER ReadTable:  the energy scale in the table read file makes no sense to me.\n" );
                        cdEXIT(EXIT_FAILURE);
                }
        }

        /* now check now the energies of the highest points agree */
        Differ = fabs(EnerLast-rfield.anu[nPoints])/rfield.anu[nPoints];
        if( Differ > 0.001 )
        {
                fprintf( ioQQQ, " DISASTER ReadTable: The energy mesh of the file read in by the TABLE READ command does not agree with this version of Cloudy.\n" );
                fprintf( ioQQQ, " DISASTER ReadTable: Was the file generated by an older version of the code?\n" );
                fprintf( ioQQQ, " DISASTER ReadTable: Did the first model do more than one iteration, but the LAST option was missing on PUNCH LAST TRANSMITTED CONTINUUM?\n");
                fprintf( ioQQQ, " DISASTER ReadTable: Number of points read in=%5ld\n", 
                  nPoints );
                fprintf( ioQQQ, " ReadTable: input, internal energies=%12.4e%12.4e\n", 
                  opac.tmn[nPoints], rfield.anu[nPoints] );
                cdEXIT(EXIT_FAILURE);
        }

        /* make sure rest of array has valid zeros */
        for( i=nPoints + 1; i < rfield.nupper; i++ )
        {
                rfield.ConTabRead[i] = 0.;
        }

        fclose(io);
        return;
}