/*HydroPhoto photoionization, recombination, radiative rates for model hydrogen atom */
#include "cddefines.h"
#include "opacpoint.h"
#include "hydrogenic.h"
#include "rfield.h"
#include "opac.h"
#include "nh.h"
#include "trace.h"
#include "hphoto.h"
#include "photrate.h"
#include "hpheat.h"
#include "hgamsc.h"
#include "heat.h"
#include "hcolrt.h"
#include "hlteok.h"
#include "printe82.h"
#include "gammas.h"

void HydroPhoto(long int ipZ)
{
	long int n, 
	  limit;

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

	/* check that we were called with valid charge */
	assert( ipZ >= 0);
	assert( ipZ < LIMELM );

	/* do photoionization rates */
	/* induced recombination; FINDUC is integral of
	 * pho rate times EXP(-hn/kt) for induc rec
	 * CIND is this times hnu-hnu0 to get ind rec cooling
	 * PhotScaleOn is 1, set to 0 with 'no photoionization' command
	 * ipSecIon points to 7.353 Ryd, lowest energy where secondary ioniz
	 * of hydrogen is possible */

	HPhoto.hgamnc[ipZ][IP1S] = GammaBn(nh.ipHn[ipZ][IP1S],
		rfield.nflux,
		OpacPoint.ipHOpac[ipZ][IP1S],
	  hydro.HLevNIonRyd[ipZ][IP1S],
	  &HPhoto.finduc[ipZ][IP1S],
	  &HPhoto.cind[ipZ][IP1S])*
	  PhotRate.PhotScaleOn;

	HPHeat.HPhotHeat[ipZ][IP1S] = heat.HeatNet*PhotRate.PhotScaleOn;

	PhotRate.PhotoRate[0][0][ipZ][ipZ] = HPhoto.hgamnc[ipZ][IP1S];
	PhotRate.PhotoRate[1][0][ipZ][ipZ] = heat.HeatLowEnr;
	PhotRate.PhotoRate[2][0][ipZ][ipZ] = heat.HeatHiEnr;

	/* option to print photoionization rates */
	if( trace.lgTrace && trace.lgHBug )
	{
		GammaPrt(nh.ipHn[ipZ][IP1S],
			rfield.nflux,
			OpacPoint.ipHOpac[ipZ][IP1S],
		  ioQQQ,
		  HPhoto.hgamnc[ipZ][IP1S],HPhoto.hgamnc[ipZ][IP1S]*0.05);
	}

	/* hgamsc is direct photioniz rate due to 
	 * bound compton scattering of very hard x-rays+Compton scat */
	HPhoto.hgamnc[ipZ][IP1S] += hgamsc.BoundCompton;

	/* all other continua, n<NHPLPHOT are actual fits to photo cross sections,
	 * NHPLPHOT is 10 in opacpoint.h */
	limit = MIN2(NHPLPHOT-1,nhlevel);
	for( n=IP2S; n <= limit; n++ )
	{
		if( n>4 && !opac.lgRedoStatic )
			break;
		if( HPhoto.lgHInducImp )
		{
			HPhoto.hgamnc[ipZ][n] = 
				GammaBn(
				nh.ipHn[ipZ][n],
			  nh.ipHn[ipZ][IP1S]-1,
			  OpacPoint.ipHOpac[ipZ][n],
			  hydro.HLevNIonRyd[ipZ][n],
			  &HPhoto.finduc[ipZ][n],
			  &HPhoto.cind[ipZ][n])*
			  PhotRate.PhotScaleOn;
		}
		else
		{
			HPhoto.hgamnc[ipZ][n] = 
				(float)(GammaK(nh.ipHn[ipZ][n],
				nh.ipHn[ipZ][IP1S]-1,
				OpacPoint.ipHOpac[ipZ][n],1.)*
				PhotRate.PhotScaleOn);
			HPhoto.finduc[ipZ][n] = 0.;
			HPhoto.cind[ipZ][n] = 0.;
		}
		HPHeat.HPhotHeat[ipZ][n] = heat.HeatNet;
	}

	/* need not redo in most cases */
	if( opac.lgRedoStatic )
	{
		/* these are nu^-3 powerlaws */
		for( n=NHPLPHOT; n <= nhlevel; n++ )
		{
			if( HPhoto.lgHInducImp )
			{
				HPhoto.hgamnc[ipZ][n] = GammaBnPL(n,ipZ,
					&HPhoto.finduc[ipZ][n],
				  &HPhoto.cind[ipZ][n])*
				  PhotRate.PhotScaleOn;
			}
			else
			{
				HPhoto.hgamnc[ipZ][n] = (float)GammaPL(n,ipZ)*PhotRate.PhotScaleOn;
				HPhoto.finduc[ipZ][n] = 0.;
				HPhoto.cind[ipZ][n] = 0.;
			}
			/* hgamnc(n,ipZ) = hgamnc(n,ipZ) */
			HPHeat.HPhotHeat[ipZ][n] = heat.HeatNet;
		}
	}

	/* now set induced recombination rate and cooling */
	if( HlteOk.lgHlteOk[ipZ] )
	{
		for( n=IP1S; n <= nhlevel; n++ )
		{
			/* induced recombination rate */
			HPhoto.rinduc[ipZ][n] = 
				(float)(HPhoto.finduc[ipZ][n]*hcolrt.hlte[ipZ][n]);

			/* cooling due to induced recombination */
			HPhoto.cinduc[ipZ][n] = 
				(float)(HPhoto.cind[ipZ][n]*hcolrt.hlte[ipZ][n]);
		}
	}
	else
	{
		for( n=IP1S; n <= nhlevel; n++ )
		{
			HPhoto.cinduc[ipZ][n] = 0.;
			HPhoto.rinduc[ipZ][n] = 0.;
		}
	}

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, "     HydroPhoto%2ld low, hi=",ipZ);
		fprintf( ioQQQ,PrintEfmt("%9.2e", HPhoto.hgamnc[ipZ][IP1S]));
		fprintf( ioQQQ,PrintEfmt("%9.2e", hgamsc.BoundCompton));
		fprintf( ioQQQ, " total=");
		fprintf( ioQQQ,PrintEfmt("%9.2e",HPhoto.hgamnc[ipZ][IP1S] ));
		fprintf( ioQQQ, "\n");
	}

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