rfield.h

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2023 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */

#ifndef RFIELD_H_
#define RFIELD_H_

/* rfield.h */
#include "energy.h"
#include "module.h"
#include "vectorize.h"
#include "mesh.h"
#include "container_classes.h"

const double WL_V_FILT = 5500.;

const double WL_B_FILT = 4400.;

const int LIMSPC = 100;

void rfield_opac_zero( long lo , long ihi );

extern bool lgRfieldAllocated;

namespace Illuminate {
        typedef enum { FORWARD , REVERSE , ISOTROPIC } IlluminationType ;
}

namespace Accumulate {
        typedef enum { INSTANTANEOUS = 0, CUMULATIVE = 1 } AccumulationType;
}


class Spectrum
{
private:
        vector<realnum> flux[2];

        double norm;

        double time_avg_flux;

        long int nflux;


public:
        // constructor
        Spectrum()
        {
                norm = 1.;
                time_avg_flux = 0.;
        }

        void resize( long int nbins );

        void accumulate_flux( const double this_cumulfac, double timestep, double total_time )
        {
                DEBUG_ENTRY( "accumulate_flux()" );

                vector<realnum> &flux_inst = flux[0];
                vector<realnum> &flux_cumul = flux[1];

                double mean = flux_cumul[ 0 ] * norm + flux_inst[ 0 ] * this_cumulfac;
        
                for( long ip = 1; ip < nflux; ip++ )
                {
                        /* >>refer      running mean    Welford, 1962, Technometrics, 4-419 */
                        double tmp_fl = flux_cumul[ ip ] * norm +
                                        flux_inst[ ip ] * this_cumulfac;
                        mean += (tmp_fl - mean) / double( ip+1 );
                }
        
                if( fp_equal( mean, 0. ) )
                       return;
        
                double mean_inv = 1. / mean;
        
                for( long ip = 0; ip < nflux; ip++ )
                {
                        double tmp_flux = flux_cumul[ ip ] * norm +
                                        flux_inst[ ip ] * this_cumulfac;
                        flux_cumul[ ip ] = realnum( tmp_flux * mean_inv );
                }
        
                norm = mean;

                if( fp_equal( time_avg_flux, 0.0 ) )
                {
                        time_avg_flux = mean;
                }
                else
                {
                        time_avg_flux += ( mean - time_avg_flux ) *
                                        timestep / ( timestep + total_time );
                }
        
                return;
        }

        double get_flux( long int accType, long int ip ) const
        {
                DEBUG_ENTRY( "get_flux" );

                ASSERT( accType == Accumulate::AccumulationType::INSTANTANEOUS ||
                        accType == Accumulate::AccumulationType::CUMULATIVE );

                ASSERT( ip >= 0 && ip < nflux );

                if( accType == Accumulate::AccumulationType::INSTANTANEOUS )
                {
                        return double( flux[0][ip] );
                }
                else
                {
                        return double( flux[1][ip] ) * norm;
                }
        }

        double get_time_avg_flux() const
        {
                return time_avg_flux;
        }

        vector<realnum>& operator[]( int accType )
        {
                ASSERT( accType == Accumulate::AccumulationType::INSTANTANEOUS ||
                        accType == Accumulate::AccumulationType::CUMULATIVE );
                return flux[ accType ];
        }
};

struct t_rfield : public module, public t_mesh {
        const char *chName() const
        {
                return "rfield";
        }
        
        void zero();
        void comment(t_warnings&) {}
        long int nflux;

        long int nflux_with_check;

        long int nPositive;

        double FluxFaint;

        vector<long> line_count;

        vector<realnum> OccNumbContEmitOut;

        Spectrum flux;

        vector<realnum> flux_isotropic;

        vector<realnum> flux_beam_time, flux_beam_const;

        vector<realnum> flux_accum;

        vector<realnum> ExtinguishFactor;
        realnum ExtinguishLeakage, 
         ExtinguishColumnDensity, 
         ExtinguishLowEnergyLimit,
         ExtinguishConvertColDen2OptDepth ,
         ExtinguishEnergyPowerLow;

        bool lgMustBlockHIon;

        /* this is set true if H-ionizing radiation is blocked with extinguish
         * command */
        bool lgBlockHIon;

        bool lgDoLineTrans;

        bool lgOpacityReevaluate;

        bool lgCMB_set;

        bool lgIonizReevaluate;

        vector<realnum> convoc;

        vector<realnum> OccNumbIncidCont;

        vector<realnum> OccNumbDiffCont;

        vector_avx<double> ContBoltz;

        vector_avx<double> ContBoltzHelp1;
        vector_avx<double> ContBoltzHelp2;

        vector_avx<double> vexp_arg;

        vector_avx<double> ContBoltzAvg;

        multi_arr<realnum,2> ConEmitLocal/* [depth][energy]*/;

        multi_arr<realnum,2> ConSourceFcnLocal/* [depth][energy]*/;

        Spectrum ConEmitReflec;

        Spectrum ConEmitOut;

        vector<realnum> ConInterOut;

        Spectrum ConRefIncid;

        vector<double> SummedCon;
        vector<realnum> SummedDif;
        vector<realnum> SummedOcc;
        vector<realnum> SummedDifSave;

        vector<realnum> ConOTS_local_photons;
        vector<realnum> ConOTS_local_OTS_rate;

        vector<realnum> DiffuseEscape;

