/home66/gary/public_html/cloudy/c08_branch/source/conv_eden_ioniz.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 */
/*ConvEdenIoniz called by ConvTempIonz, calls ConvIoniz solving for eden */
/*lgConvEden returns true if electron density is converged */
#include "cddefines.h"
#include "dense.h"
#include "trace.h"
#include "thermal.h"
#include "thirdparty.h"
#include "phycon.h"
#include "conv.h"
/* limit to how many loops we will do */
/* >>chng 04 sep 27, from 25 to 35 */
static const int LOOPMAX = 35;

/*ConvEdenIoniz called by ConvTempIonz, calls ConvIoniz solving for eden 
 * returns 0 if ok, 1 if abort */
int ConvEdenIoniz(void)
{
        long int LoopEden ,
                /* this will be LOOPMAX at first, then doubled if no oscillations occur*/
                LoopLimit ,
                LoopBig;
        bool lgFailLastTry;
        int kase = -1;
        static bool lgSlope_oscil=false;

        static bool lgEden_ever_oscil = false;
        static double eden_assumed_old ,
                eden_assumed_minus_true_old,
                eden_assumed_minus_true,
                slope_ls=0., slope_ls_lastzone=0.;
        static bool lgEvalSlop_ls=0;
        double EdenEntry;
        static bool lgMustMalloc_temp_history = true;
        static double slope=-1.;
        double 
                /* upper and lower limits to the range of the change we want to consider */
                EdenUpperLimit,
                EdenLowerLimit ,
                change_eden_rel2_tolerance ,
                PreviousChange;

        DEBUG_ENTRY( "ConvEdenIoniz()" );

        /* this routine is called by ConvTempIonz, it calls ConvIonize
         * and changes the electron density until it converges */

        /* save entry value of eden */
        EdenEntry = dense.eden;

        if( trace.lgTrace )
        {
                fprintf( ioQQQ, "  \n" );
                fprintf( ioQQQ, " ConvEdenIoniz entered \n" );
        }

        if( trace.nTrConvg>=3 )
        {
                fprintf( ioQQQ, 
                        "   ConvEdenIoniz3 called, entering eden loop using solver %sw.\n",
                        conv.chSolverEden);
        }

        /* keep track of temperature solver in this zone - 
         * conv.nTotalIoniz is zero in first call of new iteration */
        if( !conv.nTotalIoniz )
        {
                /* one time initialization of variables */
                conv.hist_temp_ntemp = -1;
                conv.hist_temp_nzone = -1;
                /*>>chng 06 jun 25, do not reset upper limit to number of
                 * points saved - only set this when malloc or remalloc */
                /*conv.hist_temp_limit = 0;*/
        }
        /* this will save the  history of density - pressure relationship
         * for the current zone */
        if( nzone!=conv.hist_temp_nzone )
        {
                /* first time in this zone - reset counters */
                conv.hist_temp_nzone = nzone;
                /* this counter will be updated after vars are saved so will be
                 * total number of saved values */
                conv.hist_temp_ntemp = 0;
        }
        /* do we need to allocate, or reallocate, memory for the vectors
         * check initial allocation first */
        /*>>chng 06 jun 25, do not test on this limit for need to malloc
         * rather set static flag in this routine - fix of memory leak
         * discovered by PvH - routine entered with conv.nTotalIoniz zero
         * during first call of each iteration so memory was malloced without
         * freeing previous memory - memory leak resulted */
        /*if( conv.hist_temp_limit==0 )*/
        if( lgMustMalloc_temp_history )
        {
                lgMustMalloc_temp_history = false;
                conv.hist_temp_limit = 200;
                conv.hist_temp_temp = (double *)MALLOC( (size_t)((unsigned long)conv.hist_temp_limit*sizeof(double) ) );
                conv.hist_temp_heat = (double *)MALLOC( (size_t)((unsigned long)conv.hist_temp_limit*sizeof(double) ) );
                conv.hist_temp_cool = (double *)MALLOC( (size_t)((unsigned long)conv.hist_temp_limit*sizeof(double) ) );
                conv.chNotConverged = (char *)MALLOC( (size_t)((unsigned long)INPUT_LINE_LENGTH*sizeof(char) ) );
                conv.chConvEden = (char *)MALLOC( (size_t)((unsigned long)INPUT_LINE_LENGTH*sizeof(char) ) );
                conv.chConvIoniz = (char *)MALLOC( (size_t)((unsigned long)INPUT_LINE_LENGTH*sizeof(char) ) );
                strcpy( conv.chNotConverged, "none" );
                strcpy( conv.chConvEden, "none" );
                strcpy( conv.chConvIoniz, "none" );
        }
        /* ran out of space - allocate more */
        if( (conv.hist_temp_ntemp+1) >= conv.hist_temp_limit )
        {
                conv.hist_temp_limit *= 3;
                conv.hist_temp_temp = (double *)REALLOC( conv.hist_temp_temp, (size_t)((unsigned long)conv.hist_temp_limit*sizeof(double)));
                conv.hist_temp_heat = (double *)REALLOC( conv.hist_temp_heat, (size_t)((unsigned long)conv.hist_temp_limit*sizeof(double)));
                conv.hist_temp_cool = (double *)REALLOC( conv.hist_temp_cool, (size_t)((unsigned long)conv.hist_temp_limit*sizeof(double)));
        }

        /* >>chng 04 feb 11, add option to remember current density and pressure */
        conv.hist_temp_temp[conv.hist_temp_ntemp] = phycon.te;
        conv.hist_temp_heat[conv.hist_temp_ntemp] = thermal.htot;
        conv.hist_temp_cool[conv.hist_temp_ntemp] = thermal.ctot;
        ++conv.hist_temp_ntemp;

        if( conv.nPres2Ioniz < 5 )
                lgEden_ever_oscil = false;

        /* this flag says wether we converged the density but failed after final tweek,
         * need to use smaller steps */
#       define PRT_FAIL_LAST_TRY        false
        lgFailLastTry = false;
        if( PRT_FAIL_LAST_TRY ) fprintf(ioQQQ,"lgFailLastTry set false\n");

        /* >>chng 03 mar 17, add this big loop, only go through one time if fully
         * converged at end, but if eden jumps off in final touchup,
         * do whole thing again */
        LoopBig = 0;
        while( LoopBig==0 || (lgFailLastTry && LoopBig==1 ) )
        {
                /* the old default solver */
                if( strcmp( conv.chSolverEden , "simple" )== 0 )
                {
                        static double damp;
                        long int nLoopOscil;
                        LoopEden = 0;
                        conv.nConvIonizFails = 0;

                        /* this will be increased by 2x if, at the end, no oscillations have occurred */
                        /* >>chng 04 aug 04, from 2x to 4x, to following through on eden front */
                        LoopLimit = LOOPMAX;

                        /* will be set true if sign of change, ever changes */
                        conv.lgEdenOscl = false;
                        nLoopOscil = 0;
#                       define PRT_EDEN_OSCIL   false
                        if( PRT_EDEN_OSCIL ) fprintf(ioQQQ," PRT_EDEN_OSCIIL setting lgEdenOscl FALSE2\n");
                        lgFailLastTry = false;
                        if( PRT_FAIL_LAST_TRY ) fprintf(ioQQQ,"lgFailLastTry set false\n");

                        /* will be used to look for oscillations in the electron density */
                        PreviousChange = 0.;

                        strcpy( conv.chConvIoniz, " NONE!!!!!" );

                        /* we have not yet  */
                        eden_assumed_old = 0.;
                        eden_assumed_minus_true_old = 0.;

                        damp = 1.;
                        /* if working with zone 0 (nzone=0) reset slope to zero */
                        if( !nzone )
                                slope = 0.;

                        /* this is electron density convergence loop */
                        do 
                        {
                                double slopenew , slopeold;
                                /* this flag will be set false below if electron density not within tolerance 
                                 * after ionization is recomputed */
                                conv.lgConvEden = true;

                                /* converge the current level of ionization, this calls eden_sum which updates EdenTrue */
                                if( ConvIoniz() || lgAbort )
                                {
                                        return 1;
                                }
                                /* this is the new error in the electron density, the difference between
                                 * the assumed and true electron density - we want this to be zero */
                                eden_assumed_minus_true = dense.eden - dense.EdenTrue;

                                {
                                        static long int limEdenHistory=1000;
                                        static bool lgMustMalloc_eden_error_history = true;
                                        bool lgHistUpdate=false;
                                        static double *eden_error_history=NULL, *eden_assumed_history=NULL ,
                                                *stdarray;
                                        static long int nEval=-1 , nzoneUsed=-1, iterUsed=-1;
                                        if( lgMustMalloc_eden_error_history )
                                        {
                                                lgMustMalloc_eden_error_history = false;
                                                lgSlope_oscil = false;
                                                eden_error_history = 
                                                        (double*)MALLOC(sizeof(double)*(unsigned)limEdenHistory );
                                                eden_assumed_history = 
                                                        (double*)MALLOC(sizeof(double)*(unsigned)limEdenHistory );
                                                stdarray = 
                                                        (double*)MALLOC(sizeof(double)*(unsigned)limEdenHistory );
                                        }
                                        if( nzoneUsed!=nzone || iterUsed!=iteration )
                                        {
                                                /* first evaluation this zone */
                                                if( nzone==0 )
                                                {
                                                        lgEvalSlop_ls = true;
                                                        slope_ls = 1.;
                                                }
                                                /* reset some vars at start of new zone */
                                                nEval = 0;
                                                nzoneUsed = nzone;
                                                iterUsed = iteration;
                                                lgSlope_oscil = false;
                                                slope_ls_lastzone = slope_ls;
                                        }
                                        /* do not evaluate this during search phase or if slope is oscillating */
                                        else if( !conv.lgSearch && !lgSlope_oscil )
                                        {
                                                /* may cycle through here with exactly same edensity before tripping if  */
                                                int ip = MAX2(0,nEval-1);
                                                if( 
                                                        /* no need to do this if already converged */
                                                        fabs(dense.eden - dense.EdenTrue)/dense.eden > conv.EdenErrorAllowed  &&
                                                        /* do it if never evaluated */
                                                        (!nEval ||
                                                        /* this test - don't want tiny corrections, within the tolerance, which
                                                         * can have noise since at the level code is converging, affecting the slope */
                                                        fabs((dense.eden - dense.EdenTrue)-eden_error_history[ip]*dense.gas_phase[ipHYDROGEN] )/dense.EdenTrue > 1e-5 ||
                                                        fabs(dense.eden-eden_assumed_history[ip]*dense.gas_phase[ipHYDROGEN])/ dense.EdenTrue > 1e-5 ))
                                                {
                                                        /* use relative fraction for constant pressure case */
                                                        eden_error_history[nEval] = (dense.eden - dense.EdenTrue)/dense.gas_phase[ipHYDROGEN];
                                                        eden_assumed_history[nEval] = dense.eden/dense.gas_phase[ipHYDROGEN];
                                                        stdarray[nEval] = 0.;
                                                        /* >>chng 05 feb 22, only update nEval if residuals are non-zero -
                                                         * this happens in constant temperature models - happened in optimization run */
                                                        if( eden_error_history[nEval] != 0. )
                                                        {
                                                                ++nEval;
                                                                lgHistUpdate = true;
                                                        }
                                                        /*fprintf(ioQQQ,"DEBUG history\t%li\t%e\t%e\t%e\n",
                                                                nEval,
                                                                dense.gas_phase[ipHYDROGEN],
                                                                eden_assumed_history[nEval],
                                                                eden_error_history[nEval]);*/
                                                }
                                                if( nEval>=limEdenHistory )
                                                {
                                                        /* used a lot of points - must create more space */
                                                        limEdenHistory *=10;
                                                        eden_error_history = (double *)REALLOC( eden_error_history, (unsigned)(limEdenHistory)*sizeof(double));
                                                        eden_assumed_history = (double *)REALLOC( eden_error_history, (unsigned)(limEdenHistory)*sizeof(double));
                                                        stdarray = (double *)REALLOC( eden_error_history, (unsigned)(limEdenHistory)*sizeof(double));
                                                }
                                                if( nEval>3 && lgHistUpdate )
                                                {
                                                        double y_intercept, stdx, stdy, slope_save;
                                                        slope_save = slope_ls;
                                                        /* enough data to do least squares */
                                                        if( linfit(nEval,
                                                                   eden_assumed_history, 
                                                                   eden_error_history, 
                                                                   y_intercept, 
                                                                   stdy, 
                                                                   slope_ls, 
                                                                   stdx) )
                                                        {
                                                                int i;
                                                                fprintf(ioQQQ," ConvEdenIoniz: linfit returns impossible condition, history follows:\n");
                                                                fprintf(ioQQQ,"eden_error_history\tstdarray\n");
                                                                for( i=0; i<nEval; ++i )
                                                                {
                                                                        fprintf(ioQQQ,"%.3e\t%.4e\n",
                                                                                eden_error_history[i] , 
                                                                                stdarray[i] );
                                                                }
                                                                fprintf(ioQQQ,"\n");
                                                                ShowMe();
                                                                cdEXIT(EXIT_FAILURE);
                                                        }
                                                        lgEvalSlop_ls = true;
                                                        /* check for slope changing sign - if does, use last zone's slope */
                                                        if( slope_ls*slope_save<0. )
                                                        {
                                                                slope_ls = slope_ls_lastzone;
                                                                lgSlope_oscil = true;
                                                        }
                                                        /*fprintf(ioQQQ,"DEBUG LS\t%.2f\t%li\t%.3e\t%.3e\t%.3e\n",
                                                                fnzone ,
                                                                nEval ,
                                                                slope_ls,
                                                                y_intercept , 
                                                                slope);*/
                                                }
                                        }
                                }
                                slopeold = slope;
                                /* may not be set below, but could print it. set to zero as flag for this case */
                                slopenew = 0.;
                                if( fabs(eden_assumed_minus_true_old) > 0. )
                                {
                                        if( fabs( eden_assumed_old-dense.eden ) > SMALLFLOAT )
                                        {
#                                               define  OLDFAC  0.75
                                                slopenew = (eden_assumed_minus_true_old - eden_assumed_minus_true) / (eden_assumed_old-dense.eden );
                                                /* >>chng 04 sep 15 do not update slope if changing sign */
#                                               if 0
                                                if( slope !=0. && slope*slopenew>=0.)
                                                        slope = OLDFAC*slope + (1.-OLDFAC)*slopenew;
                                                else if( slope==0.)
                                                        slope = slopenew;
#                                               endif
                                                if( slope !=0. )
                                                        slope = OLDFAC*slope + (1.-OLDFAC)*slopenew;
                                                else 
                                                        slope = slopenew;
#                                               undef   OLDFAC
                                        }
                                        /*else
                                                slope = 0.;*/
                                }
                                /* following block is to not let electron density change by
                                * too much
                                * conv.EdenErrorAllowed is allowed error
                                * the default value of conv.EdenErrorAllowed is 0.01 in zerologic
                                * is changed with the SET Eden Error command
                                * eden_assumed_old was incoming value of eden
                                * EdenTrue is correct value from eden_sum for current ionization
                                * new eden will be set using these, trying to prevent oscillations */

                                /* is error larger than the tolerance we are trying to beat? */
                                if( fabs(dense.eden-dense.EdenTrue)/SDIV(dense.eden) >= 
                                        conv.EdenErrorAllowed )
                                {
                                        double eden_proposed;
                                        /* difference is assumed and true electron density is larger
                                        * than tolerance, we are not yet converged, default is 0.01 */
                                        conv.lgConvEden = false;
                                        strcpy( conv.chConvIoniz, "Ne big chg" );

                                        /* these are needed since will take log below */
                                        ASSERT( dense.eden > SMALLFLOAT );

                                        /* SEARCH is true if this is only search for first solution */
                                        /* EdenTrue is allowed to go negative during search for soln - not physical,
                                         * but can happen in math search for soln */
                                        if( conv.lgSearch && dense.EdenTrue > SMALLFLOAT )
                                        {
                                                dense.eden = pow(10.,
                                                        (log10(dense.eden) + log10(dense.EdenTrue))/ 2.);
                                        }
                                        else
                                        {
                                                /* >>chng 03 jul 07, add case for grains contain significant
                                                 * fraction of electrons */
                                                /* >>chng 03 nov 28, add this test on fraction of electrons from ions,
                                                 * this branch is to identify limit where molecules and grains have
                                                 * most of the electrons */
                                                /* >>chng 04 sep 14, change from all ions to just metals */
                                                /* this was patch to fix steep slope */
                                                if( 0 && dense.eden_from_metals > 0.5 )
                                                {
                                                        static long int nzone_eval=-1;
                                                        static bool lgOscilkase10 = false;
                                                        /* this flag says be careful, if kase is this then
                                                         * small changes in assumed eden result in large changes
                                                         * in returned eden */
#                                                       define  KASE_EDEN_NOT_ION       10
                                                        if( nzone != nzone_eval )
                                                        {
                                                                nzone_eval = nzone;
                                                                lgOscilkase10 = false;
                                                        }

                                                        /* >>chng 03 dec 25, check whether oscillations have occurred */
                                                        if( PreviousChange*(dense.EdenTrue-dense.eden) < 0. && nzone_eval > 6 )
                                                                lgOscilkase10 = true;

                                                        /* few of the electrons are due to ions, most are grains and
                                                        * or molecules, be very careful */
                                                        change_eden_rel2_tolerance = 0.2;
                                                        if( lgOscilkase10 )
                                                        {
                                                                /* eden is oscillating, make changes small */
                                                                change_eden_rel2_tolerance = 0.05;
                                                        }
                                                        kase = KASE_EDEN_NOT_ION;
                                                }
                                                /* was the sign of the change in the electron density changing,
                                                 * or has it ever changed? Also use if we are close to converged? */
                                                /* >>chng 04 aug 04, rm test on eden converged, since fooled solver when
                                                 * initially very close to soln, but near eden front, stopped aggressive
                                                 * pursuit of true eden, pdr_leiden_v2 */
                                                else if( PreviousChange*(dense.EdenTrue-dense.eden) < 0. )
                                                {
                                                        /* remember that oscillations are happening  */
                                                        conv.lgEdenOscl = true;
                                                        nLoopOscil = LoopEden;
                                                        /* turn this back on 04 ep 26 */
                                                        /* >>chng 04 sep 27, from * 0.5 to * 0.75 */
                                                        damp = MAX2( 0.05, damp * 0.75 );
                                                        if( PRT_EDEN_OSCIL ) 
                                                                fprintf(ioQQQ," PRT_EDEN_OSCIIL setting lgEdenOscl true, previous change %e new change %e\n",
                                                                PreviousChange,dense.EdenTrue-dense.eden);

                                                        /* changes in electron density are changing sign - be careful,
                                                        * change_eden_rel2_tolerance multiplies the error tolerance on the electron density,
                                                        * largest change allowed is related to this */
                                                        change_eden_rel2_tolerance = 0.5;
                                                        /* >>chng 04 aug 09, from 0.5 to 0.25 */
                                                        change_eden_rel2_tolerance = 0.25;
                                                        kase = 0;
                                                }
                                                else if( lgFailLastTry )
                                                {
                                                        /* this is case where last evaluation of eden, after solution, 
                                                         * caused failure */
                                                        change_eden_rel2_tolerance = 0.5;
                                                        kase = 1;
                                                }
                                                /* >>chng 03 apr 22, add this branch 
                                                 * >>chng 04 feb 23, replace conv.lgEdenOscl with lgEden_ever_oscil */
                                                else if( LoopEden > 4 && !lgEden_ever_oscil &&
                                                        fabs( (dense.EdenTrue-dense.eden)/dense.eden) > 3.*conv.EdenErrorAllowed )
                                                {
                                                        /* we have gone a few times around this loop, the electron density is
                                                        * not oscillating, and we have a long ways to go.  Open up the 
                                                        * allowed change */
                                                        change_eden_rel2_tolerance = 3.;
                                                        kase = 2;
                                                }
                                                else if( conv.lgEdenOscl )
                                                {
                                                        /* >>chng 04 aug 14, add this to stop oscil in primal.in */
                                                        /* oscillations have occurred */
                                                        change_eden_rel2_tolerance = 0.25;
                                                        kase = 4;
                                                }
                                                else
                                                {
                                                        /* stable so far . .. */
                                                        /* change_eden_rel2_tolerance = 2.;*/
                                                        /* had been 2, change to 1 stopped oscillations from developing in 
                                                        * loc blr grid */
                                                        change_eden_rel2_tolerance = 1.;
                                                        kase = 3;
                                                }

                                                /* now remember this change for the next time through the loop  */
                                                PreviousChange = dense.EdenTrue - dense.eden;

                                                /* old difference between assumed and trye */
                                                eden_assumed_minus_true_old = eden_assumed_minus_true;

                                                /* remember electron density before we update it */
                                                eden_assumed_old = dense.eden;

                                                /* is the correct electron density - is it too different? 
                                                * default on conv.EdenErrorAllowed is 0.01, dyanmics is 0.001 */
                                                EdenUpperLimit = dense.eden * (1. + damp*conv.EdenErrorAllowed*change_eden_rel2_tolerance);
                                                EdenLowerLimit = dense.eden / (1. + damp*conv.EdenErrorAllowed*change_eden_rel2_tolerance);

#                                               define USE_NEW  true
#                                               define PRT_NEW  false
                                                if( PRT_NEW && USE_NEW )
                                                        fprintf(ioQQQ,"DEBUG slope\t%li\t%li\t%.3f\t%.3f\t%.3f\t%.4e\t%.4e\t%.4f\t%.3f\t%.3e\n",
                                                        nzone, 
                                                        conv.nPres2Ioniz,
                                                        slope , 
                                                        slopeold , 
                                                        slopenew ,
                                                        dense.eden, 
                                                        dense.EdenTrue, 
                                                        (dense.eden-dense.EdenTrue)/SDIV(dense.eden) ,
                                                        change_eden_rel2_tolerance ,
                                                        EdenLowerLimit);
                                                if( lgEvalSlop_ls )
                                                        slope = slope_ls;
                                                if( USE_NEW && slope !=0. )
                                                {
                                                        /* slope is d(ne_true) / d(ne_assumed) */
                                                        /* >>chng 04 sep 26, add damp here too */
                                                        eden_proposed = dense.eden + (dense.EdenTrue - dense.eden)/slope_ls * damp;
                                                }
                                                else
                                                {
                                                        eden_proposed = dense.EdenTrue;
                                                }
#                                               if 0
                                                {
#include "grainvar.h"
                                                fprintf(ioQQQ,"DEBUG eden grn\t%e\t%e\t%e\t%e\n",
                                                        dense.eden, eden_proposed,dense.EdenTrue, -gv.TotalEden);
                                                }
#                                               endif

                                                /* THIS IS IT !!! this is it !!! this is where eden changes. */
                                                /* get the new electron density */
                                                /*if( dense.EdenTrue > EdenUpperLimit )*/
                                                if( eden_proposed > EdenUpperLimit )
                                                {
                                                        /* this branch, proposed eden too big */
                                                        dense.eden = EdenUpperLimit;
                                                }
                                                /*else if( dense.EdenTrue < EdenLowerLimit )*/
                                                else if( eden_proposed < EdenLowerLimit )
                                                {
                                                        /* proposed eden too small */
                                                        dense.eden = EdenLowerLimit;
                                                }
                                                else
                                                {
                                                        /* eden was within fac of the correct value, use it */
                                                        /*dense.eden = dense.EdenTrue;*/
                                                        dense.eden = eden_proposed;
                                                }

                                                if( trace.lgTrace && trace.lgNeBug )
                                                {
                                                        fprintf( ioQQQ, 
                                                                "     ConvEdenIoniz, chg ne to %.3e eden_assumed_old%10.3e ratio%10.3e EdenTrue=%10.3e\n", 
                                                        dense.eden, eden_assumed_old, dense.eden/eden_assumed_old,  dense.EdenTrue );
                                                }
                                        }
                                        /* we now have the proposed new electron density */
                                }
                                /* >>chng 00 oct 21, did not have this branch before so did not make small changes to
                                 * electron density (smaller than size that determined eden not converged */
                                /* this branch electron density was pretty close, we just need to make a small update */
                                else
                                {
                                        /* this is test for whether will call ConvIoniz again - check how close old and correct
                                        * electron densities are before updating EdenTrue */
                                        if( fabs(dense.eden-dense.EdenTrue)/dense.EdenTrue  > conv.EdenErrorAllowed/10. &&
                                                /* >>chng 03 nov 29, do not update if we are in not-ion limit of eden,
                                                 * this means that molecules and grains have most of the electrons */
                                                kase !=KASE_EDEN_NOT_ION )
                                        {

                                                /* now update eden to EdenTrue, since we are so close to is */
                                                /* >>chng 04 sep 20, from just update to taking into account slope */
                                                /*dense.eden = dense.EdenTrue;*/
                                                /* max is to avoid overshoots - we don't want a large correction this late */
                                                dense.eden = dense.eden + (dense.EdenTrue-dense.eden)/MAX2(1.,slope_ls)*damp;

                                                if( trace.nTrConvg>=3 )
                                                {
                                                        fprintf( ioQQQ, 
                                                                "   ConvEdenIoniz3 converged eden, re-calling ConvIoniz with EdenTrue=%.4e (was %.4e).\n",
                                                                dense.EdenTrue , 
                                                                dense.eden );
                                                }

                                                /* we have changed the density slightly, so now must reevaluate ionization with this new value */
                                                /* converge the current level of ionization
                                                * but only do this if change was significant */
                                                /* >>chng 02 may 31, always call ConvIoniz (basically did so anyway) and remove eden_sum from here*/
                                                /* , this calls eden_sum which updates EdenTrue */
                                                if( ConvIoniz() )
                                                {
                                                        return 1;
                                                }
                                        }
                                        else if( trace.nTrConvg>=3 )
                                        {
                                                fprintf( ioQQQ, 
                                                        "   ConvEdenIoniz3 no need to re-call ConvIoniz since eden did not change much.\n");
                                        }

                                        /* >>chng 01 mar 14, check whether electron density is no longer converged 
                                        * after reevaluating ionization */
                                        /* >>chng 04 sep 26, had div by min of eden or edentrue, use eden since
                                         * EdenTrue can be negaive when not converged */
                                        if( 
                                                fabs(dense.eden-dense.EdenTrue)/SDIV(dense.eden) >= 
                                                conv.EdenErrorAllowed )
                                        {
                                                /* difference in assumed and true electron density is larger
                                                * than tolerance, we are not yet converged, default is 0.01 */
                                                conv.lgConvEden = false;
                                                /* >>chng 04 aug 04, do not set oscillations flag if not oscillating */
                                                /* make steps smaller by setting this flag
                                                conv.lgEdenOscl = true; */
                                                lgFailLastTry = true;
                                                if( PRT_FAIL_LAST_TRY ) fprintf(ioQQQ,"lgFailLastTry set true\n");
                                                /* this will stay true through this zone */
                                                lgEden_ever_oscil = true;
                                        }
                                }
                                {
                                        /*@-redef@*/
                                        /* debug print statement for change in electron density */
                                        enum {DEBUG_LOC=false};
                                        /*@+redef@*/
                                        if( DEBUG_LOC )
                                        {
                                                fprintf(ioQQQ,"edendebugg %li\t%.2e\t%.2e\t%.2e\t%.2e\n", 
                                                        nzone,dense.eden ,eden_assumed_old, (dense.eden-eden_assumed_old)/dense.eden, dense.EdenTrue);
                                        }
                                }
                                if( trace.lgTrace && trace.lgNeBug )
                                {
                                        fprintf( ioQQQ, 
                                                "     ConvEdenIoniz, chg ne to%10.3e eden_assumed_old%10.3e ratio%10.3e EdenTrue=%10.3e converge=%c\n", 
                                        dense.eden, eden_assumed_old, dense.eden/SDIV(eden_assumed_old),  dense.EdenTrue ,TorF( conv.lgConvEden));
                                }

                                if( trace.nTrConvg>=3 )
                                {
                                        /* these prints all have : delimiter to parse into excel */
                                        fprintf( ioQQQ, "   ConvEdenIoniz3 loop:%4ld ", 
                                                LoopEden );

                                        /* this is flag says whether or not eden/eden has converged */
                                        if( conv.lgConvEden )
                                        {
                                                fprintf( ioQQQ, " converged, eden rel error:%g ", 
                                                        (dense.EdenTrue-dense.eden)/dense.EdenTrue );
                                        }
                                        else if( eden_assumed_old > SMALLFLOAT )
                                        {
                                                /*                            NB - use 8.4 fmt so chng sign not misaleign chngs*/
                                                fprintf( ioQQQ, " NOT Converged! corr:%8.4f prop:%9.5f ", 
                                                (dense.EdenTrue-eden_assumed_old)/eden_assumed_old , 
                                                        (dense.eden-eden_assumed_old)/eden_assumed_old );
                                        }

                                        fprintf( ioQQQ, " new ne:%.6e true:%.6e kase:%i slope:%.3e osc:%c", 
                                                dense.eden ,  
                                                dense.EdenTrue , 
                                                kase , 
                                                slope_ls,
                                                TorF(lgSlope_oscil));

                                        if( conv.lgEdenOscl )
                                        {
                                                fprintf( ioQQQ, " EDEN OSCIL1 damp %.4f\n", damp);
                                        }
                                        else
                                        {
                                                fprintf( ioQQQ, "\n");
                                        }
                                }

                                /* loop until converged, or we give up */
                                ++LoopEden;

                                /* this is case where overshoots did occur, but no longer */
                                /* >>chng 04 sep 23, introduce this reset */
                                /* disable resetting loop unless slope is small */
                                if( LoopEden - nLoopOscil >4 && fabs(slope_ls) < 3. )
                                {
                                        conv.lgEdenOscl = false;
                                        /* turn this back on 04 ep 26 */
                                        damp = 1.;
                                }

                                /* if last iteration through here and not converged, and no oscillations, 
                                * and no ionization failures ,
                                * then increase limit by 2x */
                                if( LoopEden == (LOOPMAX-1) && !conv.lgEdenOscl && conv.nConvIonizFails==0)
                                        /* double the limit, but only one time, and only if no oscillations */
                                        /* >>chng 04 aug 04, from 2x to 4x, to follow eden changes through eden front */
                                        LoopLimit *= 4;

                        }       while( !conv.lgConvEden && LoopEden < LoopLimit && !lgAbort );
                }
                /* turned on with set eden solver new */
                else if( strcmp( conv.chSolverEden , "new" )== 0 )
                {
                        /* will be used to look for bracketing in the electron density */
                        double PreviousEdenError = 0.;
                        double CurrentEdenError = 0.;
                        double CurrentEden = 0.;
                        double ProposedEden,
                                FractionChangeEden = 0.;

                        /* new method of converging electron density */
                        LoopEden = 0;
                        conv.nConvIonizFails = 0;

                        /* this will be increased by 2x if, at the end, no oscillations have occurred */
                        LoopLimit = LOOPMAX;

                        /* will be set true if sign of change, ever changes */
                        conv.lgEdenOscl = false;
                        if( PRT_EDEN_OSCIL ) fprintf(ioQQQ," PRT_EDEN_OSCIIL setting lgEdenOscl FALSE1\n");

                        /* this is relative change in electron density that we will permit - it will become
                        * smaller each time we bracket the true electron density */
                        change_eden_rel2_tolerance = 0.02;

                        strcpy( conv.chConvIoniz, " NONE!!!!!" );

                        /* this is ionization/electron density convergence loop
                        * keep calling ionize until lgIonDone is true */
                        do 
                        {

                                /* this flag will be set false below if electron density not within tolerance 
                                * after ionization is recomputed */
                                conv.lgConvEden = true;

                                if( trace.nTrConvg>=3 )
                                        fprintf( ioQQQ, "   ConvEdenIoniz3 calling ConvIoniz with eden= %.4e\n",dense.eden);

                                /* converge the current level of ioinization, this calls eden_sum which updates EdenTrue */
                                if( ConvIoniz() )
                                {
                                        return 1;
                                }

                                /* the history of how we got here 
                                EdenAssumed3 = EdenAssumed2,
                                EdenObtained3 = EdenObtained2;
                                EdenAssumed2 = EdenAssumed1,
                                EdenObtained2 = EdenObtained1,*/

                                /* remember the old electron density for possible printout */
                                CurrentEden = dense.eden;

                                /* this is the current error in the electron density */
                                CurrentEdenError = dense.EdenTrue - dense.eden;

                                /* conv.EdenErrorAllowed is allowed error
                                * the default value of conv.EdenErrorAllowed is 0.01 in zerologic
                                * is changed with the SET Eden Error command
                                * eden_assumed_old was incoming value of eden
                                * EdenTrue is correct value from eden_sum for current ionization
                                * new eden will be set using these, trying to prevent oscillations */

                                /* is error larger than the tolerance we are trying to beat?
                                * first check is whether error is very small after very first check, since
                                * ionization soln may not have settled down yet */
                                if( 
                                        (LoopEden==0 && fabs(CurrentEdenError)/dense.EdenTrue <= conv.EdenErrorAllowed/2.) ||
                                        (LoopEden>0 && fabs(CurrentEdenError)/dense.EdenTrue <= conv.EdenErrorAllowed 
                                        && FractionChangeEden < conv.EdenErrorAllowed / 2.) )
                                        break;

                                /* diference is assumed and true electron density is larger
                                * than tolerance, we are not yet converged, default is 0.01 */
                                conv.lgConvEden = false;
                                strcpy( conv.chConvIoniz, "Ne big chg" );

                                /* SEARCH is true if this is only search for first solution */
                                if( conv.lgSearch )
                                {
                                        dense.eden = pow(10.,
                                                (log10(dense.eden) + log10(dense.EdenTrue))/ 2.);
                                }
                                else
                                {

                                        /* was the sign of the change in the electron density changing,
                                        * if so then we have bracked the target */
                                        if( PreviousEdenError*CurrentEdenError < 0. )
                                        {
                                                /* we have bracketed the correct electron density, 
                                                 * make changes smaller, also solve linear equation for zero error */
                                                slope = (PreviousEdenError - CurrentEdenError ) /
                                                        ( eden_assumed_old - CurrentEden );

                                                ProposedEden = eden_assumed_old - PreviousEdenError / slope;
                                                /* as sanity check, this must be between the two limits we examined,
                                                * since zero was between them */
                                                ASSERT( ProposedEden >= MIN2( eden_assumed_old , CurrentEden ) );
                                                ASSERT( ProposedEden <= MAX2( eden_assumed_old , CurrentEden ) );

                                                change_eden_rel2_tolerance *= 0.25;
                                                if( trace.nTrConvg>=3 )
                                                        fprintf( ioQQQ, 
                                                        "   ConvEdenIoniz3 bracketed electron density factor now %.3e error is %.4f proposed ne %.4e\n",
                                                        change_eden_rel2_tolerance,
                                                        (dense.EdenTrue-dense.eden)/dense.EdenTrue, ProposedEden);
                                        }
                                        else
                                        {
                                                /* did not bracket the target, keep moving in its direction */

                                                /* the correct electron density - is it too different? */
                                                EdenUpperLimit = dense.eden * (1. + change_eden_rel2_tolerance);
                                                EdenLowerLimit = dense.eden / (1. + change_eden_rel2_tolerance);

                                                /* get the new electron density */
                                                if( dense.EdenTrue > EdenUpperLimit )
                                                {
                                                        /* this branch, proposed eden too big */
                                                        ProposedEden = EdenUpperLimit;
                                                }
                                                else if( dense.EdenTrue < EdenLowerLimit )
                                                {
                                                        /* proposed eden too small */
                                                        ProposedEden = EdenLowerLimit;
                                                }
                                                else
                                                {
                                                        /* eden was within fac of the correct value, use it */
                                                        ProposedEden = dense.EdenTrue;
                                                }
                                                if( trace.nTrConvg>=3 )
                                                        fprintf( ioQQQ, 
                                                        "   ConvEdenIoniz3 elec dens fac %.3e err is %.4f prop ne %.4e cor ne %.4e slope=%.2e\n",
                                                        change_eden_rel2_tolerance,
                                                        (dense.EdenTrue-dense.eden)/dense.EdenTrue ,
                                                        ProposedEden,
                                                        dense.EdenTrue,slope);
                                        }

                                        /* now remember this change for the next time through the loop  */
                                        PreviousEdenError = CurrentEdenError;
                                        eden_assumed_old = CurrentEden;
                                        FractionChangeEden = fabs( dense.eden - ProposedEden ) / dense.EdenTrue;
                                        dense.eden = ProposedEden;

                                        if( trace.lgTrace && trace.lgNeBug )
                                        {
                                                fprintf( ioQQQ, 
                                                        "     ConvEdenIoniz, chg ne to%10.3e eden_assumed_old%10.3e ratio%10.3e EdenTrue=%10.3e\n", 
                                                dense.eden, eden_assumed_old, dense.eden/eden_assumed_old,  dense.EdenTrue );
                                        }
                                }
                                /* loop until we give up  --  normal exit is break in center of loop */
                                ++LoopEden;
                        }       while( LoopEden < LoopLimit &&  !lgAbort );

                        /* we now have the proposed new electron density */
                        /* we have changed the density slightly, so now must reevaluate ionization with this new value */
                        /* converge the current level of ioinization
                        * but only do this if change was significant */
                        if( trace.nTrConvg>=3 )
                                fprintf( ioQQQ, 
                                "   ConvEdenIoniz3 converged eden, current error is %.4f, calling ConvIoniz with final density=EdenTrue=%.4e\n",
                                (dense.EdenTrue - dense.eden)/dense.EdenTrue,dense.EdenTrue);
                        dense.eden = dense.EdenTrue;
                        /* , this calls eden_sum which updates EdenTrue */
                        /* >>chng 02 may 31, always call ConvIoniz (basically did so anyway) and remove eden_sum from here*/
                        if( ConvIoniz() )
                        {
                                if( trace.nTrConvg>=3 )
                                        fprintf( ioQQQ, 
                                        "   ConvEdenIoniz3 eden no longer converged eden, current error is %.4f\n",
                                        (dense.EdenTrue - dense.eden)/dense.EdenTrue);
                        }

                        /* >>chng 01 mar 14, check whether electron density is no longer converged 
                        * after reevaluating ionization */
                        if( 
                                fabs(dense.eden-dense.EdenTrue)/dense.EdenTrue > 
                                conv.EdenErrorAllowed )
                        {
                                /* diference is assumed and true electron density is larger
                                * than tolerance, we are not yet converged, default is 0.01 */
                                conv.lgConvEden = false;
                                if( trace.nTrConvg>=3 )
                                        fprintf( ioQQQ, 
                                        "   ConvEdenIoniz3 setting eden not converged, error %.4f, exiting\n",
                                        (dense.eden-dense.EdenTrue)/dense.EdenTrue );
                        }
                        else
                        {
                                conv.lgConvEden = true;
                        }

                        if( trace.lgTrace && trace.lgNeBug )
                        {
                                fprintf( ioQQQ, 
                                        "     ConvEdenIoniz, chg ne to%10.3e eden_assumed_old%10.3e ratio%10.3e EdenTrue=%10.3e converge=%c\n", 
                                dense.eden, eden_assumed_old, dense.eden/eden_assumed_old,  dense.EdenTrue ,TorF( conv.lgConvEden));
                        }

                        if( trace.nTrConvg>=3 )
                        {
                                fprintf( ioQQQ, "   ConvEdenIoniz3 %4ld new ne=%12.4e true=%12.4e", 
                                LoopEden, dense.eden ,  dense.EdenTrue );

                                /* this is flag says whether or not eden/eden has converged */
                                if( conv.lgConvEden )
                                {
                                        fprintf( ioQQQ, " converged, eden rel error is %g ", 
                                                (dense.EdenTrue-dense.eden)/dense.EdenTrue );
                                }
                                else
                                {
                                        fprintf( ioQQQ, " NOCONV corr:%7.3f prop:%7.3f ", 
                                        (dense.EdenTrue-eden_assumed_old)/eden_assumed_old , 
                                                (dense.eden-eden_assumed_old)/eden_assumed_old );
                                }
                                if( conv.lgEdenOscl )
                                        fprintf( ioQQQ, " EDEN OSCIL2 \n");
                        }

                }
                else
                {
                        fprintf( ioQQQ, "ConvEdenIoniz finds insane solver %s \n" , conv.chSolverEden );
                        ShowMe();
                }
                ++LoopBig;
        }

        if( trace.nTrConvg>=3 )
        {
                fprintf( ioQQQ, "   ConvEdenIoniz3 exits, lgConvEden=%1c entry eden %.4e -> %.4e rel change %.3f\n" , 
                        TorF(conv.lgConvEden ),
                        EdenEntry ,
                        dense.eden ,
                        (EdenEntry-dense.eden)/SDIV( dense.eden ) );
        }
        else if( trace.lgTrace )
        {
                fprintf( ioQQQ, " ConvEdenIoniz exits zn %.2f lgConvEden=%1c entry eden %.4e -> %.4e rel change %.3f\n" , 
                        fnzone,
                        TorF(conv.lgConvEden ),
                        EdenEntry ,
                        dense.eden ,
                        (EdenEntry-dense.eden)/SDIV( dense.eden ) );
        }

        return 0;
}

/* returns true if electron density is converged */
bool lgConvEden(void)
{
        bool lgRet;

        /* >>chng 04 sep 26, div by eden, not min of eden and edentrue since latter now poss neg */
        if( fabs(dense.eden-dense.EdenTrue)/SDIV(dense.eden) >= 
                conv.EdenErrorAllowed )
        {
                lgRet = false;
                conv.lgConvEden = false;
        }
        else
        {
                lgRet = true;
                conv.lgConvEden = true;
        }
        return lgRet;
}

---