/*opacrm generate photoionization cross sections from Reilman and Manson points */
#include "cddefines.h"
#define	NOP	100
#include "opac.h"
#include "rfield.h"
#include "noptot.h"
#include "opacrm.h"

void opacrm(long int low, 
  long int ihi, 
  float cross[], 
  long int ncross, 
  long int *ipop, 
  char *chLabl)
{
	long int i, 
	  ics, 
	  j, 
	  ncr;

	float cs[NOP], 
	  en[NOP], 
	  slope;

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

	/* this is the opacity entering routine designed for
	 * the Reilman and Manson tables.  It works with incident
	 * photon energy (entered in eV) and cross sections in megabarns
	 * */
	*ipop = noptot.nOpacTot + 1;
	assert( *ipop > 0 );

	ncr = ncross/2;
	if( ncr > NOP )
	{
		fprintf( ioQQQ, " Too many opacities were entered into OPACRM.  Increase the value of NOP.\n" );
		fprintf( ioQQQ, " chLabl was %4.4s\n", chLabl );
		puts( "[Stop in opacrm]" );
		exit(1);
	}

	/* the array CROSS has ordered pairs of elements.
	 * the first is the energy in eV (not Ryd)
	 * and the second is the cross section in megabarns */
	for( i=0; i < ncr; i++ )
	{
		en[i] = (float)(cross[i*2]/13.6);
		cs[i] = (float)(cross[(i+1)*2-1]*1e-18);
	}

	if( en[0] > rfield.anu[low-1] )
	{
		fprintf( ioQQQ, 
			" OPACRM: The entered opacity energy bandwidth is not large enough (low fail).\n" );
		fprintf( ioQQQ, 
			" The desired energy (Ryd) was%12.5eeV and the lowest entered in the array was%12.5e eV\n", 
		  rfield.anu[low-1]*13.606, en[0]*13.606 );

		fprintf( ioQQQ, " chLabl was %4.4s\n", chLabl );
		fprintf( ioQQQ, " The original energy (eV) and cross section (mb) arrays follow:\n" );
		fprintf( ioQQQ, " " );

		for( i=1; i <= ncross; i++ )
		{
			fprintf( ioQQQ, "%11.4e", cross[i-1] );
		}

		fprintf( ioQQQ, "\n" );
		puts( "[Stop in opacrm]" );
		exit(1);
	}

	slope = (cs[1] - cs[0])/(en[1] - en[0]);
	ics = 1;

	/* now fill in the opacities using linear interpolation */
	for( i=low-1; i < ihi; i++ )
	{
		if( rfield.anu[i] > en[ics-1] && rfield.anu[i] <= en[ics] )
		{
			opac.OpacStack[i-low+*ipop] = cs[ics-1] + slope*(rfield.anu[i] - 
			  en[ics-1]);
		}

		else
		{
			ics += 1;
			if( ics + 1 > ncr )
			{
				fprintf( ioQQQ, " OPACRM: The entered opacity energy bandwidth is not large enough (high fail).\n" );
				fprintf( ioQQQ, " The entered energy was %10.2eeV and the highest in the array was %10.2eeV\n", 
				  rfield.anu[i]*13.6, en[ncr-1]*13.6 );
				fprintf( ioQQQ, " chLabl was %4.4s\n", chLabl
				   );
				fprintf( ioQQQ, " The lowest energy enterd in the array was%10.2e eV\n", 
				  en[0]*13.65 );
				fprintf( ioQQQ, " The highest energy ever needed would be%10.2eeV\n", 
				  rfield.anu[ihi-1]*13.6 );
				fprintf( ioQQQ, " The lowest energy needed was%10.2eeV\n", 
				  rfield.anu[low-1]*13.6 );
				puts( "[Stop in opacrm]" );
				exit(1);
			}

			slope = (cs[ics] - cs[ics-1])/(en[ics] - en[ics-1]);
			if( rfield.anu[i] > en[ics-1] && rfield.anu[i] <= en[ics] )
			{
				opac.OpacStack[i-low+*ipop] = cs[ics-1] + slope*(rfield.anu[i] - 
				  en[ics-1]);
			}
			else
			{
				fprintf( ioQQQ, " Internal logical error in OPACRM.\n" );
				fprintf( ioQQQ, " The desired energy (%10.2eeV), I=%5ld, is not within the next energy bound%10.2e%10.2e\n", 
				  rfield.anu[i]*13.6, i, en[ics-1], en[ics] );

				fprintf( ioQQQ, " The previous energy (eV) was%10.2e\n", 
				  rfield.anu[i-1]*13.6 );

				fprintf( ioQQQ, " Here comes the energy array.  ICS=%4ld\n", 
				  ics );

				for( j=0; j < ncr; j++ )
				{
					fprintf( ioQQQ, "%10.2e", en[j] );
				}
				fprintf( ioQQQ, "\n" );

				fprintf( ioQQQ, " chLabl was %4.4s\n", chLabl );
				puts( "[Stop in opacrm]" );
				exit(1);
			}
		}
	}

	noptot.nOpacTot += ihi - low + 1;
	if( noptot.nOpacTot > NOPSV )
	{
		fprintf( ioQQQ, " OPACPL too many opacity points.\n" );
		fprintf( ioQQQ, " Increase parameter variable NOPSV everywhere it occurs in the code.\n" );
		fprintf( ioQQQ, " The current value is%6ld.\n", NOPSV );
		puts( "[Stop in opacrm]" );
		exit(1);
	}


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