/*
 * these are the routines that evaluate the photoionization rates, gammas,
 * throughout cloudy.  a considerable amount of time is spent in these routines,
 * so they should be compiled at the highest possible efficientcy.  
 * The routine are:
 *
 * GammaBn   evaluate photoionization rate for single shell with induced recomb
 * GammaBnPL evaluate photoionization rate for single shell with induced recomb 
 * GammaPrt  special version of gamma function to print strong contributors 
 * GammaK    evaluate photoionization rate for single shell 
 * GammaPL   evaluate photoionization rate for power law photo cross section 
 * GammaPrtRate will print photo rates for all shells of a ion and element 
 */

/*>>chng 99 apr 16, outcon had not been included in the threshold point for any
 * of these routines.  added it.  also moved thresholds above loop for a few */

/*GammaBn evaluate photoionization rate for single shell with induced recomb */
#include "cddefines.h"
#include "physconst.h"
#include "fluors.h"
#include "secondaries.h"
#include "opac.h"
#include "rfield.h"
#include "photrate.h"
#include "heat.h"
#include "linelabl.h"
#include "printe82.h"
#include "hydrophocs.h"
#include "nh.h"
#include "nsshells.h"
#include "hydrogenic.h"
#include "gammas.h"

double GammaBn(long int ipLoEnr, 
  long int ipHiEnr, 
  long int ipOpac, 
  double thresh, 
  double *ainduc, 
  double *rcool)
{
	long int i, 
	  ilo, 
	  iup, 
	  limit;
	float GamHi, 
	  bnfun_v, 
	  emin, 
	  g, 
	  phisig, 
	  RateInducRec, 
	  RateInducRecCool, 
	  prod;

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

/**********************************************************************
 *
 * special version of GAMFUN for arbitrary threshold, and induc rec
 * compute ainduc, rate of inducted recombinations, rcool, induc rec cool
 *
 **********************************************************************/

	/* special version of GAMFUN for arbitrary threshold, and induc rec
	 * compute ainduc, rate of inducted recombinations, rcool, induc rec cool */
	if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr )
	{
		bnfun_v = 0.;
		heat.HeatNet = 0.;
		heat.HeatLowEnr = 0.;
		heat.HeatHiEnr = 0.;
		*ainduc = 0.;
		*rcool = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->GammaBn()\n", debug_fp );
#		endif
		return( bnfun_v );
	}

	/* >>chng 96 oct 25, first photo elec energy may have been negative
	 * possible for first any point to be below threshold due to
	 * finite resolution of continuum mesh */
	emin = (float)MAX2(thresh,rfield.anu[ipLoEnr-1]);
	/* >>chng 99 Jun 08 heating systematically too small, change to correct
	 * threshold, and protect first cell */
	emin = (float)thresh;

	heat.HeatNet = 0.;
	g = 0.;
	RateInducRec = 0.;
	RateInducRecCool = 0.;
	iup = MIN2(ipHiEnr,rfield.nflux);

	/* do first point without otscon, which may have diffuse cont,
	 * extra minus one after ipOpac is correct since not present in i = */
	i = ipLoEnr;
	/* >>chng 99 apr 16, add outcon */
	phisig = (rfield.flux[i-1] + rfield.otslin[i-1] + rfield.outcon[i-1] )*
		opac.OpacStack[i-ipLoEnr+ipOpac-1];

	g += phisig;
	heat.HeatNet += phisig*rfield.anu[i-1];

	/* integral part of induced rec rate */
	prod = phisig*rfield.ContBoltz[i-1];
	RateInducRec += prod;
	/* incuded rec cooling */
	RateInducRecCool += prod*(rfield.anu[i-1] - emin);

	limit = MIN2(iup,Secondaries.ipSecIon-1);

	/* these is no -1 after the ipLoEnr in the following, due to c off by one counting */
	for( i=ipLoEnr; i < limit; i++ )
	{
		/* SummedCon defined in SumContinuum, 
		 * includes flux, otscon, otslin, outcon, outlin */
		phisig = rfield.SummedCon[i]*opac.OpacStack[i-ipLoEnr+ipOpac];

		g += phisig;
		heat.HeatNet += phisig*rfield.anu[i];

		/* phi-sigma times exp(-hnu/kT) */
		prod = phisig*rfield.ContBoltz[i];
		/* induced recombination rate */
		RateInducRec += prod;
		/* incuded rec cooling */
		RateInducRecCool += prod*(rfield.anu[i] - emin);
	}

	/* heating in Rydbergs, no secondary ionization */
	heat.HeatNet = heat.HeatNet - g*thresh;

	/* we will save this heat - the part that does not secondary ionize */
	heat.HeatLowEnr = heat.HeatNet;

	/* now do part where secondary ionization is possible */
	heat.HeatHiEnr = 0.;
	GamHi = 0.;

	/* make sure we don't count twice when ipSecIon = ipLoEnr */
	ilo = MAX2(ipLoEnr+1,Secondaries.ipSecIon);
	for( i=ilo-1; i < iup; i++ )
	{
		phisig = rfield.SummedCon[i]*opac.OpacStack[i-ipLoEnr+ipOpac];

		GamHi += phisig;
		heat.HeatHiEnr += phisig*rfield.anu[i];

		/* integral part of induced rec rate */
		prod = phisig*rfield.ContBoltz[i];
		RateInducRec += prod;
		/* incuded rec cooling */
		RateInducRecCool += prod*(rfield.anu[i] - emin);
	}

	/* this is total photo rate, low and high energy parts */
	bnfun_v = g + GamHi;

	/* heating in Rydbergs, uncorrected for secondary ionization */
	heat.HeatHiEnr = heat.HeatHiEnr - GamHi*thresh;

	/* add corrected high energy heat to total */
	heat.HeatNet += heat.HeatHiEnr*Secondaries.heatef;

	/* final product is heating in erg */
	heat.HeatNet *= EN1RYD;
	heat.HeatHiEnr *= EN1RYD;
	heat.HeatLowEnr *= EN1RYD;

	/* this is an option to turn off induced processes */
	if( fluors.lgFluor )
	{
		*rcool = RateInducRecCool*EN1RYD;
		*ainduc = RateInducRec;
	}
	else
	{
		/* the "no induced" command was given */
		*rcool = 0.;
		*ainduc = 0.;
	}


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

