/home66/gary/public_html/cloudy/c08_branch/source/cdspec.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 */
/*cdSPEC returns the spectrum needed for Keith Arnaud's XSPEC */
/*Spec_cont called by cdSPEC to generate actual spectrum */
#include "cddefines.h"
#include "cddrive.h"
#include "physconst.h"
#include "geometry.h"
#include "radius.h"
#include "rfield.h"
#include "opacity.h"
#include "grid.h"

/* 
 * this routine returns the spectrum needed for Keith Arnaud's XSPEC
 * X-Ray analysis code.  It should be called after cdDrive has successfully
 * computed a model.  the calling routine must ensure that the vectors
 * have enough space to store the resulting spectrum, 
 * given the bounds and energy resolution 
 */

/* array index within energy array in Cloudy grids -  this has file scope */
static long int iplo , iphi;

/* energies of lower and upper bound of XSPEC cell we are considering */
static double Elo , Ehi;

/* interpolate on cloudy's predicted spectrum to get results on XSPEC grid */
STATIC double Spec_cont( realnum cont[] );

void cdSPEC( 
        /* option - the type of spectrum to be returned
         * 1    the incident continuum 4\pi nuJ_nu, , erg cm-2 s-1
         *
         * 2    the attenuated incident continuum, same units
         * 3    the reflected continuum, same units
         *
         * 4    diffuse continuous emission outward direction
         * 5    diffuse continuous emission, reflected
         *
         * 6    collisional+recombination lines, outward
         * 7    collisional+recombination lines, reflected
         * 
         * 8    pumped lines, outward <= not implemented
         * 9    pumped lines, reflected <= not implemented
         *
         *              all lines and continuum emitted by the cloud assume full coverage of 
         *              continuum source */
        int nOption ,

        /* the energy of the lower edge of each cell 
         * (in Ryd to be consistent with all the rest of Cloudy) */
        double EnergyLow[] , 

        /* the number of cells + 1*/
        long int nEnergy ,

        /* the returned spectrum, same size is two energy spectra (see option), returns nEnergy -1 pts */
        double ReturnedSpectrum[] )

{
        /* this pointer will bet set to one of the cloudy continuum arrays */
        realnum *cont , 
                refac;
        long int ncell , j;

        /* flag saying whether we will need to free cont at the end */
        bool lgFREE;

        DEBUG_ENTRY( "cdSPEC()" );

        if( nOption == 1 )
        {
                /* this is the incident continuum, col 2 of punch continuum command */
                cont = rfield.flux_total_incident;
                lgFREE = false;
        }
        else if( nOption == 2 )
        {
                /* the attenuated transmitted continuum, no diffuse emission,
                 * col 3 of punch continuum command */
                cont = rfield.flux;
                lgFREE = false;
        }
        else if( nOption == 3 )
        {
                /* reflected incident continuum, col 6 of punch continuum command */
                lgFREE = false;
                cont = rfield.ConRefIncid;
        }
        else if( nOption == 4 )
        {
                /* diffuse continuous emission outward direction  */
                cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );

                /* need to free the vector once done */
                lgFREE = true;
                refac = (realnum)radius.r1r0sq*geometry.covgeo;
                for( j=0; j<rfield.nflux; ++j)
                {
                        cont[j] = rfield.ConEmitOut[j]*refac;
                }
        }
        else if( nOption == 5 )
        {
                /* reflected diffuse continuous emission */
                cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
                /* need to free the vector once done */
                lgFREE = true;
                refac = (realnum)radius.r1r0sq*geometry.covgeo;
                for( j=0; j<rfield.nflux; ++j)
                {
                        cont[j] = rfield.ConEmitReflec[j]*refac;
                }
        }
        else if( nOption == 6 )
        {
                /* all outward lines,   */
                cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
                /* need to free the vector once done */
                lgFREE = true;
                /* correct back to inner radius */
                refac = (realnum)radius.r1r0sq*geometry.covgeo;
                for( j=0; j<rfield.nflux; ++j)
                {
                        /* units of lines here are to cancel out with tricks applied to continuum cells
                         * when finally extracted below */
                        cont[j] = rfield.outlin[j] *rfield.widflx[j]/rfield.anu[j]*refac;
                }
        }
        else if( nOption == 7 )
        {
                /* all reflected lines */
                if( geometry.lgSphere )
                {
                        refac = 0.;
                }
                else
                {
                        refac = 1.;
                }

                cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
                /* need to free the vector once done */
                lgFREE = true;
                for( j=0; j<rfield.nflux; ++j)
                {
                        /* units of lines here are to cancel out with tricks applied to continuum cells
                         * when finally extracted below */
                        cont[j] = rfield.reflin[j] *rfield.widflx[j]/rfield.anu[j]*refac;
                }
        }
        else
        {
                fprintf(ioQQQ," cdSPEC called with impossible nOption (%i)\n", nOption);
                cdEXIT(EXIT_FAILURE);
        }

        /* this is array index used in Spec_cont */
        iplo = 0;
        iphi = 0;
        /* now generate the continua */
        for( ncell = 0; ncell < nEnergy-1; ++ncell )
        {
                /* Spec_cont will find the spectrum between these two energies */
                Elo = EnergyLow[ncell];
                Ehi = EnergyLow[ncell+1];
                ReturnedSpectrum[ncell] = Spec_cont( cont );
        }

        /* need to free the vector once done if this flag was set */
        if( lgFREE )
        {
                free(cont);
        }
        return;
}

