/home66/gary/public_html/cloudy/c08_branch/source/lines_service.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*PutLine enter local line intensity into the intensity stack for eventual printout */
/*PutExtra enter and 'extra' intensity source for some line */
/*linadd enter lines into the line storage array, called once per zone */
/*DumpLine print various information about an emission line vector, 
 * used in debugging, print to std out, ioQQQ */
/*TexcLine derive excitation temperature of line from contents of line array */
/*GetGF convert Einstein A into oscillator strength */
/*emit_frac returns fraction of populations the produce emission */
/*abscf convert gf into absorption coefficient */
/*chIonLbl use information in line array to generate a null terminated ion label in "Fe 2" */
/*chLineLbl use information in line transfer arrays to generate a line label */
/*RefIndex calculates the index of refraction of air using the line energy in wavenumbers,
 * used to convert vacuum wavelengths to air wavelengths. */
/*eina convert a gf into an Einstein A */
/*WavlenErrorGet - find difference between two wavelengths */
/*PutCS enter a collision strength into an individual line vector */
/*lindst add local line intensity to line luminosity stack */
/*PntForLine generate pointer for forbidden line */
/*TransitionZero zeros out TransitionZero */
/*EmLineJunk set all elements of transition struc to dangerous values */
/*EmLineZero set all elements of transition struc to zero */
/*LineConvRate2CS convert down coll rate back into electron cs in case other parts of code need this for reference */
/*lgTauGood returns true is we have not overrun optical depth scale */
/*OccupationNumberLine - derive the photon occupation number at line center for any line */
/*outline - adds line photons to reflin and outlin */
/*MakeCS compute collision strength by g-bar approximations */
/*gbar1 compute g-bar collision strength using Mewe approximations */
/*gbar0 compute g-bar gaunt factor for neutrals */
/*totlin sum total intensity of cooling, recombination, or intensity lines */
/*FndLineHt search through line heat arrays to find the strongest heat source */
/*ConvRate2CS convert down coll rate back into electron cs in case other parts of code need this for reference */
#include "cddefines.h"
#include "physconst.h"
#include "phycon.h"
#include "lines.h"
#include "radius.h"
#include "geometry.h"
#include "elementnames.h"
#include "rt.h"
#include "dense.h"
#include "rfield.h"
#include "opacity.h"
#include "ipoint.h"
#include "iso.h"
#include "taulines.h"
#include "hydrogenic.h"
#include "lines_service.h"

/*outline - adds line photons to reflin and outlin */
void outline( transition *t )
{
        long int ip = t->ipCont-1;
        double xInWrd = t->Emis->phots*t->Emis->FracInwd;

        DEBUG_ENTRY( "outline()" );

        ASSERT( t->Emis->phots >= 0. );
        ASSERT( t->Emis->FracInwd >= 0. );
        ASSERT( radius.BeamInIn >= 0. );
        ASSERT( radius.BeamInOut >= 0. );

        /* the reflected part */
        rfield.reflin[ip] += (realnum)(xInWrd*radius.BeamInIn);

        /* inward beam that goes out since sphere set */
        rfield.outlin[ip] += (realnum)(xInWrd*radius.BeamInOut*opac.tmn[ip]*
                t->Emis->ColOvTot);

        /* outward part */
        rfield.outlin[ip] += (realnum)(t->Emis->phots*
                (1. - t->Emis->FracInwd)*radius.BeamOutOut* t->Emis->ColOvTot);
        return;
}

/*emit_frac returns fraction of populations the produce emission */
double emit_frac( transition *t )
{
        DEBUG_ENTRY( "emit_frac()" );

        ASSERT( t->ipCont > 0 );

        /* collisional deexcitation and destruction by background opacities
         * are loss of photons without net emission */
        double deexcit_loss = t->Coll.cs * dense.cdsqte + t->Emis->Aul*t->Emis->Pdest;
        /* this is what is observed */
        double rad_deexcit = t->Emis->Aul*(t->Emis->Pelec_esc + t->Emis->Pesc);
        return rad_deexcit/(deexcit_loss + rad_deexcit);
}

/*DumpLine print various information about an emission line vector, 
 * used in debugging, print to std out, ioQQQ */
void DumpLine(transition * t)
{
        char chLbl[11];

        DEBUG_ENTRY( "DumpLine()" );

        ASSERT( t->ipCont > 0 );

        /* routine to print contents of line arrays */
        strcpy( chLbl, chLineLbl(t) );

        fprintf( ioQQQ, 
                "   %10.10s Te%.2e eden%.1e CS%.2e Aul%.1e Tex%.2e cool%.1e het%.1e conopc%.1e albdo%.2e\n", 
          chLbl, 
          phycon.te, 
          dense.eden, 
          t->Coll.cs, 
          t->Emis->Aul, 
          TexcLine(t), 
          t->Coll.cool, 
          t->Coll.heat ,
          opac.opacity_abs[t->ipCont-1],
          opac.albedo[t->ipCont-1]);

        fprintf( ioQQQ, 
                "  Tin%.1e Tout%.1e Esc%.1e eEsc%.1e DesP%.1e Pump%.1e OTS%.1e PopL,U %.1e %.1e PopOpc%.1e\n", 
          t->Emis->TauIn, 
          t->Emis->TauTot, 
          t->Emis->Pesc, 
          t->Emis->Pelec_esc, 
          t->Emis->Pdest, 
          t->Emis->pump, 
          t->Emis->ots, 
          t->Lo->Pop, 
          t->Hi->Pop ,
          t->Emis->PopOpc );
        return;
}


/*OccupationNumberLine - derive the photon occupation number at line center for any line */
double OccupationNumberLine( transition * t )
{
        double OccupationNumberLine_v;

        DEBUG_ENTRY( "OccupationNumberLine()" );

        ASSERT( t->ipCont > 0 );

        /* routine to evaluate line photon occupation number */
        if( t->Lo->Pop > SMALLFLOAT )
        {
                /* the lower population with correction for stimulated emission */
                double PopLo_corr = t->Lo->Pop / t->Lo->g - t->Hi->Pop / t->Hi->g;
                OccupationNumberLine_v = ( t->Hi->Pop / t->Hi->g )/SDIV(PopLo_corr ) *
                        (1. - t->Emis->Pesc);
        }
        else
        {
                OccupationNumberLine_v = 0.;
        }
        return( OccupationNumberLine_v );
}

