cdspec.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2017 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 */
#include "cddefines.h"
#include "cddrive.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 
 */

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 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;
        vector<realnum> cont_local;

        DEBUG_ENTRY( "cdSPEC()" );

        ASSERT( nEnergy <= rfield.nflux );

        if( nOption == 1 )
        {
                /* this is the incident continuum, col 2 of save continuum command */
                cont = rfield.flux_total_incident[0];
        }
        else if( nOption == 2 )
        {
                /* the attenuated transmitted continuum, no diffuse emission,
                 * col 3 of save continuum command */
                cont = rfield.flux[0];
        }
        else if( nOption == 3 )
        {
                /* reflected incident continuum, col 6 of save continuum command */
                cont = rfield.ConRefIncid[0];
        }
        else if( nOption == 4 )
        {
                /* diffuse continuous emission outward direction  */
                cont_local.resize(rfield.nflux_with_check);
                cont = &cont_local[0];

                /* need to free the vector once done */
                refac = (realnum)radius.r1r0sq*geometry.covgeo;
                for( j=0; j<rfield.nflux; ++j)
                {
                        cont[j] = rfield.ConEmitOut[0][j]*refac;
                }
        }
        else if( nOption == 5 )
        {
                cont_local.resize(rfield.nflux_with_check);
                /* reflected diffuse continuous emission */
                cont = &cont_local[0];
                /* need to free the vector once done */
                refac = (realnum)radius.r1r0sq*geometry.covgeo;
                for( j=0; j<rfield.nflux; ++j)
                {
                        cont[j] = rfield.ConEmitReflec[0][j]*refac;
                }
        }
        else if( nOption == 6 )
        {
                cont_local.resize(rfield.nflux_with_check);
                /* all outward lines,   */
                cont = &cont_local[0];
                /* 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[0][j] *rfield.widflx(j)/rfield.anu(j)*refac;
                }
        }
        else if( nOption == 7 )
        {
                /* all reflected lines */
                if( geometry.lgSphere )
                {
                        refac = 0.;
                }
                else
                {
                        refac = 1.;
                }

                cont_local.resize(rfield.nflux_with_check);
                cont = &cont_local[0];
                /* need to free the vector once done */
                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[0][j] *rfield.widflx(j)/rfield.anu(j)*refac;
                }
        }
        else
        {
                fprintf(ioQQQ," cdSPEC called with impossible nOption (%i)\n", nOption);
                cdEXIT(EXIT_FAILURE);
        }

        /* now generate the continua */
        for( ncell = 0; ncell < nEnergy-1; ++ncell )
        {
                ReturnedSpectrum[ncell] = cont[ncell] * EN1RYD * rfield.anu2(ncell) / rfield.widflx(ncell);
        }

        return;
}


/* 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 emission, lines + continuum, outward
         * 5    diffuse emission, lines + continuum, reflected
         *
         * 6    diffuse continuous emission outward direction
         * 7    diffuse continuous emission, 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 */
        long int nEnergy,

        /* the index of the lowest and highest energy bounds to use. */
        long ipLoEnergy,
        long ipHiEnergy,

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

