/home66/gary/public_html/cloudy/c08_branch/source/rt_ots.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 */
/*RT_OTS compute diffuse fields due to H, He atoms, ion, triplets, metal recombination,
 * called by ConvBase  */
/*RT_OTS_AddLine add local destruction of lines to ots field */
/*RT_OTS_AddCont add local destruction of continuum to ots field */
/*RT_OTS_Update sum flux, otscon, otslin, ConInterOut, outlin, to form SummeDif, SummedCon SummedOcc */
/*RT_OTS_Zero - zero out some vectors - 
 * this is only called when code initialized by ContSetIntensity */
/*RT_OTS_ChkSum sanity check confirms summed continua reflect contents of individuals */
#include "cddefines.h"
#include "taulines.h"
#include "opacity.h"
#include "dense.h"
#include "iso.h"
#include "hmi.h"
#include "h2.h"
#include "rfield.h"
#include "conv.h"
#include "rt.h"
#include "mole_co_atom.h"
#include "atomfeii.h"
#include "heavy.h"
#include "he.h"
#include "trace.h"

/* this flag may be used for debugging ots rate changes */
static int nOTS_Line_type = -1;
static int nOTS1=-1 , nOTS2=-1;
/*add local destruction of continuum to ots field */
STATIC void RT_OTS_AddCont(
        /* the ots rate itself */
        realnum ots, 
        /* pointer to continuum cell for ots, on f scale */
        long int ip);

/* =================================================================== */