/*TexcLine derive excitation temperature of line from contents of line array */
double TexcLine(transition * t)
{
        double TexcLine_v;

        DEBUG_ENTRY( "TexcLine()" );

        /* routine to evaluate line excitation temp using contents of line array
         * */
        if( t->Hi->Pop * t->Lo->Pop > 0. )
        {
                TexcLine_v = ( t->Hi->Pop / t->Hi->g )/( t->Lo->Pop / t->Lo->g );
                TexcLine_v = log(TexcLine_v);
                /* protect against infinite temp limit */
                if( fabs(TexcLine_v) > SMALLFLOAT )
                {
                        TexcLine_v = - t->EnergyK / TexcLine_v;
                }
        }
        else
        {
                TexcLine_v = 0.;
        }
        return( TexcLine_v );
}

/*eina convert a gf into an Einstein A */
double eina(double gf,
          double enercm, 
          double gup)
{
        double eina_v;

        DEBUG_ENTRY( "eina()" );

        /* derive the transition prob, given the following
         * call to function is gf, energy in cm^-1, g_up
         * gf is product of g and oscillator strength
         * eina = ( gf / 1.499e-8 ) / (wl/1e4)**2 / gup  */
        eina_v = (gf/gup)*TRANS_PROB_CONST*POW2(enercm);
        return( eina_v );
}

/*GetGF convert Einstein A into oscillator strength */
double GetGF(double trans_prob, 
          double enercm, 
          double gup)
{
        double GetGF_v;

        DEBUG_ENTRY( "GetGF()" );

        ASSERT( enercm > 0. );
        ASSERT( trans_prob > 0. );
        ASSERT( gup > 0.);

        /* derive the transition prob, given the following
         * call to function is gf, energy in cm^-1, g_up
         * gf is product of g and oscillator strength
         * trans_prob = ( GetGF/gup) / 1.499e-8 / ( 1e4/enercm )**2 */
        GetGF_v = trans_prob*gup/TRANS_PROB_CONST/POW2(enercm);
        return( GetGF_v );
}

/*abscf convert gf into absorption coefficient */
double abscf(double gf, 
          double enercm, 
          double gl)
{
        double abscf_v;

        DEBUG_ENTRY( "abscf()" );

        ASSERT(gl > 0. && enercm > 0. && gf > 0. );

        /* derive line absorption coefficient, given the following:
         * gf, enercm, g_low
         * gf is product of g and oscillator strength */
        abscf_v = 1.4974e-6*(gf/gl)*(1e4/enercm);
        return( abscf_v );
}

/*chIonLbl use information in line array to generate a null terminated ion label in "Fe 2" */
void chIonLbl(char *chIonLbl_v , transition * t )
{
        /*static char chIonLbl_v[5];*/

        DEBUG_ENTRY( "chIonLbl()" );

        /* function to use information within the line array
         * to generate an ion label, giving element and
         * ionization stage
         * */
        if( t->Hi->nelem < 0 )
        {
                /* this line is to be ignored */
                strcpy( chIonLbl_v, "Dumy" );
        }
        else if( t->Hi->nelem-1 >= LIMELM )
        {
                /* this is one of the molecules, either 12CO or 13CO */

                /* >>chng 02 may 15, from to chElementNameShort to go mole right */
                strcpy( chIonLbl_v , elementnames.chElementNameShort[t->Hi->nelem-1] );

                /* chIonStage is four char null terminated string, starting with "_1__" 
                strcat( chIonLbl_v, "CO");*/
        }

        else
        {
                ASSERT( t->Hi->nelem >= 1 );
                /* ElmntSym.chElementSym is null terminated, 2 ch + null, string giving very
                 * short form of element name */
                strcpy( chIonLbl_v , elementnames.chElementSym[t->Hi->nelem -1] );

                /* chIonStage is four char null terminated string, starting with "_1__" */
                strcat( chIonLbl_v, elementnames.chIonStage[t->Hi->IonStg-1]);
        }
        /* chIonLbl is four char null terminated string */
        return/*( chIonLbl_v )*/;
}