{
        realnum refac;

        DEBUG_ENTRY( "cdSPEC2()" );

        if( ! ( ipLoEnergy >= 0 && ipLoEnergy < ipHiEnergy && ipHiEnergy < rfield.nflux_with_check &&
                nEnergy == (ipHiEnergy-ipLoEnergy+1) && nEnergy >= 2 && nOption <= NUM_OUTPUT_TYPES ) )
        {
                fprintf( ioQQQ, "invalid parameter for cdSPEC2\n" );
                cdEXIT(EXIT_FAILURE);
        }

        const realnum *trans_coef_total = rfield.getCoarseTransCoef();

        for( long i = 0; i < nEnergy; i++ )
        {
                long j = ipLoEnergy + i;

                if( j >= rfield.nflux )
                {
                        ReturnedSpectrum[i] = SMALLFLOAT;
                        continue;
                }
                
                if( nOption == 0 )
                {
                        /* the attenuated incident continuum */
                        realnum flxatt = rfield.flux[0][j]*
                                (realnum)radius.r1r0sq * trans_coef_total[j];

                        /* the outward emitted continuum */
                        realnum conem = (rfield.ConEmitOut[0][j] + rfield.outlin[0][j])*
                                (realnum)radius.r1r0sq*geometry.covgeo;

                        /* the reflected continuum */
                        realnum flxref = rfield.ConRefIncid[0][j] + rfield.ConEmitReflec[0][j] +
                                rfield.reflin[0][j];

                        ReturnedSpectrum[i] = flxatt + conem + flxref;
                }
                else if( nOption == 1 )
                {
                        /* this is the incident continuum, col 2 of save continuum command */
                        ReturnedSpectrum[i] = rfield.flux_total_incident[0][j];
                }
                else if( nOption == 2 )
                {
                        /* the attenuated transmitted continuum, no diffuse emission,
                         * col 3 of save continuum command */
                        ReturnedSpectrum[i] = rfield.flux[0][j]*
                                (realnum)radius.r1r0sq * trans_coef_total[j];
                }
                else if( nOption == 3 )
                {
                        /* reflected incident continuum, col 6 of save continuum command */
                        ReturnedSpectrum[i] = rfield.ConRefIncid[0][j];
                }
                else if( nOption == 4 )
                {
                        /* all outward diffuse emission */
                        /* correct back to inner radius */
                        refac = (realnum)radius.r1r0sq*geometry.covgeo;
                        ReturnedSpectrum[i] = (rfield.outlin[0][j]+rfield.ConEmitOut[0][j])*refac;
                }
                else if( nOption == 5 )
                {
                        /* all reflected diffuse emission */
                        if( geometry.lgSphere )
                        {
                                refac = 0.;
                        }
                        else
                        {
                                refac = 1.;
                        }

                        ReturnedSpectrum[i] = (rfield.reflin[0][j]+rfield.ConEmitReflec[0][j])*refac;
                }
                else if( nOption == 6 )
                {
                        /* all outward line emission */
                        /* correct back to inner radius */
                        refac = (realnum)radius.r1r0sq*geometry.covgeo;
                        ReturnedSpectrum[i] = rfield.outlin[0][j]*refac;
                }
                else if( nOption == 7 )
                {
                        /* all reflected line emission */
                        if( geometry.lgSphere )
                        {
                                refac = 0.;
                        }
                        else
                        {
                                refac = 1.;
                        }

                        ReturnedSpectrum[i] = rfield.reflin[0][j]*refac;
                }
                else if( nOption == 8 )
                {
                        /* total transmitted continuum */
                        /* correct back to inner radius */
                        refac = (realnum)radius.r1r0sq*geometry.covgeo;
                        ReturnedSpectrum[i] = (rfield.ConEmitOut[0][j]+ rfield.outlin[0][j])*refac
                                + rfield.flux[0][j]*(realnum)radius.r1r0sq*trans_coef_total[j];
                }
                else if( nOption == 9 )
                {
                        /* total reflected continuum */
                        ReturnedSpectrum[i] = rfield.ConRefIncid[0][j] + rfield.ConEmitReflec[0][j] +
                                rfield.reflin[0][j];
                }
                else if( nOption == 10 )
                {
                        /* this is exp(-tau) */
                        /* This is the TOTAL attenuation in both the continuum and lines.  
                         * Jon Miller discovered that the line attenuation was missing in 07.02 */
                        ReturnedSpectrum[i] = opac.ExpmTau[j]*trans_coef_total[j];
                }
                else
                {
                        fprintf(ioQQQ," cdSPEC called with impossible nOption (%i)\n", nOption);
                        cdEXIT(EXIT_FAILURE);
                }

                ASSERT( ReturnedSpectrum[i] >=0.f );
        }

        return;
}