void RT_OTS(void)
{
        long int
                ipla,
                ipISO ,
                nelem,
                n;
        realnum 
                difflya,
                esc,
                ots;

        /* the Bowen HeII yield 
         * >>chng 06 aug 08, from 0.3 to 0.4, update with netzer */
#       define BOWEN 0.5f
        long int ipHi, 
          ipLo;

        double bwnfac;
        double ots660;
        realnum cont_phot_destroyed;
        double save_lya_dest,
          save_he2lya_dest;

        double save_he2rec_dest;

        /*static long int nCall=0;
        ++nCall;
        fprintf(ioQQQ,"debugggtos enter %li\n", nCall );*/

        DEBUG_ENTRY( "RT_OTS()" );

        /**************************************************************************
         *
         * the bowen HeII - OIII fluorescense problem
         *
         **************************************************************************/
        nOTS_Line_type = 0;
        nelem = ipHELIUM;
        if( dense.lgElmtOn[nelem] )
        {
                /* conversion per unit atom to OIII, at end of sub we divide by it,
                 * to fix lines back to proper esc/dest probs */
                bwnfac = BOWEN * MAX2(0.f,1.f- Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Pesc - 
                        Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Pelec_esc );

                /* the last factor accounts for fact that two photons are produced,
                 * and the branching ratio */
                ots660 = Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Aul*
                  StatesElem[ipH_LIKE][nelem][ipH2p].Pop*dense.xIonDense[nelem][nelem+1]*
                  /*>>chng 06 aug 08, rm 0.8 factor from below, renorm aft discuss with Netzer */
                  Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Pdest *BOWEN*2.0;

                /* now add this to the ots field */
                if( ots660 > SMALLFLOAT ) 
                        RT_OTS_AddLine(ots660 , he.ip660 );

                /* decrease the destruction prob by the amount we will add elsewhere,
                 * ok since dest probs updated on every iteration*/
                Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Pdest *= (realnum)bwnfac;
                ASSERT( Transitions[ipH_LIKE][nelem][ipH2p][ipH1s].Emis->Pdest >= 0. );
                {
                        /* debugging code for line oscillation problems 
                                * most often Lya OTS oscillations*/
                        /*@-redef@*/
                        enum {DEBUG_LOC=false};
                        /*@+redef@*/
                        if( DEBUG_LOC )
                        {
                                fprintf(ioQQQ,"DEBUG HeII Bowen nzone %li bwnfac:%.2e bwnfac esc:%.2e ots660 %.2e\n",
                                        nzone,
                                        bwnfac , 
                                        bwnfac/BOWEN , 
                                        ots660 );
                        }
                }
        }

        else
        {
                bwnfac = 1.;
        }

        /* save Lya loss rates so we can reset at end */
        save_lya_dest = Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pdest;

        /* >>chng 03 may 28, hydro.dstfe2lya should not go into the ots field since it was
         * actually lost in exciting FeII 
         * this is not ready to be used - need to recover ots with small feii atom
        fprintf(ioQQQ,"feii debug\t%.3e\t%.3e\n", 
                Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest , hydro.dstfe2lya);
        Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest = MAX2(0.f,
                Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Pdest-hydro.dstfe2lya);*/

        /* this is option to kill Lya ots rates, 
         * rfield.lgLyaOTS is usually true (1), and set false (0) with
         * no lya ots command */
        Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pdest *= rfield.lgLyaOTS;

        /* option to kill heii lya and rec continua ots */
        save_he2lya_dest = Transitions[ipH_LIKE][ipHELIUM][ipH2p][ipH1s].Emis->Pdest;
        Transitions[ipH_LIKE][ipHELIUM][ipH2p][ipH1s].Emis->Pdest *= rfield.lgHeIIOTS;

        /* option to kill heii lya and rec continua ots */
        save_he2rec_dest = iso.RadRecomb[ipH_LIKE][ipHELIUM][ipH1s][ipRecRad];
        iso.RadRecomb[ipH_LIKE][ipHELIUM][ipH1s][ipRecRad] *= rfield.lgHeIIOTS;

        nOTS_Line_type = 1;

        /* make ots fields due to lines and continua of species treated with unified 
         * isoelectronic sequence */
        /* loop over all elements */
        for( ipISO=ipH_LIKE; ipISO<NISO; ++ipISO )
        {
                for( nelem=ipISO; nelem < LIMELM; nelem++ )
                {
                        nOTS1 = ipISO;
                        nOTS2 = nelem;
                        /* do if this stage exists */
                        if( (dense.IonHigh[nelem] >= nelem+1-ipISO )  )
                        {
                                /* generate line ots rates */
                                /* now loop over all possible levels, but cannot include two photon
                                * since there is no pointer to this continuum */
                                /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */
                                for( ipHi=1; ipHi < iso.numLevels_local[ipISO][nelem]; ipHi++ )
                                {
                                        for( ipLo=0; ipLo < ipHi; ipLo++ )
                                        {
                                                /* this signifies a fake line */
                                                /* >>chng 03 may 19, DEST0 is the smallest possible
                                                 * dest prob - not a real number, don't add to ots field */
                                                if( Transitions[ipISO][nelem][ipHi][ipLo].ipCont < 1 ||
                                                        Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pdest<= DEST0 )
                                                        continue;

                                                /* ots rates, the destp prob was set in hydropesc */
                                                Transitions[ipISO][nelem][ipHi][ipLo].Emis->ots = 
                                                        Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul*
                                                        StatesElem[ipISO][nelem][ipHi].Pop*dense.xIonDense[nelem][nelem+1-ipISO]*
                                                        Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pdest;

                                                ASSERT( Transitions[ipISO][nelem][ipHi][ipLo].Emis->ots >= 0. );
                                                /* way to kill lyalpha ots rates
                                                if( nelem==ipHYDROGEN && ipHi==ipH2p && ipLo==ipH1s )
                                                {
                                                        Transitions[ipISO][nelem][ipHi][ipLo].Emis->ots = 0.;
                                                } */

                                                /* finally dump the ots rate into the stack */
                                                if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->ots > SMALLFLOAT ) 
                                                        RT_OTS_AddLine(Transitions[ipISO][nelem][ipHi][ipLo].Emis->ots,
                                                                Transitions[ipISO][nelem][ipHi][ipLo].ipCont );
                                        }
                                }
                                {
                                        /* debugging code for line oscillation problems 
                                         * most often Lya OTS oscillations*/
                                        /*@-redef@*/
                                        enum {DEBUG_LOC=false};
                                        /*@+redef@*/
                                        if( DEBUG_LOC )
                                        {
                                                long ip;
                                                if( ipISO==0 && nelem==0  && nzone>500  ) 
                                                {
                                                        ipHi = 2;
                                                        ipLo = 0;
                                                        ip = Transitions[ipISO][nelem][ipHi][ipLo].ipCont;
                                                        fprintf(ioQQQ,"DEBUG hlyaots\t%.2f\tEdenTrue\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
                                                                fnzone, 
                                                                dense.EdenTrue ,
                                                                Transitions[ipISO][nelem][ipHi][ipLo].Emis->ots,
                                                                opac.opacity_abs[ip-1],
                                                                StatesElem[ipISO][nelem][ipHi].Pop*dense.xIonDense[nelem][nelem+1-ipISO],
                                                                StatesElem[ipISO][nelem][ipHi].Pop,
                                                                dense.xIonDense[nelem][nelem+1-ipISO],
                                                                Transitions[ipISO][nelem][ipHi][ipLo].Emis->Pdest,
                                                                rfield.otslinNew[ip-1],
                                                                rfield.otslin[ip-1]);
                                                }
                                        }
                                }

                                /**************************************************************************
                                *
                                * ots recombination bound-free b-f continua continuum
                                *
                                **************************************************************************/

                                /* put in OTS continuum */
                                /* >>chng 06 aug 17, should go to numLevels_local instead of _max. */
                                for( n=0; n < iso.numLevels_local[ipISO][nelem]; n++ )
                                {
                                        cont_phot_destroyed = (realnum)(iso.RadRecomb[ipISO][nelem][n][ipRecRad]*
                                                (1. - iso.RadRecomb[ipISO][nelem][n][ipRecEsc])*dense.eden*
                                                dense.xIonDense[nelem][nelem+1-ipISO]);
                                        ASSERT( cont_phot_destroyed >= 0. );

                                        /* continuum energy index used in this routine is decremented by one there */
                                        RT_OTS_AddCont(cont_phot_destroyed,iso.ipIsoLevNIonCon[ipISO][nelem][n]);
                                        /* debugging code for rec continua */
                                        {
                                                /*@-redef@*/
                                                enum {DEBUG_LOC=false};
                                                /*@+redef@*/
                                                if( DEBUG_LOC && nzone > 400 && nelem==0 && n==2 )
                                                {
                                                        long ip = iso.ipIsoLevNIonCon[ipISO][nelem][n]-1;
                                                        fprintf(ioQQQ,"hotsdebugg\t%.3f\t%li\th con ots\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
                                                                fnzone, 
                                                                n ,
                                                                StatesElem[ipISO][ipHYDROGEN][2].Pop*dense.xIonDense[ipHYDROGEN][1],
                                                                cont_phot_destroyed,
                                                                cont_phot_destroyed/opac.opacity_abs[ip],
                                                                rfield.otsconNew[ip] ,
                                                                rfield.otscon[ip] ,
                                                                opac.opacity_abs[ip] ,
                                                                hmi.Hmolec[ipMHm] ,
                                                                hmi.HMinus_photo_rate);
                                                }
                                        }
                                }
                        }
                }
        }
        /* more debugging code for rec continua */
        {
                /*@-redef@*/
                enum {DEBUG_LOC=false};
                /*@+redef@*/
                if( DEBUG_LOC )
                {
                        nelem = 0;
                        fprintf(ioQQQ,"hotsdebugg %li \t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\t%.2e\n",
                                nzone, 
                                rfield.otsconNew[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][0]-1],
                                rfield.otsconNew[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][1]-1],
                                rfield.otsconNew[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][3]-1],
                                rfield.otsconNew[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][4]-1],
                                rfield.otsconNew[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][5]-1],
                                rfield.otsconNew[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][6]-1],
                                opac.opacity_abs[iso.ipIsoLevNIonCon[ipH_LIKE][nelem][6]-1]);
                }
        }

        /* now reset Lya dest prob in case is was clobbered */
        Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pdest = (realnum)save_lya_dest;
        Transitions[ipH_LIKE][ipHELIUM][ipH2p][ipH1s].Emis->Pdest = (realnum)save_he2lya_dest;
        iso.RadRecomb[ipH_LIKE][ipHELIUM][ipH1s][ipRecRad] = save_he2rec_dest;

        nelem = ipHELIUM;
        if( dense.lgElmtOn[nelem] && bwnfac > 0. )
        {
                /* increase the destruction prob by the amount we decreased it above */
                Transitions[ipH_LIKE][ipHELIUM][ipH2p][ipH1s].Emis->Pdest /= (realnum)bwnfac;
        }

        if( trace.lgTrace )
        {
                fprintf(ioQQQ,"     RT_OTS Pdest %.2e ots rate %.2e in otslinNew %.2e con opac %.2e\n",
                        Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->Pdest , Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].Emis->ots , 
                        rfield.otslinNew[Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1]  ,
                        opac.opacity_abs[Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1]
                        );
        }

        nOTS_Line_type = 2;
        /* recombination Lya for all elements not yet converted into std isoelectronc form */
        for( nelem=NISO; nelem < LIMELM; nelem++ )
        {
                long int ion;
                /* do not include species treated in iso-electronic fashon in the following,
                 * these were treated above */
                for( ion=0; ion < nelem+1-NISO; ion++ )
                {
                        if( dense.xIonDense[nelem][ion+1] > 0. )
                        {
                                nOTS1 = nelem;
                                nOTS2 = ion;
                                /* now do the recombination Lya */
                                ipla = Heavy.ipLyHeavy[nelem][ion];
                                ASSERT( ipla>0 );
                                esc = opac.ExpmTau[ipla-1];
                                /* xLyaHeavy is set to a fraction of total rad rec in ion_recomb, includes eden */
                                difflya = Heavy.xLyaHeavy[nelem][ion]*dense.xIonDense[nelem][ion+1];
                                /* >>chng 00 dec 22, from MIN2 to MAX2, MIN2 had effect of always
                                 * setting the ots rates to zero */
                                ots = difflya*MAX2(0.f,1.f-esc);
                                /*if( nelem==6 && ion==2 )
                                        fprintf(ioQQQ," debugggnly\t %.2f\t%.2e\n",fnzone, ots );*/
                                ASSERT( ots >= 0.);
                                /*if( iteration == 2 && nzone>290 && ipla== 2339 )
                                        fprintf(ioQQQ,"recdebugg1 %.2e %li %li %.2e %.2e \n", 
                                        ots, nelem, ion,
                                        esc , dense.xIonDense[nelem][ion+1]);*/
                                if( ots > SMALLFLOAT ) 
                                        RT_OTS_AddLine(ots,ipla);

                                /* now do the recombination balmer lines */
                                ipla = Heavy.ipBalHeavy[nelem][ion];
                                esc = opac.ExpmTau[ipla-1];
                                /* xLyaHeavy is set to a fraction of total rad rec in ion_recomb, includes eden */
                                difflya = Heavy.xLyaHeavy[nelem][ion]*dense.xIonDense[nelem][ion+1];
                                /* >>chng 00 dec 22, from MIN2 to MAX2, MIN2 had effect of always
                                 * setting the ots rates to zero */
                                ots = difflya*MAX2(0.f,1.f-esc);
                                ASSERT( ots >= 0.);
                                /*if( iteration == 2 &&nzone==294 && ipla== 2339 )
                                        fprintf(ioQQQ,"recdebugg2 %.2e %li %li\n", ots, nelem, ion );*/
                                if( ots > SMALLFLOAT ) 
                                        RT_OTS_AddLine(ots,ipla);
                        }
                }
        }

        nOTS_Line_type = 3;
        /* do OTS and outward parts of FeII lines, if large atom is turned on */
        FeII_OTS();

        nOTS_Line_type = 4;
        /* do the main set of level1 lines */
