/*boltgn evaluate Boltzmann factors for the continuum, and related variables */
#include "cddefines.h"
#include "physconst.h"
#include "nh.h"
#include "gaunts.h"
#include "trace.h"
#include "rfield.h"
#include "phycon.h"
#include "maxcor.h"
#include "tgff.h"
#include "tauff.h"
#include "gffsub.h"
#include "boltgn.h"

void boltgn(void)
{
	int lgGauntF[2];
	long int n, 
	  i, 
	  if1;
	double fac, 
	  fanu;

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

	/* correction factors for induced recombination, 
	 * also used as Boltzmann factors
	 * check for zero is because ContBoltz is zeroed out in initizationation
	 * of code, its possible this is a constant density grid of models
	 * in which the code is called as a subroutine */
	if( tgff.tgffused != phycon.te || rfield.ContBoltz[0] <= 0. )
	{
		tgff.tgffused = phycon.te;
		fac = TE1RYD/phycon.te;
		i = 0;
		fanu = fac*rfield.anu[i];
		/* NB - the 84 in the following must be kept parallel with the 84 in sexp,
		 * since level2 uses ContBoltz to see whether the rates will be significant.
		 * If the numbers did not agree then this test would be flawed, resulting in
		 * div by zero */
		while( i < rfield.nupper && fanu < 84. )
		{
			rfield.ContBoltz[i] = (float)exp(-fanu);
			++i;
			/* this is boltz factor for NEXT cell */
			fanu = fac*rfield.anu[i];
		}
		maxcor.ipMaxBolt = i;

		/* zero out remainder */
		for( i=maxcor.ipMaxBolt+1; i < rfield.nupper; i++ )
		{
			rfield.ContBoltz[i] = 0.;
		}
	}

	/* find frequency where thin to brems or plasma frequency */
	tauff();

	if( fabs(1.-tgff.tgffsued2/phycon.te) > 0.10 )
	{
		gffsub(0.,rfield.anulog,gaunts.gff,1,rfield.nflux,&if1);

		/* follwoing for ionized helium */
		gffsub(0.30103,rfield.anulog,gaunts.gffhe2,1,rfield.nflux,
		  &if1);

		tgff.tgffsued2 = phycon.te;
		lgGauntF[0] = TRUE;
	}
	else
	{
		/* this is flag that would have been set in gffsub, and
		 * printed in debug statement below.  We are not evaluating
		 * so set to -1 */
		if1 = -1;
		lgGauntF[0] = FALSE;
	}

	if( trace.lgTrace && trace.lgTrGant )
	{
		fprintf( ioQQQ, "     BOLTGN; gaunt facs?" );
		for(n=0; n < 2; n++)
			fprintf( ioQQQ, "%2c", TorF(lgGauntF[n]) );

		fprintf( ioQQQ, "%2f g 1 2=%10.2e%9.2ld flag%2f guv(1,n)%10.2e\n", 
		  gaunts.gff[0], gaunts.gff[nh.ipHn[0][2]-1], if1, gaunts.gff[nh.ipHn[0][2]], 
		  gaunts.gff[rfield.nflux-1] );
	}

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