/* interpolate on cloudy's predicted spectrum to get results on XSPEC grid */
STATIC double Spec_cont( realnum cont[] )
{
        double sum;
        long int i;

        DEBUG_ENTRY( "Spec_cont()" );

        /* iplo and iphi are set to zero before this is called to generate a full spectrum,
         * since energies are in increasing order, after the first call it should only
         * be necessary to increase ip above where it now is */
        /* return zero if continuum does not go this high */
        if( iplo >= rfield.nflux )
                return 0.;

        /* now skip out to where continuum starts, often we will already be there
         * if this is a successive call */
        /* iplo should be index where Elo is greater than anu-widflx/1. */
        while( (iplo < rfield.nflux)  &&  
                !(Elo <= (rfield.AnuOrg[iplo+1]-rfield.widflx[iplo+1]/2.) &&
                 (Elo >= (rfield.AnuOrg[iplo]-rfield.widflx[iplo]/2.)  )) )
                ++iplo;

        iphi = iplo;
        /* iphi should be index where Ehi is greater than anu-widflx/1. */
        while( (iphi < rfield.nflux)  &&  
                !(Ehi <= (rfield.AnuOrg[iphi+1]-rfield.widflx[iphi+1]/2.) &&
                 (Ehi >= (rfield.AnuOrg[iphi]-rfield.widflx[iphi]/2.)  )) )
                ++iphi;

        sum = 0.;
        if( iphi-iplo>1 )
        {
                /* this branch when more than two cells are involved - 
                 * begin by summing intermediate cells */
                for( i=iplo+1; i<iphi; ++i )
                {
                        sum += cont[i] * rfield.anu2[i];
                }
        }
        ASSERT( sum>=0.);

        /* at this point sum will often be zero, and we need to add on flux from
         * the iplo and iphi cells - there are two possibilities - iplo and iphi can
         * be the same cell, or different cells */
        if( iplo == iphi )
        {
                /* we want only a piece of the ip cell */
                sum = cont[iplo]/rfield.widflx[iplo]*(Ehi-Elo) * rfield.anu2[iplo];
                ASSERT( sum>=0.);
        }
        else
        {
                /* we want to add on a fraction of the low and high energy cells */
                double piece;

                /* the low energy cell */
                piece = cont[iplo]/rfield.widflx[iplo]*( (rfield.AnuOrg[iplo+1]-rfield.widflx[iplo+1]/2.) - Elo ) * 
                        rfield.anu2[iplo];
                ASSERT( piece >= 0.);
                sum += piece;

                /* now the high energy cell */
                piece = cont[iphi]/rfield.widflx[iphi]*(Ehi - (rfield.AnuOrg[iphi]-rfield.widflx[iphi]/2.) ) * 
                        rfield.anu2[iphi];
                ASSERT( piece >= 0.);
                sum += piece;
                ASSERT( sum>=0.);
        }

        /* now divide by total energy width */
        sum *= EN1RYD / (Ehi - Elo );

        ASSERT( sum >=0. );
        return sum;

#       if 0
        /* at this point Elo is less than the upper edge energy of the ip cell */
        /* >>chng 01 jul 23, change logic, use this branch if desired width smaller than
         * intrinsic cloudy width */
        /*if( rfield.anu[ip]+rfield.widflx[ip]/2. > Ehi )*/
        if( (Ehi-Elo) < rfield.widflx[ip] )
        {
                /* this is where we want the flux evaluated */
                double Emiddle = (Elo + Ehi)/2.;
                double f1 , f2 ,finterp;
                ASSERT( (Elo >= (rfield.anu[ip]+rfield.widflx[ip]/2.
                /* this is case where requested cell is within a single Cloudy cell -
                 * do linear interpolation */
                /* >>chng 01 jul 23, change to linear interpolation */
                /*return( cont[ip] /rfield.widflx[ip] *rfield.anu2[ip]*EN1RYD );*/
                f1 = cont[ip-1] /rfield.widflx[ip-1] *rfield.anu2[ip-1]*EN1RYD;
                f2 = cont[ip] /rfield.widflx[ip] *rfield.anu2[ip]*EN1RYD;

                finterp = f1 + (f2-f1) /rfield.widflx[ip] *
                        fabs(Emiddle-(rfield.anu[ip]-rfield.widflx[ip]/2.));
                return( finterp );
        }
        else
        {
                /* this is case where we will add on more than one Cloudy cell */
                sum = cont[ip] * rfield.anu2[ip]* EN1RYD;
                /* energy of upper edge of Cloudy cell */
                EcHi = ( rfield.anu[ip] + rfield.widflx[ip]/2. + rfield.anu[ip+1] - rfield.widflx[ip+1]/2.)/2.;
                /* mutiply this by energy width of this first cell that counts to final integral */
                sum *= (EcHi - Elo )/rfield.widflx[ip];
        }

        /* increment index since we will now work on next cell */
        ++ip;

        /* now add up total cells, this loop for cloudy cells totally within desired bins,
         * we will divide by energy width later, so factor of widflx is missing here */
        while( Ehi > rfield.anu[ip]+rfield.widflx[ip]/2. )
        {
                sum += cont[ip] * rfield.anu2[ip]*EN1RYD;
                ++ip;
        }

        /* energy of upper edge of Cloudy cell */
        EcLo = ( rfield.anu[ip] - rfield.widflx[ip]/2. + rfield.anu[ip-1] + rfield.widflx[ip-1]/2.)/2.;

        /* now finish up with last cell, whose upper bound is higher than desired bin */
        sum += cont[ip] * rfield.anu2[ip] * EN1RYD *
                (Ehi - EcLo )/rfield.widflx[ip];

        /* now divide by total energy width */
        sum /= (Ehi - Elo );

        return(sum);
#       endif
}

/* returns in units photons/cm^2/s/bin */
void cdSPEC2( 
        /* option - the type of spectrum to be returned
         * 1    the incident continuum 4\pi nuJ_nu, , erg cm-2 s-1
         *
         * 2    the attenuated incident continuum, same units
         * 3    the reflected continuum, same units
         *
         * 4    diffuse continuous emission outward direction
         * 5    diffuse continuous emission, reflected
         *
         * 6    collisional+recombination lines, outward
         * 7    collisional+recombination lines, reflected
         *
         * 8    total transmitted, lines and continuum
         * 9    total reflected, lines and continuum
         *
         *10    exp(-tau) to the illuminated face
         *
         *              all lines and continuum emitted by the cloud assume full coverage of 
         *              continuum source */
        int nOption ,

        /* the number of cells + 1*/
        long int nEnergy ,

        /* the returned spectrum, same size is two energy spectra (see option), returns nEnergy -1 pts */
        realnum ReturnedSpectrum[] )

{
        /* this pointer will bet set to one of the cloudy continuum arrays */
        realnum *cont , 
                refac;
        long int ncell , j;

        /* flag saying whether we will need to free cont at the end */
        bool lgFREE;

        DEBUG_ENTRY( "cdSPEC2()" );

        ASSERT( nOption <= NUM_OUTPUT_TYPES );

        if( nOption == 1 )
        {
                /* this is the incident continuum, col 2 of punch continuum command */
                cont = rfield.flux_total_incident;
                lgFREE = false;
        }
        else if( nOption == 2 )
        {
                /* the attenuated transmitted continuum, no diffuse emission,
                 * col 3 of punch continuum command */
                cont = rfield.flux;
                lgFREE = false;
        }
        else if( nOption == 3 )
        {
                /* reflected incident continuum, col 6 of punch continuum command */
                lgFREE = false;
                cont = rfield.ConRefIncid;
        }
        else if( nOption == 4 )
        {
                /* diffuse continuous emission outward direction  */
                cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
                /* need to free the vector once done */
                lgFREE = true;
                refac = (realnum)radius.r1r0sq*geometry.covgeo;
                for( j=0; j<rfield.nflux; ++j)
                {
                        cont[j] = rfield.ConEmitOut[j]*refac;
                }
        }
        else if( nOption == 5 )
        {
                /* reflected diffuse continuous emission */
                cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
                /* need to free the vector once done */
                lgFREE = true;
                for( j=0; j<rfield.nflux; ++j)
                {
                        cont[j] = rfield.ConEmitReflec[j];
                }
        }
        else if( nOption == 6 )
        {
                /* all outward lines,   */
                cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
                /* need to free the vector once done */
                lgFREE = true;
                /* correct back to inner radius */
                refac = (realnum)radius.r1r0sq*geometry.covgeo;
                for( j=0; j<rfield.nflux; ++j)
                {
                        cont[j] = rfield.outlin[j]*refac;
                }
        }
        else if( nOption == 7 )
        {
                /* all reflected lines */
                if( geometry.lgSphere )
                {
                        refac = 0.;
                }
                else
                {
                        refac = 1.;
                }

                cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
                /* need to free the vector once done */
                lgFREE = true;
                for( j=0; j<rfield.nflux; ++j)
                {
                        /* units of lines here are to cancel out with tricks applied to continuum cells
                         * when finally extracted below */
                        cont[j] = rfield.reflin[j]*refac;
                }
        }
        else if( nOption == 8 )
        {
                /* total transmitted continuum */
                cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
                /* need to free the vector once done */
                lgFREE = true;
                /* correct back to inner radius */
                refac = (realnum)radius.r1r0sq*geometry.covgeo;
                for( j=0; j<rfield.nflux; ++j)
                {
                        /* units of lines here are to cancel out with tricks applied to continuum cells
                         * when finally extracted below */
                        cont[j] = (rfield.ConEmitOut[j]+ rfield.outlin[j])*refac
                                          + rfield.flux[j]*(realnum)radius.r1r0sq;
                }
        }
        else if( nOption == 9 )
        {
                /* total reflected continuum */
                cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
                /* need to free the vector once done */
                lgFREE = true;
                for( j=0; j<rfield.nflux; ++j)
                {
                        /* units of lines here are to cancel out with tricks applied to continuum cells
                         * when finally extracted below */
                        cont[j] = ((rfield.ConRefIncid[j]+rfield.ConEmitReflec[j]) +
                                rfield.reflin[j]);
                }
        }
        else if( nOption == 10 )
        {
                /* this is exp(-tau) */
                cont = (realnum*)MALLOC( sizeof(realnum)*(size_t)rfield.nupper );
                /* need to free the vector once done */
                lgFREE = true;
                for( j=0; j<nEnergy; ++j)
                {
                        /* This is the TOTAL attenuation in both the continuum and lines.  
                         * Jon Miller discovered that the line attenuation was missing in 07.02 */
                        cont[j] = opac.ExpmTau[j]*rfield.trans_coef_total[j];
                }
        }
        else
        {
                fprintf(ioQQQ," cdSPEC called with impossible nOption (%i)\n", nOption);
                cdEXIT(EXIT_FAILURE);
        }

        /* now generate the continua */
        for( ncell = 0; ncell < nEnergy; ++ncell )
        {
        if( ncell >= rfield.nflux )
                {
                        ReturnedSpectrum[ncell] = SMALLFLOAT;
                }
                else
                {
                        ReturnedSpectrum[ncell] = cont[ncell];
                }
                ASSERT( ReturnedSpectrum[ncell] >=0.f );
        }

        /* need to free the vector once done if this flag was set */
        if( lgFREE )
        {
                free(cont);
        }
        return;
}

---