/*esum sum free electron density over all species, sets variable erredn.EdenTrue
 *called by ConvEdenIoniz which actually controls the electron density updates */
#include "cddefines.h"
#include "showme.h"
#include "ffsum.h"
#include "trace.h"
#include "ionfracs.h"
#include "phycon.h"
#include "hevmolec.h"
#include "hmi.h"
#include "elecden.h"

void esum(void)
{
	long int n, 
	  ion, 
	  nelem;

	double csum, 
	  edmole, 
	  edsave[LIMELM], 
	  osum;

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

	/* EdenExtra is normally zero, set with EDEN command, to add extra e-
	 * EdenSet is normally zero, set with SET EDEN command, to set e- */
	ElecDen.EdenTrue = (float)(xIonFracs[2][0] + xIonFracs[2][1] + 
		2.*xIonFracs[3][1] + ElecDen.EdenExtra);

	/* add on molecules */
	edmole = ElecDen.EdenTrue;

	ElecDen.EdenTrue += hmi.h2plus + hevmolec.hevmol[IPCHP-1] + hevmolec.hevmol[IPOHP-1] + 
	  hevmolec.hevmol[IPCOP-1] + hevmolec.hevmol[IPH2OP-1] + hevmolec.hevmol[IPO2P-1] + 
	  hevmolec.hevmol[IPC2P-1] + hevmolec.hevmol[IPH3OP-1] + hevmolec.hevmol[IPCH2P-1] + 
	  hmi.h3plus;

	/* remember electron density from molecules */
	edmole = ElecDen.EdenTrue - edmole;

	/* >>chng 97 mar 18, some extreme models have H- density very large during
	 * 666 Error! cap below should be elsewhere, this sub should produce TrueEden
	 *search phase when far from solution.  do not let eden go neg */
#if 0
	if( search.lgSearch )
	{
		ElecDen.EdenTrue -= (float)MIN2(ElecDen.EdenTrue*0.05,hmi.hminus);
	}
	else
	{
		ElecDen.EdenTrue -= (float)MIN2(ElecDen.EdenTrue/100.,hmi.hminus);
	}
#	endif
	/* correct for hminus, and flag error if we go negative */
	ElecDen.EdenTrue -= hmi.hminus;

	/* sum over all metals */
	ElecDen.EdenMet = 0.;
	/* free free sum should not include hydrogenic since already done
	 * indpenedently in coolr*/
	ffsum.EdenFFSum = 0.;
	for( nelem=2; nelem < LIMELM; nelem++ )
	{
		edsave[nelem] = 0.;
		for( ion=2; ion <= nelem; ion++ )
		{
			/* >>chng 96 oct 27, save all contributors to electron density */
			edsave[nelem] += (float)(ion-1)*xIonFracs[ion][nelem];
			ffsum.EdenFFSum += POW2( (float)(ion-1) )*xIonFracs[ion][nelem];
		}
		/* add on hydrogenic specides since now counted in free free sum */
		edsave[nelem] += (float)(nelem)*xIonFracs[nelem+1][nelem];
		ElecDen.EdenMet += (float)edsave[nelem];
	}

	ElecDen.EdenTrue += ElecDen.EdenMet;

	/* this variable is set with the set eden command, 
	 * is supposed to override physical electron density */
	if( ElecDen.EdenSet > 0. )
	{
		ElecDen.EdenTrue = ElecDen.EdenSet;
	}

	/* get all electrons not due to hydrogen */
	ElecDen.EdenNotH = ElecDen.EdenTrue - xIonFracs[2][0];

 	if( trace.lgTrace )
	{
		csum = 0.;
		for( ion=1; ion <= 6; ion++ )
		{
			csum += (float)(ion)*xIonFracs[ion+1][5];
		}

		osum = 0.;
		for( ion=1; ion <= 8; ion++ )
		{
			osum += (float)(ion)*xIonFracs[ion+1][7];
		}

		fprintf( ioQQQ, "     ESUM: TrueEden%11.4e old%11.4e ff sum%10.3e H%8.2e He%8.2e C%8.2e O%8.2e met%8.2e mol%8.2e\n", 
		  ElecDen.EdenTrue, phycon.eden, ffsum.EdenFFSum, xIonFracs[2][0], 
		  xIonFracs[2][1] + 2.*xIonFracs[3][1], csum, osum, ElecDen.EdenMet, edmole );

		if( trace.lgNeBug )
		{
			fprintf( ioQQQ, "    " );
			for(n=0; n < LIMELM; n++)
				fprintf( ioQQQ, "%8.1e", edsave[n] );
			fprintf( ioQQQ, "\n" );
		}
	}

	if( ElecDen.EdenTrue <= 0. )
	{
		fprintf( ioQQQ, "ESUM finds non-positive electron density.\n" );
		fprintf( ioQQQ, 
			"     ESUM: EdenTrue to%12.4e fm%12.4e ff sum;%10.3e Ne(H):%10.2e Ne(He):%10.2e Ne(C):\n", 
		  ElecDen.EdenTrue, 
		  phycon.eden, 
		  ffsum.EdenFFSum, 
		  xIonFracs[2][0], 
		  xIonFracs[2][1] + 2.*xIonFracs[3][1] );
		ShowMe();
		puts( "[Stop in esum]" );
		exit(1);
	}

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