/*chLineLbl use information in line transfer arrays to generate a line label */
/* this label is null terminated */
/* ContCreatePointers has test this with full range of wavelengths */
char* chLineLbl(transition * t )
{
        static char chLineLbl_v[11];

        DEBUG_ENTRY( "chLineLbl()" );


        /* function to use information within the line array
         * to generate a line label, giving element and
         * ionization stage
         * rhs are set in large block data */

        /* NB this function is profoundly slow due to sprintf statement
         * also - it cannot be evaluated within a write statement itself*/

        if( t->WLAng > (realnum)INT_MAX )
        {
                sprintf( chLineLbl_v, "%2.2s%2.2s%5i%c", 
                        elementnames.chElementSym[t->Hi->nelem -1], 
                        elementnames.chIonStage[t->Hi->IonStg-1], 
                   (int)(t->WLAng/1e8), 'c' );
        }
        else if( t->WLAng > 9999999. )
        {
                /* wavelength is very large, convert to centimeters */
                sprintf( chLineLbl_v, "%2.2s%2.2s%5.2f%c", 
                        elementnames.chElementSym[t->Hi->nelem -1], 
                        elementnames.chIonStage[t->Hi->IonStg-1], 
                        t->WLAng/1e8, 'c' );
        }
        else if( t->WLAng > 999999. )
        {
                /* wavelength is very large, convert to microns */
                sprintf( chLineLbl_v, "%2.2s%2.2s%5i%c", 
                        elementnames.chElementSym[t->Hi->nelem -1], 
                        elementnames.chIonStage[t->Hi->IonStg-1], 
                        (int)(t->WLAng/1e4), 'm' );
        }
        else if( t->WLAng > 99999. )
        {
                /* wavelength is very large, convert to microns */
                sprintf( chLineLbl_v, "%2.2s%2.2s%5.1f%c", 
                        elementnames.chElementSym[t->Hi->nelem -1], 
                        elementnames.chIonStage[t->Hi->IonStg-1], 
                        t->WLAng/1e4, 'm' );
        }
        else if( t->WLAng > 9999. )
        {
                sprintf( chLineLbl_v, "%2.2s%2.2s%5.2f%c", 
                        elementnames.chElementSym[ t->Hi->nelem -1], 
                        elementnames.chIonStage[t->Hi->IonStg-1], 
                   t->WLAng/1e4, 'm' );
        }
        else if( t->WLAng >= 100. )
        {
                sprintf( chLineLbl_v, "%2.2s%2.2s%5i%c", 
                        elementnames.chElementSym[ t->Hi->nelem -1], 
                        elementnames.chIonStage[t->Hi->IonStg-1], 
                   (int)t->WLAng, 'A' );
        }
        /* the following two formats should be changed for the
         * wavelength to get more precision */
        else if( t->WLAng >= 10. )
        {
                sprintf( chLineLbl_v, "%2.2s%2.2s%5.1f%c", 
                        elementnames.chElementSym[ t->Hi->nelem -1], 
                        elementnames.chIonStage[t->Hi->IonStg-1], 
                   t->WLAng, 'A' );
        }
        else
        {
                sprintf( chLineLbl_v, "%2.2s%2.2s%5.2f%c", 
                        elementnames.chElementSym[ t->Hi->nelem -1], 
                        elementnames.chIonStage[t->Hi->IonStg-1], 
                   t->WLAng, 'A' );
        }
        /* make sure that string ends with null character - this should
         * be redundant */
        ASSERT( chLineLbl_v[10]==0 );
        return( chLineLbl_v );
}

/*RefIndex calculates the index of refraction of air using the line energy in wavenumbers,
 * used to convert vacuum wavelengths to air wavelengths. */
double RefIndex(double EnergyWN )
{
        double RefIndex_v, 
          WaveMic, 
          xl, 
          xn;

        DEBUG_ENTRY( "RefIndex()" );

        ASSERT( EnergyWN > 0. );

        /* the wavelength in microns */
        WaveMic = 1.e4/EnergyWN;

        /* only do index of refraction if longward of 2000A */
        if( WaveMic > 0.2 )
        {
                /* longward of 2000A
                 * xl is 1/WaveMic^2 */
                xl = 1.0/WaveMic/WaveMic;
                /* use a formula from 
                 *>>refer       air     index refraction        Allen, C.W. 1973, Astrophysical quantities, 
                 *>>refercon    3rd Edition (AQ), p.124 */
                xn = 255.4/(41. - xl);
                xn += 29498.1/(146. - xl);
                xn += 64.328;
                RefIndex_v = xn/1.e6 + 1.;
        }
        else
        {
                RefIndex_v = 1.;
        }
        ASSERT( RefIndex_v > 0. );
        return( RefIndex_v );
}

/*PutCS enter a collision strength into an individual line vector */
void PutCS(double cs, 
  transition * t)
{

        DEBUG_ENTRY( "PutCS()" );

        /* collision strength must not be negative, had been test for being positive,
         * but called with zero - did not check why 98 jul 5 */
        ASSERT( cs >= 0. );

        t->Coll.cs = (realnum)cs;
        return;
}

/*WavlenErrorGet - given the real wavelength in A for a line
 * routine will find the error expected between the real 
 * wavelength and the wavelength printed in the output, with 4 sig figs,
 * function returns difference between exact and 4 sig fig wl, so 
 * we have found correct line is fabs(d wl) < return */
realnum WavlenErrorGet( realnum wavelength )
{
        double a;
        realnum errorwave;

        DEBUG_ENTRY( "WavlenErrorGet()" );

        ASSERT( LineSave.sig_figs <= 6 );

        if( wavelength > 0. )
        {
                /* normal case, positive (non zero) wavelength */
                a = log10( wavelength+FLT_EPSILON);
                a = floor(a);
        }
        else
        {
                /* might be called with wl of zero, this is that case */
                /* errorwave = 1e-4f; */
                a = 0.;
        }

        errorwave = 5.f * (realnum)pow( 10., a - (double)LineSave.sig_figs );
        return errorwave;
}

/*linadd enter lines into the line storage array, called once per zone for each line*/
void linadd(
  double xInten,        /* xInten - local emissivity per unit vol, no fill fac */
  realnum wavelength,   /* realnum wavelength */
  const char *chLab,/* string label for ion */
  char chInfo,          /* character type of entry for line - given below */
                                        /* 'c' cooling, 'h' heating, 'i' info only, 'r' recom line */
  const char *chComment )
{
        DEBUG_ENTRY( "linadd()" );

        /* main routine to actually enter lines into the line storage array
         * called at top level within routine lines
         * called series of times in routine PutLine for lines transferred
         */

        /* three values, -1 is just counting, 0 if init, 1 for calculation */
        if( LineSave.ipass > 0 )
        {
                /* not first pass, sum lines only
                 * total emission from vol */
                LineSv[LineSave.nsum].sumlin[0] += xInten*radius.dVeff;
                /* local emissivity in line */
                LineSv[LineSave.nsum].emslin[0] = xInten;
                /* no need to increment or set [1] version since this is called with no continuum
                 * index, don't know what to do */
                if( wavelength > 0 )
                {
                        /* only put informational lines, like "Q(H) 4861", in this stack
                         * many continua have a wavelength of zero and are proper intensities,
                         * it would be wrong to predict their transferred intensity */
                        LineSv[LineSave.nsum].emslin[1] = LineSv[LineSave.nsum].emslin[0];
                        LineSv[LineSave.nsum].sumlin[1] = LineSv[LineSave.nsum].sumlin[0]; 
                }
        }

        else if( LineSave.ipass == 0 )
        {
                /* first call to stuff lines in array, confirm that label is one of
                 * the four correct ones */
                /* >>chng 04 apr 24, last two had been != so this test never really happened; PvH */
                ASSERT( (chInfo == 'c') || (chInfo == 'h') || (chInfo == 'i') || (chInfo == 'r' ) );

                /* then save it into array */
                LineSv[LineSave.nsum].chSumTyp = (char)chInfo;

                /* number of lines ok, set parameters for first pass */
                LineSv[LineSave.nsum].sumlin[0] = 0.;
                LineSv[LineSave.nsum].emslin[0] = 0.;
                LineSv[LineSave.nsum].wavelength = wavelength;
                LineSv[LineSave.nsum].sumlin[1] = 0.;
                LineSv[LineSave.nsum].emslin[1] = 0.;
                LineSv[LineSave.nsum].chComment = chComment;
                /* check that null is correct, string overruns have 
                 * been a problem in the past */
                ASSERT( strlen( chLab )<5 );
                strcpy( LineSv[LineSave.nsum].chALab, chLab );
        }

        /* increment the line counter */
        ++LineSave.nsum;

        /* routine can be called with negative LineSave.ipass, in this case
         * we are just counting number of lines for current setup */
        return;
}


