/*tauff compute optical depth where cloud is thin to free-free and plasma freq */
#include "cddefines.h"
#include "opac.h"
#include "plasnu.h"
#include "rfield.h"
#include "phycon.h"
#include "tff.h"
#include "tauff.h"

void tauff(void)
{
	double fac;

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

	/* routine sets variable ntff, pointer to lowest cont cell that is thin */

	/* find frequency where continuum thin to free-free */
	while( tffCom.ntff < rfield.nflux && 
		opac.tauabs[0][tffCom.ntff-1] >= 1. )
	{
		tffCom.ntff += 1;
	}

	/* TFF will be frequency where goes thin to brems
	 * >>chng 96 May 7, had been 2, change as per Kevin Volk bug report */
	if( tffCom.ntff > 2 && opac.tauabs[0][tffCom.ntff-1] > 0.001 )
	{
		/* for case where plasma frequency important, tau can be zero */
		fac = (1. - opac.tauabs[0][tffCom.ntff-2])/(opac.tauabs[0][tffCom.ntff-1] - 
		  opac.tauabs[0][tffCom.ntff-2]);
		fac = MAX2(fac,0.);
		tffCom.tff = (float)(rfield.anu[tffCom.ntff-2] + rfield.widflx[tffCom.ntff-2]*
		  fac);
	}
	else
	{
		tffCom.tff = 0.;
	}

	/* now evaluate the plasma frequency */
	plasnu.plsfrq = (float)(2.729e-12*sqrt(phycon.eden*1.2));

	/* also remember the largest plasma frequency we encounter */
	plasnu.plsfrqmax = MAX2(plasnu.plsfrqmax, plasnu.plsfrq);

	if( plasnu.plsfrq > rfield.anu[0] )
	{
		plasnu.lgPlasNu = TRUE;
	}

	/* >>chng 96 july 15, did not include plasma freq before
	 * function returns larger of these two frequencies */
	tffCom.tff = (float)MAX2(plasnu.plsfrq,tffCom.tff);

	/* now increment ntff still further, until above plasma frequency */
	while( tffCom.ntff < rfield.nflux && 
		rfield.anu[tffCom.ntff-1] <= tffCom.tff )
	{
		tffCom.ntff += 1;
	}

	/* save max plasma frequency for case where density changes */
	if( tffCom.tff == plasnu.plsfrq )
	{
		plasnu.ipPlasmax = MAX2(plasnu.ipPlasma,tffCom.ntff);
	}


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