rfield.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 */
#include "cddefines.h"
#include "rfield.h"

#include "opacity.h"
#include "continuum.h"

t_rfield rfield;

void t_rfield::zero()
{
        DEBUG_ENTRY( "t_rfield::zero()" );
        lgHabing = false;

        /* flag to turn off Lya ots */
        lgLyaOTS = true;
        /* HeII rec and Lya ots */
        lgHeIIOTS = true;
        lgKillOTSLine = false;
        lgKillOutLine = false;
        lgKillOutCont = false;

        /* DiffPumpOn is unity unless process disabled by setting to 1
         * with no diffuse line pumping command */
        DiffPumpOn = 1.;

        /* >>chng 03 nov 28, add option to not do line transfer */
        lgDoLineTrans = true;

        /* flag saying whether to constantly reevaluated opacities -
         * set false with no opacity reevaluate command */
        lgOpacityReevaluate = true;

        /* flag saying whether to constantly reevaluated ionization -
         * set false with no ionization reevaluate command */
        lgIonizReevaluate = true;
        /* this element is default for choosing line width */
        fine_opac_nelem = ipIRON;
        /* there will be this many resolution elements in one FWHM for this element,
         * at the lowest temperature to be considered */
        fine_opac_nresolv = 1;
        /* continuum scale factor for case of time varying continuum */
        time_continuum_scale = 1.;
        /* will fine optical depths be punched? */
        lgSaveOpacityFine = false;

        /* first is set true if one of the incident continua needs to have
         * H-ionizing radiation blocked.  Second is set true is it is blocked
         * with extinguish command - want both true if first is true */
        lgMustBlockHIon = false;
        lgBlockHIon = false;
        strncpy( rfield.chCumuType, "MASS", sizeof(rfield.chCumuType));
}

const realnum *t_rfield::getCoarseTransCoef()
{
        // average opacity transmission coefficient fine to coarse
        if( opac.lgScatON && trans_coef_total_stale)
        {
                /* sum over coarse continuum */
                for( long i=0; i < nflux-1; i++ )
                {
                        // find transmission coefficient if lower and upper bounds 
                        // of coarse continuum is within boundaries of fine continuum 
                        // unity is default
                        if( ipnt_coarse_2_fine[i] && ipnt_coarse_2_fine[i+1] )
                        {
                                // first branch is normal case, where fine continuum is finer than
                                // coarse continuum.  But, when end temp is very high, fine continuum is
                                // very coarse, so may be just one cell, and following will not pass
                                if( ipnt_coarse_2_fine[i+1]>ipnt_coarse_2_fine[i] )
                                {
                                        trans_coef_total[i] = 0.;
                                        for( long j=ipnt_coarse_2_fine[i]; j<ipnt_coarse_2_fine[i+1]; ++j )
                                                trans_coef_total[i] += sexp(fine_opt_depth[j]);
                                        trans_coef_total[i] /= (ipnt_coarse_2_fine[i+1]-ipnt_coarse_2_fine[i]);
                                }
                                else
                                {
                                        // in case where fine is coarser than coarse, 
                                        // just use first cell
                                        trans_coef_total[i] = sexp(fine_opt_depth[ipnt_coarse_2_fine[i]]);
                                }
                        }
                }
                trans_coef_total_stale = false;
        }
        return trans_coef_total;
}

void t_rfield::setTrimming()
{
        DEBUG_ENTRY( "t_rfield::setTrimming()" );

        double threshold = 0.;
        if( rfield.FluxFaint > 0. )
        {
                double peak = 0.;
                for( long i=0; i < rfield.nflux; ++ i )
                        peak = max( peak, rfield.SummedCon[i]*rfield.anu(i)/rfield.widflx(i) );
                threshold = rfield.FluxFaint * peak;
        }
                
        // search for highest energy cell with positive flux, use 'i < nPositive' in loops
        nPositive = nflux;
        while( nPositive > 0 && SummedCon[nPositive-1] <= threshold )
                --nPositive;
}

/* report the total flux not including attenuated isotropic continua, if requested */
double flux_correct_isotropic( const bool lgSaveIsotr, const int nEmType, const int iflux )
{
        if( iflux < 0 || iflux >= rfield.nflux )
                return  0.;

        double this_flux = (double)rfield.flux[ nEmType ][ iflux ];
        if( nEmType == 0 &&
            ( ! continuum.lgPrtIsotropicCont || ! lgSaveIsotr ) )
        {
                        this_flux = (double)( rfield.flux_beam_const[ iflux ]
                                        + rfield.flux_beam_time[ iflux ] );
        }

        return this_flux;
}

double flux_correct_isotropic( const int nEmType, const int iflux )
{
        return flux_correct_isotropic( true, nEmType, iflux );
}