#       if 0
        {
#       include "lines_service.h"
        if( nzone>290 ) DumpLine( &TauLines[74] );
        }
#       endif
/*#     define  FRACNEW 1.0f
#       undef FRACNEW*/
        for( nOTS1=1; nOTS1 < nLevel1; nOTS1++ )
        {
                /*realnum otsold = TauLines[nOTS1].ots;*/
                TauLines[nOTS1].Emis->ots = TauLines[nOTS1].Hi->Pop * TauLines[nOTS1].Emis->Aul * TauLines[nOTS1].Emis->Pdest;
                /* * TauLines[nOTS1].ColOvTot;*/
                /*if( nzone )
                        TauLines[nOTS1].ots = TauLines[nOTS1].ots*FRACNEW + otsold*(1.f-FRACNEW);*/
                if( TauLines[nOTS1].Emis->ots > SMALLFLOAT ) 
                        RT_OTS_AddLine( TauLines[nOTS1].Emis->ots , TauLines[nOTS1].ipCont);
        }

        nOTS_Line_type = 5;
        /* do the level2 level 2 lines */
        for( nOTS1=0; nOTS1 < nWindLine; nOTS1++ )
        {
                if( TauLine2[nOTS1].Hi->IonStg < TauLine2[nOTS1].Hi->nelem+1-NISO )
                {
                        TauLine2[nOTS1].Emis->ots = TauLine2[nOTS1].Hi->Pop * TauLine2[nOTS1].Emis->Aul * TauLine2[nOTS1].Emis->Pdest;
                        if( TauLine2[nOTS1].Emis->ots > SMALLFLOAT ) 
                                RT_OTS_AddLine( TauLine2[nOTS1].Emis->ots , TauLine2[nOTS1].ipCont);
                }
        }
        /*fprintf(ioQQQ,"debuggne1 \t%.2e\t%.2e\t%.2e\t%.2e\n", 
                TauLine2[743].ots/SDIV(opac.opacity_abs[743]), 
                TauLine2[744].ots/SDIV(opac.opacity_abs[744]) ,
                TauLine2[743].ots,opac.opacity_abs[743]);*/
        /*CO_OTS - add CO lines to ots fields, called by RT_OTS */

        nOTS_Line_type = 6;
        CO_OTS();

        /* the large H2 molecule */
        nOTS_Line_type = 7;
        H2_RT_OTS();
        /*fprintf(ioQQQ,"DEBUG ggtos exit\n");*/
        return;
}

