/*dtempr compute grains temperature, heating */
#include <float.h>
#include "cddefines.h"
#include "physconst.h"
#include "grainvar.h"
#include "trace.h"
#include "zonecnt.h"
#include "rfield.h"
#include "dhla.h"
#include "heating.h"
#include "dheton.h"
#include "phycon.h"
#include "dclgas.h"
#include "splint.h"
#include "dcolid.h"
#include "qheat1.h"
#include "taulines.h"
#include "dtempr.h"

void dtempr(void)
{
	long int i, 
	  nd;
	double 
	  dtherm, 
	  eyhat[NDUST][NCELL], 
	  facmax ,
	  hcon, 
	  hla, 
	  hots, 
	  phi, 
	  qtmax0, 
	  sum1, 
	  ted1, 
	  tqhout, 
	  x, 
	  y;
	float dccool, 
	  dcheat;

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

	ted1 = -DBL_MAX;
	tqhout = -DBL_MAX;

	dclgas.DustCoolGas = 0.;
	for( nd=0; nd < NDUST; nd++ )
	{
		if( GrainVar.lgDustOn1[nd] )
		{
			/* find heating
			 * add on fractional bit */
			hcon = 0.;
			hots = 0.;

			/* integrate over incident continuum for non-ionizing energies */
			for( i=0; i < (GrainVar.ipDustPot[nd] - 1); i++ )
			{
				/*  direct heating by incident continuum */
				hcon += rfield.flux[i]*rfield.anu[i]*GrainVar.dstab1[nd][i];
				hots += rfield.SummedDif[i]*rfield.anu[i]*GrainVar.dstab1[nd][i];
			}

			/* hcon is radiative heating by NON-ionizing radiation */
			hcon *= (float)(EN1RYD*phycon.hden*GrainVar.dstAbund[nd]);
			hots *= GrainVar.dstAbund[nd];

			/* integrate over ionizing continuum; energy goes to dust and gas
			 * dheat is what heats the gas */
			GrainVar.dheat[nd] = 0.;

			/*  dtherm is part of h.nu-IP which heats dust before ejection */
			dtherm = 0.;

			for( i=GrainVar.ipDustPot[nd]-1; i < rfield.nflux; i++ )
			{	
				facmax = POW2(GrainVar.DustPhotoE[nd][i])-POW2(GrainVar.dstpot[nd]);
				facmax = MAX2(0.,facmax);

				eyhat[nd][i] = 
					0.5*GrainVar.yn[nd][i]*
					MIN2(GrainVar.DustPhotoE[nd][i],facmax/GrainVar.DustPhotoE[nd][i]);

				/*  this is inciident continuum heating */
				phi = rfield.flux[i]*GrainVar.dstab1[nd][i]*GrainVar.dstAbund[nd];

				/*  dtherm is part of h.nu-IP which heats dust before ejection */
				dtherm += phi*(rfield.anu[i] - eyhat[nd][i]);

				/*  this is heating by all diffuse fields:
				 *  SummedDif has all continua and lines */
				sum1 = rfield.SummedDif[i]*GrainVar.dstab1[nd][i]*
				  GrainVar.dstAbund[nd];

				hots += sum1*(rfield.anu[i] - eyhat[nd][i]);

				/*  dheat is rate grain photoionization heats the gas */
				facmax = MAX2(0.,1.-GrainVar.dstpot[nd]/GrainVar.DustPhotoE[nd][i]);
				GrainVar.yhat[nd][i] = (float)(GrainVar.yn[nd][i]*MIN2(1.,facmax));

				GrainVar.dheat[nd] += 
					(float)(rfield.SummedCon[i]*GrainVar.dstab1[nd][i]*
				  GrainVar.dstAbund[nd]*(eyhat[nd][i] - GrainVar.dstpot[nd]*
				  GrainVar.yhat[nd][i]));
			}

			/* hots is total heating of the grain by diffuse fields */
			hots *= (float)(EN1RYD*phycon.hden);

			/* dtherm is part of h.nu-IP which heats dust before ejection */
			dtherm *= (float)(phycon.hden*EN1RYD);

			/* dheat is what heats the gas */
			if( dheton.lgDHetOn )
			{
				GrainVar.dheat[nd] *= (float)(phycon.hden*EN1RYD);
				HeatingCom.heating[13][0] += GrainVar.dheat[nd];
			}
			else
			{
				HeatingCom.heating[13][0] = 0.;
				GrainVar.dheat[nd] = 0.;
			}

			/*  find power absorbed by dust and resulting temperature
			 *  DustHeatLya - total heating by LA in this zone, erg cm-3 s-1, only here
			 *   for eventual printout, hots is total ots line heating
			 *  DTOTLA - integrated destroyed LA, erg s-1 */

			/*  heating by Ly A on dust in this zone, and total LA heating
			 *  only used for printout; Ly-a is already in OTS fields */
			hla = 
				rfield.otslin[HydroLines[0][IP2P][IP1S].ipCont-1]*
				GrainVar.dstab1[nd][HydroLines[0][IP2P][IP1S].ipCont-1]*
				0.75*EN1RYD*phycon.hden*GrainVar.dstAbund[nd];

			dhla.DustHeatLya += hla;

			/*  heating by thermal collisions with gas does work
			 *  factor of 1e8 is SQRT(8k/pi m_H)2k/4*1e20
			 *  DCHEAT is grain collisional heating by gas
			 *  DCCOOL is gas cooling due to collisions with grains
			 *  they should wqual each other */
			dcolid(&dcheat,&dccool,nd );

			/*  this will be total collisional heating, for printing in lines */
			dhla.dhcol += dcheat;

			/*  dust often hotter than gas during initial TE search */
			if( ZoneCnt.nzone <= 2 )
				dcheat = MAX2(0.f,dcheat);

			/* hcon is heating from incident continuum by non-ionizing radiation
			 * hots is heating from ots continua and lines
			 * dcheat is grain collisional heating by gas
			 * dtherm is part of h.nu-IP which heats dust before ejection,
			 * this is only direct heating by ionizing radiation
			 *
			 * dstsum is total heating of this grain type, to be balanced by cooling
			 * dstsum is the grain heat input - this is the actual number entering
			 * the energy balance */
			dhla.dstsum = hcon + hots + dcheat + dtherm;

			/* if( dstsum.gt.0. )write(66,'(i4,1p,10e10.2)') nzone,
			 * 1   hcon/dstsum , hots/dstsum , dcheat/dstsum , dtherm/dstsum
			 *
			 * following used to keep track of heating agents in printout
			 * no physics done with dincon
			 * dust heating by incident continuum, and elec friction before ejection */
			dhla.dincon += hcon + dtherm;
			dhla.HDustDif += hots;

			/* dstotl is total heating of all grain types in this zone,
			 * will be carried by total cooling, only used in lines to print tot heat
			 * printed as entry "GraT    0 " */
			dhla.dstotl += dhla.dstsum;

			/* remember total heating by diffuse fields, for printout */
			dhla.HDustDif += dcheat;

			/* >>chng 96 nov 30, add Kevin Volk's quantum heating
			 * (K Volk) added calculation of QTCOL here */
			if( GrainVar.lgQHeat[nd] )
			{
				if( dcheat > 0. )
				{
					GrainVar.qcu1 = (float)(dcheat*4.5911e10*4.);

					/* QCU1 is the total grain collisional heating rate
					 * in rydberg/s/H atom
					 * DCHEAT is in erg/s/H atom; factor 4 for the entire surface, factor
					 * 4.5911e+10 to convert to rydberg energy units
					 * change, 97 july 17, div by dstAbund due to variable abundances */
					x = log10(dcheat/phycon.hden/GrainVar.dstAbund[nd]);

					if( x < GrainVar.dstems[nd][0] )
					{
						GrainVar.qtcol = (float)pow(10.,GrainVar.dsttmp[0]);
						GrainVar.qcu1 = (float)(phycon.hden*(pow(10.,GrainVar.dstems[nd][0]))*
						  4.5911e10*4.*GrainVar.dstAbund[nd]);
						/* correct qtcol and qcu1 to the minimum values
						 * if the collisional heating rate is very low;
						 * otherwise the quantum heating routines fail since some of the
						 * temperature bins are out of range of the dstems array.
						 * This will produce erroneous results if both the
						 * collisional and quantum heating rates are very low. */
					}
					else if( x > GrainVar.dsttmp[NDEMS-1] )
					{
						GrainVar.qtcol = (float)(pow(10.,GrainVar.dsttmp[NDEMS-1]));
					}
					else
					{
						splint(&GrainVar.dstems[nd][0],GrainVar.dsttmp,
						  &GrainVar.dstslp[nd][0],NDEMS,x,&y);
						GrainVar.qtcol = (float)pow(10.,y);
					}
				}
				else
				{
					GrainVar.qtcol = 1.0;
				}
			}

			/*  End of changes here (K Volk) */

			/*  now find temperature, dstsum is sum of total heating of grain
			 *  >>chng 97 july 17, divide by abundance here */
			x = log10(MAX2(1e-37,dhla.dstsum/phycon.hden/GrainVar.dstAbund[nd]));

			if( x < GrainVar.dstems[nd][0] )
			{
				GrainVar.tedust[nd] = (float)pow(10.,GrainVar.dsttmp[0]);
				/* >>chng 96 apr 27, as per Peter van Hoof comment */
			}
			else if( x > GrainVar.dstems[nd][NDEMS-1] )
			{
				GrainVar.tedust[nd] = (float)pow(10.,GrainVar.dsttmp[NDEMS-1]);
			}
			else
			{
				splint(&GrainVar.dstems[nd][0],GrainVar.dsttmp,&GrainVar.dstslp[nd][0],
				  NDEMS,x,&y);
				GrainVar.tedust[nd] = (float)pow(10.,y);
			}

			/* >>chng 96 nov 30, add Kevin Volk quantum heating
			 * Save the equilibrium temperature for possible print-out when the
			 * quantum heating is asked for. (K Volk) */
			GrainVar.tequil1[nd] = GrainVar.tedust[nd];

			/*  if qheat is set do the quantum heating calculation.
			 *  test qtmax first, and adjust if needed. */
			if( GrainVar.lgQHeat[nd] )
			{
				ted1 = GrainVar.tedust[nd];

				/* use this value if something strange happens in qheat1 */
				if( GrainVar.qtmax < 5.*GrainVar.tedust[nd] )
				{
					qtmax0 = GrainVar.qtmax;
					GrainVar.qtmax = (float)MIN2(20.*GrainVar.tedust[nd],
					  5000.);
					qheat1(nd,&tqhout,ted1);
					GrainVar.qtmax = (float)qtmax0;
				}
				else
				{
					qheat1(nd,&tqhout,ted1);
				}
				GrainVar.tedust[nd] = (float)tqhout;
			}

			/*  end of changes here (K Volk)
			 *  save for latter possible printout */
			GrainVar.TeGrainMax[nd] = (float)MAX2(GrainVar.TeGrainMax[nd],
			  GrainVar.tedust[nd]);

			/*  rate dust cools gas by collisions */
			if( dheton.lgDColOn )
			{
				/*  dccool is gas cooling due to collisions with grains */
				dclgas.DustCoolGas += dccool;
			}

			if( trace.lgTrace && trace.lgDustBug )
			{
				fprintf( ioQQQ, " DTEMPR finds%2ld Tdst%10.2e totHeat%10.2e Cont;%10.2e OTS;%10.2e HIIcol:%10.2e pho ej;%10.2e dheat%10.2e\n", 
				  nd+1, GrainVar.tedust[nd], dhla.dstsum, hcon/dhla.dstsum, 
				  hots/dhla.dstsum, dcheat/dhla.dstsum, dtherm/dhla.dstsum, 
				  GrainVar.dheat[nd] );

				if( GrainVar.lgQHeat[nd] )
				{
					fprintf( ioQQQ, " DTEMPR: QHEAT1 found tqhout, ted1=%10.2e%10.2e\n", 
					  tqhout, ted1 );
				}
			}
		}
	}

	/* if grains hotter than gas then collisions with gas act
	 * to heat the gas, add this in here
	 * a symmetric statement appears in COOLR, where cooling is added on */
	if( dheton.lgDHetOn )
	{
		HeatingCom.heating[14][0] = MAX2(0.,-dclgas.DustCoolGas);
	}
	else
	{
		HeatingCom.heating[14][0] = 0.;
	}


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