/*dintg compute total radiative cooling due to large grains */
#include "cddefines.h"
#include "physconst.h"
#include "grainvar.h"
#include "rfield.h"
#include "trace.h"
#include "powi.h"
#include "dintg.h"

double dintg(double tdust, 
	/* pointer to type of grain */
  long int nd, 
  long int ip)
{
	long int i;

	double Planck1, 
	  arg, 
	  dintg_v, 
	  TDustRyg, 
	  sum ,
	  ExpM1 ;

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

	/******************************************************************
	 *
	 * >>>chng 99 mar 12, this sub rewritten following Peter van Hoof
	 * comments.  Oringinal coding was in single precision, and for
	 * very low temperature the exponential was set to zero.  As 
	 * a result Q was far too large for grain temperatures below 10K
	 *
	 ******************************************************************/

	/* Boltzmann factors for Planck integration */
	TDustRyg = TE1RYD/tdust;

	dintg_v = 0.;
	sum = 0.;

	for( i=0; i < rfield.nupper; i++ )
	{
		/* this is hnu/kT for grain at this temp and photon energy */
		arg = TDustRyg*rfield.anu[i];

		/* want the number exp(hnu/kT) - 1, two expansions */
		if( arg < 1e-6 )
		{
			/* for small arg expand exp(hnu/kT) - 1 = hnu/kT */
			ExpM1 = arg;
		}
		else
		{
			/* for large arg, evaluate the full expansion */
			ExpM1 = exp(MIN2(700.,arg)) - 1.;
		}

		Planck1 = 2.1675e16*rfield.anu3[i]/ExpM1*rfield.widflx[i];
		sum += Planck1;
		dintg_v += Planck1*GrainVar.dstab1[nd][i];
	}

	/* this is an option to print out every few steps, when 'trace grains' is set */
	if( (trace.lgDustBug && trace.lgTrace) && 
		((ip - 1)/5)*5 == (ip - 1) )
	{
		fprintf( ioQQQ, "  %4ld %11.4e %11.4e %11.4e\n", nd, tdust, 
		  dintg_v, sum/4./5.6696e-5/powi(tdust,4) );
	}


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