/*GrainRateDr called by nextdr to find grains heating rate dr */
#include "cddefines.h"
#include "showme.h"
#include "hydrogenic.h"
#include "opac.h"
#include "rfield.h"
#include "heavy.h"
#include "tff.h"
#include "heat.h"
#include "grainratedr.h"

void GrainRateDr(double *grfreqm, 
  double *gropacm)
{
	long int i, 
	  iplow;
	double xMax;

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

	/* in all following changed from anu2 to anu  july 25 95
	 *
	 * find maximum continuum interaction rate */

	/* 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];
	}

	xMax = 0.;
	/* integrate up to heii edge */
	for( i=iplow-1; i < Heavy.ipHeavy[1][1]; i++ )
	{
		if( rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*opac.opac[i] > xMax )
		{
			xMax = rfield.anu[i]*rfield.flux[i]/rfield.widflx[i]*
			  opac.opac[i];
			*grfreqm = rfield.anu[i];
			*gropacm = opac.opac[i];
		}
	}

	/*begin sanity check */
	if( *gropacm <= 0. || *grfreqm <= 0. )
	{
		fprintf( ioQQQ, " GrainRateDr found non positive interaction, or energy.\n" );
		fprintf( ioQQQ, " This is impossible.\n" );
		ShowMe();
		puts( "[Stop in grainratedr]" );
		exit(1);
	}
	/*end sanity check */

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