/*PhotoIonize fill array PhotoRate with photoionization rates for heavy elements */
#include "cddefines.h"
#include "opacpoint.h"
#include "heat.h"
#include "photrate.h"
#include "yield.h"
#include "showme.h"
#include "nsshells.h"
#include "opac.h"
#include "ionrange.h"
#include "elementnames.h"
#include "shell.h"
#include "gammas.h"
#include "photoionize.h"

void PhotoIonize(
	/* nlem is atomic number on C scale, 0 for H */
	long int nelem , 
	/* debugging flag to turn on print */
	int lgPrintIt )
{
	long int nion, 
	  iphi, 
	  iplow, 
	  ipop, 
	  limit, 
	  ns;

#	ifdef DEBUG_FUN
	fputs( "<+>PhotoIonize()\n", debug_fp );
#	endif
	
	/* IonLow(nelem) and IonHigh(nelem) are bounds for evaluation*/

	/* this puts it onto the c scale, which is used below */

	/*begin sanity check */
	if( (IonRange.IonLow[nelem] <= 0 || IonRange.IonLow[nelem] > nelem+1) || 
		IonRange.IonHigh[nelem] > nelem + 2 )
	{
		fprintf( ioQQQ, " PhotoIonize called with insane low, high:nelem=%4ld low, high=%10ld%10ld\n", 
		  nelem, IonRange.IonLow[nelem], IonRange.IonHigh[nelem] );
		ShowMe();
		puts( "[Stop in photoionize]" );
		exit(1);
	}
	/*end sanity check */

	/* NB - in following, nelem is on c scale, so is 29 for Zn */

	/* min since hydrogenic rates produced in hydrgn routine.
	 * IonHigh and IonLow range from 1 to nelem+1, but for photo rates
	 * we want 1 to nelem since cannot destroy fully stripped ion.  
	 * since hydrogenic done elsewhere, want to actually do IonHigh-2.
	 * with limit below, and with following for loop only going up to nion < limit,
	 * this works out ok.  for nelem=30-1 (Zn), IonHigh can be 31, and then
	 * limit is 29.  for loop will then go up to 28. With logic below do not
	 * need to fill in PhotoRate[0][1][nelemm][nelem] */
	limit = MIN2(nelem,IonRange.IonHigh[nelem]-1);

	for( nion=IonRange.IonLow[nelem]-1; nion < limit; nion++ )
	{
		for( ns=0; ns < nsShellsCom.nsShells[nion][nelem]; ns++ )
		{
			/* always reevaluate the outer shell, and all shells if lgRedoStatic is set */
			if( (ns==(nsShellsCom.nsShells[nion][nelem]-1) || opac.lgRedoStatic) )
			{
				/* option to redo the rates only on occasion */
				iplow = OpacPoint.ipElement[0][ns][nion][nelem];
				iphi = OpacPoint.ipElement[1][ns][nion][nelem];
				ipop = OpacPoint.ipElement[2][ns][nion][nelem];

				/* compute the photoionization rate,PhotScaleOn is 1, set 0
				 * with "no photoionization" command */
				PhotRate.PhotoRate[0][ns][nion][nelem] = 
					GammaK(iplow,iphi,
				  ipop,yield.vyield[0][ns][nion][nelem])*PhotRate.PhotScaleOn;

				/* these three lines must be kept parallel with the lines
				 * in GammaK nion*/

				/* the heating rate */
				PhotRate.PhotoRate[1][ns][nion][nelem] = heat.HeatLowEnr;
				PhotRate.PhotoRate[2][ns][nion][nelem] = heat.HeatHiEnr;

			}
		}
	}

	/* option to print information about these rates for this element */
	if( lgPrintIt )
	{
		/* option to print rates for particular shell */
		ns = 5;
		nion = 1;
		GammaPrt(
		  OpacPoint.ipElement[0][ns][nion][nelem], 
		  OpacPoint.ipElement[1][ns][nion][nelem], 
		  OpacPoint.ipElement[2][ns][nion][nelem], 
		  ioQQQ, /* io unit we will write to */
		  PhotRate.PhotoRate[0][ns][nion][nelem], 
		  0.05);

		/* outer loop is from K to most number of shells present in atom */
		for( ns=0; ns < nsShellsCom.nsShells[0][nelem]; ns++ )
		{
			fprintf( ioQQQ, "\n %s", elementnames.chElementNameShort[nelem] );
			fprintf( ioQQQ, " %s" , Shell.chShell[ns]);
			/* MB hydrogenic photo rate may not be included in beow */
			for( nion=0; nion < (IonRange.IonHigh[nelem] - 1); nion++ )
			{
				if( nsShellsCom.nsShells[nion][nelem] > ns )
				{
					fprintf( ioQQQ, " %8.1e", PhotRate.PhotoRate[0][ns][nion][nelem] );
				}
				else
				{
					break;
				}
			}
		}
	}


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