iso_level.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*iso_level solve for iso-sequence level populations */
#include "cddefines.h"
#include "atmdat.h"
#include "continuum.h"
#include "conv.h"
#include "dense.h"
#include "dynamics.h"
#include "elementnames.h"
#include "grainvar.h"
#include "he.h"
#include "heavy.h"
#include "helike.h"
#include "hmi.h"
#include "hydrogenic.h"
#include "ionbal.h"
#include "iso.h"
#include "mole.h"
#include "opacity.h"
#include "phycon.h"
#include "physconst.h"
#include "rfield.h"
#include "secondaries.h"
#include "taulines.h"
#include "thirdparty.h"
#include "trace.h"
#include "yield.h"
 
/* solve for level populations  */
void iso_level( const long int ipISO, const long int nelem, double &renorm)
{
        long int ipHi,
                ipLo,
                i,
                level,
                level_error;
 
        t_iso_sp* sp = &iso_sp[ipISO][nelem]; 
        const long int numlevels_local = sp->numLevels_local;
 
        double BigError;
 
        int32 nerror;
        double HighestPColOut = 0.;
        bool lgNegPop=false;
        valarray<int32> ipiv(numlevels_local); 
        /* this block of variables will be obtained and freed here */
        double source=0.,
                sink=0.;
        valarray<double> PopPerN(sp->n_HighestResolved_local+1);
 
        DEBUG_ENTRY( "iso_level()" );
        renorm = 1.;
 
        /* check that we were called with valid charge */
        ASSERT( nelem >= ipISO );
        ASSERT( nelem < LIMELM );
 
        /* these two collision rates must be the same or we are in big trouble,
         * since used interchangeably */
        ASSERT( ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0]< SMALLFLOAT ||
                fabs( (sp->fb[0].ColIoniz* dense.EdenHCorr) /
                SDIV(ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0] ) - 1.) < 0.001 );
 
        if( (dense.IonHigh[nelem] >= nelem - ipISO) &&
                 (dense.IonLow[nelem] <= nelem - ipISO) &&
                 ionbal.RateRecomIso[nelem][ipISO] > 0. )
        {
                /* get simple estimate of atom to ion ratio */
                sp->xIonSimple = ionbal.RateIonizTot(nelem,nelem-ipISO)/ionbal.RateRecomIso[nelem][ipISO];
        }
        else
        {
                sp->xIonSimple = 0.;
        }
 
        /* which case atom to solve??? */
        if( dense.xIonDense[nelem][nelem+1-ipISO] < SMALLFLOAT || sp->xIonSimple < 1e-35 )
        {
                /* don't bother if no ionizing radiation */
                strcpy( sp->chTypeAtomUsed, "zero " );
                if( trace.lgTrace && (nelem == trace.ipIsoTrace[ipISO]) )
                {
                        fprintf( ioQQQ, "     iso_level iso=%2ld nelem=%2ld simple II/I=%10.2e so not doing equilibrium, doing %s.\n", 
                                ipISO, nelem, sp->xIonSimple, sp->chTypeAtomUsed );
                }
 
                /* total ionization will just be the ground state */
                lgNegPop = false;
                sp->st[0].Pop() = dense.xIonDense[nelem][nelem-ipISO];
                for( long n=1; n < numlevels_local; n++ )
                {
                        sp->st[n].Pop() =  0.;
                }
                sp->qTot2S = 0.;
        }
        else
        {
                /* fill in recombination vector - values were set in iso_ionize_recombine.cpp */
                valarray<double>
                        creation(numlevels_local),
                        error(numlevels_local),
                        work(numlevels_local),
                        Save_creation(numlevels_local);
 
                ASSERT( dense.xIonDense[nelem][nelem+1-ipISO] >= 0.f );
                for( level=0; level < numlevels_local; level++ )
                {
                        /* total recombination from once more ionized [cm-3 s-1] */
                        creation[level] = sp->fb[level].RateCont2Level * dense.xIonDense[nelem][nelem+1-ipISO];
                }
                
                double ctsource=0.0, ctsink=0.0, ctrec=0.0;
                /* now charge transfer - all into/from ground, two cases, H and not H */
                if( nelem==ipHYDROGEN )
                {
                        /* charge transfer, hydrogen onto everything else */
                        /* charge exchange recombination */
                        ctrec += atmdat.CharExcRecTotal[nelem];
                        ctsink += atmdat.CharExcIonTotal[nelem];
                }
                else if( nelem==ipHELIUM && ipISO==ipHE_LIKE )
                {
                        /* this is recom of He due to ct with all other gas constituents */
                        ctrec += atmdat.CharExcRecTotal[nelem];
                        ctsink += atmdat.CharExcIonTotal[nelem];
                }
                else
                {
                        for (long nelem1=ipHYDROGEN; nelem1 < t_atmdat::NCX; ++nelem1)
                        {
                                long ipISO=nelem1;
                                ctrec += atmdat.CharExcRecTo[nelem1][nelem][nelem-ipISO]*iso_sp[ipISO][nelem1].st[0].Pop();
                                ctsink += atmdat.CharExcIonOf[nelem1][nelem][nelem-ipISO]*dense.xIonDense[nelem1][1];
                        }                                       
                }
                ctsource += ctrec*dense.xIonDense[nelem][nelem+1-ipISO];
                
                if ( nelem > ipISO )
                {
                        double ction=0.0;
                        if( nelem==ipHELIUM )
                        {
                                ctsink += atmdat.CharExcRecTotal[nelem];
                                ction += atmdat.CharExcIonTotal[nelem];
                        }
                        else
                        {
                                for (long nelem1=ipHYDROGEN; nelem1 < t_atmdat::NCX; ++nelem1)
                                {
                                        long ipISO1=nelem1;
                                        ctsink += atmdat.CharExcRecTo[nelem1][nelem][nelem-(ipISO+1)]*iso_sp[ipISO1][nelem1].st[0].Pop();
                                        ction += atmdat.CharExcIonOf[nelem1][nelem][nelem-(ipISO+1)]*dense.xIonDense[nelem1][1];
                                }
                        }
                        ctsource += ction * dense.xIonDense[nelem][nelem-(ipISO+1)];
                }
                
                /* now do the 2D array */
                multi_arr<double,2,C_TYPE> z, SaveZ;
                z.alloc(numlevels_local,numlevels_local);
                z.zero();
                vector<double> Boltzmann(numlevels_local-1);
                for (unsigned lev = 0; lev < Boltzmann.size(); ++lev)
                        Boltzmann[lev] = 1.0;
 
                /* this branch is main solver, full level populations 
                 * assert since this code must change if NISO ever increased */
                ASSERT( NISO == 2 );
                long SpinUsed[NISO] = {2,1};
                long indexNmaxP =
                        sp->QuantumNumbers2Index[ sp->n_HighestResolved_local ][1][SpinUsed[ipISO]];
 
                strcpy( sp->chTypeAtomUsed, "Popul" );
                if( trace.lgTrace && (nelem == trace.ipIsoTrace[ipISO]) )
                {
                        fprintf( ioQQQ, "     iso_level iso=%2ld nelem=%2ld doing regular matrix inversion, %s\n", 
                                ipISO, nelem, sp->chTypeAtomUsed );
                }
 
                //qList::const_iterator StElm = StatesElemNEW[nelem][nelem-ipISO].begin();
                qList::const_iterator StElm = sp->st.begin();
                
                /* master balance equation, use when significant population */
                for( level=0; level < numlevels_local; level++ )
                {
                        /* all process depopulating level and placing into the continuum
                         * this does NOT include grain charge transfer ionization, added below */
                        z[level][level] += sp->fb[level].RateLevel2Cont;
 
                        if( level == 1 )
                                /* >>chng 05 dec 21, rm eden to make into rate coefficient */
                                sp->qTot2S = sp->fb[level].ColIoniz;
 
                        //TransitionList::iterator Trans = sp->tr.begin() + sp->ipTrans[level][0];
 
                        if (level != 0)
                        {
                                double arg = (StElm[level].energy().K()-StElm[level-1].energy().K()) /
                                        phycon.te;
                                arg = MAX2( -38., arg); fixit(); // Wouldn't need to mask this out if levels were in order
                                double bstep =  sexp( arg );
                                for ( ipLo = 0; ipLo < level; ++ipLo )
                                        Boltzmann[ipLo] *= bstep;
                        }
 
                        /* all processes populating level from below */
                        for( ipLo=0; ipLo < level; ipLo++ )
                        {
                                double coll_down = sp->trans(level,ipLo).Coll().ColUL( colliders );
 
                                double RadDecay = MAX2( iso_ctrl.SmallA, sp->trans(level,ipLo).Emis().Aul()*
                                                        (sp->trans(level,ipLo).Emis().Pesc() + 
                                                         sp->trans(level,ipLo).Emis().Pelec_esc() + 
                                                         sp->trans(level,ipLo).Emis().Pdest())*
                                                        KILL_BELOW_PLASMA(sp->trans(level,ipLo).EnergyRyd()) );
 
                                double pump = MAX2( iso_ctrl.SmallA, sp->trans(level,ipLo).Emis().pump() *
                                                    KILL_BELOW_PLASMA(sp->trans(level,ipLo).EnergyRyd()) );
 
                                if( iso_ctrl.lgRandErrGen[ipISO] )
                                {
                                        coll_down *= sp->ex[level][ipLo].ErrorFactor[IPCOLLIS];
                                        RadDecay *= sp->ex[level][ipLo].ErrorFactor[IPRAD];
                                        pump *= sp->ex[level][ipLo].ErrorFactor[IPRAD];
                                }
 
                                double coll_up = coll_down * 
                                        (double)StElm[level].g()/
                                        (double)StElm[ipLo].g()*
                                        Boltzmann[ipLo];
 
                                z[ipLo][ipLo] += coll_up + pump ;
                                z[ipLo][level] -= coll_up + pump ;
 
                                double pump_down = pump *
                                        (double)StElm[ipLo].g()/
                                        (double)StElm[level].g();
 
                                z[level][level] += RadDecay + pump_down + coll_down;
                                z[level][ipLo] -= RadDecay + pump_down + coll_down;
 
                                if( level == indexNmaxP )
                                {
                                        HighestPColOut += coll_down;
                                }                                
                                if( ipLo == indexNmaxP )
                                {
                                        HighestPColOut += coll_up;
                                }
 
                                /* find total collisions out of 2^3S to singlets. */
                                if( (level == 1) && (ipLo==0) )
                                {
                                        sp->qTot2S += coll_down/dense.EdenHCorr;
                                }
                                if( (ipLo == 1) && (ipISO==ipH_LIKE || StElm[level].S()==1) )
                                {
                                        sp->qTot2S += coll_up/dense.EdenHCorr;
                                }
                        }
                }
                
                if( ipISO == nelem )
                {
                        /* iso_ctrl.lgCritDensLMix[ipISO] is a flag used to print warning if density is
                        * too low for first collapsed level to be l-mixed.  Check is if l-mixing collisions
                        * out of highest resolved singlet P are greater than sum of transition probs out.       */
                        if( HighestPColOut < 1./sp->st[indexNmaxP].lifetime() )
                        {
                                iso_ctrl.lgCritDensLMix[ipISO] = false;
                        }
                }
 
                ASSERT( ipISO <= ipHE_LIKE );
                for( vector<two_photon>::iterator tnu = sp->TwoNu.begin(); tnu != sp->TwoNu.end(); ++tnu )
                {
                        // induced two photon emission - special because upward and downward are
                        // not related by ratio of statistical weights 
                        // iso.lgInd2nu_On is controlled with SET IND2 ON/OFF command 
 
                        fixit(); // need Pesc or the equivalent to multiply AulTotal?
                        // downward rate
                        z[tnu->ipHi][tnu->ipLo] -= tnu->AulTotal;
                        z[tnu->ipHi][tnu->ipHi] += tnu->AulTotal; 
 
                }
                
                if (iso_ctrl.lgInd2nu_On)
                {
                        for( vector<two_photon>::iterator tnu = sp->TwoNu.begin(); tnu != sp->TwoNu.end(); ++tnu )
                        {
                                // induced two photon emission - special because upward and downward are
                                // not related by ratio of statistical weights 
                                // iso.lgInd2nu_On is controlled with SET IND2 ON/OFF command 
                                
                                z[tnu->ipHi][tnu->ipLo] -= tnu->induc_dn;
                                z[tnu->ipHi][tnu->ipHi] += tnu->induc_dn;
                                
                                // upward rate
                                z[tnu->ipLo][tnu->ipHi] -= tnu->induc_up;
                                z[tnu->ipLo][tnu->ipLo] += tnu->induc_up;
                        }
                }
                
                /* grain charge transfer recombination and ionization to ALL other stages */
                for( long ion=dense.IonLow[nelem]; ion<=dense.IonHigh[nelem]; ++ion )
                {
                        if( ion!=nelem-ipISO )
                        {
                                source += gv.GrainChTrRate[nelem][ion][nelem-ipISO] *
                                        dense.xIonDense[nelem][ion];
                                sink += gv.GrainChTrRate[nelem][nelem-ipISO][ion];
 
                                source += mole.xMoleChTrRate[nelem][ion][nelem-ipISO] * 
                                        dense.xIonDense[nelem][ion] * atmdat.lgCTOn;
                                sink += mole.xMoleChTrRate[nelem][nelem-ipISO][ion] * atmdat.lgCTOn;
 
                        }
                }
                
                /* add in source and sink terms from molecular network. */
                source += mole.source[nelem][nelem-ipISO];
                sink += mole.sink[nelem][nelem-ipISO];
 
#if     1
                /* >>chng 02 Sep 06 rjrw -- all elements have these terms */
                /*>>>chng 02 oct 01, only include if lgAdvection is set */
                if( iteration > dynamics.n_initial_relax+1 && dynamics.lgAdvection && 
                                dynamics.Rate != 0.0 &&
                        !dynamics.lgEquilibrium && dynamics.lgISO[ipISO])
                {
                        /* add in advection - these terms normally zero */
                        source += dynamics.Source[nelem][nelem-ipISO];
                        /* >>chng 02 Sep 06 rjrw -- advective term not recombination */
                        sink += dynamics.Rate;
                }
#else
                /*>>>chng 02 oct 01, only include if lgAdvection is set */
                if( iteration > dynamics.n_initial_relax+1 && dynamics.lgAdvection &&
                                dynamics.Rate != 0.0 &&
                        !dynamics.lgEquilibrium && dynamics.lgISO[ipISO])
                {
                        for( level=0; level < numlevels_local; level++ )
                        {
                                creation[level] += dynamics.StatesElem[nelem][nelem-ipISO][level];
                                z[level][level] += dynamics.Rate;
                        }
                }
#endif
 
                /* ionization from/recombination from lower ionization stages */
                for(long ion_from=dense.IonLow[nelem]; ion_from < MIN2( dense.IonHigh[nelem], nelem-ipISO ) ; ion_from++ )
                {
                        if( ionbal.RateIoniz[nelem][ion_from][nelem-ipISO] >= 0. )
                        {                         
                                /* ionization from lower ionization stages, cm-3 s-1 */
                                source += ionbal.RateIoniz[nelem][ion_from][nelem-ipISO] * dense.xIonDense[nelem][ion_from];
                        }
                        /* recombination to next lower ionization stage, s-1 */
                        if( ion_from == nelem-1-ipISO )
                                sink += ionbal.RateRecomTot[nelem][ion_from];
                }
 
                ASSERT( source >= 0.f );
                if (0)
                {
                        creation[0] += source;
                        for( level=0; level < numlevels_local; level++ )
                        {
                                z[level][level] += sink;
                        }
                }
                else
                {
                        // Try Boltzmann weighting to capture LTE limit correctly
                        t_iso_sp* sp = &iso_sp[ipISO][nelem];
                        double partfun=0.0;
                        for ( level = 0; level < numlevels_local; level++ )
                        {
                                partfun += sp->st[level].Boltzmann()*sp->st[level].g();
                        }
                        source /= partfun;
                        for( level=0; level < numlevels_local; level++ )
                        {
                                creation[level] += source*
                                        sp->st[level].Boltzmann()*sp->st[level].g();
                                z[level][level] += sink;
                        }
                }
                creation[0] += ctsource;
                z[0][0] += ctsink;
 
                /* >>chng 04 nov 30, atom XX-like collisions off turns this off too */
                if( sp->trans(iso_ctrl.nLyaLevel[ipISO],0).Coll().rate_lu_nontherm() * iso_ctrl.lgColl_excite[ipISO] > 0. )
                {
                        /* now add on supra thermal excitation */
                        for( level=1; level < numlevels_local; level++ )
                        {
                                double RateUp , RateDown;
 
                                RateUp = sp->trans(level,0).Coll().rate_lu_nontherm();
                                RateDown = RateUp * (double)sp->st[ipH1s].g() /
                                        (double)sp->st[level].g();
 
                                /* total rate out of lower level */
                                z[ipH1s][ipH1s] += RateUp;
 
                                /* rate from the upper level to ground */
                                z[level][ipH1s] -= RateDown;
 
                                /* rate from ground to upper level */
                                z[ipH1s][level] -= RateUp;
 
                                z[level][level] += RateDown;  
                        }
                }
 
                /* =================================================================== 
                 *
                 * at this point all matrix elements have been established 
                 *
                 * ==================================================================== */
                /* save matrix, this allocates SaveZ */
                SaveZ = z;
 
                for( ipLo=0; ipLo < numlevels_local; ipLo++ )
                        Save_creation[ipLo] = creation[ipLo];
 
                if( trace.lgTrace && trace.lgIsoTraceFull[ipISO] && (nelem == trace.ipIsoTrace[ipISO]) )
                {
                        const long numlevels_print = numlevels_local;
                        fprintf( ioQQQ, "  pop level     others => (iso_level)\n" );
                        for( long n=0; n < numlevels_print; n++ )
                        {
                                fprintf( ioQQQ, "  %s %s %2ld", iso_ctrl.chISO[ipISO], elementnames.chElementNameShort[nelem], n );
                                for( long j=0; j < numlevels_print; j++ )
                                {
                                        fprintf( ioQQQ,"\t%.9e", z[j][n] );
                                }
                                fprintf( ioQQQ, "\n" );
                        }
                        fprintf(ioQQQ," recomb ct %.2e co %.2e hectr %.2e hctr %.2e\n",
                                          atmdat.CharExcRecTotal[ipHELIUM],
                                          findspecieslocal("CO")->den ,
                                          atmdat.CharExcRecTo[ipHELIUM][nelem][nelem-ipISO]*iso_sp[ipHE_LIKE][ipHELIUM].st[0].Pop() ,
                                          atmdat.CharExcRecTo[ipHYDROGEN][nelem][nelem-ipISO]*iso_sp[ipH_LIKE][ipHYDROGEN].st[0].Pop() );
                        fprintf(ioQQQ," recomb          ");
                        for( long n=0; n < numlevels_print; n++ )
                        {
                                fprintf( ioQQQ,"\t%.9e", creation[n] );
                        }
                        fprintf( ioQQQ, "\n" );
                }
 
                nerror = 0;
 
                getrf_wrapper(numlevels_local,numlevels_local,
                                                  z.data(),numlevels_local,&ipiv[0],&nerror);
                
                getrs_wrapper('N',numlevels_local,1,z.data(),numlevels_local,&ipiv[0],
                                                  &creation[0],numlevels_local,&nerror);
 
                if( nerror != 0 )
                {
                        fprintf( ioQQQ, " iso_level: dgetrs finds singular or ill-conditioned matrix\n" );
                        cdEXIT(EXIT_FAILURE);
                }
 
                /* check whether solution is valid */
                /* >>chng 06 aug 28, both of these from numLevels_max to _local. */
                for( level=ipH1s; level < numlevels_local; level++ )
                {
                        double qn = 0., qx = 0.;
                        error[level] = 0.;
                        for( long n=ipH1s; n < numlevels_local; n++ )
                        {
                                double q = SaveZ[n][level]*creation[n];
                                
                                /* remember the largest size of element in sum to div by below */
                                if ( q > qx )
                                        qx = q;
                                else if (q < qn)
                                        qn = q;
 
                                error[level] += q;
                        }
                        
                        if (-qn > qx)
                                qx = -qn;
 
                        if( qx > 0. )
                        {
                                error[level] = (error[level] - Save_creation[level])/qx;
                        }
                        else
                        {
                                error[level] = 0.;
                        }
                }
 
                /* remember largest residual in matrix inversion */
                BigError = -1.;
                level_error = -1;
                /* >>chng 06 aug 28, from numLevels_max to _local. */
                for( level=ipH1s; level < numlevels_local; level++ )
                {
                        double abserror;
                        abserror = fabs( error[level]);
                        /* this will be the largest residual in the matrix inversion */
                        if( abserror > BigError )
                        {
                                BigError = abserror;
                                level_error = level;
                        }
                }
 
                /* matrix inversion should be nearly as good as the accuracy of a double,
                 * but demand that it is better than epsilon for a float */
                if( BigError > FLT_EPSILON ) 
                {
                        if( !conv.lgSearch )
                                fprintf(ioQQQ," PROBLEM" );
 
                        fprintf(ioQQQ,
                                " iso_level zone %.2f - largest residual in iso=%li %s nelem=%li matrix inversion is %g "
                                "level was %li Search?%c \n", 
                                fnzone,
                                ipISO,
                                elementnames.chElementName[nelem],
                                nelem , 
                                BigError , 
                                level_error,
                                TorF(conv.lgSearch) );
                }
                
                // Force level balance to LTE
                if ( iso_ctrl.lgLTE_levels[ipISO] )
                {
                        t_iso_sp* sp = &iso_sp[ipISO][nelem];
                        double partfun=0.0;
                        for ( level = 0; level < numlevels_local; level++ )
                        {
                                partfun += sp->st[level].Boltzmann()*sp->st[level].g();
                        }
                        double scale = dense.xIonDense[nelem][nelem-ipISO]/partfun;
                        for ( level = 0; level < numlevels_local; level++ )
                        {
                                creation[level] = sp->st[level].Boltzmann()*sp->st[level].g()*scale;
                        }
                }
 
                for( level = numlevels_local-1; level > 0; --level )
                {
                        /* check for negative populations */
                        if( creation[level] < 0. )
                                lgNegPop = true;
                }
 
                if( lgNegPop && dense.lgSetIoniz[nelem] )
                {
                        //  simulation can become unphysical if ionization is fixed.
                        //  in this case, just put everything in ground.
                        //  It's really the best we can do.
                        for( level = 1; level < numlevels_local; ++level )
                                creation[level] = 0.;
                        creation[0] = dense.xIonDense[nelem][nelem-ipISO];
                        // flip flag back
                        lgNegPop = false;
                }
 
                /* put level populations into master array */
                for( level=0; level < numlevels_local; level++ )
                {
                        ASSERT( creation[level] >= 0. );
                        sp->st[level].Pop() = creation[level];
 
                        if( sp->st[level].Pop() <= 0 && !conv.lgSearch )
                        {
                                fprintf(ioQQQ,
                                        "PROBLEM non-positive level pop for iso = %li, nelem = "
                                        "%li = %s, level=%li val=%.3e nTotalIoniz %li\n", 
                                        ipISO,
                                        nelem , 
                                        elementnames.chElementSym[nelem],
                                        level,
                                        sp->st[level].Pop() ,
                                        conv.nTotalIoniz);
                        }
                }
 
                /* zero populations of unused levels. */
                for( level=numlevels_local; level < sp->numLevels_max; level++ )
                {
                        sp->st[level].Pop() = 0.;
                        /* >>chng 06 jul 25, no need to zero this out, fix limit to 3-body heating elsewhere. */
                        /* sp->st[level].PopLTE = 0.; */
                }
 
                /* TotalPopExcited is sum of excited level pops */
                /* renormalize the populations to agree with ion solver */
                iso_renorm( nelem, ipISO, renorm );
 
                double TotalPopExcited = 0.;
                /* create sum of populations */
                for( level=1; level < numlevels_local; level++ )
                        TotalPopExcited += sp->st[level].Pop();
                ASSERT( TotalPopExcited >= 0. );
                double TotalPop = TotalPopExcited + sp->st[0].Pop();
 
                /* option to force ionization */
                if( dense.lgSetIoniz[nelem] )
                {
                        if( !fp_equal( TotalPop, (double)dense.xIonDense[nelem][nelem-ipISO] ) )
                        {
                                if( TotalPopExcited >= dense.xIonDense[nelem][nelem-ipISO] )
                                {
                                        for( level=0; level < numlevels_local; level++ )
                                                sp->st[level].Pop() *=
                                                        dense.xIonDense[nelem][nelem-ipISO] / TotalPop;
                                }
                                else
                                {
                                        sp->st[0].Pop() = 
                                                MAX2( 1e-30 * dense.xIonDense[nelem][nelem-ipISO], 
                                                dense.xIonDense[nelem][nelem-ipISO] - TotalPopExcited );
                                }
                                sp->lgPopsRescaled = true;
                        }       
                        ASSERT( sp->st[0].Pop() >= 0. );
                }
        }
        /* all solvers end up here */
 
        /* check on the sum of the populations */
        if( lgNegPop || dense.xIonDense[nelem][nelem-ipISO] < 0. )
        {
                fprintf( ioQQQ, 
                        " DISASTER iso_level finds negative ion fraction for iso=%2ld nelem=%2ld "\
                        "%s using solver %s, IonFrac=%.3e simple=%.3e TE=%.3e ZONE=%4ld\n", 
                        ipISO,
                        nelem,
                        elementnames.chElementSym[nelem],
                        sp->chTypeAtomUsed,
                        dense.xIonDense[nelem][nelem+1-ipISO]/SDIV(dense.xIonDense[nelem][nelem-ipISO]),
                        sp->xIonSimple,
                        phycon.te,
                        nzone );
 
                fprintf( ioQQQ, " level pop are:\n" );
                for( i=0; i < numlevels_local; i++ )
                {
                        fprintf( ioQQQ,PrintEfmt("%8.1e", sp->st[i].Pop() ));
                        if( (i!=0) && !(i%10) ) fprintf( ioQQQ,"\n" );
                }
                fprintf( ioQQQ, "\n" );
                ContNegative();
                ShowMe();
                cdEXIT(EXIT_FAILURE);
        }
 
        for( ipHi=1; ipHi<numlevels_local; ++ipHi )
        {
                for( ipLo=0; ipLo<ipHi; ++ipLo )
                {
                        if( sp->trans(ipHi,ipLo).Emis().Aul() <= iso_ctrl.SmallA )
                                continue;
 
                        /* population of lower level, corrected for stimulated emission */
                        sp->trans(ipHi,ipLo).Emis().PopOpc() = 
                                sp->st[ipLo].Pop() - sp->st[ipHi].Pop()*
                                sp->st[ipLo].g()/sp->st[ipHi].g();
 
                        // don't allow masers from collapsed levels or in Case B
                        if( N_(ipHi) > sp->n_HighestResolved_local || opac.lgCaseB )
                                sp->trans(ipHi,ipLo).Emis().PopOpc() = MAX2( 0., sp->trans(ipHi,ipLo).Emis().PopOpc() );
                }
        }
 
        return;
}
 