/* this is energy in Ryd as set by call to PntForLine */
static double EnergyRyd;

/*emergent_line find absorption due to continuous opacity */
double emergent_line( 
        /* emissivity [erg cm-3 s-1] in inward direction */
        double emissivity_in , 
        /* emissivity [erg cm-3 s-1] in outward direction */
        double emissivity_out , 
        /* array index for continuum frequency */
        long int ipCont )
{

        double emergent_in , emergent_out;
        long int i = ipCont-1;

        DEBUG_ENTRY( "emergent_line()" );

        ASSERT( i >= 0 && i < rfield.nupper-1 );

        /* do nothing if first iteration since we do not know the outward-looking
         * optical depths.  In version C07.02.00 we assumed an infinite optical
         * depth in the outward direction, which would be appropriate for a 
         * HII region on the surface of a molecular cloud.  This converged onto
         * the correct solution in later iterations, but on the first iteration
         * this underestimated total emission if the infinite cloud were not
         * present.  With C07.02.01 we make no assuptions about what is in the
         * outward direction and simply use the local emission. 
         * Behavior is unchanged on later iterations */
        if( iteration == 1 )
        {
                /* first iteration - do not know outer optical depths so assume very large 
                 * optical depths */
                emergent_in = emissivity_in*opac.E2TauAbsFace[i];
                emergent_out = emissivity_out;
        }
        else
        {
                if( geometry.lgSphere )
                {
                        /* second or later iteration in closed or spherical geometry */
                        /* inwardly directed emission must get to central hole then across entire
                         * far side of shell */
                        emergent_in = emissivity_in  * opac.E2TauAbsFace[i] *opac.E2TauAbsTotal[i];

                        /* E2 is outwardly directed emission to get to outer edge of cloud */
                        emergent_out = emissivity_out * opac.E2TauAbsOut[i];
                }
                else
                {
                        /* open geometry in second or later iteration, outer optical depths are known 
                         * this is light emitted into the outer direction and backscattered
                         * into the inner */
                        double reflected = emissivity_out * opac.albedo[i] * (1.-opac.E2TauAbsOut[i]);
                        /* E2 is to get to central hole */
                        emergent_in = (emissivity_in + reflected) * opac.E2TauAbsFace[i];
                        /* E2 is to get to outer edge */
                        emergent_out = (emissivity_out - reflected) * opac.E2TauAbsOut[i];
                }
        }
        /* return the net emissivity times a vol element */
        return( (emergent_in + emergent_out)  );
}

/*lindst add line with destruction and outward */
void lindst(
  // xInten - local emissivity per unit vol
  double xInten, 
  // wavelength of line in Angstroms
  realnum wavelength, 
  // *chLab string label for ion
  const char *chLab, 
  // ipnt offset of line in continuum mesh
  long int ipnt, 
  // chInfo character type of entry for line - 'c' cooling, 'h' heating, 'i' info only, 'r' recom line
  char chInfo, 
  // lgOutToo should line be included in outward beam?
  bool lgOutToo,
  // *chComment string explaining line 
  const char *chComment )
{
        double saveemis;
        DEBUG_ENTRY( "lindst()" );

        /* >>chng 06 feb 08, add test on xInten positive, no need to evaluate
         * for majority of zero */
        if( LineSave.ipass > 0 && xInten > 0. )
        {
                /* LineSave.ipass > 0, integration across simulation, sum lines only 
                 * emissivity, emission per unit vol, for this zone */
                LineSv[LineSave.nsum].emslin[0] = xInten;
                /* integrated intensity or luminosity, the emissivity times the volume */
                LineSv[LineSave.nsum].sumlin[0] += xInten*radius.dVeff;

                /* add line to outward beam 
                 * there are lots of lines that are sums of other lines, or
                 * just for info of some sort.  These have flag lgOutToo false.
                 * Note that the EnergyRyd variable only has a rational
                 * value if PntForLine was called just before this routine - in all
                 * cases where this did not happen the flag is false. */
                if( lgOutToo )
                {
                        rfield.outlin[ipnt-1] += (realnum)(xInten/(rfield.anu[ipnt-1]*EN1RYD)*
                          radius.dVolOutwrd*opac.ExpZone[ipnt-1]);
                        rfield.reflin[ipnt-1] += (realnum)(xInten/(rfield.anu[ipnt-1]*EN1RYD)*
                          radius.dVolReflec);
                }
                /* main production pass, sum lines only */
                if( ipnt <= rfield.nflux )
                {
                        /* emergent_line accounts for destruction by absorption outside
                         * the line-forming region */
                        saveemis = emergent_line( 
                                xInten*rt.fracin , xInten*(1.-rt.fracin) , ipnt );
                        LineSv[LineSave.nsum].emslin[1] = saveemis;
                        LineSv[LineSave.nsum].sumlin[1] += saveemis*radius.dVeff; 
                }
        }
        else if( LineSave.ipass == 0 )
        {
                ASSERT( ipnt > 0 );

                /* first call to stuff lines in array, confirm that label is one of
                 * the four correct ones */
                /* >>chng 04 apr 24, last two had been != so this test never really happened; PvH */
                ASSERT( (chInfo == 'c') || (chInfo == 'h') || (chInfo == 'i') || (chInfo == 'r' ) );
                LineSv[LineSave.nsum].chSumTyp = (char)chInfo;
                /* number of lines ok, set parameters for first pass */
                LineSv[LineSave.nsum].sumlin[0] = 0.;
                LineSv[LineSave.nsum].emslin[0] = 0.;
                LineSv[LineSave.nsum].wavelength = fabs(wavelength);
                LineSv[LineSave.nsum].emslin[1] = 0.;
                LineSv[LineSave.nsum].sumlin[1] = 0.;
                LineSv[LineSave.nsum].chComment = chComment;

                /* check that null is correct, string overruns have 
                * been a problem in the past */
                ASSERT( strlen( chLab )<5 );
                strcpy( LineSv[LineSave.nsum].chALab, chLab );

                // check that line wavelength and continuum index agree to some extent
                // this check cannot be very precise because some lines have 
                // "wavelengths" that are set by common usage rather than the correct
                // wavelength derived from energy and index of refraction of air
                double error = MAX2(0.1*rfield.AnuOrg[ipnt-1] , rfield.widflx[ipnt-1] );
                ASSERT( wavelength<=0 ||
                        fabs( rfield.AnuOrg[ipnt-1] - RYDLAM / wavelength) < error );
        }

        /* increment the line counter */
        ++LineSave.nsum;

        EnergyRyd = 0.;
        return;
}

