prt_lines.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 */
/*lines main routine to put emission line intensities into line stack,
 * calls lineset1, 2, 3, 4 */
/*Drive_cdLine do the drive cdLine command */
/*FindStrongestLineLabels find strongest lines contributing to point in continuum energy mesh, output in some save commands */
#include "cddefines.h"
#include "physconst.h"
#include "taulines.h"
#include "thermal.h"
#include "yield.h"
#include "ipoint.h"
#include "ionbal.h"
#include "cddrive.h"
#include "iterations.h"
#include "trace.h"
#include "dense.h"
#include "prt.h"
#include "rt.h"
#include "coolheavy.h"
#include "rfield.h"
#include "phycon.h"
#include "elementnames.h"
#include "iso.h"
#include "hyperfine.h"
#include "hydrogenic.h"
#include "lines_service.h"
#include "atmdat.h"
#include "lines.h"
#include "radius.h"
 
STATIC void Drive_cdLine( void );
 
void lines(void)
{
        char chLabel[5];
        long int i, 
          ipnt,
          nelem;
        double BigstExtra, 
          ExtraCool,  
          f2, sum; 
 
        DEBUG_ENTRY( "lines()" );
 
        /* LineSave.ipass
         * -1 - space not yet allocated - just count number of lines entered into stack
         *  0 - first call with space allocated - must create labels and add in wavelengths
         * +1 - later calls - just add intensity 
         */
 
        /* major routines used here:
         *
         * PutLine( tarray )
         * this uses information in tarray to give various
         * contributions to lines, and their intensities
         *
         * PutExtra( extra )
         * extra is some extra intensity to add to the line
         * it will go into the totl contribution put out by PutLine,
         * and this contribution should be indicated by independent
         * call to linadd
         * */
 
        if( trace.lgTrace )
        {
                fprintf( ioQQQ, " lines called\n" );
        }
 
        /* the drive cdline command, which checks that all lines can be pulled out by cdLine */
        if( trace.lgDrv_cdLine  && LineSave.ipass > 0 )
                Drive_cdLine();
 
        /* total luminosity radiated by this model - will be compared with energy in incident
         * continuum when calculation is complete */
        thermal.power += thermal.htot*radius.dVeffAper;
 
        /* remember the total free-free heating */
        fixit(); // get rid of brems_heat_total entirely (only net heating is added into stack, so use net here to avoid nonsense ratios elsewhere)  
        //thermal.FreeFreeTotHeat += CoolHeavy.brems_heat_total*radius.dVeffAper;
        thermal.FreeFreeTotHeat += thermal.heating[0][11]*radius.dVeffAper;
 
        /* total Compton heating - cooling */
        rfield.comtot += rfield.cmheat*radius.dVeffAper;
        thermal.totcol += thermal.ctot*radius.dVeffAper;
 
        /* up up induced recombination cooling */
        for( nelem=0; nelem<LIMELM; ++nelem )
        {
                hydro.cintot += iso_sp[ipH_LIKE][nelem].RecomInducCool_Rate*radius.dVeffAper;
        }
 
        /* nsum is line pointer within large stack of line intensities */
        LineSave.nsum = 0;
        LineSave.nComment = 0;
 
        /* this is used by lindst to proportion inward and outward.  should be 50% for
         * optically thin line.  putline sets this to actual value for particular line
         * and calls lindst then rests to 50% */
        rt.fracin = 0.5;
 
        /* last arg in call to lindst and linadd is info on line
         * info is char variable indicating type of line this is
         * 'c' cooling
         * 'h' heating
         * 'i' information only
         * 'r' recombination line
         *
         * all components of lines are entered into the main line stack here
         * when printing, filters exist to not print Inwd component */
 
        /* initialize routine that can generate pointers for forbidden lines,
         * these are lines that are not transferred otherwise,
         * in following routines there will be pairs of calls, first to
         * PntForLine to get pointer, then lindst to add to stack */
        PntForLine(0.,"FILL",&i);
 
        /* evaluate rec coefficient for rec lines of C, N, O
         * some will be used in LineSet2 and then zeroed out,
         * others left alone and used below */
        t_ADfA::Inst().rec_lines(phycon.te,LineSave.RecCoefCNO);
 
        /* initialize ExtraInten with a zero */
        PutExtra(0.);
 
        /* put in something impossible in element 0 of line stack */
        linadd(0.f,0,"zero",'i' , "null placeholder");
 
        /* this is compared with true volume in final.  The number can't
         * actually be unity since this would overflow on a 32 bit machine */
        /* integrate the volume as a sanity check */
        linadd( 1.e-10 , 1 , "Unit" , 'i' , "unit integration placeholder");
        static long int ipOneAng=-1;
        if( LineSave.ipass<0 )
                ipOneAng = ipoint( RYDLAM );
        lindst( 1.e-10 , 1. , "UntD" , ipOneAng , 'i' , false,"unit integration placeholder");
 
        /* initial set of general properties */
        lines_general();
 
        /* do all continua */
        lines_continuum();
 
        /* information about grains */
        lines_grains();
 
        /* update all satellite lines */
        for( long nelem=ipHYDROGEN; nelem < LIMELM; ++nelem )
                iso_satellite_update(nelem);
 
        /* do all hydrogenic ions */
        lines_hydro();
 
        /* enter He-iso sequence lines */
        lines_helium();
 
#if     0
        /* This is Ryan's code for dumping lots of Helium lines according to
         * quantum number rather than wavelength, principally for comparison with Rob
         * Bauman's code. */
        if( iteration > 1 )
        {
                fprintf( ioQQQ,"ipHi\tipLo\tnu\tlu\tsu\tnl\tll\tsl\tWL\tintens\n" );
                for( long ipHi=5; ipHi<= iso_sp[ipHE_LIKE][ipHELIUM].numLevels_local - iso_sp[ipHE_LIKE][ipHELIUM].nCollapsed_local; ipHi++ )
                {
                        for( long ipLo=0; ipLo<ipHi; ipLo++ )
                        {
                                if( iso_sp[ipHE_LIKE][ipHELIUM].trans(ipHi,ipLo).ipCont() > 0 )
                                {
                                        double relint, absint;
 
                                        if( cdLine("He 1", 
                                                iso_sp[ipHE_LIKE][ipHELIUM].trans(ipHi,ipLo).WLAng(),
                                                &relint, &absint ) )
                                        {
                                                //iso_sp[ipHE_LIKE][ipHELIUM].trans(ipHi,ipLo).Hi()->chLabel
 
                                                if( iso_sp[ipHE_LIKE][ipHELIUM].trans(ipHi,ipLo).WLAng() < 1.1E4 &&
                                                        iso_sp[ipHE_LIKE][ipHELIUM].trans(ipHi,ipLo).WLAng() > 3.59E3 &&
                                                        ipLo!=3 && ipLo!=4 && relint >= 0.0009 )
                                                {
                                                        fprintf( ioQQQ,"%li\t%li\t%li\t%li\t%li\t%li\t%li\t%li\t%e\t%e\n",
                                                                ipHi,
                                                                ipLo,
                                                                iso_sp[ipHE_LIKE][ipHELIUM].st[ipHi].n(),
                                                                iso_sp[ipHE_LIKE][ipHELIUM].st[ipHi].l(),
                                                                iso_sp[ipHE_LIKE][ipHELIUM].st[ipHi].S(),
                                                                iso_sp[ipHE_LIKE][ipHELIUM].st[ipLo].n(),
                                                                iso_sp[ipHE_LIKE][ipHELIUM].st[ipLo].l(),
                                                                iso_sp[ipHE_LIKE][ipHELIUM].st[ipLo].S(),
                                                                iso_sp[ipHE_LIKE][ipHELIUM].trans(ipHi,ipLo).WLAng(),
                                                                relint );
                                                }
                                        }
                                }
                        }
                }
        }
#endif
 
        // these must come before the old level 1 or 2 lines.  we now have the option
        // to use the external database by default, or turn it off (set chianti off)
        // and fall back to the old line treatment.  When database is off those lines
        // do not exist.  But when database is on the level 1 lines are still evaluated
        // but with the ion abundance set to zero.  If the level 1 lines were entered
        // into the stack then searches for the line with cdLine would turn up the level 1
        // line, which would be zero.
        // The long term goal is to have all lines be external databases & rm level 1&2 liens
 
        /* external database lines */
        i = StuffComment( "database lines" );
        linadd( 0., (realnum)i , "####", 'i' ,
                "database lines ");
        for (int ipSpecies=0; ipSpecies < nSpecies; ++ipSpecies)
        {
                for( EmissionList::iterator em=dBaseTrans[ipSpecies].Emis().begin();
                          em != dBaseTrans[ipSpecies].Emis().end(); ++em)
                {
                        /* \todo 2 say which database in the comment */
                        if( (*em).Tran().ipCont() > 0)
                        {
                                PutLine((*em).Tran(), "lines from third-party databases", (*(*em).Tran().Hi()).chLabel());
                        }
                }
        }
 
        /* do heavies, lithium through neon */
        lines_lv1_li_ne();
 
        /* do heavies, sodium through argon */
        lines_lv1_na_ar();
 
        /* do heavies, potassium through zinc */
        lines_lv1_k_zn();
 
        /* add up line intensities for certain set of lines */
        sum = PrtLineSum();
        /* zero out the location that will receive this information, 
         * remember that memory not allocated until ipass >= 0 */
        if( LineSave.ipass > 0 )
        {
                LineSv[LineSave.nsum].SumLine[0] = 0.;
                LineSv[LineSave.nsum].SumLine[1] = 0.;
        }
        /* optional sum of certain emission lines, set with "print sum" */
        linadd(sum/radius.dVeffAper,0,"Stoy",'i' , 
                "Stoy method energy sum ");
 
        /* next come some recombination lines */
        i = StuffComment( "recombination" );
        linadd( 0., (realnum)i , "####", 'i' ,
                "recombination lines");
 
        /***********************************************************************
         * large number of C, N, and O recombination lines                     *
         *************************************************************************/
 
        for( i=0; i < 471; i++ )
        {
                /* generate label for the line if ipass is -1 or 0 - saved in arrays
                 * so no need to do this during production */
                if( LineSave.ipass <= 0 )
                {
                        /* generate label for the line */
                        strcpy( chLabel, elementnames.chElementSym[(long)(LineSave.RecCoefCNO[0][i])-1] );
                        strcat( chLabel, elementnames.chIonStage[(long)(LineSave.RecCoefCNO[0][i]-
                                LineSave.RecCoefCNO[1][i]+1.01)-1] );
                }
                else
                        chLabel[0] = ' ';
 
                /* number of rec per unit vol
                 * do not predict these simple reccombination intensities at high densities
                 * since lines neglect collisional deexciation and line optical depth effects.
                 * They were not intended for high densities or column densities.
                 * As a result they become unphysically bright at high densities and
                 * violate the black body limit.  There would be major
                 * energy conservation problems if they were added in the outward beam in
                 * dense simulations.
                 * */
                if( dense.eden < 1e8  )
                {
                        nelem = (long)LineSave.RecCoefCNO[0][i]-1;
                        long int ion = (long)(LineSave.RecCoefCNO[0][i]-LineSave.RecCoefCNO[1][i]+2)-1;
                        f2 = LineSave.RecCoefCNO[3][i]*dense.eden*
                                dense.xIonDense[nelem][ion];
 
                        /* convert to intensity */
                        f2 = f2*1.99e-8/LineSave.RecCoefCNO[2][i];
                }
                else
                {
                        f2 = 0.;
                }
                /* stuff it into the stack */
                PntForLine(LineSave.RecCoefCNO[2][i],chLabel,&ipnt);
                lindst(f2,LineSave.RecCoefCNO[2][i],chLabel,ipnt, 'r',true ,
                        "recombination line");
        }
 
        /* next come the atom_level2 lines */
        i = StuffComment( "level2 lines" );
        linadd( 0., (realnum)i , "####", 'i' ,
                "level2 lines");
 
        /* add in all the other level 2 wind lines
         * Dima's 6k lines */
        ExtraCool = 0.;
        BigstExtra = 0.;
        for( i=0; i < nWindLine; i++ )
        {
                if( (*TauLine2[i].Hi()).IonStg() < (*TauLine2[i].Hi()).nelem()+1-NISO )
                {
                        PutLine(TauLine2[i],"level 2 line");
                        if( TauLine2[i].Coll().cool() > BigstExtra )
                        {
                                BigstExtra = TauLine2[i].Coll().cool();
                                thermal.ipMaxExtra = i+1;
                        }
                        ExtraCool += TauLine2[i].Coll().cool();
                }
        }
        /* keep track of how important this is */
        thermal.GBarMax = MAX2(thermal.GBarMax,(realnum)(ExtraCool/thermal.ctot));
 
        /* next come the hyperfine structure lines */
        i = StuffComment( "hyperfine structure" );
        linadd( 0., (realnum)i , "####", 'i' ,
                "hyperfine structure lines ");
 
        /* this is total cooling due to all HF lines */
        linadd( hyperfine.cooling_total, 0., "hfin", 'i' ,
                "total cooling all hyperfine structure lines ");
 
        /* remember largest local cooling for possible printout in comments */
        hyperfine.cooling_max = (realnum)MAX2(hyperfine.cooling_max,hyperfine.cooling_total/thermal.ctot);
 
        /* the hyperfine lines */
        for( i=0; i < nHFLines; i++ )
        {
                PutLine(HFLines[i],
                        "hyperfine structure line");
        }
 
        /* next come the inner shell fluorescent lines */
        i = StuffComment( "inner shell" );
        linadd( 0., (realnum)i , "####", 'i' ,
                "inner shell lines");
 
        /* the group of inner shell fluorescent lines */
        for( i=0; i < t_yield::Inst().nlines(); ++i )
        {
                double xInten = 
                        /* density of parent ion, cm-3 */
                        dense.xIonDense[t_yield::Inst().nelem(i)][t_yield::Inst().ion(i)] *
                        /* photo rate per atom per second, s-1 */
                        ionbal.PhotoRate_Shell[t_yield::Inst().nelem(i)][t_yield::Inst().ion(i)][t_yield::Inst().nshell(i)][0]*
                        /* fluor yield - dimensionless */
                        t_yield::Inst().yield(i) *
                        /* photon energy - ryd, converted into ergs */
                        t_yield::Inst().energy(i) * EN1RYD;
 
                /* create label if initializing line stack */
                if( LineSave.ipass == 0 )
                {
                        /* only generate the line label if it is going to be used */
                        strcpy( chLabel , elementnames.chElementSym[t_yield::Inst().nelem(i)] );
                        strcat( chLabel , elementnames.chIonStage[t_yield::Inst().ion_emit(i)] );
#                       if 0
                        /* only print yields for atoms */
                        if( t_yield::Inst().ion(i) == 0 && t_yield::Inst().nelem()(i) == ipIRON )
                        fprintf(ioQQQ,"DEBUGyeild\t%s\t%.3f\t%.3e\n",
                                /* line designation, energy in eV, yield */
                                chLabel , t_yield::Inst().energy()(i)*EVRYD, t_yield::Inst().yield(i) );
#                       endif
                }
 
                /* the group of inner shell fluorescent lines */
                lindst(
                        /* intensity of line */
                        xInten,
                        /* wavelength of line in Angstroms */
                        (realnum)RYDLAM / t_yield::Inst().energy(i),
                        /* label */
                        chLabel ,
                        /* continuum array offset for line as set in ipoint */
                        t_yield::Inst().ipoint(i), 
                        /* type of line - count as a recombination line */
                        'r',
                        /* include line in continuum? */
                        true ,
                        "inner shell line");
        }
 
        /* >>chng 06 jan 03, confirm that number of lines never changes once we have
         * created the labels */
        {
                static long nLineSave=-1 , ndLineSave=-1;
                if( LineSave.ipass == 0 )
                {
                        nLineSave = LineSave.nsum;
                        ndLineSave = LineSave.nsum;
                }
                else if( LineSave.ipass > 0 )
                {
                        /* this can't happen */
                        if( nLineSave<= 0 || ndLineSave < 0 )
                                TotalInsanity();
 
                        /* now make sure that we have the same number of lines as we had previously
                         * created labels.  This would not pass if we did not add in exactly the same
                         * number of lines on each pass */
                        if( nLineSave != LineSave.nsum )
                        {
                                fprintf( ioQQQ, "DISASTER number of lines in LineSave.nsum changed between pass 0 and 1 - this is impossible\n" );
                                fprintf( ioQQQ, "DISASTER LineSave.nsum is %li and nLineSave is %li\n",
                                        LineSave.nsum , 
                                        nLineSave);
                                ShowMe();
                                cdEXIT(EXIT_FAILURE);
                        }
                        if( ndLineSave != LineSave.nsum )
                        {
                                fprintf( ioQQQ, "DISASTER number of lines in LineSave.nsum changed between pass 0 and 1 - this is impossible\n" );
                                fprintf( ioQQQ, "DISASTER LineSave.nsum is %li and ndLineSave is %li\n",
                                        LineSave.nsum , 
                                        ndLineSave);
                                ShowMe();
                                cdEXIT(EXIT_FAILURE);
                        }
                }
        }
 
        /* now do all molecules - do last since so many H2 lines */
        lines_molecules();
 
        if( trace.lgTrace )
        {
                fprintf( ioQQQ, " lines returns\n" );
        }
        return;
}
 
/*Drive_cdLine do the drive cdLine command */
STATIC void Drive_cdLine( void )
{
        long int j;
        bool lgMustPrintHeader = true;
        double absval , rel;
 
        DEBUG_ENTRY( "Drive_cdLine()" );
 
        for( j=1; j < LineSave.nsum; j++ )
        {
                if( cdLine( LineSv[j].chALab , LineSv[j].wavelength , &absval , &rel ) <= 0 )
                {
                        /* print header if first miss */
                        if( lgMustPrintHeader )
                        {
                                fprintf(ioQQQ,"n\tlab\twl\n");
                                lgMustPrintHeader = false;
                        }
 
                        fprintf(ioQQQ,"%li\t%s\t%f\n", j, LineSv[j].chALab , LineSv[j].wavelength );
                }
        }
        fprintf( ioQQQ, " Thanks for checking on the cdLine routine!\n" );
        cdEXIT(EXIT_FAILURE);
}