void iso_set_ion_rates( long ipISO, long nelem)
{
        DEBUG_ENTRY("iso_set_ion_rates()");
        t_iso_sp* sp = &iso_sp[ipISO][nelem]; 
        const long int numlevels_local = sp->numLevels_local;
        /* this is total ionization rate, s-1, of this species referenced to
         * the total abundance */
        double TotalPop = 0.;
        ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] = 0.;
        for( long level=0; level < numlevels_local; level++ )
        {
                /* sum of all ionization processes from this atom to ion, cm-3 s-1 now,
                 * but is divided by TotalPop below to become s-1 */
                ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] += 
                        sp->st[level].Pop() * sp->fb[level].RateLevel2Cont;
                TotalPop += sp->st[level].Pop();
        }
        
        if( ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] > BIGDOUBLE )
        {
                fprintf( ioQQQ, "DISASTER RateIonizTot for Z=%li, ion %li is larger than BIGDOUBLE.  This is a big problem.",
                                        nelem+1, nelem-ipISO);
                cdEXIT(EXIT_FAILURE);
        }
 
        if (TotalPop <= SMALLFLOAT)
                ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] = sp->fb[0].RateLevel2Cont;
        else
                ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] /= TotalPop;
 
        if( ionbal.RateRecomIso[nelem][ipISO] > 0. )
                sp->xIonSimple = ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1]/ionbal.RateRecomIso[nelem][ipISO];
        else
                sp->xIonSimple = 0.;
 
        ASSERT( ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] >= 0. );
 
        if( ipISO == ipHE_LIKE && nelem==ipHELIUM && nzone>0 )
        {
                /* find fraction of He0 destructions due to photoionization of 2^3S */
                double ratio;
                double rateOutOf2TripS = sp->st[ipHe2s3S].Pop() * sp->fb[ipHe2s3S].RateLevel2Cont;
                if( rateOutOf2TripS > SMALLFLOAT )
                {
                        ratio = rateOutOf2TripS /
                                ( rateOutOf2TripS + sp->st[ipHe1s1S].Pop()*ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] );
                }
                else
                {
                        ratio = 0.;
                }
                if( ratio > he.frac_he0dest_23S )
                {
                        /* remember zone where this happended and fraction, and frac due to photoionization */
                        he.nzone = nzone;
                        he.frac_he0dest_23S = ratio;
                        rateOutOf2TripS = sp->st[ipHe2s3S].Pop() *sp->fb[ipHe2s3S].gamnc;
                        if( rateOutOf2TripS > SMALLFLOAT )
                        {
                                he.frac_he0dest_23S_photo = rateOutOf2TripS /
                                        ( rateOutOf2TripS + sp->st[ipHe1s1S].Pop()*ionbal.RateIoniz[nelem][nelem-ipISO][nelem-ipISO+1] );
                        }
                        else
                        {
                                he.frac_he0dest_23S_photo = 0.;
                        }
                }
        }
}