/*PntForLine generate pointer for forbidden line */
void PntForLine(
  /* wavelength of transition in Angstroms */
  double wavelength, 
  /* label for this line */
  const char *chLabel,
  /* this is array index on the f, not c scale,
   * for the continuum cell holding the line */
  long int *ipnt)
{
        /* 
         * maximum number of forbidden lines - this is a good bet since
         * new lines do not go into this group, and lines are slowly 
         * moving to level 1 
         */
#       define  MAXFORLIN       1000
        static long int ipForLin[MAXFORLIN]={0};

        /* number of forbidden lines entered into continuum array */
        static long int nForLin;

        DEBUG_ENTRY( "PntForLine()" );

        /* must be 0 or greater */
        ASSERT( wavelength >= 0. );

        if( wavelength == 0. )
        {
                /* zero is special flag to initialize */
                nForLin = 0;
                /* ipLineEnergy will only put in line label if nothing already there */
                EnergyRyd = 0.;
        }
        else
        {

                if( LineSave.ipass > 0 )
                {
                        /* not first pass, sum lines only */
                        *ipnt = ipForLin[nForLin];
                }
                else if( LineSave.ipass == 0 )
                {
                        /* check if number of lines in arrays exceeded */
                        if( nForLin >= MAXFORLIN )
                        {
                                fprintf( ioQQQ, "PROBLEM %5ld lines is too many for PntForLine.\n", 
                                  nForLin );
                                fprintf( ioQQQ, " Increase the value of maxForLine everywhere in the code.\n" );
                                cdEXIT(EXIT_FAILURE);
                        }

                        /* ipLineEnergy will only put in line label if nothing already there */
                        EnergyRyd = RYDLAM/wavelength;
                        ipForLin[nForLin] = ipLineEnergy(EnergyRyd,chLabel , 0);
                        *ipnt = ipForLin[nForLin];
                }
                else
                {
                        /* this is case where we are only counting lines */
                        *ipnt = 0;
                }
                ++nForLin;
        }
        return;
}

static realnum ExtraInten;

/*PutLine enter local line intensity into the intensity stack for eventual printout */
void PutLine(transition * t, const char *chComment, const char *chLabelTemp)
{
        char chLabel[5];
        double xIntensity,
                other,
                xIntensity_in;

        DEBUG_ENTRY( "PutLine()" );

        /* routine to use line array data to generate input
         * for emission line array */
        ASSERT( t->ipCont > 0. );

        strcpy( chLabel, chLabelTemp );
        chLabel[4] = '\0';

        /* if ipass=0 then we must generate label info since first pass
         * gt.0 then only need line intensity data */
        if( LineSave.ipass == 0 )
        {
                xIntensity = 0.;
        }
        else
        {
                /* both the counting and integrating modes comes here */
                /* not actually used so set to impossible values */
                chLabel[0] = '\n';
                /* total line intensity or luminosity 
                 * these may not be defined in initial calls so define here */
                xIntensity = t->Emis->xIntensity + ExtraInten;
        }
        /* initial counting case, where ipass == -1, just ignored above, call linadd below */

        /* ExtraInten is option that allows extra intensity (i.e., recomb)
         * to be added to this line  with Call PutExtra( exta )
         * in main lines this extra
         * contribution must be identified explicitly */
        ExtraInten = 0.;
        /*linadd(xIntensity,wl,chLabel,'i');*/
        /*lindst add line with destruction and outward */
        rt.fracin = t->Emis->FracInwd;
        lindst(xIntensity, 
          t->WLAng, 
          chLabel, 
          t->ipCont, 
          /* this is information only - has been counted in cooling already */
          'i', 
          /* do not add to outward beam, also done separately */
          false,
          chComment);
        rt.fracin = 0.5;

        /* inward part of line - do not move this away from previous lines
         * since xIntensity is used here */
        xIntensity_in = xIntensity*t->Emis->FracInwd;
        linadd(xIntensity_in,t->WLAng,"Inwd",'i',chComment);

        /* cooling part of line */
        other = t->Coll.cool;
        linadd(other,t->WLAng,"Coll",'i',chComment);

        /* fluorescent excited part of line */
        other = t->Emis->PopOpc * t->Emis->pump * (1.-t->Emis->ColOvTot) * t->EnergyErg;
        linadd(other,t->WLAng,"Pump",'i',chComment);

        /* heating part of line */
        other = t->Coll.heat;
        linadd(other,t->WLAng,"Heat",'i',chComment);
        return;
}

