rt_tau_inc.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*RT_tau_inc increment optical depths once per zone, called after radius_increment */
#include "cddefines.h"
#include "taulines.h"
#include "iso.h"
#include "rfield.h"
#include "trace.h"
#include "dense.h"
#include "hyperfine.h"
#include "wind.h"
#include "prt.h"
#include "conv.h"
#include "h2.h"
#include "mole.h"
#include "hmi.h"
#include "opacity.h"
#include "cooling.h"
#include "thermal.h"
#include "radius.h"
#include "atomfeii.h"
#include "rt.h"
#include "doppvel.h"
#include "mole.h"

/*RT_tau_inc increment optical depths once per zone, called after radius_increment */
void RT_tau_inc(void)
{

        long int i, 
          ipHi, 
          ipLo;

        DEBUG_ENTRY( "RT_tau_inc()" );

        if( trace.lgTrace )
        {
                fprintf( ioQQQ, " RT_tau_inc called.\n" );
        }

        /* call RT_line_all one last time in this zone, to get fine opacities defined */
        ASSERT( !conv.lgFirstSweepThisZone );
        conv.lgLastSweepThisZone = true;
        RT_line_all( );

        /* rfield.lgOpacityFine flag set false with no fine opacities command 
         * tests show that always evaluating this changes fast run of
         * pn_paris from 26.7 sec to 35.1 sec 
         * but must always update fine opacities since used for transmission */
        if( rfield.lgOpacityFine )
        {
                /* increment the fine opacity array */
                for( long i=0; i<rfield.nfine; ++i )
                {
                        realnum tauzone = rfield.fine_opac_zone[i]*(realnum)radius.drad_x_fillfac;
                        rfield.fine_opt_depth[i] += tauzone;
                }
                rfield.trans_coef_total_stale = true;
        }

        /* this may have updated some escape/destruction rates - force update
         * to all cooling lines */
        CoolEvaluate( &thermal.ctot );

        if( nzone <=1 )
        {
                opac.telec = (realnum)(radius.drad_x_fillfac*dense.eden*6.65e-25);
                opac.thmin = (realnum)(radius.drad_x_fillfac*findspecieslocal("H-")->den*3.9e-17*
                        (1. - rfield.ContBoltz[hmi.iphmin-1]/ hmi.hmidep));
        }
        else
        {
                opac.telec += (realnum)(radius.drad_x_fillfac*dense.eden*6.65e-25);
                opac.thmin += (realnum)(radius.drad_x_fillfac*findspecieslocal("H-")->den*3.9e-17*
                        (1. - rfield.ContBoltz[hmi.iphmin-1]/ hmi.hmidep));
        }

        /* prevent maser runaway */
        rt.dTauMase = 0;
        rt.mas_species = 0;
        rt.mas_ion = 0;
        rt.mas_hi = 0;
        rt.mas_lo = 0;

        /* all lines in iso sequences */
        for( long ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
        {
                for( long nelem=ipISO; nelem < LIMELM; nelem++ )
                {
                        /* this is the parent ion, for HI lines, is 1, 
                         * for element He is 1 for He-like (HeI) and 2 for H-like (HeII) */
                        int ion = nelem+1-ipISO;
                        /* do not evaluate in case where trivial parent ion */
                        if( ion <=dense.IonHigh[nelem] && dense.xIonDense[nelem][ion] > dense.density_low_limit )
                        {
                                if( iso_ctrl.lgDielRecom[ipISO] )
                                {
                                        // SatelliteLines are indexed by lower level
                                        for( ipLo=0; ipLo < iso_sp[ipISO][nelem].numLevels_local; ipLo++ )
                                        {
                                                RT_line_one_tauinc(SatelliteLines[ipISO][nelem][ipSatelliteLines[ipISO][nelem][ipLo]], ipISO, nelem, -1, ipLo, 
                                                        GetDopplerWidth(dense.AtomicWeight[nelem]) );
                                        }
                                }

                                for( ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
                                {
                                        for( ipLo=0; ipLo < ipHi; ipLo++ )
                                        {
                                                if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() <= 0 )
                                                        continue;

                                                /* actually do the work */
                                                RT_line_one_tauinc(iso_sp[ipISO][nelem].trans(ipHi,ipLo), ipISO, nelem, ipHi, ipLo, 
                                                        GetDopplerWidth(dense.AtomicWeight[nelem]) );
                                        }
                                }
                                ipLo = 0;
                                /* these are the extra Lyman lines */
                                for( ipHi=iso_sp[ipISO][nelem].st[iso_sp[ipISO][nelem].numLevels_local-1].n()+1; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ )
                                {
                                        TransitionList::iterator tr = ExtraLymanLines[ipISO][nelem].begin()+ipExtraLymanLines[ipISO][nelem][ipHi];
                                        (*tr).Emis().PopOpc() = iso_sp[ipISO][nelem].st[0].Pop();

                                        /* actually do the work */
                                        RT_line_one_tauinc(*tr, -1 ,ipISO, nelem, ipHi,
                                                GetDopplerWidth(dense.AtomicWeight[nelem]) );
                                }
                        }
                }
        }

        /* increment optical depths for all heavy element lines
         * same routine does wind and static,
         * does not start from 0 since first line is dummy */
        for( i=1; i <= nLevel1; i++ )
        {
                RT_line_one_tauinc(TauLines[i], -2, -2, -2, i, GetDopplerWidth(dense.AtomicWeight[(*TauLines[i].Hi()).nelem()-1]) );
        }

        /* all lines in cooling with g-bar */
        for( i=0; i < nWindLine; i++ )
        {
                /* do not include H-like or He-like in the level two lines since
                 * these are already counted in iso sequences */
                if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
                {
                        RT_line_one_tauinc(TauLine2[i], -3, -3, -3, i, GetDopplerWidth(dense.AtomicWeight[(*TauLine2[i].Hi()).nelem()-1]) );
                }
        }

        /* the block of inner shell lines */
        for( i=0; i < nUTA; i++ )
        {
                /* populations have not been set */
                UTALines[i].Emis().PopOpc() = dense.xIonDense[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1];
                (*UTALines[i].Lo()).Pop() = dense.xIonDense[(*UTALines[i].Hi()).nelem()-1][(*UTALines[i].Hi()).IonStg()-1];
                (*UTALines[i].Hi()).Pop() = 0.;
                RT_line_one_tauinc(UTALines[i], -4 , -4 , -4 , i, GetDopplerWidth(dense.AtomicWeight[(*UTALines[i].Hi()).nelem()-1]) );
        }

        /* all hyper fine structure lines  */
        for( i=0; i < nHFLines; i++ )
        {
                /* remember current gas-phase abundances */
                realnum save = dense.xIonDense[(*HFLines[i].Hi()).nelem()-1][(*HFLines[i].Hi()).IonStg()-1];

                /* bail if no abundance */
                if( save<=0. ) continue;

                /* set gas-phase abundance to total times isotope ratio */
                dense.xIonDense[(*HFLines[i].Hi()).nelem()-1][(*HFLines[i].Hi()).IonStg()-1] *= hyperfine.HFLabundance[i];

                RT_line_one_tauinc(HFLines[i] , -5 , -5 , -5 , i, GetDopplerWidth(dense.AtomicWeight[(*HFLines[i].Hi()).nelem()-1]) );

                /* put the correct gas-phase abundance back in the array */
                dense.xIonDense[(*HFLines[i].Hi()).nelem()-1][(*HFLines[i].Hi()).IonStg()-1] = save;
        }

        /* do large FeII atom if this is enabled */
        FeII_RT_TauInc();

        /* increment optical depth for the H2 molecule */
        for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
                (*diatom)->H2_RT_tau_inc();

        /* database Lines*/
        for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
        {
                if( dBaseSpecies[ipSpecies].lgActive )
                {
                        realnum DopplerWidth = GetDopplerWidth( dBaseSpecies[ipSpecies].fmolweight );
                        for (TransitionList::iterator tr=dBaseTrans[ipSpecies].begin(); 
                                  tr != dBaseTrans[ipSpecies].end(); ++tr)
                        {       
                                int ipHi = (*tr).ipHi();
                                if (ipHi >= dBaseSpecies[ipSpecies].numLevels_local || (*tr).ipCont() <= 0)
                                        continue;
                                int ipLo = (*tr).ipLo();

                                RT_line_one_tauinc( *tr, -10, ipSpecies, ipHi, ipLo, DopplerWidth );
                        }
                }
        }
        
        /* following is for static atmosphere */
        if( wind.lgStatic() )
        {
                /* iron fe feii fe2  - overlapping feii lines */
                t_fe2ovr_la::Inst().tau_inc();
        }

        if( trace.lgTrace && trace.lgOptcBug )
        {
                fprintf( ioQQQ, " RT_tau_inc updated optical depths:\n" );
                prtmet();
        }

        if( trace.lgTrace )
                fprintf( ioQQQ, " RT_tau_inc returns.\n" );

        return;
}