/home66/gary/public_html/cloudy/c08_branch/source/eden_sum.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 */
/*eden_sum sum free electron density over all species, sets variable erredn.EdenTrue
 *called by ConvBase - ConvEdenIoniz actually updates the electron density 
 * returns 0 if all is ok, 1 if need to abort calc */
#include "cddefines.h"
#include "hmi.h"
#include "trace.h"
#include "grainvar.h"
#include "rfield.h"
#include "mole.h"
#include "dense.h"
#include "conv.h"

/*eden_sum sum free electron density over all species, sets variable erredn.EdenTrue
 *called by ConvBase - ConvEdenIoniz actually updates the electron density 
 * returns 0 if all is ok, 1 if need to abort calc */
int eden_sum(void)
{
        long int i,
          ion, 
          nelem;

        realnum sum_all_ions ,
                sum_metals ,
                hmole_eden,
                eden_ions[LIMELM];

        DEBUG_ENTRY( "eden_sum()" );

        /* EdenExtra is normally zero, set with EDEN command, to add extra e- */
        dense.EdenTrue = dense.EdenExtra;

        /* sum over all ions */
        sum_all_ions = 0.;
        sum_metals = 0.;
        for( nelem=ipHYDROGEN; nelem < LIMELM; nelem++ )
        {
                if( nelem==ipLITHIUM )
                        sum_metals = 0.;
                eden_ions[nelem] = dense.xIonDense[nelem][1];
                /* >>chng 02 jul 15, upper limit now includes H-like, so all species
                 * in this one loop (except hydrogen  and helium ) */
                /* >>chng 02 apr 26, upper limit had been nelem, not nelem+1, so H-like
                 * metals were missed - ok for H and He since not part of this loop before */
                for( ion=2; ion <= nelem+1; ion++ )
                {
                        /* >>chng 96 oct 27, save all contributors to electron density */
                        eden_ions[nelem] += ion*dense.xIonDense[nelem][ion];
                }

                sum_all_ions += eden_ions[nelem];
                sum_metals += eden_ions[nelem];
        }
        dense.EdenTrue += sum_all_ions;
        ASSERT( dense.EdenTrue >= 0. );

        /* add on molecules */
        co.comole_eden = 0.;
        for( i=0; i < mole.num_comole_calc; i++ )
        {
                if(COmole[i]->n_nuclei != 1)
                {
                        co.comole_eden += COmole[i]->hevmol*COmole[i]->nElec;
                }
        }

        /* electrons contributed by molecules */
        dense.EdenTrue += co.comole_eden;
        ASSERT( dense.EdenTrue >= 0. );

        /* >>chng 03 nov 28, loop over all H molecules, had just been H- */
        hmole_eden = 0.;
        for(i=0;i<N_H_MOLEC;i++)
        {
                /* hmi.nElectron is zero for H+ since counted as an ion, -1 for H-, etc */
                hmole_eden += hmi.Hmolec[i]*hmi.nElectron[i];
        }
        /* >>chng 01 jan 08, following logic added to stop H- from forcing
         * negative electron density during very first search, when H- is badly off */
        /* >>chng 03 nov 28, add test for lgSearch, had been absent */
        /* >>chng 04 feb 20, recoil_ion.in shows this important,
         * update test so prevent neg eden */
        if( (-hmole_eden) > dense.EdenTrue/4. && conv.lgSearch )
        {
                /*dense.EdenTrue += MIN2(dense.EdenTrue*0.05,hmole_eden);*/
                dense.EdenTrue *= 0.9;
        }
        else if( (-hmole_eden) > dense.EdenTrue )
        {
                /* >>chng 04 mar 14, add this branch to prevent
                 * negative electron density.  occurred in pdr
                 * with large h2, after hmole failure */
                fprintf(ioQQQ," PROBLEM sum eden from hmole too neg, set to  limt.  EdenTrue:%.3e hmole_eden:%.3e \n",
                        dense.EdenTrue, 
                        hmole_eden);
                dense.EdenTrue = dense.EdenTrue/2.;
        }
        else
        {
                /* account for electrons on all H-bearing molecules */
                dense.EdenTrue += hmole_eden;
        }
        ASSERT(dense.EdenTrue >= 0. );

        /* this variable is set with the set eden command, 
         * is supposed to override physical electron density */
        if( dense.EdenSet > 0. )
        {
                dense.EdenTrue = dense.EdenSet;
                dense.eden_from_metals = 1.;
        }
        else
        {
                /* fraction of electrons from si, s, mg, al */
                /*dense.eden_from_metals = sum_all_ions / SDIV(dense.EdenTrue);*/
                dense.eden_from_metals = sum_metals / SDIV(dense.EdenTrue);
                /*fprintf(ioQQQ," debuggg ionfrac %.2f %.3f\n",fnzone,dense.eden_from_metals );*/
        }
        /* dense.eden itself is actually changed in ConvEdenIoniz */

        /* >>chng 00 dec 19, include electrons on grains in total sum */
        /* >>chng 02 jul 15, only add grain electrons when we have stable soln -
         * everett.in had very negative grain elec total, so drove total eden negative,
         * but due to bad inital electron density - do not add grain electrons until we
         * are close to a solution */
        /* >>chng 04 sep 26, allow nEdenTrue to become -ve  */
        if( dense.EdenTrue+gv.TotalEden*gv.lgGrainElectrons >= 0. )
        {
                /* gv.lgGrainElectrons - should grain electron source/sink be included in overall electron sum?
                 * default is true, set false with no grain electrons command */
                dense.EdenTrue += gv.TotalEden*gv.lgGrainElectrons;
        }
        else
        {
                /* >>chng 05 jan 24, only print if not in search phase */
                if( !conv.lgSearch )
                        fprintf(ioQQQ,
                        " PROBLEM eden grain neg limt %.2f eden %.4e new eden bef grn %.4e"
                        "grain eden %.4e  set edentrue to %.4e Search?%c\n",
                        fnzone,
                        dense.eden, 
                        dense.EdenTrue, 
                        gv.TotalEden,
                        dense.eden*0.9,
                        TorF(conv.lgSearch));

                /*dense.EdenTrue = dense.eden*0.9;*/
                dense.EdenTrue += gv.TotalEden*gv.lgGrainElectrons;
        }

        if( trace.lgTrace || trace.lgESOURCE )
        {
                fprintf( ioQQQ, 
                        "     eden_sum zn:%.2f current:%.4e new true:%.4e ions:%.4e comole:%.4e hmole:%.4e grain:%.4e extra:%.4e sum:%.4e LaOTS:%.4e\n",
                  fnzone ,
                  dense.eden , 
                  dense.EdenTrue , 
                  sum_all_ions ,
                  co.comole_eden ,
                  hmole_eden ,
                  gv.TotalEden*gv.lgGrainElectrons,
                  dense.EdenExtra ,
                  sum_all_ions + co.comole_eden + hmole_eden + gv.TotalEden*gv.lgGrainElectrons + dense.EdenExtra,
                  rfield.otslin[2182] );

                if( trace.lgNeBug )
                {
                        for(nelem=ipHYDROGEN; nelem < LIMELM; nelem++)
                        {
                                if( nelem==0 )
                                {
                                        fprintf( ioQQQ, "      eden_sum H -Ne:" );
                                }
                                else if( nelem==10 )
                                {
                                        fprintf( ioQQQ, "\n      eden_sum Na-Ca:" );
                                }
                                else if( nelem==20 )
                                {
                                        fprintf( ioQQQ, "\n      eden_sum Sc-Zn:" );
                                }
                                fprintf( ioQQQ, " %.4e", eden_ions[nelem] );
                                if( nelem==29 )
                                {
                                        fprintf( ioQQQ, "\n" );
                                }
                                /*if( nelem==9 || nelem==19 || nelem==29 )
                                        fprintf( ioQQQ, "\n      " );*/
                        }
                }
        }

        /* abort if negative electron density */
        /*>>chng 04 sep 25, allow negative true elec density so solver can work  */
        if( 0 && dense.EdenTrue < 0. )
        {
                fprintf( ioQQQ, "eden_sum finds non-positive electron density.\n" );
                fprintf( ioQQQ, 
                        "     eden_sum: EdenTrue to%12.4e fm%12.4e Ne(H):%10.2e Ne(He):%10.2e Ne(C):\n", 
                  dense.EdenTrue, 
                  dense.eden, 
                  dense.xIonDense[ipHYDROGEN][1], 
                  dense.xIonDense[ipHELIUM][1] + 2.*dense.xIonDense[ipHELIUM][2] );
                ShowMe();
                cdEXIT(EXIT_FAILURE);
        }
        else if( dense.EdenTrue == 0. )
        {
                fprintf( ioQQQ, "\nDISASTER PROBLEM eden_sum finds an electron density of zero, this is unphysical.\n" );
                fprintf( ioQQQ, "There appears to be no source of ionization.\n");
                fprintf( ioQQQ, "Consider adding background cosmic rays with COSMIC RAY BACKGROUND and BACKGROUND commands.\n\n");
                lgAbort = true;

                return 1;
        }

        {
                /*@-redef@*/
                enum {DEBUG_LOC=false};
                /*@+redef@*/
                if( DEBUG_LOC )
                {
                        fprintf(ioQQQ,"esumdebugg\t%li\t%.2e\t%.2e\n",
                                nzone,
                          dense.eden, dense.EdenTrue);
                }
        }

        /* >>chng 05 jan 05, don't let elec den be zero - logs are taken */
        dense.eden = MAX2( SMALLFLOAT , dense.eden );

        return 0;
}

---