/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * 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;

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

#	ifdef DEBUG_FUN
	fputs( "<+>eden_sum()\n", debug_fp );
#	endif

	/* 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 < NUM_COMOLE_CALC; i++ )
	{
		co.comole_eden += co.hevmol[i]*co.nElec[i];
	}

	/* 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;

		/*conv.lgConvEden = FALSE;
		strcpy( conv.chConvEden, "EdnNegGrn" );*/
	}
	/* >>chng 04 sep 26, chng to allow neg true eden */
	/*ASSERT( dense.EdenTrue != 0. );*/
	/*ASSERT( dense.EdenTrue != 0. );*/

 	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( 1 || trace.lgNeBug )
		{
			/*lint -e771 conceivably not initialized */
			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      " );*/
			}
			/*lint +e771 conceivably not initialized */
		}
	}

	/* 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();
		puts( "[Stop in eden_sum]" );
		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;

#		ifdef DEBUG_FUN
		fputs( " <->eden_sum()\n", debug_fp );
#		endif

		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 );

#	ifdef DEBUG_FUN
	fputs( " <->eden_sum()\n", debug_fp );
#	endif

	return 0 ;
}

