/home66/gary/public_html/cloudy/c08_branch/source/prt_lines_continuum.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 */
/*lines_continuum put energetics, H, and He lines into line intensity stack */
#include "cddefines.h"
#include "taulines.h"
#include "physconst.h"
#include "iso.h"
#include "geometry.h"
#include "dense.h"
#include "prt.h"
#include "opacity.h"
#include "coolheavy.h"
#include "phycon.h"
#include "rfield.h"
#include "predcont.h"
#include "lines_service.h"
#include "radius.h"
#include "continuum.h"
#include "lines.h"

void lines_continuum(void)
{

        double f1, 
                f2 , 
                bac , 
                flow;
        long i,nBand;

        DEBUG_ENTRY( "lines_continuum()" );

        /* code has all local emissivities zeroed out with cryptic comment about being
         * situation dependent.  Why?  this is option to turn back on */
        const bool KILL_CONT = false;

        i = StuffComment( "continua" );
        linadd( 0., (realnum)i , "####", 'i',
                " start continua");

        /***********************************************************************
         * stuff in Bac ratio - continuum above the Balmer Jump 
         * this is trick, zeroing out saved continuum integrated so far,
         * and adding the current version, so that the line array gives the 
         * value in the final continuum 
         *
         * reflected continuum is different from others since relative to inner
         * radius, others for for this radius 
         *************************************************************************/

        /***************************************************************************
         * "Bac " , 3646,  this is residual continuum at peak of Balmer Jump
         * flux below - flux above                                   
         ***************************************************************************/
        /* >>chng 00 dec 02, remove opac.tmn */
        /* >>chng 00 dec 19, remove / radius.GeoDil */
        /* extrapolated continuum above head */
        /* >>chng 01 jul 13, from ConInterOut to ConEmitOut */
        f1 = (rfield.ConEmitOut[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1] + 
                rfield.ConEmitReflec[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1]/radius.r1r0sq )/
          rfield.widflx[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1]/*/radius.GeoDil /opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1]*/;

        /* extrapolated continuum below head */
        /* >>chng 00 dec 19, remove / radius.GeoDil */
        f2 = (rfield.ConEmitOut[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2]+ 
                rfield.ConEmitReflec[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2]/radius.r1r0sq )/
          rfield.widflx[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2]/*/radius.GeoDil /opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2]*/;

        /* convert to nuFnu units */
        f1 = f1*0.250*0.250*EN1RYD*radius.r1r0sq;
        f2 = f2*0.250*0.250*EN1RYD*radius.r1r0sq;
        bac = (f1 - f2);

        /* memory not allocated until ipass >= 0 
         * clear summed intrinsic and emergent intensity of following
         * entry - following call to linadd will enter the total and
         * keep entering the total but is done for each zone hence need to
         * keep resetting to zero*/
        if( LineSave.ipass > 0 )
        {
                LineSv[LineSave.nsum].sumlin[0] = 0.;
                LineSv[LineSave.nsum].sumlin[1] = 0.;
        }

        linadd(MAX2(0.,bac)/radius.dVeff,3646,"Bac ",'i',
                "residual flux at head of Balmer continuum, nuFnu ");
        /* >>chng 03 feb 06, set to zero */
        /* emslin saves the per unit vol emissivity of a line, which is normally 
         * what goes into linadd.  We zero this unit emissivity which was set
         * FOR THE PREVIOUS LINE since it is so situation dependent */
        if( KILL_CONT && LineSave.ipass > 0 )
        {
                LineSv[LineSave.nsum-1].emslin[0] = 0.;
                LineSv[LineSave.nsum-1].emslin[1] = 0.;
        }

        /* memory not allocated until ipass >= 0 */
        if( LineSave.ipass > 0 )
        {
                LineSv[LineSave.nsum].sumlin[0] = 0.;
                LineSv[LineSave.nsum].sumlin[1] = 0.;
        }

        linadd(f1/radius.dVeff,3645,"nFnu",'i',
                "total flux above head of Balmer continuum, nuFnu ");
        /* >>chng 03 feb 06, set to zero */
        /* emslin saves the per unit vol emissivity of a line, which is normally 
         * what goes into linadd.  We zero this unit emissivity which was set
         * FOR THE PREVIOUS LINE since it is so situation dependent */
        if( KILL_CONT && LineSave.ipass > 0 )
        {
                LineSv[LineSave.nsum-1].emslin[0] = 0.;
                LineSv[LineSave.nsum-1].emslin[1] = 0.;
        }

        /* memory not allocated until ipass >= 0 */
        if( LineSave.ipass > 0 )
        {
                LineSv[LineSave.nsum].sumlin[0] = 0.;
                LineSv[LineSave.nsum].sumlin[1] = 0.;
        }

        linadd(f2/radius.dVeff,3647,"nFnu",'i',
                "total flux above head of Balmer continuum, nuFnu ");
        /* >>chng 03 feb 06, set to zero */
        /* emslin saves the per unit vol emissivity of a line, which is normally 
         * what goes into linadd.  We zero this unit emissivity which was set
         * FOR THE PREVIOUS LINE since it is so situation dependent */
        if( KILL_CONT && LineSave.ipass > 0 )
        {
                LineSv[LineSave.nsum-1].emslin[0] = 0.;
                LineSv[LineSave.nsum-1].emslin[1] = 0.;
        }

        /******************************************************************************
         * "cout" , 3646,  this is outward residual continuum at peak of Balmer Jump  *
         * equal to total in spherical geometry, half in opt thin open geometry       *
         ******************************************************************************/
        /* >>chng 00 dec 02, remove opac.tmn */
        /* >>chng 00 dec 19, remove / radius.GeoDil */
        f1 = rfield.ConEmitOut[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1]/
                rfield.widflx[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1]/*/radius.GeoDil /opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1]*/;

        /* >>chng 00 dec 19, remove / radius.GeoDil */
        f2 = rfield.ConEmitOut[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2]/
                rfield.widflx[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2]/*/radius.GeoDil /opac.tmn[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2]*/;

        /* net Balmer jump */
        bac = (f1 - f2)*0.250*0.250*EN1RYD*radius.r1r0sq;

        /* memory not allocated until ipass >= 0 */
        if( LineSave.ipass > 0 )
        {
                LineSv[LineSave.nsum].sumlin[0] = 0.;
                LineSv[LineSave.nsum].sumlin[1] = 0.;
        }

        linadd(MAX2(0.,bac)/radius.dVeff,3646,"cout",'i',
                "residual flux in Balmer continuum, nuFnu ");
        /* >>chng 03 feb 06, set to zero */
        /* emslin saves the per unit vol emissivity of a line, which is normally 
         * what goes into linadd.  We zero this unit emissivity which was set
         * FOR THE PREVIOUS LINE since it is so situation dependent */
        if( KILL_CONT && LineSave.ipass > 0 )
        {
                LineSv[LineSave.nsum-1].emslin[0] = 0.;
                LineSv[LineSave.nsum-1].emslin[1] = 0.;
        }

        /*********************************************************************
         * "cref" , 3646,  this is reflected continuum at peak of Balmer Jump*
         * equal to zero in spherical geometry, half of total in op thin opn *
         *********************************************************************/
        /* >>chng 00 dec 02, remove opac.tmn */
        /* >>chng 00 dec 19, remove / radius.GeoDil */
        f1 = rfield.ConEmitReflec[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1]/
                rfield.widflx[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1];

        f2 = rfield.ConEmitReflec[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2]/
                rfield.widflx[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2];

        /* net Balmer jump */
        bac = (f1 - f2)*0.250*0.250*EN1RYD;

        /* memory not allocated until ipass >= 0 */
        if( LineSave.ipass > 0 )
        {
                LineSv[LineSave.nsum].sumlin[0] = 0.;
                LineSv[LineSave.nsum].sumlin[1] = 0.;
        }

        linadd(MAX2(0.,bac)/radius.dVeff,3646,"cref",'i',
                "residual flux in Balmer continuum, nuFnu ");
        /* >>chng 03 feb 06, set to zero */
        /* emslin saves the per unit vol emissivity of a line, which is normally 
         * what goes into linadd.  We zero this unit emissivity which was set
         * FOR THE PREVIOUS LINE since it is so situation dependent */
        if( KILL_CONT && LineSave.ipass > 0 )
        {
                LineSv[LineSave.nsum-1].emslin[0] = 0.;
                LineSv[LineSave.nsum-1].emslin[1] = 0.;
        }

        /*********************************************************************
         * "thin" , 3646, tot optically thin continuum at peak of Balmer Jump*/
        if( nzone > 0 )
        {
                /* rfield.ConEmitLocal is not defined initially, only evaluate when into model */
                f1 = rfield.ConEmitLocal[nzone][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1]/
                        rfield.widflx[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1];
                f2 = rfield.ConEmitLocal[nzone][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2]/
                        rfield.widflx[iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-2];
                bac = (f1 - f2)*0.250*0.250*EN1RYD;
        }
        else
        {
                f1 = 0.;
                f2 = 0.;
                bac = 0.;
        }

        linadd(MAX2(0.,bac),3646,"thin",'i',
                "residual flux in Balmer continuum, nuFnu ");

        continuum.cn4861 /= (realnum)radius.dVeff;
        linadd(continuum.cn4861,4860,"Inci",'i',
                "incident continuum nu*f_nu at H-beta, at illuminated face of cloud ");
        /* >>chng 07 jun 13, at this point nsum is the number of lines in the stack
         * so nsum-1 is the entry created by the above call to linadd
         * we want to set the emergent incident continuum to the intrinsic
         * incident continuum - there are reports of the incident continuum 
         * striking the cloud and the emergent - incident distinction does
         * not apply - code had left emergent intensity at zero, bug reported
         * by K Korista on discussion group 2007 jun 14 */
        if( LineSave.ipass > 0 )
        {
                LineSv[LineSave.nsum-1].sumlin[1] = LineSv[LineSave.nsum-1].sumlin[0];
        }

        continuum.cn1216 /= (realnum)radius.dVeff;
        linadd(continuum.cn1216,1215,"Inci",'i',
                "incident continuum nu*f_nu near Ly-alpha, at illuminated face of cloud");
        /* see comment concerning parallel code immediately above */
        if( LineSave.ipass > 0 )
        {
                LineSv[LineSave.nsum-1].sumlin[1] = LineSv[LineSave.nsum-1].sumlin[0];
        }

        continuum.cn1216 *= (realnum)radius.dVeff;
        continuum.cn4861 *= (realnum)radius.dVeff;

        if( LineSave.ipass > 0 )
        {
                continuum.cn4861 = 0.;
                continuum.cn1216 = 0.;
        }

        flow = (iso.RadRecomb[ipH_LIKE][ipHYDROGEN][ipH2p][ipRecRad] + 
                iso.RadRecomb[ipH_LIKE][ipHYDROGEN][ipH2s][ipRecRad])*
          iso.RadRecomb[ipH_LIKE][ipHYDROGEN][ipH2p][ipRecEsc]*
          dense.eden*dense.xIonDense[ipHYDROGEN][1]* 5.45e-12;
        linadd(flow,0,"Ba C",'i',
                "integrated Balmer continuum emission");

        if( iso.n_HighestResolved_max[ipH_LIKE][ipHYDROGEN] >= 3 )
        {
                flow = ( iso.RadRecomb[ipH_LIKE][ipHYDROGEN][ipH3s][ipRecRad]*
                        iso.RadRecomb[ipH_LIKE][ipHYDROGEN][ipH3s][ipRecEsc] +
                        iso.RadRecomb[ipH_LIKE][ipHYDROGEN][ipH3p][ipRecRad]*
                        iso.RadRecomb[ipH_LIKE][ipHYDROGEN][ipH3p][ipRecEsc] +
                        iso.RadRecomb[ipH_LIKE][ipHYDROGEN][ipH3d][ipRecRad]*
                        iso.RadRecomb[ipH_LIKE][ipHYDROGEN][ipH3d][ipRecEsc] ) * 
                        dense.eden*dense.xIonDense[ipHYDROGEN][1]*3.53e-12;
        }
        else
        {
                flow = iso.RadRecomb[ipH_LIKE][ipHYDROGEN][3][ipRecRad]*
                        iso.RadRecomb[ipH_LIKE][ipHYDROGEN][3][ipRecEsc]*
                  dense.eden*dense.xIonDense[ipHYDROGEN][1]*3.53e-12;
        }
        linadd(flow,0,"PA C",'i',
                "Paschen continuum emission ");

        /* >>chng 05 oct 16, add this output option, these are a
         * series of continuum bands defined in the file bands_continuum.dat
         * this makes it possible to enter any integrated total emission
         * into the emission-line stack */
        for( nBand=0; nBand<continuum.nContBand; ++nBand )
        {
                double total = 0.;
                /* find total emission over band - units will be erg cm-2 s-1 */
                for( i=continuum.ipContBandLow[nBand]; i<continuum.ipContBandHi[nBand]; ++i )
                {
                        double xIntenOut = 
                                /* the attenuated incident continuum */
                                rfield.flux[i-1] +

                                /* the outward emitted continuum radiation field 
                                 * >>chng 06 aug 07, add geomery.covgeo here */
                                (rfield.ConEmitOut[i-1] + 

                                /* outward emitted lines */
                                rfield.outlin[i-1])*geometry.covgeo;
                        /* div by opac.E2TauAbsFace[i] because ConEmitReflec has already
                         * been corrected for this - emergent_line will introduce a 
                         * further correction so this will cancel the extra factor */
                        double xIntenIn = 
                                ((double)rfield.ConEmitReflec[i-1]/SDIV((double)opac.E2TauAbsFace[i-1]) + 
                                 /* outward emitted lines */
                                 rfield.reflin[i-1] )*geometry.covgeo;

                        /* the fraction of this that gets out */
                        total += rfield.anu[i-1] * 
                                emergent_line( xIntenIn , xIntenOut , i )
                                /* >>chng 06 aug 07, add opac.tmn here */
                                / SDIV(opac.tmn[i-1]);
                }
                /* we will call lindst with an energy half way between the two
                 * ends of the band.  This will make an extinction correction that
                 * has already been applied above by multiplying by emergent_line.
                 * Find this factor and remove it before the call */
                double corr = emergent_line( 0.5 , 0.5 , 
                        (continuum.ipContBandLow[nBand]+continuum.ipContBandHi[nBand])/2 );
                if( corr < SMALLFLOAT )
                        total = 0.;
                else
                        total /= corr;

                /* convert to physical units */
                total *= EN1RYD*radius.r1r0sq/radius.dVeff;
                /* memory not allocated until ipass >= 0 */
                if( LineSave.ipass > 0 )
                {
                        LineSv[LineSave.nsum].sumlin[0] = 0.;
                        LineSv[LineSave.nsum].sumlin[1] = 0.;
                }

                lindst( total , 
                        // the negative wavelength is a sentinel that the wavelength
                        // may be bogus - very often the "center" of the band, as defined
                        // by observers, is far to the blue end of the range.  use the
                        // observed definition of the wavelength.  This introduces an error
                        // since the wavelength is used to determine the transfer of the
                        // band against background opacities
                        -continuum.ContBandWavelength[nBand] , 
                        continuum.chContBandLabels[nBand],
                        (continuum.ipContBandLow[nBand]+continuum.ipContBandHi[nBand])/2 ,
                        'i' , false,
                        "continuum bands defined in bands_continuum.dat");
        }

        linadd(MAX2(0.,CoolHeavy.brems_cool_net),0,"HFFc",'c',
                "net free-free cooling, ALL species, free-free heating subtracted, so nearly cancels with cooling in LTE ");

        linadd(MAX2(0.,-CoolHeavy.brems_cool_net),0,"HFFh",'h',
                "net free-free heating, nearly cancels with cooling in LTE ");

        linadd(CoolHeavy.brems_cool_h,0,"H FF",'i',
                " H brems (free-free) cooling ");

        linadd(CoolHeavy.brems_heat_total,0,"FF H",'i',
                "total free-free heating ");

        linadd(CoolHeavy.brems_cool_he,0,"HeFF",'i',
                "He brems emission ");

        linadd(CoolHeavy.herec,0,"HeFB",'c',
                "He recombination cooling ");

        linadd(CoolHeavy.heavfb,0,"MeFB",'c',
                "heavy element recombination cooling ");

        linadd(CoolHeavy.brems_cool_metals,0,"MeFF",'i',
                "heavy elements (metals) brems cooling, heat not subtracted ");

        linadd(CoolHeavy.brems_cool_h+CoolHeavy.brems_cool_he+CoolHeavy.brems_cool_metals,0,"ToFF",'i',
                "total brems emission - total cooling but not minus heating ");

        linadd((CoolHeavy.brems_cool_h+CoolHeavy.brems_cool_he)*sexp(5.8e6/phycon.te),0,"FF X",'i',
                "part of H brems, in x-ray beyond 0.5KeV ");

        linadd(CoolHeavy.eebrm,0,"eeff",'c',
                "electron - electron brems ");

        /* predict emitted continuum at series of continuum points */
        /* variables are located in predcont.h, 
         * NPREDCONT number of continuum points,
         * next two arrays are defined in zerologic.c
         * EnrPredCont - energy (Ryd) where we want to predict the continuum,
         *               these are set in zerologic.c,
         * ipPreCont -   pointers to energy within continuum array, 
         *               these are set in ContCreatePointers
         *
         * the array of continuum energies where this will be done is
         * stored as EnrPredCont which lives in predcont.h
         * and is initialized in zerologic.c
         *
         * the entry nFnu will only be printed if the command
         * print diffuse continuum
         * is entered - 
         *
         * this code should be kept parallel with that in dopunch, where
         * punch continuum is produced, since two must agree */
        for( i=0; i < NPREDCONT; i++ )
        {
                double SourceTransmitted , Cont_nInu;
                double SourceReflected, 
                  DiffuseOutward,
                  DiffuseInward;
                double renorm;

                /* put wavelength in Angstroms into dummy structure, so that we can use iWavLen
                 * to get a proper wavelength with units, continuum energies are stored in EnrPredCont */
                TauDummy.WLAng = (realnum)(RYDLAM/EnrPredCont[i]);
                /*lambda = iWavLen(&TauDummy , &chUnits , &chShift );*/

                /* >>chng 00 dec 02, there were three occurrences of /opac.tmn which had the
                 * effect of raising the summed continuum by the local opacity correction factor.
                 * in the case of the Lyman continuum this raised the reported value by orders
                 * of magnitude.  There have been commented out in the following for now. */
                /* reflected total continuum (diff+incident emitted inward direction) */

                /* >>chng 00 dec 08, implement the "set nFnu [SOURCE_REFLECTED] ... command, PvH */
                /* >>chng 00 dec 19, remove / radius.GeoDil */
                renorm = rfield.anu2[ipPredCont[i]]*EN1RYD/rfield.widflx[ipPredCont[i]]/*/radius.GeoDil*/;

                /* this is the reflected diffuse continuum */
                if( prt.lgDiffuseInward )
                {
                        DiffuseInward = rfield.ConEmitReflec[ipPredCont[i]]*renorm;
                        /* /opac.tmn[ipPredCont[i]]*/
                }
                else
                {
                        DiffuseInward = 0.;
                }

                /* the outward diffuse continuum */
                if( prt.lgDiffuseOutward )
                {
                        DiffuseOutward = rfield.ConEmitOut[ipPredCont[i]]*renorm*radius.r1r0sq;
                }
                else
                {
                        DiffuseOutward = 0.;
                }

                /* reflected part of INCIDENT continuum (only incident, not diffuse, which was above) */
                if( prt.lgSourceReflected )
                {
                        SourceReflected =  rfield.ConRefIncid[ipPredCont[i]]*renorm;
                        /* /opac.tmn[ipPredCont[i]]*/
                }
                else
                {
                        SourceReflected =  0.;
                }

                /* the attenuated incident continuum */
                if( prt.lgSourceTransmitted )
                {
                        SourceTransmitted = rfield.flux[ipPredCont[i]]*renorm*radius.r1r0sq;
                }
                else
                {
                        SourceTransmitted = 0.;
                }

                /* memory has not been allocated until ipass >= 0, so must not access this element,
                 * this element will be used to save the following quantity */
                if( LineSave.ipass > 0 )
                {
                        LineSv[LineSave.nsum].sumlin[0] = 0.;
                        LineSv[LineSave.nsum].sumlin[1] = 0.;
                }

                linadd((DiffuseInward+SourceReflected+DiffuseOutward+SourceTransmitted)/radius.dVeff,
                        TauDummy.WLAng,"nFnu",'i',
                        "total continuum at selected energy points " );

                /* emslin saves the per unit vol emissivity of a line, which is normally 
                 * what goes into linadd.  We zero this unit emissivity which was set
                 * FOR THE PREVIOUS LINE since it is so situation dependent */
                if( KILL_CONT && LineSave.ipass > 0 )
                {
                        LineSv[LineSave.nsum-1].emslin[0] = 0.;
                        LineSv[LineSave.nsum-1].emslin[1] = 0.;
                }

                /* this is the normal set to zero to trick the NEXT line into going in properly */
                if( LineSave.ipass > 0 )
                {
                        LineSv[LineSave.nsum].sumlin[0] = 0.;
                        LineSv[LineSave.nsum].sumlin[1] = 0.;
                }

                /* the nsum-1 -- emslin and nsum -- sumlin is not a bug, look above - they do
                 * different things to different saves */
                Cont_nInu = rfield.flux[ipPredCont[i]]*renorm*radius.r1r0sq +
                        rfield.ConRefIncid[ipPredCont[i]]*renorm;

#               if 0
                /* this code can be used to create assert statements for the continuum shape */
                if( !i )
                        fprintf(ioQQQ,"\n");
                char chWL[1000];
                sprt_wl( chWL , TauDummy.WLAng );
                fprintf( ioQQQ,"assert line luminosity \"nInu\" %s  %.3f\n",
                        chWL , 
                        log10(SDIV(Cont_nInu/radius.dVeff)) + radius.Conv2PrtInten  );
#               endif

                linadd( Cont_nInu/radius.dVeff,TauDummy.WLAng,"nInu",'i',
                        "transmitted and reflected incident continuum at selected energy points " );

                /* emslin saves the per unit volume emissivity of a line, which is normally 
                 * what goes into linadd.  We zero this unit emissivity since it is so situation dependent */
                if( KILL_CONT && LineSave.ipass > 0 )
                {
                        LineSv[LineSave.nsum-1].emslin[0] = 0.;
                        LineSv[LineSave.nsum-1].emslin[1] = 0.;
                }

                /* memory has not been allocated until ipass >= 0 */
                if( LineSave.ipass > 0 )
                {
                        LineSv[LineSave.nsum].sumlin[0] = 0.;
                        LineSv[LineSave.nsum].sumlin[1] = 0.;
                }

                linadd( (DiffuseInward+SourceReflected)/radius.dVeff,TauDummy.WLAng,"InwT",'i',
                        "total reflected continuum, total inward emission plus reflected (XXdiffuseXX) total continuum ");

                if( KILL_CONT && LineSave.ipass > 0 )
                {
                        LineSv[LineSave.nsum-1].emslin[0] = 0.;
                        LineSv[LineSave.nsum-1].emslin[1] = 0.;
                }

                /* memory has not been allocated until ipass >= 0 */
                if( LineSave.ipass > 0 )
                {
                        LineSv[LineSave.nsum].sumlin[0] = 0.;
                        LineSv[LineSave.nsum].sumlin[1] = 0.;
                }

                linadd(SourceReflected/radius.dVeff,TauDummy.WLAng,"InwC",'i',
                        "reflected incident continuum (only incident) ");

                if( KILL_CONT && LineSave.ipass > 0 )
                {
                        LineSv[LineSave.nsum-1].emslin[0] = 0.;
                        LineSv[LineSave.nsum-1].emslin[1] = 0.;
                }
        }
        return;
}

---