prt_linepres.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 */
/*PrtLinePres print line radiation pressures for current conditions */
#include "cddefines.h"
#include "pressure.h"
#include "taulines.h"
#include "iso.h"
#include "dense.h"
#include "lines_service.h"
#include "h2.h"
#include "prt.h"
#include "mole.h"
/* faintest pressure we will bother with */
#define THRESH  0.05
 
void PrtLinePres( FILE *ioPRESSURE )
{
        long int i, 
          ip,
          ipLo,
          ipHi,
          nelem;
        int ier;
        double RadPres1;
 
        // this is limit to number of lines we can store at one time
        const int nLine = 100;
        /* labels, wavelengths, and fraction of total pressure, for important lines */
        char chLab[nLine][5];
        realnum wl[nLine] , frac[nLine];
        long int iperm[nLine];
 
        /* will be used to check on size of opacity, was capped at this value */
        realnum smallfloat=SMALLFLOAT*10.f;
 
        DEBUG_ENTRY( "PrtLinePres()" );
 
        /* this routine only called if printout of contributors to line
         * radiation pressure is desired */
 
        /* compute radiation pressure due to iso lines */
        ip = 0;
 
        for(i=0; i<nLine; ++i)
        {
                frac[i] = FLT_MAX;
                wl[i] = FLT_MAX;
        }
 
        if( pressure.pres_radiation_lines_curr > 1e-30 )
        {
                for( long ipISO = 0; ipISO<NISO; ipISO++ )
                {
                        for( nelem=ipISO; nelem < LIMELM; nelem++ )
                        {
                                /* does this ion stage exist? */
                                if( dense.IonHigh[nelem] >= nelem + 1 - ipISO )
                                {
                                        /* do not include highest levels since maser can occur due to topoff,
                                         * and pops are set to small number in this case */
                                        for( ipHi=1; ipHi <iso_sp[ipISO][nelem].numLevels_local - iso_sp[ipISO][nelem].nCollapsed_local; ipHi++ )
                                        {
                                                for( ipLo=0; ipLo < ipHi; ipLo++ ) 
                                                {
                                                        if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).ipCont() <= 0 )
                                                                continue;
 
                                                        ASSERT( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().Aul() > iso_ctrl.SmallA );
 
                                                        /*NB - this code must be kept parallel with code in pressure_total */
                                                        if( iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().PopOpc() > smallfloat &&
                                                                /* >>chng 01 nov 01, add test that have not overrun optical depth scale */
                                                                ( (iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauTot() - iso_sp[ipISO][nelem].trans(ipHi,ipLo).Emis().TauIn()) > smallfloat ) )
                                                        {
                                                                RadPres1 =  PressureRadiationLine( iso_sp[ipISO][nelem].trans(ipHi,ipLo), GetDopplerWidth(dense.AtomicWeight[nelem]) );
 
                                                                if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH )
                                                                {
                                                                        wl[ip] = iso_sp[ipISO][nelem].trans(ipHi,ipLo).WLAng();
 
                                                                        /* put null terminated line label into chLab using line array*/
                                                                        chIonLbl(chLab[ip], iso_sp[ipISO][nelem].trans(ipHi,ipLo));
                                                                        frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr);
 
                                                                        ip = MIN2((long)nLine-1,ip+1);
                                                                }
                                                        }
                                                }
                                        }
                                }
                        }
                }
 
                /* line radiation pressure from large set of level 1 lines */
                for( i=1; i <= nLevel1; i++ )
                {
                        if( (*TauLines[i].Hi()).Pop() > 1e-30 )
                        {
                                RadPres1 =  PressureRadiationLine( TauLines[i], GetDopplerWidth(dense.AtomicWeight[(*TauLines[i].Hi()).nelem()-1]) );
                        }
                        else
                        {
                                RadPres1 = 0.;
                        }
 
                        if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH )
                        {
                                wl[ip] = TauLines[i].WLAng();
 
                                /* put null terminated line label into chLab using line array*/
                                chIonLbl(chLab[ip], TauLines[i]);
                                frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr);
 
                                ip = MIN2((long)nLine-1,ip+1);
                        }
                }
 
                for( i=0; i < nWindLine; i++ )
                {
                        if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
                        {
                                if( (*TauLine2[i].Hi()).Pop() > 1e-30 )
                                {
                                        RadPres1 =  PressureRadiationLine( TauLine2[i], GetDopplerWidth(dense.AtomicWeight[(*TauLine2[i].Hi()).nelem()-1]) );
                                }
                                else
                                {
                                        RadPres1 = 0.;
                                }
 
                                if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH )
                                {
                                        wl[ip] = TauLine2[i].WLAng();
 
                                        /* put null terminated line label into chLab using line array*/
                                        chIonLbl(chLab[ip], TauLine2[i]);
                                        frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr);
 
                                        ip = MIN2((long)nLine-1,ip+1);
                                }
                        }
                }
 
                for( i=0; i < nHFLines; i++ )
                {
                        if( (*HFLines[i].Hi()).Pop() > 1e-30 )
                        {
                                RadPres1 =  PressureRadiationLine( HFLines[i], GetDopplerWidth(dense.AtomicWeight[(*HFLines[i].Hi()).nelem()-1]) ); 
                        }
                        else
                        {
                                RadPres1 = 0.;
                        }
 
                        if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH )
                        {
                                wl[ip] = HFLines[i].WLAng();
 
                                /* put null terminated line label into chLab using line array*/
                                chIonLbl(chLab[ip], HFLines[i]);
                                frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr);
 
                                ip = MIN2((long)nLine-1,ip+1);
                        }
                }
 
                /* lines from external databases */
                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;
                                        if( (*(*tr).Hi()).Pop() > 1e-30 )
                                        {
                                                RadPres1 =  PressureRadiationLine( *tr, DopplerWidth );
                                        }
                                        else
                                        {
                                                RadPres1 = 0.;
                                        }
                                        
                                        if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH )
                                        {
                                                wl[ip] = (*tr).WLAng();
                                                
                                                /* put null terminated line label into chLab using line array*/
                                                chIonLbl(chLab[ip], *tr );
                                                frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr);
 
                                                ip = MIN2((long)nLine-1,ip+1);
                                        }
                                }       
                        }       
                }       
 
                /* radiation pressure due to H2 */
                for( diatom_iter diatom = diatoms.begin(); diatom != diatoms.end(); ++diatom )
                {
                        RadPres1 = (*diatom)->H2_RadPress();
                        if( RadPres1 > pressure.pres_radiation_lines_curr*THRESH )
                        {
                                wl[ip] = 0;
 
                                /* put null terminated 4 char line label into chLab using line array*/
                                strcpy(chLab[ip], "H2  ");
                                frac[ip] = (realnum)(RadPres1/pressure.pres_radiation_lines_curr);
 
                                ip = MIN2((long)nLine-1,ip+1);
                        }
                }
 
                /* return if no significant contributors to radiation pressure found */
                if( ip<= 0 )
                {
                        fprintf( ioPRESSURE, "\n" );
                        return;
                }
 
                /* sort from strong to weak, then only print strongest */
                spsort(
                        /* input array to be sorted */
                        frac, 
                        /* number of values in x */
                        ip, 
                        /* permutation output array */
                        iperm, 
                        /* flag saying what to do - 1 sorts into increasing order, not changing
                        * the original routine */
                        -1, 
                        /* error condition, should be 0 */
                        &ier);
 
                /* now print up to ten strongest contributors to radiation pressure */
                fprintf( ioPRESSURE, " P(Lines):" );
                for( i=0; i < MIN2(10,ip); i++ )
                {
                        int ipline = iperm[i];
                        fprintf( ioPRESSURE, "(%4.4s ", chLab[ipline]);
                        prt_wl(ioPRESSURE, wl[ipline] );
                        fprintf(ioPRESSURE," %.2f) ",frac[ipline]);
                }
 
                /* finally end the line */
                fprintf( ioPRESSURE, "\n" );
        }
        return;
}