/*he1gma evaluate photoionization rates for helium singlets */
#include "cddefines.h"
#include "nhe1lvl.h"
#include "opacpoint.h"
#include "he1crt.h"
#include "he1ind.h"
#include "rfield.h"
#include "phe1lv.h"
#include "heat.h"
#include "ip2gam.h"
#include "nh.h"
#include "secondaries.h"
#include "trace.h"
#include "heopfr.h"
#include "nhe1.h"
#include "hhe1.h"
#include "he1nionryd.h"
#include "gammas.h"
#include "photrate.h"
#include "he1.h"

void he1gma( void )
{
	long int i, 
	  n;
	float g2elec, 
	  rectot;
	static double ghe1[9], 
	  he1cn[NHE1LVL], 
	  p1, 
	  p2;

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

	/* photoionization rates
	 * He I ground state ionization
	 *
	 * photoionization He */
	/* with induced recombination, PhotScaleOn is 1, set to
	 * zero with 'no photoionization' command */
	ghe1[0] = GammaBn(nhe1Com.nhe1[0],ip2gam.ip2ePhoto-1,OpacPoint.iophe1[0],
	  1.808,&he1ind.he1qin[0],&he1cn[0])*PhotRate.PhotScaleOn;

	hhe1Com.hhe1[0] = heat.HeatNet;
	PhotRate.PhotoRate[0][0][0][1] = ghe1[0];
	PhotRate.PhotoRate[1][0][0][1] = heat.HeatLowEnr;
	PhotRate.PhotoRate[2][0][0][1] = heat.HeatHiEnr;

	/* now add on rest of integral; these energies can double photoionize */
	g2elec = (float)(GammaBn(ip2gam.ip2ePhoto,rfield.nflux,OpacPoint.iophe1[0]+
	  ip2gam.ip2ePhoto-nhe1Com.nhe1[0],1.808,&p1,&p2)*PhotRate.PhotScaleOn);

	hhe1Com.hhe1[0] += heat.HeatNet;
	PhotRate.PhotoRate[0][0][0][1] += ghe1[0];
	PhotRate.PhotoRate[1][0][0][1] += heat.HeatLowEnr;
	PhotRate.PhotoRate[2][0][0][1] += heat.HeatHiEnr;

	if( trace.lgTrace && trace.lgHeBug )
	{
		GammaPrt(nhe1Com.nhe1[0],rfield.nflux,OpacPoint.iophe1[0],
		  ioQQQ,ghe1[0]+g2elec,(ghe1[0]+g2elec)*0.05);
	}

	he1ind.he1qin[0] += p1;
	he1cn[0] += p2;
	ghe1[0] += g2elec;

	/* gamma => 2 elec cross section - Brown (1970) */
	g2elec *= 0.03f;

	/* HeII LA and diffuse continua */
	phe1lv.he1gam[0] = ghe1[0];

	/* He I n>=2 state ionization */
	/* F2 is related to induced rec rate, HE1CN(N) is related to cooling */
	for( n=1; n < NHE1LVL; n++ )
	{
		ghe1[n] = GammaBn(nhe1Com.nhe1[n],nhe1Com.nhe1[0],OpacPoint.iophe1[n],
		  He1NIonRyd.He1IonRyd[n],&he1ind.he1qin[n],&he1cn[n])*
		  PhotRate.PhotScaleOn;
		hhe1Com.hhe1[n] = heat.HeatNet;
		phe1lv.he1gam[n] = ghe1[n];
	}

	/*  induced recombination rate */
	he1ind.he1rin[0] = he1ind.he1qin[0]*he1crt.he1lte[0];
	he1ind.he1rin[1] = he1ind.he1qin[1]*(he1crt.he1lte[1] + he1crt.he1lte[9]);

	for( n=2; n < NHE1LVL; n++ )
	{
		/* induced recombination rate */
		he1ind.he1rin[n] = he1crt.he1lte[n]*he1ind.he1qin[n];
	}

	if( trace.lgTrace && trace.lgHeBug )
	{
		rectot = phe1lv.he1rec[1][0]*phe1lv.he1rec[0][0];
		for( i=1; i < NHE1LVL; i++ )
		{
			rectot += phe1lv.he1rec[1][i]*phe1lv.he1rec[0][i];
		}

		fprintf( ioQQQ, 
			"     He sin Gam:%10.2e RadRec:%10.2e G eff%10.2e GndDrate%10.2e FrOp(other)%10.2e csup%10.2e\n", 
		  phe1lv.he1gam[0], rectot, phe1lv.he1rec[1][0], ghe1[0], 
		  heopfr.ophe1f[0], Secondaries.csupra );
	}

	if( trace.lgTrace && (trace.lgHe1Bug || trace.lgHeBug) )
	{
		if( phe1lv.he1gam[0] > 0. )
		{
			fprintf( ioQQQ, 
				"     HE1GAMgr ate%10.2e pht frac Prim%10.2e ots228%10.2e\n", 
			  phe1lv.he1gam[0], ghe1[0]/phe1lv.he1gam[0], 
			  rfield.otscon[nh.ipHn[1][IP1S]-1]*1.75e-18/phe1lv.he1gam[0] );
		}
		else
		{
			fprintf( ioQQQ, "     HE1GAMgr pht frac = none.\n" );
		}
	}

	if( trace.lgTrace && trace.lgHe1Bug )
	{
		fprintf( ioQQQ, "     HE1GMA finds photo rates:" );
		for( i=1; i <= 9; i++ )
		{
			fprintf( ioQQQ, "%2ld%10.2e", i, phe1lv.he1gam[i-1] );
		}
		fprintf( ioQQQ, "\n" );
	}


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