/* =================================================================== */

void RT_OTS_AddLine(double ots, 
        /* pointer on the f scale */
  long int ip )
{

        DEBUG_ENTRY( "RT_OTS_AddLine()" );

        /* add ots due to line destruction to radiation field */

        /* return if outside bounds of this continuum source, ip > rfield.nflux
         * first case ip==0 happens when called with dummy line */
        if( ip==0 || ip > rfield.nflux )
        { 
                return;
        }

        /*the local ots rate must be non-negative */
        ASSERT( ots >= 0. );
        /* continuum pointer must be positive */
        ASSERT( ip > 0 );

        /* add locally destroyed flux of photons to line OTS array */
        /* check whether local gas opacity  (units cm-1) is positive, if so
         * convert line destruction rate into ots rate by dividing by it */
        if( opac.opacity_abs[ip-1] > 0. )
        {
                rfield.otslinNew[ip-1] += (realnum)(ots/opac.opacity_abs[ip-1]);
        }
        /* first iteration is 1, second is two */
        {
                /*@-redef@*/
                enum {DEBUG_LOC=false};
                /*@+redef@*/
                if( DEBUG_LOC && /*iteration == 2 && nzone>294 &&*/ ip== 2363 /*&& ots > 1e16*/)
                {
                        fprintf(ioQQQ,"DEBUG ots, opc, otsr %.3e\t%.3e\t%.3e\t",
                                ots ,
                                opac.opacity_abs[ip-1],
                                ots/opac.opacity_abs[ip-1] );
                        fprintf(ioQQQ,"iteration %li type %i %i %i \n", 
                                iteration, 
                                nOTS_Line_type,
                                nOTS1,nOTS2 );
                }
        }
        return;
}

/* =================================================================== */

/*add local destruction of continuum to ots field */
STATIC void RT_OTS_AddCont(
                                        /* the ots rate itself */
                                        realnum ots, 
                                        /* pointer to continuum cell for ots, on f scale */
                                        long int ip)
{

        DEBUG_ENTRY( "RT_OTS_AddCont()" );

        /* 
         * routine called to add ots due to continuum destruction to
         * radiation field
         */

        /* check if outside bounds of this continuum source */
        if( ip > rfield.nflux )
        { 
                return;
        }

        ASSERT( ip> 0 );
        ASSERT( ots >= 0. );
        ASSERT( ip <= rfield.nupper );

        /* add locally destroyed flux of photons to continuum OTS array */
        /* check whether local gas opacity  (units cm-1) is positive, if so
         * convert continuum destruction rate into ots rate by dividing by it */
        if( opac.opacity_abs[ip-1] > 0. )
        {
                rfield.otsconNew[ip-1] += (realnum)(ots/opac.opacity_abs[ip-1]);
                /*if( ots>0. ) fprintf(ioQQQ,
                        "buggg ip %li ots %.2e opac %.2e cont_phot_destroyed %.2e\n",
                        ip, ots , opac.opacity_abs[ip-1] ,rfield.otsconNew[ip-1]);*/
        }
        return;
}

#if 0
/* check whether Lya OTS rate is oscillating */
STATIC bool lgLyaOTS_oscil( realnum ots_new )
{
        static realnum old , older;
        bool lgReturn = false;

        if( !conv.nTotalIoniz || conv.lgSearch )
        {
                /* this is very first call on this iteration, zero everything out */
                old = 0.;
                older = 0.;
        }

        /* is sign of change oscillating? */
        if( conv.nPres2Ioniz>1 && (ots_new-old) * (old-older) < 0.)
                lgReturn = true;

        older = old;
        old = ots_new;
        return( lgReturn );
}

        /* >>chng 04 jan 26, check whether Lya is oscillating */
        if( lgLyaOTS_oscil( rfield.otslinNew[Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1] ) )
        {
                fprintf(ioQQQ,"DEBUG lya ots %.2f %.3e\n", 
                        fnzone, 
                        rfield.otslinNew[Transitions[ipH_LIKE][ipHYDROGEN][ipH2p][ipH1s].ipCont-1] );
        }
#endif
/* =================================================================== */

