/*MakeOpacity generate ionic subshell opacities by calling phfit */
#include "cddefines.h"
#include "physconst.h"
#include "opacpoint.h"
#include "rfield.h"
#include "ph1com.h"
#include "noptot.h"
#include "punpoint.h"
#include "nsshells.h"
#include "opac.h"
#include "showme.h"
#include "kshllenr.h"
#include "phfiton.h"
#include "phfit.h"
#include "makeopacity.h"

void MakeOpacity(
					  /* atomic number on the C scale, lowest ever called will be Li=2 */
					  long int nelem)
{
	int lgOK;
	long int ihi, 
	  ip, 
	  ipop, 
	  limit, 
	  low, 
	  need, 
	  nelec, 
	  nion, 
	  nshell;
	float cs; 
	double energy;

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

	/* if this flag ever false then we have problems */
	lgOK = TRUE;

	/* confirm range of validity of atomic number, Li=2 should be the lightest */
	assert( nelem >= 2 );
	assert( nelem < LIMELM );

	/*>>chng 99 jan 27, no longer redo hydrogenic opacity here */
	/* for( nion=0; nion <= nelem; nion++ )*/
	for( nion=0; nion < nelem; nion++ )
	{

		/* will be used for a sanity check on number of hits in a cell*/
		for( ip=0; ip < NCELL; ip++ )
		{
			opac.opac[ip] = 0.;
		}

		/* number of bound electrons */
		nelec = nelem+1 - nion;

		for( nshell=0; nshell < nsShellsCom.nsShells[nion][nelem]; nshell++ )
		{
			OpacPoint.ipElement[2][nshell][nion][nelem] = noptot.nOpacTot +  1;

			/* this is continuum pointer upper limit for energy range of this shell */
			limit = OpacPoint.ipElement[1][nshell][nion][nelem];

			/* this is number of cells in continuum needed to store opacity */
			need = limit - OpacPoint.ipElement[0][nshell][nion][nelem] + 1;

			/* check that opac will have enough frequeny cells */
			if( noptot.nOpacTot + need > NOPSV )
			{
				fprintf( ioQQQ, " MakeOpacity needs more opacity cells than NOPSV,  increase it.\n" );
				fprintf( ioQQQ, " element, Ion, are%3ld%3ld\n", 
				  nelem, nion );
				puts( "[Stop in makeopacity]" );
				exit(1);
			}

			/* set lower and upper limits to this range */
			low = OpacPoint.ipElement[0][nshell][nion][nelem];
			ihi = OpacPoint.ipElement[1][nshell][nion][nelem];
			ipop = OpacPoint.ipElement[2][nshell][nion][nelem];

			/* begin sanity check */
			if( low > ihi && low > 5 )
			{
				fprintf( ioQQQ, " MakeOpacity has insane pointers.\n" );
				fprintf( ioQQQ, " nelem,nion,shell%10ld%10ld%10ld%10ld%10ld\n", 
				  nelem, nion, nshell, low, ihi );
				ShowMe();
				puts( "[Stop in makeopacity]" );
				exit(1);
			}
			/*end sanity check */

			for( ip=low-1; ip < ihi; ip++ )
			{
				/* >>chng 96 oct 8 add min to ip, so never eval below threshold */
				energy = MAX2(rfield.AnuOrg[ip]*EVRYD , 
					PH1COM.PH1[nshell][nelec-1][nelem][0]);

				phfit(nelem+1,nelec,nshell+1,energy,&cs,PhFitOn.lgPhFit);
				opac.OpacStack[ip-low+ipop] = (float)(cs*1e-18);

				/*begin sanity check */
				opac.opac[ip] += cs;
				/*end sanity check */
			}

			noptot.nOpacTot += ihi - low + 1;

			/* punch pointers option */
			if( PunPoint.lgPunPoint )
			{
				fprintf( PunPoint.ipPoint, "%3ld%3ld%3ld%10.2e%10.2e%10.2e%10.2e\n", 
				  nelem, nion, nshell, rfield.anu[low-1], rfield.anu[ihi-1], 
				  opac.OpacStack[ipop-1], opac.OpacStack[ihi-low+ipop-1] );
			}
		}

		/*begin sanity check  */
		for( 
			ip=OpacPoint.ipElement[0][nsShellsCom.nsShells[nion][nelem]-1][nion][nelem]-1;
		   ip < KshllEnr.KshellLimit; ip++ )
		{
			if( opac.opac[ip] <= 0. )
			{
				lgOK = FALSE;
				fprintf( ioQQQ, 
					" MakeOpacity finds insane opacity, nelem, nu=%3ld%3ld%10.2e%10.2e\n", 
				  nelem, nion, rfield.anu[ip], opac.opac[ip] );
			}
		}
		/*end sanity check */

	}

	if( !lgOK )
	{
		fprintf( ioQQQ, " MakeOpacity: insanity has occurred.\n" );
		ShowMe();
		puts( "[Stop in makeopacity]" );
		exit(1);
	}


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