/*opacpl generate array of cross sections using a simple power law fit */
#include "cddefines.h"
#include "opac.h"
#include "rfield.h"
#include "noptot.h"
#include "opacpl.h"

void opacpl(
	/* lower energy limit on continuum mesh */
	long int ilo, 
	/* upper energy limit on continuum mesh */
	long int ihi, 
	/* threshold cross section */
	double cross, 
	/* power law index */
	double s, 
	/* pointer to opacity offset where this starts */
	long int *ip)
{
	long int i;
	double thres;

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

	/* non-positive cross section is unphysical */
	assert( cross > 0. );

	/* place in the opacity stack where we will stuff cross sections */
	*ip = noptot.nOpacTot + 1;
	assert( *ip > 0 );
	assert( ilo > 0 );
	thres = rfield.anu[ilo-1];

	for( i=ilo-1; i < ihi; i++ )
	{
		opac.OpacStack[i-ilo+*ip] = (float)(cross*pow(rfield.anu[i]/thres,-s));
	}

	noptot.nOpacTot += ihi - ilo + 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 opacpl]" );
		exit(1);
	}


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