/*PutLine enter local line intensity into the intensity stack for eventual printout */
void PutLine(transition * t,
  const char *chComment)
{
        char chLabel[5];
        double xIntensity,
                other,
                xIntensity_in;

        DEBUG_ENTRY( "PutLine()" );

        /* routine to use line array data to generate input
         * for emission line array */
        ASSERT( t->ipCont > 0. );

        /* if ipass=0 then we must generate label info since first pass
         * gt.0 then only need line intensity data */
        if( LineSave.ipass == 0 )
        {
                /* these variables not used by linadd if ipass ne 0 */
                /* chIonLbl is function that generates a null terminated 4 char string, of form "C  2" */
                chIonLbl(chLabel, t);

                xIntensity = 0.;
        }
        else
        {
                /* both the counting and integrating modes comes here */
                /* not actually used so set to impossible values */
                /*strcpy( chLabel, "    " );*/
                chLabel[0] = '\n';
                /* total line intensity or luminosity 
                 * these may not be defined in initial calls so define here */
                xIntensity = t->Emis->xIntensity + ExtraInten;
        }
        /* initial counting case, where ipass == -1, just ignored above, call linadd below */

        /* ExtraInten is option that allows extra intensity (i.e., recomb)
         * to be added to this line  with Call PutExtra( exta )
         * in main lines this extra
         * contribution must be identified explicitly */
        ExtraInten = 0.;
        /*lindst add line with destruction and outward */
        rt.fracin = t->Emis->FracInwd;
        lindst(xIntensity, 
          t->WLAng, 
          chLabel, 
          t->ipCont, 
          /* this is information only - has been counted in cooling already */
          'i', 
          /* do not add to outward beam, also done separately */
          false,
          chComment);
        rt.fracin = 0.5;

        /* inward part of line - do not move this away from previous lines
         * since xIntensity is used here */
        xIntensity_in = xIntensity*t->Emis->FracInwd;
        linadd(xIntensity_in,t->WLAng,"Inwd",'i',chComment);

        /* cooling part of line */
        other = t->Coll.cool;
        linadd(other,t->WLAng,"Coll",'i',chComment);

        /* fluorescent excited part of line */
        other = t->Emis->PopOpc * t->Emis->pump * (1.-t->Emis->ColOvTot) * t->EnergyErg;
        linadd(other,t->WLAng,"Pump",'i',chComment);

        /* heating part of line */
        other = t->Coll.heat;
        linadd(other,t->WLAng,"Heat",'i',chComment);
        return;
}

/* ==================================================================== */
/*PutExtra enter and 'extra' intensity source for some line */
void PutExtra(double Extra)
{

        DEBUG_ENTRY( "PutExtra()" );

        ExtraInten = (realnum)Extra;
        return;
}

void TransitionJunk( transition *t )
{

        DEBUG_ENTRY( "TransitionJunk()" );

                /* wavelength, usually in A, used for printout */
        t->WLAng = -FLT_MAX;

                /* transition energy in degrees kelvin*/
        t->EnergyK = -FLT_MAX;
                /* transition energy in ergs */
        t->EnergyErg = -FLT_MAX;
                /* transition energy in wavenumbers */
        t->EnergyWN = -FLT_MAX;

        /* array offset for radiative transition within continuum array 
         * is negative if transition is non-radiative. */
        t->ipCont = -10000;

        CollisionJunk( &t->Coll );

        /* set these equal to NULL first. That will cause the code to crash if
         * the variables are ever used before being deliberately set. */
        t->Emis = &DummyEmis;
        
        t->Lo = NULL;
        t->Hi = NULL;
        return;
}


/*EmLineJunk set all elements of transition struc to dangerous values */
void EmLineJunk( emission *t )
{

        DEBUG_ENTRY( "EmLineJunk()" );

        /* optical depth in continuum to ill face */
        t->TauCon = -FLT_MAX;

        /* inward and total line optical depths */
        t->TauIn = -FLT_MAX;
        t->TauTot = -FLT_MAX;

        /* type of redistribution function, */
        t->iRedisFun = INT_MIN;

        /* array offset for line within fine opacity */
        t->ipFine = -10000;

        /* inward fraction */
        t->FracInwd = -FLT_MAX;

        /* continuum pumping rate */
        t->pump = -FLT_MAX;

        /* line intensity */
        t->xIntensity = -FLT_MAX;

        /* number of photons emitted per sec in the line */
        t->phots = -FLT_MAX;

        /* gf value */
        t->gf = -FLT_MAX;

        /* escape and destruction probs */
        t->Pesc = -FLT_MAX;
        t->Pdest = -FLT_MAX;
        t->Pelec_esc = -FLT_MAX;

        /* damping constant, and number related to it */
        t->dampXvel = -FLT_MAX;
        t->damp = -FLT_MAX;

        /* ratio of collisional to radiative excitation*/
        t->ColOvTot = -FLT_MAX;

        /* line opacity */
        t->opacity = -FLT_MAX;

        t->PopOpc = -FLT_MAX;

        /* transition prob, Einstein A upper to lower */
        t->Aul = -FLT_MAX;

        /* ots rate */
        t->ots = -FLT_MAX;
        return;
}

/*CollisionJunk set all elements of transition struc to dangerous values */
void CollisionJunk( collision * t )
{

        DEBUG_ENTRY( "CollisionJunk()" );

        t->ColUL = -FLT_MAX;

        /* Coll->cooling and Coll->heating due to collisional excitation */
        t->cool = -FLT_MAX;
        t->heat = -FLT_MAX;

        /* collision strengths for transition */
        t->cs = -FLT_MAX;
        for( long i=0; i<ipNCOLLIDER; i++ )
                t->csi[i] = -FLT_MAX;
        return;
}

/*StateJunk set all elements of transition struc to dangerous values */
void StateJunk( quantumState * t )
{

        DEBUG_ENTRY( "StateJunk()" );

        /* t->chLabel = ''; */

        t->g = -FLT_MAX;

        t->Pop = -FLT_MAX;

        t->IonStg = -10000;

        t->nelem = -10000;
        return;
}