/**********************************************************************
 *
 * GammaPrt special version of gamma function to print strong contributors
 * this only prints info - does not update heating rates properly since
 * this is only a diagnostic
 *
 **********************************************************************/

void GammaPrt(
	  /* first three par are identical to GammaK */
	  long int ipLoEnr, 
	  long int ipHiEnr, 
	  long int ipOpac, 
	  /* io unit we will write to */
	  FILE * ioFILE, 
	  /* total photo rate from previous call to GammaK */
	  double total, 
	  /* we will print contributors that are more than this rate */
	  double threshold)
{
	long int i, 
	  j, 
	  k;
	float flxcor, 
	  phisig;

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

	/* special special version of GAMFUN to punch step-by-step results */
	/* >>chng 99 apr 02, test for lower greater than higher (when shell
	 * does not exist) added.  This caused incorrect photo rates for
	 * non-existant shells */
	if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr )
	{ 
#		ifdef DEBUG_FUN
		fputs( " <->GammaPrt()\n", debug_fp );
#		endif
		return;
	}

	fprintf( ioFILE, " GammaPrt from ");
	fprintf(ioFILE,PrintEfmt("%9.2e",rfield.anu[ipLoEnr-1]));
	fprintf( ioFILE, " to ");
	fprintf(ioFILE,PrintEfmt("%9.2e",rfield.anu[ipHiEnr-1]));
	fprintf( ioFILE, "R rates >");
	fprintf(ioFILE,PrintEfmt("%9.2e",threshold));
	fprintf( ioFILE, " of total=");
	fprintf(ioFILE,PrintEfmt("%9.2e",total));
	fprintf( ioFILE, " (frac inc, otslin, otscon, outcon+lin) chL, C\n");

	if( threshold <= 0. || total <= 0. )
	{ 
#		ifdef DEBUG_FUN
		fputs( " <->GammaPrt()\n", debug_fp );
#		endif
		return;
	}

	k = ipLoEnr;
	j = MIN2(ipHiEnr,rfield.nflux);

	/* do theshold special, do not pick up otscon */
	i = k-1;

	phisig = (rfield.flux[i] + rfield.otslin[i]+ rfield.outcon[i] )*
		opac.OpacStack[i-k+ipOpac];
	if( phisig > threshold )
	{
		flxcor = rfield.flux[i] + rfield.otslin[i] + rfield.outcon[i];

		fprintf( ioFILE, "%5ld", i );
		fprintf(ioFILE,PrintEfmt("%9.2e",rfield.anu[i]));
		fprintf(ioFILE,PrintEfmt("%9.2e",phisig/total));
		fprintf( ioFILE, "%5.2f%5.2f%5.2f%5.2f %4.4s %4.4s \n", 
		  rfield.flux[i]/flxcor, 
		  rfield.otslin[i]/flxcor, 
		  rfield.otscon[i]/flxcor, 
		  (rfield.outcon[i] + rfield.outlin[i])/flxcor, 
		  LineLabl.chLineLabel[i], 
		  LineLabl.chContLabel[i] );
	}

	for( i=k; i < j; i++ )
	{
		phisig = rfield.SummedCon[i]*opac.OpacStack[i-k+ipOpac];
		if( phisig > threshold )
		{
			flxcor = rfield.flux[i] + rfield.otslin[i] + rfield.otscon[i] + 
			  rfield.outlin[i] + rfield.outcon[i];

			fprintf( ioFILE, "%5ld", i+1 );
			fprintf(ioFILE,PrintEfmt("%9.2e",rfield.anu[i]));
			fprintf(ioFILE,PrintEfmt("%9.2e",phisig/total));
			fprintf( ioFILE, "%5.2f%5.2f%5.2f%5.2f %4.4s %4.4s \n", 
			  rfield.flux[i]/flxcor, 
			  rfield.otslin[i]/flxcor, 
			  rfield.otscon[i]/flxcor, 
			  (rfield.outcon[i] + rfield.outlin[i])/flxcor, 
			  LineLabl.chLineLabel[i], 
			  LineLabl.chContLabel[i] );
		}
	}


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

