/*he3gma compute photoionization rates for helium triplets */
#include "cddefines.h"
#include "nhe3lvl.h"
#include "nhe1lvl.h"
#include "opacpoint.h"
#include "heat.h"
#include "he3lte.h"
#include "rfield.h"
#include "he.h"
#include "he3hsv.h"
#include "tepowers.h"
#include "photrate.h"
#include "phycon.h"
#include "he3rr.h"
#include "he3gams.h"
#include "phe1lv.h"
#include "trace.h"
#include "nhe1.h"
#include "herec.h"
#include "csphot.h"
#include "taulines.h"
#include "gammas.h"
#include "he3.h"

void he3gma(void)
{
	double hirate, 
	  r;
	double fac;
	static double f1, 
	  ghe32s, 
	  he3h;

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

	if( trace.lgTrace && trace.lgHe3Bug )
	{
		fprintf( ioQQQ, " HE3GMA called.\n" );
	}

	/* following are corrected radiative recombination coefficient
	 * not actual ones
	 * >>chng 96 oct 25, removed eden from edensqte, so not cancel of eden below
	 * some evidence that eden and edensqte not updated together */
	he3rr.r2s = (float)(1.31e-12/TePowers.sqrte*TePowers.te10);

	he3rr.r2p = (float)(3.09e-11/(TePowers.te10*TePowers.te03)/TePowers.sqrte);

	he3rr.r3s = (float)(7.1e-13/TePowers.sqrte);

	he3rr.r3p = (float)(2.24e-11/(TePowers.te10*TePowers.te10*TePowers.te03)/
	  TePowers.sqrte);

	he3rr.r3d = (float)(1.25e-9/phycon.te/TePowers.te10);

	herecCom.rectri = he3rr.r2s + he3rr.r2p + he3rr.r3s + he3rr.r3p + 
	  he3rr.r3d;

	/* with induced recombination */
	ghe32s = GammaBn(he.nhei3,nhe1Com.nhe1[0],OpacPoint.ioptri,
	  0.3506,&f1,&r)*PhotRate.PhotScaleOn;

	he3h = heat.HeatNet;
	if( trace.lgTrace && trace.lgHeBug )
	{
		GammaPrt(he.nhei3,nhe1Com.nhe1[0],OpacPoint.ioptri,ioQQQ,
		  ghe32s,ghe32s*0.05);
	}

	/* induced recombination rates */
	he3rr.r2s += (float)(f1*he3lteCom.he3lte[0]/phycon.eden);
	he3rr.r2p += (float)(f1*he3lteCom.he3lte[1]/phycon.eden);
	he3hsvCom.he3hsv[0] = he3h;

	/* put in rest of flux integral */
	if( nhe1Com.nhe1[0] < rfield.nflux )
	{
		/* ratio of trip cross to singlet ground */
		fac = csphot(nhe1Com.nhe1[0],nhe1Com.nhe1[0],  OpacPoint.iophe1[0]);
		hirate = (float)(phe1lv.he1gam[0]*csphot(nhe1Com.nhe1[0],he.nhei3,
		  OpacPoint.ioptri)/MAX2(1e-37,fac));
	}
	else
	{
		hirate = 0.;
	}

	/* now save total rates for 2^3S and 2^3P */
	he3gams.he3gam[0] = ghe32s + hirate;
	he3gams.he3gam[1] = (float)(he3gams.he3gam[0]*3.);

	he3gams.He3FrLya = (float)(
		rfield.otslin[HydroLines[0][IP2P][IP1S].ipCont-1]*
		csphot(HydroLines[0][IP2P][IP1S].ipCont,
	  he.nhei3,OpacPoint.ioptri));

	if( he3gams.he3gam[0] > 1e-30 )
	{
		he3gams.He3FrLya /= (float)he3gams.he3gam[0];
	}
	else
	{
		he3gams.He3FrLya = 0.;
	}

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