/*RT_OTS_Update update ots rates, called in ConvBase */
void RT_OTS_Update(
        /* summed ots rates */
        double *SumOTS , 
        /* array index for constituent that changed the most
        long int *ipOTSchange,  */
        /* limit on how large a change to allow, 0 for no limit */
        double BigFrac)
{
        long int i;

        static bool lgNeedSpace=true;
        static double *old_OTS_line_x_opac , *old_OTS_cont_x_opac;
        realnum FacBig , FacSml;

        DEBUG_ENTRY( "RT_OTS_Update()" );

        if( BigFrac <= 0. )
                BigFrac = 10.;
        FacBig = (realnum)(1. + BigFrac);
        FacSml = 1.f/FacBig;

        /* create space to save arrays if needed, one time per coreload */
        if( lgNeedSpace )
        {
                old_OTS_line_x_opac = (double*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(double)) );
                old_OTS_cont_x_opac = (double*)MALLOC((size_t)((unsigned)rfield.nupper*sizeof(double)) );
        }
        lgNeedSpace = false;

        *SumOTS = 0.;

        /* option to kill ots rates with no ots lines command */
        if( rfield.lgKillOTSLine )
        {
                for( i=0; i < rfield.nflux; i++ )
                {
                        rfield.otslinNew[i] = 0.;
                }
        }

        /* during search phase set ots to current values */
        if( nzone==0 )
        {
                for( i=0; i < rfield.nflux; i++ )
                {
                        old_OTS_line_x_opac[i] = rfield.otslinNew[i] * opac.opacity_abs[i];
                        old_OTS_cont_x_opac[i] = rfield.otsconNew[i] * opac.opacity_abs[i];
                }
        }

        /* remember largest change in ots rates */
        *SumOTS = 0.;
        /* now update new ots rates */
        for( i=0; i < rfield.nflux; ++i )
        {
                double OTS_line_x_opac = rfield.otslinNew[i] * opac.opacity_abs[i];
                double OTS_cont_x_opac = rfield.otsconNew[i] * opac.opacity_abs[i];
                double CurrentInverseOpacity = 1./MAX2( SMALLDOUBLE , opac.opacity_abs[i] );

                /* get new ots line rates for each cell, but don't let change be too big */
                if( OTS_line_x_opac > old_OTS_line_x_opac[i]*FacBig )
                {
                        rfield.otslin[i] = rfield.otslin[i]*FacBig;
                }
                else if( OTS_line_x_opac < old_OTS_line_x_opac[i]*FacSml )
                {
                        rfield.otslin[i] = rfield.otslin[i]*FacSml;
                }
                else 
                {
                        rfield.otslin[i] = rfield.otslinNew[i];
                }

                /* get new ots continuum rates for each cell */
                if( OTS_cont_x_opac > old_OTS_cont_x_opac[i]*FacBig )
                {
                        rfield.otscon[i] = rfield.otscon[i]*FacBig;
                }
                else if( OTS_cont_x_opac < old_OTS_cont_x_opac[i]*FacSml )
                {
                        rfield.otscon[i] = rfield.otscon[i]*FacSml;
                }
                else 
                {
                        rfield.otscon[i] = rfield.otsconNew[i];
                }

                /* zero out the new otscon vector*/
                rfield.otsconNew[i] = 0.;

                /* zero out the new otslin vector*/
                rfield.otslinNew[i] = 0.;

                /* now update to new values */
                old_OTS_line_x_opac[i] = rfield.otslin[i] * opac.opacity_abs[i];
                old_OTS_cont_x_opac[i] = rfield.otscon[i] * opac.opacity_abs[i];

                /* this is local ots continuum created by destroyed diffuse continua, 
                 * currently only two-photon */
                rfield.ConOTS_local_OTS_rate[i] = (realnum)((double)rfield.ConOTS_local_photons[i]*CurrentInverseOpacity);

                /*if( i==14 )
                {
                        fprintf(ioQQQ,"DEBUG OTS rate %.2e %.2e %.2e %.2e \n", *SumOTS ,
                                rfield.otscon[14] , rfield.otslin[14],opac.opacity_abs[14]);
                                fflush(ioQQQ);
                }*/

                /* remember sum of ots rates for convergence criteria */
                *SumOTS += (rfield.otscon[i] + rfield.otslin[i])*opac.opacity_abs[i];

                rfield.SummedDif[i] = rfield.otscon[i] + rfield.otslin[i] + rfield.outlin_noplot[i]+ 
                  rfield.ConInterOut[i]*rfield.lgOutOnly + rfield.outlin[i] +
                  rfield.ConOTS_local_OTS_rate[i];

                rfield.SummedCon[i] = rfield.flux[i] + rfield.SummedDif[i];
                rfield.SummedOcc[i] = rfield.SummedCon[i]*rfield.convoc[i];

        }


        /* >>chng 02 oct 18, add this */
        /* sum of accumulated flux from particular frequency to infinity */
        rfield.flux_accum[rfield.nflux-1] = 0;
        for( i=1; i < rfield.nflux; i++ )
        {
                rfield.flux_accum[rfield.nflux-i-1] = rfield.flux_accum[rfield.nflux-i] +
                        rfield.SummedCon[rfield.nflux-i-1];
        }
        /* this has to be positive since is sum of all photons in SummedCon */
        ASSERT( rfield.flux_accum[0]> 0. );

        /* >>chng 02 jul 23, set to black body at local temp if in optically thick continuum,
         * between plasma frequency and energy where brems is thin */
        ASSERT(rfield.ipPlasma>0 );

        /* >>chng 02 jul 25, set all radiation fields to zero below plasma frequency */
        for( i=0; i < rfield.ipPlasma-1; i++ )
        {
                rfield.otscon[i] = 0.;
                rfield.ConOTS_local_OTS_rate[i] = 0.;
                rfield.ConOTS_local_photons[i] = 0.;
                rfield.otslin[i] = 0.;
                rfield.SummedDif[i] = 0.;
                rfield.OccNumbBremsCont[i] = 0.;
                rfield.SummedCon[i] = 0.;
                rfield.SummedOcc[i] = 0.;
                rfield.ConInterOut[i] = 0.;
        }
        /* this loop evaluates occupation number for brems continuum,
         * only used for induced two photon emission */
        if( rfield.ipEnergyBremsThin > 0 )
        {
                for( i=rfield.ipPlasma-1; i < rfield.nflux; i++ )
                {
                        /* this corrects for opacity / optical depth in brems - brems opacity goes as
                         * energy squared. */
                        /* when totally optically thin to brems rfield.ipEnergyBremsThin is zero,
                         * so need this max */
                        realnum factor = MIN2(1.f,rfield.anu2[MAX2(0,rfield.ipEnergyBremsThin-1)] / rfield.anu2[i]);

                        /* occupation number for black body is 1/ (exp hn/kT) -1) */
                        rfield.OccNumbBremsCont[i] = (realnum)(1./(1./SDIV(rfield.ContBoltz[i]) - 1.)) * factor;
                }
        }
        return;