/*GammaK evaluate photoionization rate for single shell */

/* this routine is a major pace setter for the code
 * carefully check anything put into the following do loop */

double GammaK(
	/* ipLoEnr and ipHiEnr are pointers within frequency array to upper and
	 * lower bounds of evaluation */
	long int ipLoEnr, 
	long int ipHiEnr, 
	/* ipOpac is offset to map onto OPSV*/
	long int ipOpac, 
	/* yield1 is fraction of ionizations that emit 1 electron only,
	 * only used for energy balance */
	double yield1)
{
	long int i, 
	  ilo, 
	  ipOffset, 
	  iup, 
	  limit;
	float GamHi, 
	  eauger; 
	  
	float gamk_v ,
		phisig; 

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

	/* evaluate photoioinzation rate and photoelectric heating
	 * returns photoionization rate as GAMK, heating is H	 */

	/* >>chng 99 apr 02, test for lower greater than higher (when shell
	 * does not exist) added.  This caused incorrect photo rates for
	 * non-existant shells */
	if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr)
	{
		gamk_v = 0.;
		heat.HeatNet = 0.;
		heat.HeatHiEnr = 0.;
		heat.HeatLowEnr = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->GammaK()\n", debug_fp );
#		endif
		return( gamk_v );
	}

	iup = MIN2(ipHiEnr,rfield.nflux);

	/* anu(i) is threshold, assume each auger electron has this energy
	 * less threshold energy, IP lost in initial photoionizaiton
	 * yield1 is the fraction of photos that emit 1 electron */
	eauger = rfield.anu[ipLoEnr-1]*(float)yield1;

	/* low energies where secondary ionization cannot occur
	 * will do threshold point, ipLoEnr, later */
	gamk_v = 0.;
	heat.HeatNet = 0.;

	/* set up total offset for pointer addition not in loop */
	ipOffset = ipOpac - ipLoEnr;

	/* >>>chng 99 apr 16, this had followed the loop below, moved here to 
	 * be like rest of gamma functions */
	/* first do the threshold point
	 * do not include otscon, which may contain diffuse continuum */
	phisig = (rfield.flux[ipLoEnr-1] + rfield.otslin[ipLoEnr-1]+ rfield.outcon[ipLoEnr-1])*
		opac.OpacStack[ipOpac-1];
	gamk_v += phisig;
	heat.HeatNet += phisig*rfield.anu[ipLoEnr-1];

	/* this loop only executed for energies than cannot sec ioniz */
	limit = MIN2(iup,Secondaries.ipSecIon-1);
	for( i=ipLoEnr; i < limit; i++ )
	{
		phisig = rfield.SummedCon[i]*opac.OpacStack[i+ipOffset];
		gamk_v += phisig;
		heat.HeatNet += phisig*rfield.anu[i];
	}

	/* correct heating for work function */
	heat.HeatNet += -gamk_v*eauger;

	/* this is the total low-energy heating, in ryd, that cannot secondary ionize */
	heat.HeatLowEnr = heat.HeatNet;

	/* now do part where secondary ionization possible */
	heat.HeatHiEnr = 0.;
	GamHi = 0.;
	/* make sure we don't count twice when ipSecIon = ipLoEnr */
	ilo = MAX2(ipLoEnr+1,Secondaries.ipSecIon);
	for( i=ilo-1; i < iup; i++ )
	{
		/* produce of flux and cross section */
		phisig = rfield.SummedCon[i]*opac.OpacStack[i+ipOffset];
		GamHi += phisig;
		heat.HeatHiEnr += phisig*rfield.anu[i];
	}

	/* add on the high energy part */
	gamk_v += GamHi;
	/* correct for work function */
	heat.HeatHiEnr -= GamHi*eauger;

	/* net heating include high energy heat, with correction for
	 * secondary ionization */
	heat.HeatNet += heat.HeatHiEnr*Secondaries.heatef;

	/* final product is heating in erg */
	heat.HeatNet *= EN1RYD;
	heat.HeatLowEnr *= EN1RYD;
	heat.HeatHiEnr *= EN1RYD;

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


