/*helium solve ionization balance helium ion, helium singlets, helium triplets */
#include <string.h>
#include "cddefines.h"
#include "taulines.h"
#include "taumin.h"
#include "nhe3lvl.h"
#include "nhe1lvl.h"
#include "phe1lv.h"
#include "herec.h"
#include "trace.h"
#include "ionfracs.h"
#include "heating.h"
#include "he.h"
#include "he3tau.h"
#include "he1bn.h"
#include "he3n.h"
#include "he3bn.h"
#include "hhe1.h"
#include "hydrogenic.h"
#include "he3hsv.h"
#include "elmton.h"
#include "ioreco.h"
#include "hionfraccom.h"
#include "hhe2phtots.h"
#include "itercnt.h"
#include "ionrange.h"
#include "converge.h"
#include "heopfr.h"
#include "he1nxt.h"
#include "he1tau.h"
#include "he1.h"
#include "he3.h"
#include "helium.h"

void helium( void )
{
	long int i;
	long int n;
	float SimpHePlusOvSing;
	double dtri, 
	  he3old, 
	  he3new ,
	  SimpHePlusPlusOvPlus, 
	  ssav[7], 
	  sum,
	  HePlusOvSing;

	static double 
	  HePlusOvTrip, 
	  HePlusPlusOvPlus, 
	  r1ov2;

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

	/* option to "turn off" helium */
	if( !elmton.lgElmtOn[1] )
	{
		he.hei3 = 0.;
		he.heii = 0.;
		he.heiii = 0.;
		he.hei1 = 0.;
		he.hei = 0.;
		xIonFracs[1][1] = 0.;
		xIonFracs[2][1] = 0.;
		xIonFracs[3][1] = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->helium()\n", debug_fp );
#		endif
		return;
	}

	/* if this is the very first call to the helium routines in this model, but not the
	 * first in this coreload, then photo rates from previous soln are still in he1gam,
	 * which will be used below to add onto the triplet rates.  We need to zero them
	 * in this case, so that second models act like inital models in the coreload */

	if( !conv.nTotalIoniz && IterCnt.iter==0 )
	{
		/*these will be set in RTDiffuse, are two photon diffuse fields
		 * deep in cloud */
		hhe2phtots.H12phtOTS = 0.;
		hhe2phtots.He12phtOTS = 0.;
		hhe2phtots.He22phtOTS = 0.;

		for( i=0; i < (NHE1LVL - 1); i++ )
		{
			phe1lv.he1gam[i] = 0.;
			heopfr.ophe1f[i] = 0.;
		}

		phe1lv.he12s = 1e-30f;
		phe1lv.he12p = 1e-30f;
		for( i=0; i < NHE3TAU; i++ )
		{
			he3tau[i].TauIn = tauminCom.taumin;
			he3tau[i].TauTot = 1e20f;
			he3tau[i].Pesc = 1.;
			he3tau[i].FracInwd = 1.;
		}
		/* zero out inner optical depths */
		for( i=0; i < (NHE1LVL - 1); i++ )
		{
			for( n=0; n < NHE1LVL; n++ )
			{
				he1nxtCOM.he1nxt[i][n] = tauminCom.taumin;
				he1tauCOM.he1tau[i][n] = tauminCom.taumin;
			}
		}
		t206.TauIn = tauminCom.taumin;
		t206.TauTot = tauminCom.taumin;
		t206.Pesc = 1.;
		t206.FracInwd = 0.5;
		t206.Aul = 1.976e6;
	}

	/* remember old triplet ionization so we can test for changes */
	if( xIonFracs[2][ipHE]/xIonFracs[0][ipHE] > 0.05 )
	{
		/* helium escape and destruction probs were done in RTMake */
		he3old = he.hei3 / xIonFracs[2][ipHE];
	}
	else
	{
		he3old = 0.;
	}

	/* we do not generate a He++/He+ ratio here 
	 * since this was done in hydrogenic routines.  They defined the following
	 * which will be imposed on the solution here -  these are ratios (total bound)/ion */
	HePlusPlusOvPlus = HIonFraccom.HIonFrac[ipHE];

	/* this is the simple estimate from ground state photoionization 
	 * and radiative recombination*/
	SimpHePlusPlusOvPlus = HIonFraccom.HIonSimple[ipHE];

	/* do trip before singlets since need rec coef to singlets from trip
	 * R2OVT is HeII/He(trip) */
	if( IonRange.IonHigh[1] >= 2 )
	{
		he3col();
		he3rad();
		he3gma();

		he3lev(&HePlusOvTrip);

		/* do singlets; HE2OV1 is HeII/He(sin) */
		he1col();
		he1rad();
		he1gma();
		he1lev(&HePlusOvSing,&SimpHePlusOvSing);
	}
	else 
	{
		HePlusOvSing = 0.;
		SimpHePlusOvSing = 0.;
		HePlusOvTrip = 0.;
	}

	if( trace.lgTrace && trace.lgHeBug )
	{
		fprintf( ioQQQ, 
			"     helium I/2%10.3e simple=%10.3e 2/1%10.3e simple=%10.3e 2/3%10.3e\n", 
		  HePlusPlusOvPlus, SimpHePlusPlusOvPlus, 
		  HePlusOvSing, SimpHePlusOvSing, HePlusOvTrip );
	}

	/* populations */
	if( HePlusOvSing == 0.0 || HePlusOvTrip == 0.0 )
	{
		/* gas is completely neutral (although trace He++ possible) */
		he.hei3 = 0.;
		he.heii = 0.;
		he.heiii = 0.;
		he.hei1 = xIonFracs[0][ipHE];
		he.hei = xIonFracs[0][ipHE];
	}
	else if( HePlusPlusOvPlus == 0. && HePlusOvSing != 0. )
	{
		/* gas is singly ionized, but not doubly */
		he.heiii = 0.;
		/* ratio of total neutral to ionized helium */
		r1ov2 = 1.e0/HePlusOvTrip + 1.e0/HePlusOvSing;
		he.heii = (float)(xIonFracs[0][ipHE]/(1.e0 + r1ov2));
		he.hei = (float)(xIonFracs[0][ipHE]/(1.e0 + 1.e0/r1ov2));
		he.hei1 = (float)(he.hei/(1.e0 + HePlusOvSing/HePlusOvTrip));
		he.hei3 = (float)(he.hei/(1.e0 + HePlusOvTrip/HePlusOvSing));
	}
	else
	{
		/* all three stages of ionization populated
		 * ratio of total neutral to ionized helium */
		r1ov2 = 1.e0/HePlusOvTrip + 1.e0/HePlusOvSing;
		he.heiii = (float)(xIonFracs[0][ipHE]/(1.e0 + (1.e0 + r1ov2)/HePlusPlusOvPlus));
		he.heii = (float)(xIonFracs[0][ipHE]/(1.e0 + r1ov2 + HePlusPlusOvPlus));
		he.hei = (float)(xIonFracs[0][ipHE]/(1.e0 + (1.e0 + HePlusPlusOvPlus)/r1ov2));
		he.hei1 = (float)(he.hei/(1.e0 + HePlusOvSing/HePlusOvTrip));
		he.hei3 = (float)(he.hei/(1.e0 + HePlusOvTrip/HePlusOvSing));
	}

	/* chng>>> 99 apr 16, to damp out HeI 3 - Lya oscillations in pphheonly.in */
	{
		float FacNew=0.5f;
		he.hei3 = FacNew*he.hei3 + (float)( (1.-FacNew)*he3old*xIonFracs[2][ipHE] );
	}

	/* save into master ionization array */
	xIonFracs[1][1] = he.hei;
	xIonFracs[2][1] = he.heii;
	xIonFracs[3][1] = he.heiii;

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, "     HE; 23S:%10.2e HEo:%10.2e HE+:%10.2e HE++:%10.2e\n", 
		  he.hei3, he.hei, he.heii, he.heiii );
	}

	/* July 93 logic added for partially neutral soln, not yet converged */
	/* >>chng 99 feb 27, declare non-convergence if changes relative to ion,
	 * only if ion is present */
	if( xIonFracs[2][1]/xIonFracs[0][1] > 0.05 )
	{
		he3new = he.hei3 / xIonFracs[2][1];
		dtri = fabs(he3new-he3old)/he3new;
		if( he.hei3 > 0. && dtri > 0.02 )
		{
			conv.lgConvIoniz = FALSE;
			strcpy( conv.chConvIoniz, "He trip ch" );
			conv.BadConvIoniz[0] = he3old;
			conv.BadConvIoniz[1] = he3new;
		}
	}
	else
	{
		dtri = 0.;
	}

	/* get total heating for helium singlets */
	ssav[0] = hhe1Com.hhe1[0]*phe1lv.he1n[0];
	sum = ssav[0];
	for( n=1; n < 7; n++ )
	{
		ssav[n] = hhe1Com.hhe1[n]*phe1lv.he1n[n];
		sum += ssav[n];
	}
	HeatingCom.heating[0][1] = sum*xIonFracs[2][1];

	/* heating(2,3) is helium triplets */
	HeatingCom.heating[2][1] = he3hsvCom.he3hsv[0]*he3nCom.he3n[0]*
	  he.hei3;

	/* option to punch recombination coefficients, set with punch recom coef command */
	if( lgioRecom )
	{
		fprintf( ioRecom, "Helium, singlets+triplets= %g\n", 
			herecCom.reci + herecCom.rectri );
	}

	if( trace.lgTrace && trace.lgHeBug )
	{
		fprintf( ioQQQ, "     HE; bn(2tripS)=%10.2e bn(HeI)=%10.2e bn(HeII)=%10.2e and helium returns.\n", 
		  he3bnCom.he3bn[0], he1bnCOM.he1bn[0], hydro.hbn[1][IP1S] );
	}


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