        vector<realnum> TotDiff2Pht;

        vector<realnum> otslin;
        vector<realnum> otscon;
        multi_arr<realnum,2> otssav;

        Spectrum outlin;
        vector<realnum> outlin_noplot;

        vector<realnum> DiffuseLineEmission;

        Spectrum reflin;

        Spectrum flux_total_incident;
        vector<realnum> flux_beam_const_save, flux_time_beam_save, flux_isotropic_save;

        realnum time_continuum_scale;

        realnum DiffPumpOn;

        vector<string> chLineLabel/*[rfield.nflux_with_check][5]*/;

        vector<string> chContLabel/*[rfield.nflux_with_check][5]*/;

        char chDffTrns[4];

        bool lgOutOnly;

        long int ipEnergyBremsThin;
        realnum EnergyBremsThin;

        long int ipMaxBolt;

        bool lgInducProcess;

        vector<double> comup, comdn;

        long int ipB_filter , ipV_filter;

        long int ipG0_TH85_lo , ipG0_TH85_hi;

        long int ipG0_DB96_lo , ipG0_DB96_hi;

        long int ipG0_spec_lo , ipG0_spec_hi;

        long int ip1000A;

        double extin_mag_B_point , extin_mag_V_point;

        double extin_mag_B_extended , extin_mag_V_extended;

        double opac_mag_B_point, opac_mag_V_point, opac_mag_B_extended , opac_mag_V_extended;

        vector<realnum> csigh, csigc;

        vector<realnum> eeBremsDif;

        double comtot, 
          cmheat, 
          cmcool, 
          cinrat;
        bool lgComptonOn;

        bool lgComUndr;

        double totpow[LIMSPC], 
          slope[LIMSPC], 
          cutoff[LIMSPC][3], 
          spfac[LIMSPC];

        bool lgTimeVary[LIMSPC];

        /* beamed or isotropic continuum?  if isotropic then does not vary
         * with time */
        bool lgBeamed[LIMSPC];

        realnum OpticalDepthScaleFactor[LIMSPC];

        Illuminate::IlluminationType Illumination[LIMSPC];

        long int nShape, 
          ipSpec;

        vector<Energy> tNu[LIMSPC];
        vector<realnum> tslop[LIMSPC];
        vector<realnum> tFluxLog[LIMSPC];
        double TableRadius[LIMSPC];
        bool lgSphericalDilution[LIMSPC];
        
        long ncont[LIMSPC];

        double RSFCheck[LIMSPC];

        double range[LIMSPC][2];

        char chSpNorm[LIMSPC][5], 
                chRSpec[LIMSPC][5], 
                chSpType[LIMSPC][6],
                chCumuType[5];

        realnum qhtot, 
          qhe, 
          qheii, 
          qbal, 
          qrad, 
          qtot;

        realnum uh;

        realnum uheii;

        bool lgUSphON;
        realnum rstrom;

        bool lgHabing;

        long int fine_opac_nelem;

        long int fine_opac_nresolv;

        realnum fine_opac_velocity_width;

private:
        vector<realnum> trans_coef_total;

public:
        bool trans_coef_total_stale;

        vector<long> ipnt_coarse_2_fine;

        realnum fine_ener_lo, fine_ener_hi;
        long nfine;
        double fine_resol;
        vector<realnum> fine_opac_zone;
        vector<realnum> fine_opt_depth;
        vector<realnum> fine_anu;

        unordered_map<long, vector<int> > fine_lstack;

        long int ipFineConVelShift;

        bool lgOpacityFine;
        bool lgSaveOpacityFine;

        bool lgLyaOTS;

        bool lgHeIIOTS;

        bool lgKillOutLine;

        bool lgKillOutCont;

        bool lgKillOTSLine;

        bool lgPlasNu;
        realnum plsfrq, 
                plsfrqmax;
        // the zone where the plasma frequency is evaluated 
        long int nZonePlsFrqEval;
        long int ipPlasma, 
          ipPlasmax;

        bool lgMMok, 
          lgHPhtOK, 
          lgXRayOK, 
          lgGamrOK;

        realnum EnerGammaRay;
        long int ipEnerGammaRay;

        bool lgHionRad;

        realnum occmax, 
          occmnu, 
          tbrmax, 
          tbrmnu, 
          tbr4nu, 
          occ1nu;

        bool lgOcc1Hi;

        realnum EnergyIncidCont ,
                EnergyDiffCont;

        // constructor
        t_rfield( )
        {
                nZonePlsFrqEval = -1;

                // the constant that multiplies the column density to get optical depth at 1 Ryd
                ExtinguishConvertColDen2OptDepth = (realnum)6.22e-18;
                // the power on the energy for the extinction
                ExtinguishEnergyPowerLow = (realnum)-2.43;
        }
        
        const realnum *getCoarseTransCoef();

        void setCoarseTransCoefPtr(size_t size)
        {
                trans_coef_total.resize(size);
        }
        void resetCoarseTransCoef()
        {
                for (long i=0; i<nflux_with_check; ++i)
                        trans_coef_total[i] = 1.0;
                trans_coef_total_stale = true;
        }
        void setTrimming();
};
extern t_rfield rfield;

double flux_correct_isotropic( const bool lgSaveIsotrp, const int nEmType, const int j );

double flux_correct_isotropic( const int nEmType, const int j );

#endif /* RFIELD_H_ */