/*GammaPL evaluate photoionization rate for power law photoionization cross section */

double GammaPL(
	/* this is prin quantum number for level */
	long int n, 
	/* ipZ = 0 for H */
	long int ipZ)
{
	long int i, 
	  ilo, 
	  iup, 
	  limit, 
	  ipLoEnr, 
	  ipHiEnr;
	float GamHi, 
	  csThresh, 
	  gamPL_v, 
	  phisig;

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

	/* evaluate photoioinzation rate and photoelectric heating
	 * returns photoionization rate as GammaPL, heating is H
	 * n and ipZ and principal quantum number and charge
	 */

	/* validate the input */
	assert( ipZ >= 0);
	assert( ipZ < LIMELM);
	assert( n >= 0);

	/* get pointer to photo threshold */
	ipLoEnr = nh.ipHn[ipZ][n];
	assert( ipLoEnr>0 );

	if( ipLoEnr >= rfield.nflux )
	{
		gamPL_v = 0.;
		heat.HeatNet = 0.;
		heat.HeatHiEnr = 0. ;
		heat.HeatLowEnr = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->GammaPL()\n", debug_fp );
#		endif
		return( gamPL_v );
	}

	ipHiEnr = nh.ipHn[ipZ][IP2S];
	iup = MIN2(ipHiEnr,rfield.nflux);

	csThresh = HydroPhoCS.STH[n]*rfield.anu3[ipLoEnr-1]/ POW2(ipZ+1.f);

	/* low energies where secondary ionization cannot occur
	 * will do threshold point, ipLoEnr, later */
	gamPL_v = 0.;
	heat.HeatNet = 0.;

	/* first do the threshold point
	 * do not include otscon, which may contain diffuse continuum */
	phisig = (rfield.flux[ipLoEnr-1] + rfield.otslin[ipLoEnr-1] + 
		rfield.outcon[ipLoEnr-1] ) /rfield.anu3[ipLoEnr-1];
	gamPL_v += phisig;

	/* this loop only executed for energies than cannot sec ioniz */
	limit = MIN2(iup,Secondaries.ipSecIon-1);
	for( i=ipLoEnr; i < limit; i++ )
	{
		phisig = rfield.SummedCon[i]/rfield.anu3[i];
		gamPL_v += phisig;
		heat.HeatNet += phisig*rfield.anu[i];
	}

	/* convert photo rate into proper units with CS at threshold */
	gamPL_v *= csThresh;

	/* correct heating for CS at threshold and work function */
	heat.HeatNet *= csThresh;
	heat.HeatNet += -gamPL_v*rfield.anu[ipLoEnr-1];
	heat.HeatLowEnr = heat.HeatNet;

	/* now do part where secondary ionization possible */
	heat.HeatHiEnr = 0.;
	GamHi = 0.;

	/* make sure we don't count twice when ipSecIon = ipLoEnr */
	ilo = MAX2(ipLoEnr+1,Secondaries.ipSecIon);
	for( i=ilo-1; i < iup; i++ )
	{
		/* produce of flux and cross section */
		phisig = rfield.SummedCon[i]/rfield.anu3[i];
		GamHi += phisig;
		heat.HeatHiEnr += phisig*rfield.anu[i];
	}

	/* add on the high energy part */
	gamPL_v += GamHi*csThresh;

	/* correct heating for work function and convert to true cross section */
	heat.HeatHiEnr = (heat.HeatHiEnr - GamHi*rfield.anu[ipLoEnr-1])*csThresh;

	/* net heating includes high energy heating corrected for secondary ionizations */
	heat.HeatNet += heat.HeatHiEnr*Secondaries.heatef;

	/* final product is heating in erg */
	heat.HeatNet *= EN1RYD;
	heat.HeatHiEnr *= EN1RYD;
	heat.HeatLowEnr *= EN1RYD;

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

/*GammaBnPL evaluate hydrogenic photoionization rate for single level with induced recomb */

double GammaBnPL(long int n, 
  long int ipZ, /* 0 for H, etc */
  double *ainduc, 
  double *rcool)
{
	long int i, 
	  ilo, 
	  iup, 
	  limit, 
	  ipLoEnr, 
	  ipHiEnr;
	float GamHi, 
	  bnfunPL_v, 
	  csThresh, 
	  emin, 
	  g, 
	  phisig, 
	  RateInducRec, 
	  RateInducRecCool, 
	  prod, 
	  thresh;

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

	/* validate the incoming charge */
	assert( ipZ >= 0);
	assert( ipZ < LIMELM );

	/* special version of GAMFUN for arbitrary threshold, and induc rec
	 * compute ainduc, rate of inducted recombinations, rcool, induc rec cool */

	/* these are the upper and lower energy bounds, which we know from iso-sequence */
	ipLoEnr = nh.ipHn[ipZ][n];
	ipHiEnr = nh.ipHn[ipZ][IP1S];

	ilo = ipLoEnr;
	iup = MIN2(ipHiEnr,rfield.nflux);
	if( ipLoEnr >= rfield.nflux || ipLoEnr >= ipHiEnr )
	{
		bnfunPL_v = 0.;
		heat.HeatNet = 0.;
		heat.HeatHiEnr = 0. ;
		heat.HeatLowEnr = 0.;
		*ainduc = 0.;
		*rcool = 0.;
		
#		ifdef DEBUG_FUN
		fputs( " <->GammaBnPL()\n", debug_fp );
#		endif
		return( bnfunPL_v );
	}

	/* >>chng 96 oct 25, first photo elec energy may have been negative
	 * possible for first any point to be below threshold due to
	 * finite resolution of continuum mesh */
	thresh = hydro.HLevNIonRyd[ipZ][n];
	emin = (float)MAX2(thresh,rfield.anu[ipLoEnr-1]);
	/* >>chng 99 Jun 08 heating systematically too small, change to correct
	 * threshold, and protect first cell */
	emin = thresh;

	heat.HeatNet = 0.;
	g = 0.;
	RateInducRec = 0.;
	RateInducRecCool = 0.;

	/*csThresh = HydroPhoCS.STH[n]*POW3((float)thresh) / POW2(ipZ+1.f);*/
	csThresh = HydroPhoCS.STH[n]*rfield.anu3[ipLoEnr-1] / POW2(ipZ+1.f);

	/* do first point without otscon, which may have diffuse cont */
	phisig = (rfield.flux[ipLoEnr-1] + rfield.otslin[ipLoEnr-1]+ rfield.outcon[ipLoEnr-1])/
		rfield.anu3[ipLoEnr-1];

	g += phisig;
	heat.HeatNet += phisig*rfield.anu[ipLoEnr-1];

	/* integral part of induced rec rate */
	prod = phisig*rfield.ContBoltz[ipLoEnr-1];
	RateInducRec += prod;
	/* induced rec cooling */
	RateInducRecCool += prod*(rfield.anu[ipLoEnr-1] - emin);

	limit = MIN2(iup,Secondaries.ipSecIon-1);
	for( i=ipLoEnr; i < limit; i++ )
	{
		phisig = rfield.SummedCon[i]/rfield.anu3[i];

		g += phisig;
		heat.HeatNet += phisig*rfield.anu[i];

		/* integral part of induced rec rate */
		prod = phisig*rfield.ContBoltz[i];
		RateInducRec += prod;
		/* incuded rec cooling */
		RateInducRecCool += prod*(rfield.anu[i] - emin);
	}

	/* heating in Rydbergs, no secondary ionization */
	heat.HeatNet = (heat.HeatNet - g*thresh)*csThresh;
	heat.HeatLowEnr = heat.HeatNet;

	/* now do part where secondary ionization is possible */
	heat.HeatHiEnr = 0.;
	GamHi = 0.;
	/* make sure we don't count twice when ipSecIon = ipLoEnr */
	ilo = MAX2(ipLoEnr+1,Secondaries.ipSecIon);
	for( i=ilo-1; i < iup; i++ )
	{
		phisig = rfield.SummedCon[i]/rfield.anu3[i];
		GamHi += phisig;
		heat.HeatHiEnr += phisig*rfield.anu[i];

		/* integral part of induced rec rate */
		prod = phisig*rfield.ContBoltz[i];
		RateInducRec += prod;
		/* incuded rec cooling */
		RateInducRecCool += prod*(rfield.anu[i] - emin);
	}
	RateInducRec *= csThresh;
	RateInducRecCool *= csThresh;

	/* this is total photo rate, low and high energy parts */
	bnfunPL_v = (g + GamHi)*csThresh;

	/* heating in Rydbergs, uncorrected for secondary ionization */
	heat.HeatHiEnr = (heat.HeatHiEnr - GamHi*thresh)*csThresh;

	/* net heating includes high energy heat corrected for secondary ionization */
	heat.HeatNet += heat.HeatHiEnr*Secondaries.heatef;

	/* final product is heating in erg */
	heat.HeatNet *= EN1RYD;
	heat.HeatHiEnr *= EN1RYD;
	heat.HeatLowEnr *= EN1RYD;

	/* this is an option to turn off induced processes */
	if( fluors.lgFluor )
	{
		*rcool = RateInducRecCool*EN1RYD;
		*ainduc = RateInducRec;
	}
	else
	{
		/* the "no induced" command was given */
		*rcool = 0.;
		*ainduc = 0.;
	}

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

/* GammaPrtRate will print resulting rates for ion and element */
void GammaPrtRate(
	/* io unit we will write to */
	FILE * ioFILE, 
	/* stage of ionization on C scale, 0 for atom */
	long int ion ,
	/* 0 for H, etc */
	long int ipZ)
{
	long int i , ns;

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

	/* number of shells for this species */
	ns = nsShellsCom.nsShells[ion][ipZ];

	/* now print subshell photo rate */
	fprintf(ioFILE , "GammaPrtRate: %li %li",ion , ipZ );
	for( i=0; i<ns ; ++i )
	{
		fprintf(ioFILE , " %.2e" ,  PhotRate.PhotoRate[0][i][ion][ipZ]);
	}
	fprintf(ioFILE , "\n");

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