/*TransitionZero zeros out TransitionZero */
void TransitionZero( transition *t )
{

        DEBUG_ENTRY( "TransitionZero()" );

        CollisionZero( &t->Coll );

        StateZero( t->Lo );
        StateZero( t->Hi );
        EmLineZero( t->Emis );

        return;
}

/*EmLineZero zeros out the emission line structure */
void EmLineZero( emission * t )
{

        DEBUG_ENTRY( "EmLineZero()" );

        /* total optical depth in all overlapping lines to illuminated face,
         * used for pumping */
        t->TauCon = opac.taumin;

        /* inward and total line optical depths */
        /* >>chng 03 feb 14, from 0 to opac.taumin */
        t->TauIn = opac.taumin;

        /* total optical depths */
        t->TauTot = 1e20f;

        /* inward fraction */
        /* >>chng 03 feb 14, from 0 to 1 */
        t->FracInwd = 1.;

        /* continuum pumping rate */
        t->pump = 0.;

        /* line intensity */
        t->xIntensity = 0.;

        /* number of photons emitted per sec in the line */
        t->phots = 0.;

        /* escape and destruction probs */
        /* >>chng 03 feb 14, change from 0 to 1 */
        t->Pesc = 1.;
        t->Pdest = 0.;
        t->Pelec_esc = 0.;

        /* ratio of collisional to radiative excitation*/
        t->ColOvTot = 0.;

        /* pop that enters net opacity */
        t->PopOpc = 0.;

        /* ots rate */
        t->ots = 0.;
        return;
}

/*CollisionZero zeros out the structure */
void CollisionZero( collision * t )
{

        DEBUG_ENTRY( "CollisionZero()" );

        /* Coll->cooling and Coll->heating due to collisional excitation */
        t->cool = 0.;
        t->heat = 0.;
        return;
}

/*StateZero zeros out the structure */
void StateZero( quantumState * t )
{

        DEBUG_ENTRY( "StateZero()" );

        t->Pop = 0.;
        return;
}

/*LineConvRate2CS convert down coll rate back into electron cs in case other parts of code need this for reference */
void LineConvRate2CS( transition* t , realnum rate )
{

        DEBUG_ENTRY( "LineConvRate2CS()" );

        /* return is collision strength, convert from collision rate from 
         * upper to lower, this assumes pure electron collisions, but that will
         * also be assumed by anything that uses cs, for self-consistency */
        t->Coll.cs = rate * t->Hi->g / (realnum)dense.cdsqte;

        /* change assert to non-negative - there can be cases (Iin H2) where cs has
         * underflowed to 0 on some platforms */
        ASSERT( t->Coll.cs >= 0. );
        return;
}

/*ConvRate2CS convert down coll rate back into electron cs in case other parts of code need this for reference */
double ConvRate2CS( realnum gHi , realnum rate )
{

        double cs;

        DEBUG_ENTRY( "ConvRate2CS()" );

        /* return is collision strength, convert from collision rate from 
         * upper to lower, this assumes pure electron collisions, but that will
         * also be assumed by anything that uses cs, for self-consistency */
        cs = rate * gHi / dense.cdsqte;

        /* change assert to non-negative - there can be cases (Iin H2) where cs has
         * underflowed to 0 on some platforms */
        ASSERT( cs >= 0. );
        return cs;
}

/* returns true is we have not overrun optical depth scale */
bool lgTauGood( transition * t)
{
        bool lgGoodTau;

        DEBUG_ENTRY( "lgTauGood()" );

        if( (t->Emis->TauTot*0.9 - t->Emis->TauIn < 0. && t->Emis->TauIn > 0.) && 
            ! fp_equal( t->Emis->TauTot, opac.taumin ) )
        {
                /* do not do anything if we have overrun the scale, */
                lgGoodTau = false;
        }
        else
        {
                lgGoodTau = true;
        }
        return lgGoodTau;
}

/*gbar0 compute g-bar gaunt factor for neutrals */
STATIC void gbar0(double ex, 
          realnum *g)
{
        double a, 
          b, 
          c, 
          d, 
          y;

        DEBUG_ENTRY( "gbar0()" );

        /* written by Dima Verner
         *
         * Calculation of the effective Gaunt-factor by use of 
         * >>refer      gbar    cs      Van Regemorter, H., 1962, ApJ 136, 906
         * fits for neutrals
         *  Input parameters: 
         * ex - energy ryd - now K
         * t  - temperature in K
         *  Output parameter:
         * g  - effective Gaunt factor
         * */

        /* y = ex*157813.7/te */
        y = ex/phycon.te;
        if( y < 0.01 )
        {
                *g = (realnum)(0.29*(log(1.0+1.0/y) - 0.4/POW2(y + 1.0))/exp(y));
        }
        else
        {
                if( y > 10.0 )
                {
                        *g = (realnum)(0.066/sqrt(y));
                }
                else
                {
                        a = 1.5819068e-02;
                        b = 1.3018207e00;
                        c = 2.6896230e-03;
                        d = 2.5486007e00;
                        d = log(y/c)/d;
                        *g = (realnum)(a + b*exp(-0.5*POW2(d)));
                }
        }
        return;
}

/*gbar1 compute g-bar collision strength using Mewe approximations */
STATIC void gbar1(double ex, 
          realnum *g)
{
        double y;

        DEBUG_ENTRY( "gbar1()" );

        /*      *written by Dima Verner
         *
         * Calculation of the effective Gaunt-factor by use of 
         * >>refer      gbar    cs      Mewe,R., 1972, A&A 20, 215
         * fits for permitted transitions in ions MgII, CaII, FeII (delta n = 0)
         * Input parameters: 
         * ex - excitation energy in Ryd - now K
         * te  - temperature in K
         * Output parameter:
         * g  - effective Gaunt factor
         */

        /* y = ex*157813.7/te */
        y = ex/phycon.te;
        *g = (realnum)(0.6 + 0.28*(log(1.0+1.0/y) - 0.4/POW2(y + 1.0)));
        return;
}

