/home66/gary/public_html/cloudy/c08_branch/source/iso_level.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 */
/*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 "helike.h"
#include "hydrogenic.h"
#include "ionbal.h"
#include "iso.h"
#include "mole.h"
#include "phycon.h"
#include "physconst.h"
#include "rfield.h"
#include "secondaries.h"
#include "taulines.h"
#include "thirdparty.h"
#include "trace.h"

/* solve for level populations  */
void iso_level( const long int ipISO, const long int nelem)
{
        long int ipHi,
                ipLo,
                i,
                level,
                level_error;

  const long int numlevels_local = iso.numLevels_local[ipISO][nelem];

        double BigError;

        int32 nerror;
        double HighestPColOut = 0.,
                collider,
                sum_popn_ov_ion,
                TooSmall;
        bool lgNegPop=false;
        valarray<int32> ipiv(numlevels_local); 
        /* this block of variables will be obtained and freed here */
        valarray<double>
                creation(numlevels_local),
                error(numlevels_local),
                work(numlevels_local),
                Save_creation(numlevels_local);
        double source=0.,
                sink=0.;
        valarray<double> PopPerN(iso.n_HighestResolved_local[ipISO][nelem]+1);

        multi_arr<double,2,C_TYPE> z, SaveZ;

        DEBUG_ENTRY( "iso_level()" );

        /* check that we were called with valid charge */
        ASSERT( nelem >= ipISO );
        ASSERT( nelem < LIMELM );

        /* now do the 2D array */
        z.alloc(numlevels_local,numlevels_local);

        /* >>chng 05 aug 17, must use real electron density for collisional ionization of H
         * since in Leiden v4 pdr with its artificial high temperature coll ion can be important
         * H on H is homonuclear and scaling laws for other elements does not apply 
         * next two were multiplied by dense.EdenHCorr and changed to dense.eden for H,
         * EdenHCorr for rest */
        if( ipISO == ipH_LIKE && nelem==ipHYDROGEN )
        {
                /* special version for H0 onto H0 */
                collider = dense.EdenHontoHCorr;
        }
        else
        {
                collider = dense.EdenHCorr;
        }

        /* fill in recombination vector - values were set in iso_ionize_recombine.cpp */
        for( level=0; level < numlevels_local; level++ )
        {
                /* total recombination to level [s-1] */
                creation[level] = iso.RateCont2Level[ipISO][nelem][level];
        }

        /* these two collision rates must be the same or we are in big trouble,
         * since used interchangably */
        ASSERT( ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0]< SMALLFLOAT ||
                fabs( iso.ColIoniz[ipISO][nelem][0]* collider /
                SDIV(ionbal.CollIonRate_Ground[nelem][nelem-ipISO][0] ) - 1.) < 0.001 );


        if( ipISO == ipH_LIKE )
        {
                if( nelem==ipHYDROGEN )
                        TooSmall = 1e-28;
                else if( nelem==ipHELIUM )
                        TooSmall = 1e-18;
                else
                        TooSmall = 2e-18;
        }
        else if( ipISO == ipHE_LIKE )
        {
                /* >>chng 03 apr 30, different limit for He itself,
                 * since so important in molecular regions */
                if( nelem==ipHELIUM )
                        TooSmall = 1e-20;
                else
                        TooSmall = 1e-10;
        }
        else
                TotalInsanity();

        /* which case atom to solve??? */
        /* >>chng 03 apr 11, from ae-30 bail, to ae-35, will catch it in hydro 2 low */
        if( ipISO == ipH_LIKE &&
                ( iso.xIonSimple[ipISO][nelem] < 1e-35 ) )
        {
                /* don't bother if no ionizing radiation */
                strcpy( iso.chTypeAtomUsed[ipISO][nelem], "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, iso.xIonSimple[ipISO][nelem], iso.chTypeAtomUsed[ipISO][nelem] );
                }

                /* okay to leave as iso.numLevels_max */
                for( i=0; i < iso.numLevels_max[ipISO][nelem]; i++ )
                {
                        iso.DepartCoef[ipISO][nelem][i] = 1.;
                        StatesElem[ipISO][nelem][i].Pop = 0.;
                }

                ionbal.RateRecomTot[nelem][nelem-ipISO] = 0.;
                iso.xIonSimple[ipISO][nelem] = 0.;
                iso.pop_ion_ov_neut[ipISO][nelem] = 0.;
                iso.qTot2S[ipISO][nelem] = 0.;
        }
        /* >>chng 99 nov 23, very dense model close to lte had very low simple
         * ionization fraction but moderate ionization since all pops in
         * excited states - added test for density 
         * second change - for all cases use this routine if ionization is
         * very low indeed */
        /* >>chng 00 aug 26, added matrix force option,
         * logic is to override default options, "none was set in zero */

        /* NB - this test is mostly bogus since chTypeAtom is "POPU" in zero */

        /* >>chng 02 mar 13, iso.chTypeAtomSet had been set to POPU in zero, look for
         * comment with this date.  This had the effect of killing this following test.
         * This also meant that the code had been happily inverting these impossible
         * matrices for quite some time.  This test changed to stringest one, in
         * light of this */
        else if( ipISO == ipH_LIKE &&
                ((strcmp( iso.chTypeAtomSet[ipISO] , "LOWT" )==0) ||
                (iso.xIonSimple[ipISO][nelem] < TooSmall) ) )
        {
                strcpy( iso.chTypeAtomUsed[ipISO][nelem], "Low T" );
                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, iso.xIonSimple[ipISO][nelem], iso.chTypeAtomUsed[ipISO][nelem] );
                }

                /* the ionization ratio WILL be equal to iso.xIonSimple[ipISO][nelem] */
                iso.pop_ion_ov_neut[ipISO][nelem] = iso.xIonSimple[ipISO][nelem];

                /* do simple cascade, should always work, 
                 * but not very well at high densities */
                HydroT2Low( ipISO, nelem );
                iso.qTot2S[ipISO][nelem] = 0.;
        }
        /* branch for very little ionization, use simple estimate but do not do level populations
         * never set to zero, use simple two-level system, since He+ important for chemistry */
        else if( ipISO == ipHE_LIKE &&
                ((strcmp( iso.chTypeAtomSet[ipISO] , "LOWT" )==0) ||
                (iso.xIonSimple[ipISO][nelem] < TooSmall) ) )
        {
                strcpy( iso.chTypeAtomUsed[ipISO][nelem], "Low T" );
                /* total ionization will just be the ground state */
                ionbal.RateIonizTot[nelem][nelem-ipISO] = iso.RateLevel2Cont[ipISO][nelem][0];
                /* >>chng 01 nov 23, increase upper limit by 2 */
                lgNegPop = false;
                /* >>chng 02 apr 10, had set to zero,
                 * we will leave the level populations at zero but define an ion ratio,
                 * since He+ abundance cannot go to zero in a molecular cloud,
                 * - He+ charge transfer is the main CO destruction mechanism */
                iso.pop_ion_ov_neut[ipISO][nelem] = iso.xIonSimple[ipISO][nelem];
                StatesElem[ipISO][nelem][0].Pop =  1./SDIV(iso.pop_ion_ov_neut[ipISO][nelem]);
                for( long n=1; n < numlevels_local; n++ )
                {
                        StatesElem[ipISO][nelem][n].Pop =  0.;
                }
                iso.qTot2S[ipISO][nelem] = 0.;
        }
        else
        {
                ASSERT( NISO == 2 );
                long SpinUsed[NISO] = {2,1};
                long indexNmaxP =
                        iso.QuantumNumbers2Index[ipISO][nelem][ iso.n_HighestResolved_local[ipISO][nelem] ][1][SpinUsed[ipISO]];

                strcpy( iso.chTypeAtomUsed[ipISO][nelem], "Popul" );
                if( trace.lgTrace && (nelem == trace.ipIsoTrace[ipISO]) )
                {
                        fprintf( ioQQQ, "     iso_level iso=%2ld nelem=%2ld doing regular matrix inversion, %s\n", 
                                ipISO, nelem, iso.chTypeAtomUsed[ipISO][nelem] );
                }

                multi_arr<quantumState,3>::const_iterator StElm = StatesElem.begin(ipISO,nelem);

                /* 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] = iso.RateLevel2Cont[ipISO][nelem][level];
                        
                        if( level == 1 )
                                /* >>chng 05 dec 21, rm eden to make into rate coefficient */
                                iso.qTot2S[ipISO][nelem] = iso.ColIoniz[ipISO][nelem][level];

                        multi_arr<transition,4>::const_iterator Trans = Transitions.begin(ipISO,nelem,level);
                        md4ci Boltz = iso.Boltzmann.begin(ipISO,nelem,level);

                        /* all processes populating level from below */
                        for( ipLo=0; ipLo < level; ipLo++ )
                        {
                                double coll_down = Trans[ipLo].Coll.ColUL * collider;

                                double RadDecay = MAX2( iso.SmallA, Trans[ipLo].Emis->Aul*
                                                        (Trans[ipLo].Emis->Pesc + 
                                                         Trans[ipLo].Emis->Pelec_esc + 
                                                         Trans[ipLo].Emis->Pdest)*
                                                        KILL_BELOW_PLASMA(Trans[ipLo].EnergyWN*WAVNRYD) );

                                double pump = MAX2( iso.SmallA, Trans[ipLo].Emis->pump *
                                                    KILL_BELOW_PLASMA(Trans[ipLo].EnergyWN*WAVNRYD) );

                                if( iso.lgRandErrGen[ipISO] )
                                {
                                        coll_down *= iso.ErrorFactor[ipISO][nelem][level][ipLo][IPCOLLIS];
                                        RadDecay *= iso.ErrorFactor[ipISO][nelem][level][ipLo][IPRAD];
                                        pump *= iso.ErrorFactor[ipISO][nelem][level][ipLo][IPRAD];
                                }

                                double coll_up = coll_down * 
                                        (double)StElm[level].g/
                                        (double)StElm[ipLo].g*
                                        Boltz[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) )
                                {
                                        iso.qTot2S[ipISO][nelem] += coll_down/dense.EdenHCorr;
                                }
                                if( (ipLo == 1) && (ipISO==ipH_LIKE || StElm[level].S==1) )
                                {
                                        iso.qTot2S[ipISO][nelem] += coll_up/dense.EdenHCorr;
                                }
                        }
                }

                if( ipISO == nelem )
                {
                        /* iso.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./StatesElem[ipISO][nelem][indexNmaxP].lifetime )
                        {
                                iso.lgCritDensLMix[ipISO] = false;
                        }
                }

                ASSERT( ipISO <= ipHE_LIKE );

                /* 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[1+ipISO][0] -= iso.TwoNu_induc_dn[ipISO][nelem]*iso.lgInd2nu_On;
                z[0][1+ipISO] -= iso.TwoNu_induc_up[ipISO][nelem]*iso.lgInd2nu_On;

                /* rates out of 1s, and out of 2s */
                z[0][0] += iso.TwoNu_induc_up[ipISO][nelem]*iso.lgInd2nu_On;
                z[1+ipISO][1+ipISO] += iso.TwoNu_induc_dn[ipISO][nelem]*iso.lgInd2nu_On;

                /* grain charge transfer recombination and ionization to ALL other stages */
                for( long ion=0; ion<=nelem+1; ++ion )
                {
                        if( ion!=nelem-ipISO )
                        {
                                /* recombination must be multiplied by a ratio of densities to get proper rate. */
                                source += gv.GrainChTrRate[nelem][ion][nelem-ipISO] * 
                                        dense.xIonDense[nelem][ion] ;
                                /* take ionization out of every level. */
                                sink += gv.GrainChTrRate[nelem][nelem-ipISO][ion];
                        }
                }

                /* >>chng 02 Sep 06 rjrw -- all elements have these terms */
                /*>>>chng 02 oct 01, only include if lgAdvection is set */
                if( dynamics.lgAdvection && dynamics.lgISO[ipISO])
                {
                        /* add in advection - these terms normally zero */
                        /* assume for present that all advection is into ground state */
                        source += dynamics.Source[nelem][nelem-ipISO];
                        /* >>chng 02 Sep 06 rjrw -- advective term not recombination */
                        /* can sink from all components (must do, for conservation) */
                        sink += dynamics.Rate;
                }

#if     0
                /* add in source and sink terms from molecular network. */
                CodeReview(); /* Check... */
                if( nelem == ipISO && conv.nTotalIoniz )
                {
                        /* these are the external source and sink terms */
                        /* source first */
                        source += mole.source[nelem][nelem-ipISO];
                        sink += mole.sink[nelem][nelem-ipISO];
                }
#endif

                creation[0] += source/SDIV(dense.xIonDense[nelem][nelem+1-ipISO]);
                for( level=0; level < numlevels_local; level++ )
                {
                        z[level][level] += sink;
                }

                /* >>chng 04 nov 30, atom XX-like collisions off turns this off too */
                if( secondaries.Hx12[ipISO][nelem][iso.nLyaLevel[ipISO]] * iso.lgColl_excite[ipISO] > 0. )
                {
                        /* now add on supra thermal excitation */
                        for( level=1; level < numlevels_local; level++ )
                        {
                                double RateUp , RateDown;

                                RateUp = secondaries.Hx12[ipISO][nelem][level];
                                RateDown = RateUp * (double)StatesElem[ipISO][nelem][ipH1s].g /
                                        (double)StatesElem[ipISO][nelem][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]) )
                {
                        fprintf( ioQQQ, "  pop level     others => (iso_level)\n" );
                        for( long n=0; n < numlevels_local; n++ )
                        {
                                fprintf( ioQQQ, "  %s %s %2ld", iso.chISO[ipISO], elementnames.chElementNameShort[nelem], n );
                                for( long j=0; j < numlevels_local; j++ )
                                {
                                        fprintf( ioQQQ,"\t%.9e", z[j][n] );
                                }
                                fprintf( ioQQQ, "\n" );
                        }
                        fprintf(ioQQQ," recomb          ");
                        for( long n=0; n < numlevels_local; n++ )
                        {
                                fprintf( ioQQQ,"\t%.9e", creation[n] );
                        }
                        fprintf( ioQQQ, "\n" );
                        fprintf(ioQQQ," recomb ct %.2e co %.2e hectr %.2e hctr %.2e\n",
                                atmdat.HeCharExcRecTotal,
                                findspecies("CO")->hevmol ,
                                atmdat.HeCharExcRecTo[nelem][nelem-ipISO]*dense.xIonDense[ipHELIUM][0] ,
                                atmdat.HCharExcRecTo[nelem][nelem-ipISO]*dense.xIonDense[ipHYDROGEN][0] );
                }

                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 BigSoln= 0.;
                        error[level] = 0.;
                        for( long n=ipH1s; n < numlevels_local; n++ )
                        {
                                /* remember the largest value of the soln matrix to div by below */
                                if( fabs(SaveZ[n][level] ) > BigSoln )
                                        BigSoln = fabs(SaveZ[n][level]);

                                error[level] += SaveZ[n][level]*creation[n];
                        }

                        if( BigSoln > 0. )
                        {
                                error[level] = (error[level] - Save_creation[level])/ BigSoln;
                        }
                        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) );
                }

                /* put departure coefficients and level populations into master array */
                for( level=0; level < numlevels_local; level++ )
                {
                        StatesElem[ipISO][nelem][level].Pop = creation[level];

                        /* check for negative level populations */
                        if( StatesElem[ipISO][nelem][level].Pop < 0. )
                                lgNegPop = true;

                        if( StatesElem[ipISO][nelem][level].Pop <= 0 )
                        {
                                fprintf(ioQQQ," non-positive level pop for iso = %li, nelem = %li = %s, level=%li val=%.3e\n", 
                                        ipISO,
                                        nelem , 
                                        elementnames.chElementSym[nelem],
                                        level,
                                        StatesElem[ipISO][nelem][level].Pop );
                        }

                        if( iso.PopLTE[ipISO][nelem][level] > 0. )
                        {
                                /* the LTE departure coefficients */
                                iso.DepartCoef[ipISO][nelem][level] = 
                                        StatesElem[ipISO][nelem][level].Pop/
                                        (iso.PopLTE[ipISO][nelem][level]* dense.eden);
                        }
                        else
                        {
                                iso.DepartCoef[ipISO][nelem][level] = 0.;
                        }
                }

                /* zero populations of unused levels. */
                for( level=numlevels_local; level < iso.numLevels_max[ipISO][nelem]; level++ )
                {
                        StatesElem[ipISO][nelem][level].Pop = 0.;
                        /* >>chng 06 jul 25, no need to zero this out, fix limit to 3-body heating elsewhere. */
                        /* iso.PopLTE[ipISO][nelem][level] = 0.; */
                }
        }

        /* all solvers end up here */
        /* sum_popn_ov_ion will become the ratio of iso to parent
         * species, create sum of level pops per ion first */
        sum_popn_ov_ion = 0.;

        /* this is total ionization of this species referenced to the ground state */
        ionbal.RateIonizTot[nelem][nelem-ipISO] = 0.;

        for( level=0; level < numlevels_local; level++ )
        {
                /* sum of all ionization processes from this atom to ion */
                ionbal.RateIonizTot[nelem][nelem-ipISO] += 
                        StatesElem[ipISO][nelem][level].Pop * iso.RateLevel2Cont[ipISO][nelem][level];

                /* create sum of populations relative to ion */
                sum_popn_ov_ion += StatesElem[ipISO][nelem][level].Pop;
        }

        /* convert back to scaled from ground */
        if( ( ionbal.RateIonizTot[nelem][nelem-ipISO]/MAX2(SMALLDOUBLE , sum_popn_ov_ion) ) > 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);
        }
        else
                ionbal.RateIonizTot[nelem][nelem-ipISO] /= MAX2(SMALLDOUBLE , sum_popn_ov_ion);

        /* this is sum of all populations in model, div by ion 
         * create ratio of ion pops to total atom */
        if( sum_popn_ov_ion >= 1e-30 )
        {
                iso.pop_ion_ov_neut[ipISO][nelem] = 1./sum_popn_ov_ion;
        }
        else
        {
                iso.pop_ion_ov_neut[ipISO][nelem] = 0.;

                /* >>chng 03 aug 16, zero out parent density - had not been done */
                dense.xIonDense[nelem][nelem+1-ipISO] = 0.;

                /* reset pointer to one lower stage of ionization so this not
                 * considered again, hydrogenic considered if IonHigh is nelem+2 */
                dense.IonHigh[nelem] = nelem;
                ionbal.RateIonizTot[nelem][nelem-ipISO] = 0.;

                /* now zero this out */
                /* okay to leave as iso.numLevels_max */
                for( ipHi=0; ipHi < iso.numLevels_max[ipISO][nelem]; ipHi++ )
                {
                        /* population of level relative to ion */
                        StatesElem[ipISO][nelem][ipHi].Pop = 0.;
                }
        }

        ASSERT( ionbal.RateIonizTot[nelem][nelem-ipISO] >= 0. );

        /* check on the sum of the populations */
        if( lgNegPop || iso.pop_ion_ov_neut[ipISO][nelem] < 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],
                        iso.chTypeAtomUsed[ipISO][nelem],
                        iso.pop_ion_ov_neut[ipISO][nelem],
                        iso.xIonSimple[ipISO][nelem],
                        phycon.te,
                        nzone );

                fprintf( ioQQQ, " level pop are:\n" );
                for( i=0; i < numlevels_local; i++ )
                {
                        fprintf( ioQQQ,PrintEfmt("%8.1e", StatesElem[ipISO][nelem][i].Pop ));
                        if( (i!=0) && !(i%10) ) fprintf( ioQQQ,"\n" );
                }
                fprintf( ioQQQ, "\n" );
                ContNegative();
                ShowMe();
                cdEXIT(EXIT_FAILURE);
        }

        if( ipISO == ipHE_LIKE && nelem==ipHELIUM && nzone>0 )
        {
                /* find fraction of He0 destructions due to photoionization of 2^3S */
                double ratio;
                double rateOutOf2TripS = StatesElem[ipISO][nelem][ipHe2s3S].Pop * iso.RateLevel2Cont[ipISO][nelem][ipHe2s3S];
                if( rateOutOf2TripS > SMALLFLOAT )
                {
                        ratio = rateOutOf2TripS /
                                ( rateOutOf2TripS + StatesElem[ipISO][nelem][ipHe1s1S].Pop*ionbal.RateIonizTot[nelem][nelem-ipISO] );
                }
                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 = StatesElem[ipISO][nelem][ipHe2s3S].Pop *iso.gamnc[ipISO][nelem][ipHe2s3S];
                        if( rateOutOf2TripS > SMALLFLOAT )
                        {
                                he.frac_he0dest_23S_photo = rateOutOf2TripS /
                                        ( rateOutOf2TripS + StatesElem[ipISO][nelem][ipHe1s1S].Pop*ionbal.RateIonizTot[nelem][nelem-ipISO] );
                        }
                        else
                        {
                                he.frac_he0dest_23S_photo = 0.;
                        }
                }
        }

        for( ipHi=1; ipHi<numlevels_local; ++ipHi )
        {
                for( ipLo=0; ipLo<ipHi; ++ipLo )
                {
                        if( Transitions[ipISO][nelem][ipHi][ipLo].Emis->Aul <= iso.SmallA )
                                continue;

                        /* population of lower level rel to ion, corrected for stim em */
                        Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc = 
                                StatesElem[ipISO][nelem][ipLo].Pop - StatesElem[ipISO][nelem][ipHi].Pop*
                                StatesElem[ipISO][nelem][ipLo].g/StatesElem[ipISO][nelem][ipHi].g;

                        /* don't allow masers from collapsed levels */
                        if( N_(ipHi) > iso.n_HighestResolved_local[ipISO][nelem] )
                                Transitions[ipISO][nelem][ipHi][ipLo].Emis->PopOpc = StatesElem[ipISO][nelem][ipLo].Pop;
                }
        }

        return;
}

---