#       if 0
        OTSOld = 0.;
        OTSNew = 0.;
        BigOTSNew = 0.;
        /* find summed ots rates, old and new */
        for( i=0; i < rfield.nflux; i++ )
        {
                OTSOld += rfield.otscon[i] + rfield.otslin[i];
                OTSNew += rfield.otsconNew[i] + rfield.otslinNew[i];
                if( BigOTSNew < rfield.otsconNew[i] + rfield.otslinNew[i] )
                {
                        BigOTSNew = rfield.otsconNew[i] + rfield.otslinNew[i];
                }
        }

        /* we now have old and new rates, what is the ratio, and by how much will
         * we allow this to change? */
        if( BigFrac == 0. || conv.lgSearch || OTSOld < SMALLFLOAT )
        {
                /* this branch, want a clear update of ots rates, either because
                 * requested with BigFrac == 0, or we are in search phase */
                FracOld = 0.;
        }
        else
        {
                if( OTSNew > OTSOld )
                {
                        chng = fabs(1. - OTSOld/SDIV(OTSNew) );
                }
                else
                {
                        chng = fabs(1. - OTSNew/SDIV(OTSOld) );
                }

                if( chng < BigFrac )
                {
                        /* this branch, old and new do not differ by much, so use new */
                        FracOld = 0.;
                }
                else
                {
                        /* this branch, too large a difference between old and new, cap it to BigFrac */
                        FracOld = (1. - BigFrac / chng);
                        ASSERT( FracOld >= 0. );
                        FracOld = MIN2( 0.25 , FracOld );
                }
        }

        /* fraction old and new ots rates */
        FracNew = 1. - FracOld;
        fprintf(ioQQQ," DEBUG FracNew\t%.2e\n", FracNew);

        /* remember largest change in ots rates */
        changeOTS = 0.;
        /*fprintf(ioQQQ," sumcontinuum zone%li 1168=%e\n", nzone,rfield.otslin[1167]);*/

        for( i=0; i < rfield.nflux; i++ )
        {
                /* >>chng 01 feb 01, define inverse opacity in safe manner */
                double CurrentInverseOpacity = 1./MAX2( SMALLDOUBLE , opac.opacity_abs[i] );

                /* remember which energy had largest change in ots rates */
                if( fabs( rfield.otscon[i]+rfield.otslin[i]-rfield.otsconNew[i]-rfield.otslinNew[i])> changeOTS)
                {
                        changeOTS = 
                          fabs( rfield.otscon[i]+rfield.otslin[i]-rfield.otsconNew[i]-rfield.otslinNew[i]);
                }

                /* >>chng 01 apr 10, only change ots with means if FracOld not zero,
                 * this is to avoid taking means when this routine called by StartEndIter 
                 * to reset sums of rates */
                if( BigFrac > 0. && conv.nTotalIoniz > 0 )
                {
                        /* the New vectors are the ones updated by the AddOTS routines,
                         * and have the current OTS rates.  The otscon and otslin vectors
                         * have the rates from the previous refresh of the New vectors*/
                        /* here is the refresh.  all were initially set to zero in call to 
                         * RT_OTS_Zero (below) from zero */
                        rfield.otscon[i] = (realnum)((rfield.otscon[i]*FracOld +rfield.otsconNew[i]*FracNew));

                        /* this is local ots continuum created by destroyed diffuse continua, 
                         * currently only two-photon */
                        /* >>chng 00 dec 27, define this continuum and do not count 2-photon in old var */
                        rfield.ConOTS_local_OTS_rate[i] = (realnum)(rfield.ConOTS_local_photons[i]*CurrentInverseOpacity);

                        /* current ots will be average of old and new */
                        rfield.otslin[i] = (realnum)((rfield.otslin[i]*FracOld + rfield.otslinNew[i]*FracNew));
                }

                /* zero out the new otscon vector*/
                rfield.otsconNew[i] = 0.;

                /* zero out the new otslin vector*/
                rfield.otslinNew[i] = 0.;

                /* remember sum of ots rates for convergence criteria */
                *SumOTS += (rfield.otscon[i] + rfield.otslin[i])*opac.opacity_abs[i];
                {
                        /* following should be set true to print strongest ots contributors */
                        /*@-redef@*/
                        enum {DEBUG_LOC=false};
                        /*@+redef@*/
                        if( DEBUG_LOC && (nzone==378) )
                        {
                                if( conv.nPres2Ioniz > 3 ) 
                                        cdEXIT( EXIT_FAILURE );
                                fprintf(ioQQQ,"rtotsbugggg\t%li\t%.3e\t%.3e\t%.3e\t%.3e\n",
                                        conv.nPres2Ioniz,
                                        rfield.anu[i],
                                        opac.opacity_abs[i],
                                        rfield.otscon[i],
                                        rfield.otslin[i]);
                        }
                }

                /* >>chng 00 dec 27, include ConOTS_local_OTS_rate */
                /* >>chng 01 jul 04, add *rfield.lgOutOnly logic */
                rfield.SummedDif[i] = rfield.otscon[i] + rfield.otslin[i] + rfield.outlin_noplot[i]+ 
                  rfield.ConInterOut[i]*rfield.lgOutOnly + rfield.outlin[i] +
                  rfield.ConOTS_local_OTS_rate[i];

                rfield.SummedCon[i] = rfield.flux[i] + rfield.SummedDif[i];
                rfield.SummedOcc[i] = rfield.SummedCon[i]*rfield.convoc[i];
        }

        ASSERT(*SumOTS >= 0. );

        /* >>chng 02 oct 18, add this */
        /* sum of accumulated flux from particular frequency to infinity */
        rfield.flux_accum[rfield.nflux-1] = 0;
        for( i=1; i < rfield.nflux; i++ )
        {
                rfield.flux_accum[rfield.nflux-i-1] = rfield.flux_accum[rfield.nflux-i] +
                        rfield.SummedCon[rfield.nflux-i-1];
        }
        /* this has to be positive since is sum of all photons in SummedCon */
        ASSERT( rfield.flux_accum[0]> 0. );

        /* >>chng 02 jul 23, set to black body at local temp if in optically thick continuum,
         * between plasma frequency and energy where brems is thin */
        ASSERT(rfield.ipPlasma>0 );

        /* >>chng 02 jul 25, set all radiation fields to zero below plasma frequency */
        for( i=0; i < rfield.ipPlasma-1; i++ )
        {
                rfield.otscon[i] = 0.;
                rfield.ConOTS_local_OTS_rate[i] = 0.;
                rfield.ConOTS_local_photons[i] = 0.;
                rfield.otslin[i] = 0.;
                rfield.SummedDif[i] = 0.;
                rfield.OccNumbBremsCont[i] = 0.;
                rfield.SummedCon[i] = 0.;
                rfield.SummedOcc[i] = 0.;
                rfield.ConInterOut[i] = 0.;
        }
        /* this loop evaluates occupation number for brems continuum,
         * only used for induced two photon emission */
        if( rfield.ipEnergyBremsThin > 0 )
        {
                for( i=rfield.ipPlasma-1; i < rfield.nflux; i++ )
                {
                        /* this corrects for opacity / optical depth in brems - brems opacity goes as
                         * energy squared. */
                        /* when totally optically thin to brems rfield.ipEnergyBremsThin is zero,
                         * so need this max */
                        realnum factor = MIN2(1.f,rfield.anu2[MAX2(0,rfield.ipEnergyBremsThin-1)] / rfield.anu2[i]);

                        /* occupation number for black body is 1/ (exp hn/kT) -1) */
                        rfield.OccNumbBremsCont[i] = (realnum)(1./(1./SDIV(rfield.ContBoltz[i]) - 1.)) * factor;
                }
        }

        {
                /* following should be set true to print strongest ots contributors */
                /*@-redef@*/
                enum {DEBUG_LOC=false};
                /*@+redef@*/
                if( DEBUG_LOC && (nzone>0) )
                {
                        double BigOTS;
                        int ipOTS=0;
                        BigOTS = -1.;
                        /* find the biggest ots contributor */
                        /*for( i=iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][0]; i < rfield.nflux; i++ )*/
                        for( i=0; i < rfield.nflux; i++ )
                        {
                                if( (rfield.otscon[i] + rfield.otslin[i])*opac.opacity_abs[i] > BigOTS )
                                {
                                        BigOTS = (rfield.otscon[i] + rfield.otslin[i])*opac.opacity_abs[i];
                                        ipOTS = (int)i;
                                }
                        }
                        fprintf(ioQQQ,
                                " sumcontinuunots zone %li SumOTS %.2e %.2eRyd opc:%.2e OTS%.2e lin%.2e con:%.2e %s %s cell%i FracNew %.2f \n",
                                nzone,
                                *SumOTS,
                                rfield.anu[ipOTS] , 
                                opac.opacity_abs[ipOTS],
                                BigOTS ,
                                rfield.otslin[ipOTS],
                                rfield.otscon[ipOTS] , 
                                /* line label */
                                rfield.chLineLabel[ipOTS] ,
                                /* cont label*/
                                rfield.chContLabel[ipOTS] ,
                                ipOTS,
                                FracNew);
                }
        }
        return;
