/*ContRate called by nextdr to find energy of maximum continuum-gas interaction */
#include "cddefines.h"
#include "showme.h"
#include "hydrogenic.h"
#include "grainvar.h"
#include "nh.h"
#include "opac.h"
#include "rfield.h"
#include "heavy.h"
#include "phycon.h"
#include "tff.h"
#include "heat.h"
#include "contrate.h"

void ContRate(double *freqm, 
  double *opacm)
{
	long int i, 
	  iplow, 
	  limit;
	double FreqH, 
	  FreqSub, 
	  OpacH, 
	  OpacSub, 
	  xMaxH, 
	  xMaxSub;

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

	/* 
	 * find maximum continuum interaction rate,
	 * these should be reset in following logic without exception,
	 * if they are still zero at the end we have a logical error 
	 */
	xMaxSub = 0.;
	FreqSub = 0.;
	OpacSub = 0.;

	/* do up to carbon photo edge */
	for( i=1; i < (Heavy.ipHeavy[0][5] - 1); i++ )
	{
		/* this does not have grain opacity since grains totally passive
		 * at energies below CI edge */
		if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*(opac.opac[i] - 
		  GrainVar.dstab[i]*phycon.hden) > xMaxSub )
		{
			xMaxSub = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
			  (opac.opac[i] - GrainVar.dstab[i]*phycon.hden);
			FreqSub = rfield.anu[i];
			OpacSub = opac.opac[i] - GrainVar.dstab[i]*phycon.hden;
		}
	}

	/* not every continuum extends beyond C edge-this whole loop can add to zero
	 * use total opacity here
	 * test is to put in fir continuum if free free heating is important */
	if( hydro.FreeFreeHeat/heat.htot > 0.05 )
	{
		/* this is pointer to energy where cloud free free optical depth is unity */
		iplow = tffCom.ntff;
	}
	else
	{
		iplow = Heavy.ipHeavy[0][5];
	}

	limit = MIN2(rfield.nflux,nh.ipHn[0][IP1S]-1);
	for( i=iplow-1; i < limit; i++ )
	{
		if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*(opac.opac[i] - 
		  GrainVar.dstab[i]*phycon.hden) > xMaxSub )
		{
			xMaxSub = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
			  (opac.opac[i] - GrainVar.dstab[i]*phycon.hden);
			FreqSub = rfield.anu[i];
			OpacSub = opac.opac[i] - GrainVar.dstab[i]*phycon.hden;
		}
	}

	/* variables to check continuum interactions over lyman continuum */
	xMaxH = 0.;
	OpacH = 0.;
	FreqH = 0.;

	/* not every continuum extends beyond 1 Ryd-this whole loop can add to zero */
	for( i=nh.ipHn[0][IP1S]-1; i < rfield.nflux; i++ )
	{
		if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*(opac.opac[i] - 
		  GrainVar.dstab[i]*phycon.hden) > xMaxH )
		{
			/* xMaxH = anu(i)*flux(i)/widflx(i)*opac(i) */
			xMaxH = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
			  (opac.opac[i] - GrainVar.dstab[i]*phycon.hden);
			FreqH = rfield.anu[i];
			OpacH = opac.opac[i] - GrainVar.dstab[i]*phycon.hden;
		}
	}

	/* use lyman continuum if its opacity is larger than non-h ion */
	if( xMaxSub < 1e-30 )
	{
		/* this happens for laser source - use Lyman continuum */
		*opacm = OpacH;
		*freqm = FreqH;
	}
	else if( OpacH > OpacSub && xMaxH/xMaxSub > 1e-10 )
	{
		/* use Lyman continuum */
		*opacm = OpacH;
		*freqm = FreqH;
	}
	else
	{
		/* not much rate in the lyman continuum, stick with low energy */
		*opacm = OpacSub;
		*freqm = FreqSub;
	}

	if( xMaxH > xMaxSub )
	{
		*opacm = OpacH;
		*freqm = FreqH;
	}
	else
	{
		*opacm = OpacSub;
		*freqm = FreqSub;
	}

	/*begin sanity check
	 * these were set to zero at start, and must have been reset if one of the
	 * two loops.  Logic error if still zero. 
	 */
	if( *opacm <= 0. || *freqm <= 0. )
	{
		fprintf( ioQQQ, " ContRate found non positive interaction, or energy.\n" );
		fprintf( ioQQQ, " This is impossible.\n" );
		ShowMe();
		puts( "[Stop in contrate]" );
		exit(1);
	}
	/*end sanity check*/

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