/*freeht evaluate free-free heating due to incident continuum */
#include "cddefines.h"
#include "physconst.h"
#include "hydrogenic.h"
#include "opacpoint.h"
#include "freeon.h"
#include "gaunts.h"
#include "opac.h"
#include "phycon.h"
#include "tepowers.h"
#include "he.h"
#include "rfield.h"
#include "h.h"
#include "heat.h"
#include "ffsum.h"
#include "zonecnt.h"
#include "freeht.h"

void freeht(void)
{
	long int i;

	double damper, 
	  fac, 
	  fac1, 
	  fachmi, 
	  fachmi1;
	double dfactor;
	static double oldfac = 0.;

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

	if( !freeon.lgFreeOn )
	{
		hydro.FreeFreeHeat = 0.;
#		ifdef DEBUG_FUN
		fputs( " <->freeht()\n", debug_fp );
#		endif
		return;
	}

	/* must constantly reevaluate if free free heating is important */
	if( opac.lgRedoStatic || hydro.FreeFreeHeat/heat.htot > 0.05 )
	{
		if( ZoneCnt.nzone < 1 )
		{
			damper = 1.;
		}
		else
		{
			damper = 0.5;
		}

		/* OPSV is missing factor of 1E-20 to avoid underflow */
		fac = (phycon.eden/1e10)*((h.hii + he.heii + 4.*he.heiii + 
		  ffsum.EdenFFSum)/1e10)/TePowers.sqrte;
		fac1 = fac*TE1RYD/phycon.te;

		/* >>chng 96 oct 30, put following after fac1 defn */
		fac = damper*fac + oldfac*(1. - damper);
		fachmi = (phycon.eden/1e10)*(h.hii*hydro.hn[0][IP1S]/1e10)/
		  TePowers.sqrte;
		fachmi1 = fachmi*TE1RYD/phycon.te;

		for( i=0; i < rfield.nflux; i++ )
		{
			if( rfield.ContBoltz[i] < 0.995 )
			{
				/* use 1-exp(hn/kT)) only if exp is small
				 * last term is scaled h minus brems free-free absorption */
				dfactor = opac.OpacStack[i-1+OpacPoint.ipBrems]*(1. - 
				  rfield.ContBoltz[i])*gaunts.gff[i];
				opac.FreeFreeOpacity[i] = dfactor*fac;
				/* this is the H- part, opac here is ratio of H- to ff */
				opac.FreeFreeOpacity[i] += dfactor*fachmi*opac.OpacStack[i-1+OpacPoint.iphmra];
			}
			else
			{
				dfactor = opac.OpacStack[i-1+OpacPoint.ipBrems]*rfield.anu[i]*
				  gaunts.gff[i];
				opac.FreeFreeOpacity[i] = dfactor*fac1;
				/* h- free free absorption */
				opac.FreeFreeOpacity[i] += dfactor*fachmi1*opac.OpacStack[i-1+OpacPoint.iphmra];
			}
		}
		oldfac = fac;

		/* now integrate over this */
		hydro.FreeFreeHeat = 0.;
		for( i=0; i < rfield.nflux; i++ )
		{
			hydro.FreeFreeHeat += opac.FreeFreeOpacity[i]*rfield.flux[i]*
			  rfield.anu[i];
		}
		hydro.FreeFreeHeat *= EN1RYD;
	}

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