rt_line_driving.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_line_driving derive radiative acceleration due to line absorption of incident continuum,
 * return value is line radiative acceleration */
#include "cddefines.h"
#include "physconst.h"
#include "iso.h"
#include "dense.h"
#include "taulines.h"
#include "h2.h"
#include "atomfeii.h"
#include "rt.h"
 
/*RT_line_driving derive radiative acceleration due to line absorption of incident continuum,
 * return value is line radiative acceleration */
double RT_line_driving(void)
{
        long int i, 
          ipHi, 
          nelem, 
          ipLo,
          ipISO;
 
        double AllHeavy, 
          AllRest, 
          OneLine, 
          fe2drive, 
          forlin_v, 
          h2drive,
          accel_iso[NISO];
 
        /* following used for debugging */
        /* double 
          RestMax, 
          HeavMax, 
          hydromax;
          long int 
          ipRestMax, 
          ihmax; */
 
        DEBUG_ENTRY( "RT_line_driving()" );
 
        /* this function finds the total rate the gas absorbs energy
         * this result is divided by the calling routine to find the
         * momentum absorbed by the gas, and eventually the radiative acceleration
         *
         * the energy absorbed by the line is
         * Abundance * energy * A *(g_up/g_lo) * occnum * escape prob
         * where occnum is the photon occupation number, and the g's are
         * the ratios of statistical weights */
 
        /* total energy absorbed in this zone, per cubic cm
         * do hydrogen first */
 
        for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
        {
                accel_iso[ipISO] = 0;
                for( nelem=ipISO; nelem < LIMELM; nelem++ )
                {
                        if( (dense.IonHigh[nelem] >= nelem + 1-ipISO)  )
                        {
                                for( ipHi=1; ipHi < iso_sp[ipISO][nelem].numLevels_local; ipHi++ )
                                {
                                        /* do not put in highest level since its not real */
                                        for( ipLo=0; ipLo < ipHi - 1; ipLo++ )
                                        {
                                                /* do not include bogus lines */
                                                if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() > 0 )
                                                {
                                                        OneLine = iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().pump()*
                                                        iso_sp[ipISO][nelem].trans(ipHi,ipLo).EnergyErg()*
                                                        iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().PopOpc();
 
                                                        accel_iso[ipISO] += OneLine;
                                                }
                                        }
                                }
                        
                                if( iso_ctrl.lgDielRecom[ipISO] )
                                {       
                                        // SatelliteLines are indexed by lower level, summed over satellite levels
                                        for( ipLo=0; ipLo < iso_sp[ipISO][nelem].numLevels_local; ipLo++ )
                                        {
                                                /* do not include bogus lines */
                                                TransitionList::iterator tr = SatelliteLines[ipISO][nelem].begin()+ipSatelliteLines[ipISO][nelem][ipLo];
                                                if((*tr).ipCont() > 0 )
                                                {
                                                        OneLine = (*tr).Emis().pump()*
                                                                (*tr).EnergyErg()*
                                                                (*tr).Emis().PopOpc();
 
                                                        accel_iso[ipISO] += OneLine;
                                                }
                                        }
                                }
 
                                for( ipHi=iso_sp[ipISO][nelem].st[iso_sp[ipISO][nelem].numLevels_local-1].n()+1; ipHi < iso_ctrl.nLyman[ipISO]; ipHi++ )
                                {
                                        /* do not include bogus lines */
                                        TransitionList::iterator tr = ExtraLymanLines[ipISO][nelem].begin()+ipExtraLymanLines[ipISO][nelem][ipHi];
                                        if( (*tr).ipCont() > 0 )
                                        {
                                                OneLine = (*tr).Emis().pump()*
                                                        (*tr).EnergyErg()*
                                                        (*tr).Emis().PopOpc();
 
                                                accel_iso[ipISO] += OneLine;
                                        }
 
                                }
                        }
                }
        }
 
        /* all heavy element lines in calculation of cooling 
         * these are the level 1 lines */
        AllHeavy = 0.;
        for( i=1; i <= nLevel1; i++ )
        {
                OneLine = 
                  TauLines[i].Emis().pump()*
                  TauLines[i].EnergyErg()*
                  TauLines[i].Emis().PopOpc();
                AllHeavy += OneLine;
        }
 
        /* all heavy element lines treated with g-bar 
         * these are the level 2 lines, f should be ok */
        AllRest = 0.;
        for( i=0; i < nWindLine; i++ )
        {
                OneLine = 
                        TauLine2[i].Emis().pump()*
                        TauLine2[i].EnergyErg()*
                        TauLine2[i].Emis().PopOpc();
                AllRest += OneLine;
        }
        for( i=0; i < nUTA; i++ )
        {
                OneLine = 
                        UTALines[i].Emis().pump()*
                        UTALines[i].EnergyErg()*
                        UTALines[i].Emis().PopOpc();
                AllRest += OneLine;
        }
        for( i=0; i < nHFLines; i++ )
        {
                OneLine = 
                        HFLines[i].Emis().pump()*
                        HFLines[i].EnergyErg()*
                        HFLines[i].Emis().PopOpc();
                AllRest += OneLine;
        }
        
        for( long ipSpecies=0; ipSpecies<nSpecies; ipSpecies++ )
        {
                if( dBaseSpecies[ipSpecies].lgActive )
                {
                        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;
                                OneLine = (*tr).EnergyErg()*(*tr).Emis().pump()*(*tr).Emis().PopOpc();
                                AllRest += OneLine;
                        }
                }
        }
 
        /* the H2 molecule */
        h2drive = 0.;
        for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
                h2drive += (*diatom)->H2_Accel();
 
        /* The large model FeII atom */
        fe2drive = 0.;
        FeIIAccel(&fe2drive);
 
        forlin_v = AllHeavy + accel_iso[ipH_LIKE] + accel_iso[ipHE_LIKE] + 
                fe2drive + h2drive + AllRest;
 
        /*fprintf(ioQQQ," wind te %e %e %e %e %e\n",    
                AllHeavy , HydroAccel , fe2drive , he1l , AllRest );*/
        return( forlin_v );
}