/*dcolid compute grains collisional heating cooling */
#include "cddefines.h"
#include "physconst.h"
#include "grainvar.h"
#include "tepowers.h"
#include "phycon.h"
#include "doppvel.h"
#include "sexp.h"
#include "ionfracs.h"
#include "dcolid.h"

void dcolid(float *dcheat, 
  float *dccool, 
  /* pointer to grain species */
  long int nd)
{
	long int i, 
	  ion, 
	  nelem;

	double cool, 
	  eta, 
	  heat, 
	  psi, 
	  vnz[4], 
	  xi, 
	  z;

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


	/* this subroutine evaluates the gas heating-cooling rate
	 * (erg cm^-3 s^-1) due to grain gas collisions.
	 * the net effect can be positive or negative,
	 * depending on whether the grains or gas are hotter
	 *
	 * this array is set of charges striking the grain surface
	 * vnz(1) are neutrals, vnz(2) are single ions, etc
	 * each population is weighted by the PROJECTED velocity
	 * */
	for( ion=0; ion < 4; ion++ )
	{
		vnz[ion] = 0.;
		for( nelem=0; nelem < LIMELM; nelem++ )
		{
			vnz[ion] += xIonFracs[ion+1][nelem]*DoppVel.doppler[nelem];
		}
	}

	/* following loop sums up collisions between these nuclei
	 * electrons are done later
	 * HEAT will be rate collisions heat the grain
	 * COOL will be rate collisions cool the gas */
	heat = 0.;
	cool = 0.;
	for( i=0; i < 4; i++ )
	{
		z = i;
		psi = z*GrainVar.dstpot[nd]*EVRYD/(8.617e-5*phycon.te);
		if( psi <= 0. )
		{
			eta = 1. - psi;
			xi = 1. - psi/2.;
		}
		else
		{
			eta = sexp(psi);
			xi = (1. + psi/2.)*eta;
		}

		/*  factor of two makes VNZ(I) arrival rate, not projected velocity */
		heat += vnz[i]*2.*BOLTZMANN*(2.*phycon.te*xi - eta*(z*GrainVar.dstpot[nd]*
		  EVRYD*1.6022e-12 - z*z*EN1RYD + 2.*GrainVar.tedust[nd]));

		/*  rate kinetic energy lost from gas - gas cooling */
		cool += vnz[i]*2.*BOLTZMANN*(2.*phycon.te*xi - eta*2.*GrainVar.tedust[nd]);
	}

	/* heating due to electrons */
	z = -1.;
	psi = z*GrainVar.dstpot[nd]*EVRYD/(8.617e-5*phycon.te);
	if( psi <= 0. )
	{
		eta = 1. - psi;
		xi = 1. - psi/2.;
	}
	else
	{
		eta = sexp(psi);
		xi = (1. + psi/2.)*eta;
	}

	/* electron arrival rate */
	vnz[0] = TePowers.sqrte*6.2124e5*phycon.eden;
	heat += vnz[0]*2.*1.38062e-16*phycon.te*xi;

	/* heat = heat + vnz(1)* ( 2.*1.38062e-16*te*xi -
	 *  1eta*z *dstpot(nd)*evryd*1.6022e-12 )
	 * */
	cool += vnz[0]*2.*1.38062e-16*phycon.te*xi;

	*dcheat = (float)(heat*(phycon.hden*GrainVar.darea[nd]*
		GrainVar.dstAbund[nd]/4.));
	/* cooling is per projected area */
	*dccool = (float)(cool*(phycon.hden*GrainVar.darea[nd]*
		GrainVar.dstAbund[nd]/ 4.));


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