/*ipShells assign continuum energy pointers to shells for all atoms,
 * called by CreatePoint */
#include <string.h>
#include <float.h>
#include "cddefines.h"
#include "physconst.h"
#include "opacpoint.h"
#include "rfield.h"
#include "heavy.h"
#include "nh.h"
#include "ph1com.h"
#include "trace.h"
#include "kshllenr.h"
#include "elementnames.h"
#include "nsshells.h"
#include "shell.h"
#include "ipoint.h"
#include "ipconsafe.h"
#include "iplinsafe.h"
#include "ipshells.h"
#include "limitsh.h"

void ipShells(
				  /* nelem is the atomic number on the C scale, Li is 2 */
				  long int nelem)
{
	char chLab[5];
	long int 
	  imax, 
	  ion, 
	  nelec, 
	  ns, 
	  nshell;
	/* following value cannot be used - will be set to proper threshold */
	double thresh=-DBL_MAX;

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

	assert( nelem >= 2);
	assert( nelem < LIMELM );

	/* fills in pointers to valence shell ionization threshold
	 * PH1(a,b,c,d)
	 * a=1 => thresh, others fitting parameters
	 * b atomic number
	 * c number of electrons
	 * d shell number 7-1 */

	/* threshold in Ryd
	 * ion=0 for atom, up to nelem-1 for helium like, hydrogenic is elsewhere */
	for( ion=0; ion < nelem; ion++ )
	{
		/* number of bound electrons */
		nelec = nelem+1 - ion ;

		/* nsShells(nelem,ion) is the number of shells for ion with nelec electrons,
		 * physical not c scale */
		imax = nsShellsCom.nsShells[ion][nelem];

		/* loop on all inner shells, valence shell */
		for( nshell=0; nshell < imax; nshell++ )
		{
			/* ionization potential of this shell in rydbergs */
			thresh = (double)(PH1COM.PH1[nshell][nelec-1][nelem][0]/EVRYD* 0.9998787);
			if( thresh <= 0.1 )
			{
				/* negative ip shell does not exist, set upper limit
				 * to less than lower limit so this never looped upon
				 * these are used as flags by LimitSh to check whether
				 * this is a real shell - if 1 or 2 is changed - change LimitSh!! */
				OpacPoint.ipElement[0][nshell][ion][nelem] = 2;
				OpacPoint.ipElement[1][nshell][ion][nelem] = 1;
			}
			else
			{
				/* this is lower limit to range */
				OpacPoint.ipElement[0][nshell][ion][nelem] = ipoint(thresh);

				/* this is upper limit to range
				 * LimitSh is an integer function, returns pointer
				 * to threshold of next major shell.  For k-shell it
				 * returns the values KshellLimit, default=7.35e4
				 * >>chng 96 sep 26, had been below, result zero cross sec at 
				 * many energies where opacity project did not produce state specific 
				 * cross section */
				OpacPoint.ipElement[1][nshell][ion][nelem] = 
					LimitSh(ion+1,  nshell+1,nelem+1);
			}
		}

		/* this is valence pointer
		 * generate label for ionization stage */
		/* this is short form of element name */
		strcpy( chLab, elementnames.chElementSym[nelem] );

		/* this is a number between 1 and 31 */
		strcat( chLab, elementnames.chIonStage[ion] );

		/* [0] is pointer to threshold in energy array */
		OpacPoint.ipElement[0][imax-1][ion][nelem] = 
			ipConSafe(thresh, chLab);

		/* pointer to valence electron ionization potential */
		Heavy.ipHeavy[ion][nelem] = OpacPoint.ipElement[0][imax-1][ion][nelem];

		/* this is set of 3/4 of valence shell IP, this is important
		 * source of ots deep in cloud */
		Heavy.ipLyHeavy[ion][nelem] = 
			ipLinSafe(thresh*0.75,chLab , 0);
	}

	/* above loop did up to hydrogenic, now do hydrogenic - 
	 * hydrogenic is special since arrays already set up */
	nsShellsCom.nsShells[nelem][nelem] = 1;

	/* this is lower limit to range */
	/* hydrogenic photoionization set to special hydro array 
	 * this is pointer to threshold energy */
	OpacPoint.ipElement[0][0][nelem][nelem] = nh.ipHn[nelem][IP1S];

	/* this is the high energy limit */
	OpacPoint.ipElement[1][0][nelem][nelem] = KshllEnr.KshellLimit;

	Heavy.ipHeavy[nelem][nelem] = OpacPoint.ipElement[0][0][nelem][nelem];

	/* this is for backwards compatability with cambridge code */
	if( trace.lgTrace && trace.lgPointBug )
	{
		for( ion=0; ion < (nelem+1); ion++ )
		{
			fprintf( ioQQQ, "Ion:%3ld%3ld %2.2s%2.2s total shells:%3ld\n", 
			  nelem, ion+1, elementnames.chElementSym[nelem], elementnames.chIonStage[ion]
			  , nsShellsCom.nsShells[ion][nelem] );
			for( ns=0; ns < nsShellsCom.nsShells[ion][nelem]; ns++ )
			{
				fprintf( ioQQQ, " shell%3ld %2.2s range eV%10.2e-%8.2e\n", 
				  ns+1, Shell.chShell[ns], rfield.anu[OpacPoint.ipElement[0][ns][ion][nelem]-1]*
				  EVRYD, rfield.anu[OpacPoint.ipElement[1][ns][ion][nelem]-1]*EVRYD );
			}
		}
	}

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