/*grain main routine to converge grains thermal solution */
#include <string.h>
#include "cddefines.h"
#include "grainvar.h"
#include "heating.h"
#include "rfield.h"
#include "dhla.h"
#include "dclgas.h"
#include "grnvrydpth.h"
#include "dcharg.h"
#include "dtempr.h"
#include "ddrift.h"
#include "grain.h"

void grain(void)
{
	long int i, 
	  nd;

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

	/* this is the main routine that drives the grain physics */

	/* these will sum heating agents over grain populations */
	dhla.DustHeatLya = 0.;
	dhla.dhcol = 0.;
	dhla.dincon = 0.;

	/* dstotl will be total cooling by grains here */
	dhla.dstotl = 0.;

	/* total heating due to absorption of diffuse fields, not Lya */
	dhla.HDustDif = 0.;

	/* grains are not on, set all observables to zero */
	if( !GrainVar.lgDustOn )
	{
		dhla.dstsum = 0.;
		HeatingCom.heating[13][0] = 0.;
		dclgas.DustCoolGas = 0.;
		for( nd=0; nd < NDUST; nd++ )
		{
			GrainVar.dheat[nd] = 0.;
			GrainVar.tedust[nd] = 0.;
			GrainVar.DustDftVel[nd] = 0.;
		}
		
#		ifdef DEBUG_FUN
		fputs( " <->grain()\n", debug_fp );
#		endif
		return;
	}

	/* option to have grain abundance depend on depth */
	for( i=0; i < rfield.nupper; i++ )
	{
		/* these are total absorption and scattering cross sections */
		GrainVar.dstab[i] = 0.;
		GrainVar.dstsc[i] = 0.;
	}

	for( nd=0; nd < NDUST; nd++ )
	{
		/* refresh grain opacities if variable grain abundances enabled */
		if( GrainVar.lgDustOn1[nd] && GrainVar.lgDustVary[nd] )
		{
			GrainVar.dstAbund[nd] = (float)(GrainVar.dstfactor[nd]*GrnVryDpth(nd));
		}

		for( i=0; i < rfield.nupper; i++ )
		{
			/* these are total absorption and scattering cross sections */
			GrainVar.dstab[i] += GrainVar.dstab1[nd][i]*GrainVar.dstAbund[nd];
			GrainVar.dstsc[i] += GrainVar.dstsc1[nd][i]*GrainVar.dstAbund[nd];
		}
	}

	/* must redo grain solution all the time since they
	 * generally have non-trivial contributions to heating-cooling */
	HeatingCom.heating[13][0] = 0.;

	/* find dust charge */
	dcharg();

	/* find heating (gas and dust) and grain temperature */
	dtempr();

	/* find dust drift velocity */
	ddrift();

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