/*MakeCS compute collision strength by g-bar approximations */
void MakeCS(transition * t)
{
        long int ion;
        double Abun, 
          cs;
        realnum
          gbar;

        DEBUG_ENTRY( "MakeCS()" );

        /* routine to get cs from various approximations */

        /* check if abundance greater than 0 */
        ion = t->Hi->IonStg;

        Abun = dense.xIonDense[ t->Hi->nelem -1 ][ ion-1 ];
        if( Abun <= 0. )
        {
                gbar = 1.;
        }
        else
        {
                /* check if neutral or ion */
                if( ion == 1 )
                {
                        /* neutral - compute gbar for eventual cs */
                        gbar0(t->EnergyK,&gbar);
                }
                else
                {
                        /* ion - compute gbar for eventual cs */
                        gbar1(t->EnergyK,&gbar);
                }
        }

        /* above was g-bar, convert to cs */
        cs = gbar*(14.5104/WAVNRYD)*t->Emis->gf/t->EnergyWN;

        /* stuff the cs in place */
        t->Coll.cs = (realnum)cs;
        return;
}

/*totlin sum total intensity of cooling, recombination, or intensity lines */
double totlin(
        /* chInfor is 1 char, 
        'i' information, 
        'r' recombination or 
        'c' collision */
        int chInfo)
{
        long int i;
        double totlin_v;

        DEBUG_ENTRY( "totlin()" );

        /* routine goes through set of entered line
         * intensities and picks out those which have
         * types agreeing with chInfo.  Valid types are
         * 'c', 'r', and 'i'
         *begin sanity check */
        if( (chInfo != 'i' && chInfo != 'r') && chInfo != 'c' )
        {
                fprintf( ioQQQ, " TOTLIN does not understand chInfo=%c\n", 
                  chInfo );
                cdEXIT(EXIT_FAILURE);
        }
        /*end sanity check */

        /* now find sum of lines of type chInfo */
        totlin_v = 0.;
        for( i=0; i < LineSave.nsum; i++ )
        {
                if( LineSv[i].chSumTyp == chInfo )
                {
                        totlin_v += LineSv[i].sumlin[0];
                }
        }
        return( totlin_v );
}


/*FndLineHt search through line heat arrays to find the strongest heat source */
void FndLineHt(long int *level, 
  /* this is the index of the strongest line in the array on the c scale */
  long int *ipStrong, 
  double *Strong)
{
        long int i; 

        DEBUG_ENTRY( "FndLineHt()" );

        *Strong = 0.;
        *level = 0;

        /* do the level 1 lines, 0 is dummy line, <=nLevel1 is correct for c scale */
        for( i=1; i <= nLevel1; i++ )
        {
                /* check if a line was the major heat agent */
                if( TauLines[i].Coll.heat > *Strong )
                {
                        *ipStrong = i;
                        *level = 1;
                        *Strong = TauLines[i].Coll.heat;
                }
        }

        /* now do the level 2 lines */
        for( i=0; i < nWindLine; i++ )
        {
                if( TauLine2[i].Hi->IonStg < TauLine2[i].Hi->nelem+1-NISO )
                {
                        /* check if a line was the major heat agent */
                        if( TauLine2[i].Coll.heat > *Strong )
                        {
                                *ipStrong = i;
                                *level = 2;
                                *Strong = TauLine2[i].Coll.heat;
                        }
                }
        }

        /* now do co carbon monoxide lines */
        for( i=0; i < nCORotate; i++ )
        {
                /* check if a line was the major heat agent */
                if( C12O16Rotate[i].Coll.heat > *Strong )
                {
                        *ipStrong = i;
                        *level = 3;
                        *Strong = C12O16Rotate[i].Coll.heat;
                }
        }
        for( i=0; i < nCORotate; i++ )
        {
                /* check if a line was the major heat agent */
                if( C13O16Rotate[i].Coll.heat > *Strong )
                {
                        *ipStrong = i;
                        *level = 4;
                        *Strong = C13O16Rotate[i].Coll.heat;
                }
        }

        /* now do the hyperfine structure lines */
        for( i=0; i < nHFLines; i++ )
        {
                /* check if a line was the major heat agent */
                if( HFLines[i].Coll.heat > *Strong )
                {
                        *ipStrong = i;
                        *level = 5;
                        *Strong = HFLines[i].Coll.heat;
                }
        }

        /*Chianti & Leiden Atoms & Molecules */
        for( i=0; i <linesAdded2; i++)
        {
                /* check if a line was the major heat agent */
                if( atmolEmis[i].tran->Coll.heat > *Strong )
                {
                        *ipStrong = i;
                        *level = 6;
                        *Strong = atmolEmis[i].tran->Coll.heat;
                }
        }
        return;
}

quantumState *AddState2Stack( void )
{
        DEBUG_ENTRY( "AddState2Stack()" );

        ASSERT( !lgStatesAdded );

        currentState = new quantumState;

        StateJunk( currentState );

        if( statesAdded == 0 )
        {
                GenericStates = currentState;
                GenericStates->next = NULL;
                lastState = GenericStates;
        }
        else
        {
                StateZero( currentState );
                lastState->next = currentState;
                lastState = lastState->next;
        }

        statesAdded++;

        return currentState;
}

emission *AddLine2Stack( bool lgRadiativeTrans )
{
        DEBUG_ENTRY( "AddLine2Stack()" );

        if( !lgRadiativeTrans )
        {
                return &DummyEmis;
        }
        else
        {
                ASSERT( lgLinesAdded == false );

                currentLine = new emission;

                EmLineJunk( currentLine );

                if( linesAdded == 0 )
                {
                        GenericLines = currentLine;
                        GenericLines->next = NULL;
                        lastLine = GenericLines;
                }
                else
                {
                        /* 
                        \todo 2 Does doing EmLineZero here defeat the purpose of EmLineJunk? 
                        * maybe we should pass full set of Emis components, fill everything in 
                        * here, and THEN use EmLineZero?  */
                        EmLineZero( currentLine );

                        lastLine->next = currentLine;
                        lastLine = lastLine->next;
                }

                linesAdded++;
                return currentLine;
        }
}

---