#       endif
}

/* =================================================================== */

/*RT_OTS_Zero zero out some vectors - this is only called when code 
 * initialized by ContSetIntensity */
void RT_OTS_Zero( void )
{
        long int i;

        DEBUG_ENTRY( "RT_OTS_Zero()" );

        /* this loop goes up to nflux itself (<=) since the highest cell
         * will be used to pass unity through the code to verify integrations */
        for( i=0; i <= rfield.nflux; i++ )
        {
                rfield.SummedDif[i] = 0.;
                /* the four main ots vectors */
                rfield.otscon[i] = 0.;
                rfield.otslin[i] = 0.;
                rfield.otslinNew[i] = 0.;
                rfield.otsconNew[i] = 0.;
                /* >>chng 05 mar 03, add this */
                rfield.trans_coef_zone[i] = 1.f;
                rfield.trans_coef_total[i] = 1.f;

                rfield.ConInterOut[i] = 0.;
                rfield.outlin[i] = 0.;
                rfield.outlin_noplot[i] = 0.;
                rfield.SummedDif[i] = 0.;
                /* "zero" for the summed con will be just the incident radiation field */
                rfield.SummedCon[i] = rfield.flux[i];
                rfield.SummedOcc[i] = rfield.SummedCon[i]*rfield.convoc[i];
                rfield.ConOTS_local_photons[i] = 0.;
                rfield.ConOTS_local_OTS_rate[i] = 0.;
        }
        return;
}

/* =================================================================== */

/*RT_OTS_ChkSum sanity check confirms summed continua reflect contents of individuals */
void RT_OTS_ChkSum(long int ipPnt)
{
        static long int nInsane=0;
        long int i;
        double phisig;
#       define LIM_INSAME_PRT   30

        DEBUG_ENTRY( "RT_OTS_ChkSum()" );

        /* this entire sub is a sanity check */
        /* >>chng 02 jul 23, lower bound from 0 to rfield.ipEnergyBremsThin - since now
         * set radiation field to black body below this energy */
        for( i=rfield.ipEnergyBremsThin; i < rfield.nflux; i++ )
        {
                phisig = rfield.otscon[i] + rfield.otslin[i] + rfield.ConInterOut[i]*rfield.lgOutOnly + 
                  rfield.outlin[i]+
                  rfield.outlin_noplot[i]+
                  rfield.ConOTS_local_OTS_rate[i];
                /* >>chng 02 sep 19, add sec test on SummedDif since it can be zero whild 
                 * phisig is just above small float */
                if( phisig > 0. && rfield.SummedDif[i]> 0.)
                {
                        if( fabs(rfield.SummedDif[i]/phisig-1.) > 1e-3 )
                        {
                                ++nInsane;
                                /* limit amount of printout - in many failures would print entire
                                 * continuum array */
                                if( nInsane < LIM_INSAME_PRT )
                                {
                                        fprintf( ioQQQ, " PROBLEM RT_OTS_ChkSum insane SummedDif at energy %.5e error= %.2e i=%4ld\n", 
                                        rfield.anu[i], rfield.SummedDif[i]/phisig - 1., i );
                                        fprintf( ioQQQ, " SummedDif, sum are%11.4e%11.4e\n", 
                                        rfield.SummedDif[i], phisig );
                                        fprintf( ioQQQ, " otscon otslin ConInterOut outlin are%11.4e%11.4e%11.4e%11.4e\n", 
                                        rfield.otscon[i], rfield.otslin[i]+rfield.outlin_noplot[i], rfield.ConInterOut[i], 
                                        rfield.outlin[i]+rfield.outlin_noplot[i] );
                                        fprintf( ioQQQ, " line continuum here are %4.4s %4.4s\n", 
                                        rfield.chLineLabel[i], rfield.chContLabel[i] );
                                }
                        }
                }

                phisig += rfield.flux[i];
                /* >>chng 02 sep 19, add sec test on SummedDif since it can be zero when 
                 * phisig is just above small float */
                if( phisig > 0. && rfield.SummedDif[i]> 0. )
                {
                        if( fabs(rfield.SummedCon[i]/phisig-1.) > 1e-3 )
                        {
                                ++nInsane;
                                /* limit amount of printout - in many failures would print entire
                                 * continuum array */
                                if( nInsane < LIM_INSAME_PRT )
                                {
                                        fprintf( ioQQQ, " PROBLEM RT_OTS_ChkSum %3ld, insane SummedCon at energy %.5e error=%.2e i=%ld\n", 
                                        ipPnt, rfield.anu[i], rfield.SummedCon[i]/phisig - 1., i );
                                        fprintf( ioQQQ, " SummedCon, sum are %.4e %.4e\n", 
                                        rfield.SummedCon[i], phisig );
                                        fprintf( ioQQQ, " otscon otslin ConInterOut outlin flux are%.4e %.4e %.4e %.4e %.4e\n", 
                                        rfield.otscon[i], rfield.otslin[i]+rfield.outlin_noplot[i], rfield.ConInterOut[i], 
                                        rfield.outlin[i]+rfield.outlin_noplot[i], rfield.flux[i] );
                                        fprintf( ioQQQ, " line continuum here are %s %s\n", 
                                        rfield.chLineLabel[i], rfield.chContLabel[i]
                                        );
                                }
                        }
                }
        }

        if( nInsane > 0 )
        {
                fprintf( ioQQQ, " PROBLEM RT_OTS_ChkSum too much insanity to continue.\n");
                /* TotalInsanity exits after announcing a problem */
                TotalInsanity();
        }
        return;
}

/* =================================================================== */

/*RT_OTS_PrtRate print continuum and line ots rates when trace ots is on */
void RT_OTS_PrtRate(
          /* weakest rate to print */
          double weak ,
          /* flag, 'c' continuum, 'l' line, 'b' both */
          int chFlag )
{
        long int i;

        DEBUG_ENTRY( "RT_OTS_PrtRate()" );

        /* arg must be one of these three */
        ASSERT( chFlag=='l' || chFlag=='c' || chFlag=='b' );

        /* 
         * both printouts have cell number (on C array scale)
         * energy in ryd
         * the actual value of the ots rate 
         * the ots rate relative to the continuum at that energy 
         * rate times opacity 
         * all are only printed if greater than weak 
         */

        /*===================================================================*/
        /* first print ots continua                                          */
        /*===================================================================*/
        if( chFlag == 'c' || chFlag == 'b' )
        {
                fprintf( ioQQQ, "     DEBUG OTSCON array, anu, otscon, opac, OTS*opac limit:%.2e zone:%.2f IonConv?%c\n",
                        weak,fnzone ,TorF(conv.lgConvIoniz) );

                for( i=0; i < rfield.nupper; i++ )
                {
                        if( rfield.otscon[i]*opac.opacity_abs[i] > weak )
                        {
                                fprintf( ioQQQ, "     %4ld%12.4e%12.4e%12.4e%12.4e %s \n", 
                                  i, 
                                  rfield.anu[i], 
                                  rfield.otscon[i], 
                                  opac.opacity_abs[i], 
                                  rfield.otscon[i]*opac.opacity_abs[i], 
                                  rfield.chContLabel[i]);

                        }
                }
        }

        /*===================================================================*/
        /* second print ots line rates                                       */
        /*===================================================================*/
        if( chFlag == 'l' || chFlag == 'b' )
        {
                fprintf( ioQQQ, "     DEBUG OTSLIN array, anu, otslin, opac, OTS*opac Lab nLine limit:%.2e zone:%.2f IonConv?%c\n",
                        weak,fnzone,TorF(conv.lgConvIoniz) );

                for( i=0; i < rfield.nupper; i++ )
                {
                        if( rfield.otslin[i]*opac.opacity_abs[i] > weak )
                        {
                                fprintf( ioQQQ, "     %4ld%12.4e%12.4e%12.4e%12.4e %s %3li\n", 
                                  i, 
                                  rfield.anu[i], 
                                  rfield.otslin[i], 
                                  opac.opacity_abs[i],
                                  rfield.otslin[i]*opac.opacity_abs[i],
                                  rfield.chLineLabel[i] ,
                                  rfield.line_count[i] );
                        }
                }
        }
        return;
}

---