/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
/*tfidle update some temperature dependent variables */
/*tauff compute optical depth where cloud is thin to free-free and plasma freq */
/*velset set thermal velocities for all particles in gas */
#include "cddefines.h"
#include "physconst.h"
#include "thermal.h"
#include "opacity.h"
#include "iso.h"
#include "dense.h"
#include "phycon.h"
#include "gaunt.h"
#include "trace.h"
#include "rfield.h"
#include "doppvel.h"
#include "radius.h"
#include "wind.h"
#include "path.h"
#include "tfidle.h"

/*tauff compute optical depth where cloud is thin to free-free and plasma freq */
static void tauff(void);
/*velset set thermal velocities for all particles in gas */
static void velset(void);
/* On first run, fill GauntFF with gaunt factors	*/
static void FillGFF(void);
/* Interpolate on GauntFF to calc gaunt at current temp, phycon.te	*/
static float InterpolateGff( long charge , double ERyd );
static int LinterpTable(float **t, float *v, long int lta, long int ltb, float x, float *a, long int *pipx);
static int LinterpVector(float **t, float *v, long lta , long ltb, float *yy , long ny, float **a);
void fhunt(float *xx, long int n, float x, long int *j);

static long lgGffNotFilled = TRUE;

#define N_TE_GFF		21L
#define N_PHOTON_GFF	145L	/* log(photon energy) = -8 to 10 in one-eighth dec steps	*/
static float ***GauntFF;
static float **GauntFF_T;
/* the array of logs of temperatures at which GauntFF is defined */
static float TeGFF[N_TE_GFF];
/* the array of logs of u at which GauntFF is defined	*/
static float PhoGFF[N_PHOTON_GFF];

void tfidle(
	/* option to force update of all variables */
	int lgForceUpdate)
{
	static float tgffused=-1.f, 
	  tgffsued2=-1.f;
	static long int nff_defined=-1;
	static long maxion = 0, oldmaxion = 0;
	static float ttused = 0.f;
	static double edused = 0.;
	static int lgZLogSet = FALSE;
	int lgGauntF[2];
	long int ion;
	long int n, 
	  i,
	  nelem,
	  if1,
		ipTe,
		ret;
	double fac,  
	  fanu;

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

	/* called with lgForceUpdate true in zero.c, when we must update everything */
	if( lgForceUpdate )
	{
		ttused = -1.f;
		edused = 0.;
	}

	/* check that eden not negative */
	if( dense.eden <= 0. )
	{
		fprintf( ioQQQ, " I found a zero or negative electron density,%10.2e\n", 
		  dense.eden );
		puts( "[Stop in tfidle]" );
		ShowMe();
		puts( "[Stop in tfidle]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* check that temperature not negative */
	if( phycon.te <= 0. )
	{
		fprintf( ioQQQ, " I found a zero or negative electron temperature,%10.2e\n", 
		  phycon.te );
		puts( "[Stop in tfidle]" );
		ShowMe();
		puts( "[Stop in tfidle]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* only time only, set up array of logs of charge squared */
	if( !lgZLogSet )
	{
		for( nelem=0; nelem<LIMELM; ++nelem )
		{
			/* this array is used to modify the log temperature array
			 * defined below, for hydrogenic species of charge nelem+1 */
			phycon.sqlogz[nelem] = (float)log10( POW2(nelem+1.) );
		}
		lgZLogSet = TRUE;
	}

	/*lint -e777 test float equality */
	if( phycon.te != ttused )
	{
		ttused = phycon.te;
		thermal.te_update = phycon.te;
		/* current temperature in various units */
		phycon.te_eV = phycon.te/(float)EVDEGK;
		phycon.te_ryd = phycon.te/(float)TE1RYD;
		phycon.te_wn = phycon.te / (float)T1CM;

		phycon.tesqrd = POW2(phycon.te);
		phycon.sqrte = (float)sqrt(phycon.te);
		thermal.halfte = (float)(0.5/phycon.te);
		thermal.tsq1 = 1.f/phycon.tesqrd;
		phycon.te32 = phycon.te*phycon.sqrte;
		phycon.te70 = (float)pow(phycon.te,0.70);
		phycon.te30 = (float)pow(phycon.te,0.30);
		phycon.te20 = (float)pow(phycon.te,0.20);
		phycon.te90 =  phycon.te70 * phycon.te20;
		phycon.te10 = (float)pow(phycon.te,0.10);
		phycon.te07 = (float)pow(phycon.te,0.07);
		phycon.te05 = (float)pow(phycon.te,0.05);
		phycon.te03 = (float)pow(phycon.te,0.03);
		phycon.te02 = (float)pow(phycon.te,0.02);
		phycon.te04 = POW2( phycon.te02 );
		phycon.te01 = (float)pow(phycon.te,0.01);
		phycon.te005 = (float)pow(phycon.te,0.005);
		phycon.te003 = (float)pow(phycon.te,0.003);
		phycon.te001 = (float)pow(phycon.te,0.001);
		phycon.teinv = 1.f/phycon.te;
		phycon.alogte = (float)log10(phycon.te);
		phycon.telogn[0] = phycon.alogte;

		for( i=1; i < 7; i++ )
		{
			phycon.telogn[i] = phycon.telogn[i-1]*phycon.telogn[0];
		}
	}
	/*lint +e777 test float equality */

	/* term with hi added June 4, 93, to account for warm pdr */
	phycon.edensqte = (float)((dense.eden + dense.xIonDense[ipHYDROGEN][0]/1e4)/phycon.sqrte);
	phycon.cdsqte = (float)(phycon.edensqte*COLL_CONST);

	/* >>>chng 99 nov 23, removed this line, so back to old method of h coll */
	/* used for hydrogenic collisions */
	/* 
	 * following electron density has approximate correction for neutrals
	 * corr of hi*1.7e-4 accounts for col ion by HI; 
	 * >>refer	H0	correction for collisional contribution		Drawin Zs Phys 225, 483.
	 * used EdenHCorr instead
	 * edhi = eden + hi * 1.7e-4
	 */
	dense.EdenHCorr = dense.eden + dense.xIonDense[ipHYDROGEN][0]*1.7e-4;

	if( fabs(1.-edused/phycon.te)>=0.00001 )
	{
		edused = dense.eden;
		dense.SqrtEden = (float)sqrt(dense.eden);
	}

	/* finally reset velocities */
	/* find line widths for thermal and turbulent motion
	 * CLOUDY uses line center optical depths, =0.015 F / DELTA NU */
	velset();

	/* rest have to do with radiation field which may not be defined yet */
	if( !lgRfieldMalloced )
	{
#		ifdef DEBUG_FUN
		fputs( " <->tfidle()\n", debug_fp );
#		endif
		return;
	}

	/* 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 */
	/* >>chng 01 aug 21, must also test on size of continuum nflux because 
	 * conintitemp can increase nflux then call this routine, although 
	 * temp may not have changed */
	if( fabs(1.-tgffused/phycon.te)>=0.00001 /* tgffused != phycon.te */ 
		|| rfield.ContBoltz[0] <= 0. || nff_defined<rfield.nflux )
	{
		tgffused = phycon.te;
		fac = TE1RYD/phycon.te;
		i = 0;
		fanu = fac*rfield.anu[i];
		/* SEXP_LIMIT is 84 in cddefines.h - it is the -ln of smallest number
		 * that sexp can handle, and is used elsewhere in the code.  
		 * atom_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 < SEXP_LIMIT )
		{
			rfield.ContBoltz[i] = exp(-fanu);
			++i;
			/* this is boltz factor for NEXT cell */
			fanu = fac*rfield.anu[i];
		}
		/* ipMaxBolt is number of cells defined, so defined up through ipMaxBolt-1 */
		rfield.ipMaxBolt = i;

		/* zero out remainder */
		/* >>chng 01 apr 14, upper limit has been ipMaxBolt+1, off by one */
		for( i=rfield.ipMaxBolt; i < rfield.nupper; i++ )
		{
			rfield.ContBoltz[i] = 0.;
		}
	}

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

	oldmaxion = maxion;
	maxion = 0;

	/* Find highest maximum stage of ionization	*/
	for( nelem = 0; nelem < LIMELM; nelem++ )
	{
		maxion = MAX2( maxion, dense.IonHigh[nelem] );
	}

	/* reevaluate if temperature or number of cells has changed */
	/* >>chng 03 sep 26, following had been 0.1, cauing jumps when evaluated,
	 * changed to 0.02 */
#define LIM 0.02
	if( fabs(1.-tgffsued2/phycon.te) > LIM || 
		/* this test - reevaluate if upper bound of defined values is
		 * above nflux, the highest point.  This only triggers in
		 * large grids when continuum source gets harder */
		nff_defined<rfield.nflux 
		/* this occurs when now have more stages of ionization than in previous defn */
		|| maxion>oldmaxion )
	{
		static long lgFirstRunDone = FALSE;
		long lowion;
		/* >>chng 02 jan 10, only reevaluate needed states of ionization */
		if( fabs(1.-tgffsued2/phycon.te) <= LIM && nff_defined>=rfield.nflux && 
			maxion>oldmaxion )
		{
			/* this case temperture did not change by much, but upper
			 * stage of ionization increased.  only need to evaluate
			 * stages that have not been already done */
			lowion = oldmaxion ;
		}
		else
		{
			/* temperature changed - do them all */
			lowion = 1;
		}

		/* if1 will certainly be set to some positive value in gffsub */
		if1 = 1;

		/* chng 02 may 16, by Ryan...one gaunt factor array for all charges	*/
		/* First index is EFFECTIVE CHARGE!	*/
		/* highest stage of ionization is LIMELM, 
         * index goes from 1 to LIMELM */
		
		nff_defined = rfield.nflux;

		ASSERT( if1 >= 0 );

		tgffsued2 = phycon.te;
		lgGauntF[0] = TRUE;

		/* new gaunt factors	*/
		if( lgGffNotFilled )
		{
			FillGFF();
		}
		
		if( lgFirstRunDone == FALSE )
		{
			maxion = LIMELM;
			lgFirstRunDone = TRUE;
		}

		/* >> chng 03 jan 23, rjrw -- move a couple of loops down into
		 * subroutines, and make those subroutines generic
		 * (i.e. dependences only on arguments, may be useful elsewhere...) */

		/* Make Gaunt table for new temperature */
		ipTe = -1;
		for( ion=1; ion<=LIMELM; ion++ )
		{
			/* Interpolate for all tabulated photon energies at this temperature */
			ret = LinterpTable(GauntFF[ion], TeGFF, N_PHOTON_GFF, N_TE_GFF, phycon.alogte, GauntFF_T[ion], &ipTe);
			if (ret == -1) 
			{
				fprintf(ioQQQ," LinterpTable for GffTable called with te out of bounds \n");
				puts( "[Stop in Tfidle]" );
				cdEXIT(EXIT_FAILURE);
			}
		}

		/* Interpolate for all ions at required photon energies -- similar
		 * to LinterpTable, but not quite similar enough... */
		/* >>chng 04 jun 30, change nflux to nupper */
		ret = LinterpVector(GauntFF_T+lowion, PhoGFF, maxion-lowion+1, N_PHOTON_GFF,
			rfield.anulog, rfield.nupper,/*rfield.nflux,*/ rfield.gff + lowion); 
		if (ret == -1) 
		{
			fprintf(ioQQQ," LinterpVector for GffTable called with photon energy out of bounds \n");
			puts( "[Stop in Tfidle]" );
			cdEXIT(EXIT_FAILURE);
		}
	}
	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, "     tfidle; 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", 
		  rfield.gff[1][0], rfield.gff[1][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]-1],
		  if1, rfield.gff[1][iso.ipIsoLevNIonCon[ipH_LIKE][ipHYDROGEN][2]], 
		  rfield.gff[1][rfield.nflux-1] );
	}

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

/*tauff compute optical depth where cloud is thin to free-free and plasma freq */
static void tauff(void)
{
	double fac;

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

	/* simply return if space not yet allocated */
	if( !lgOpacMalloced )
		return;

	/* routine sets variable ipEnergyBremsThin, index for lowest cont cell that is optically thin */
	/* find frequency where continuum thin to free-free */
	while( rfield.ipEnergyBremsThin < rfield.nflux && 
		opac.TauAbsGeo[0][rfield.ipEnergyBremsThin] >= 1. )
	{
		++rfield.ipEnergyBremsThin;
	}

	/* TFF will be frequency where cloud becomes optically thin to brems
	 * >>chng 96 may 7, had been 2, change as per Kevin Volk bug report */
	if( rfield.ipEnergyBremsThin > 1 && opac.TauAbsGeo[0][rfield.ipEnergyBremsThin] > 0.001 )
	{
		/* tau can be zero when plasma frequency is within energy grid, */
		fac = (1. - opac.TauAbsGeo[0][rfield.ipEnergyBremsThin-1])/(opac.TauAbsGeo[0][rfield.ipEnergyBremsThin] - 
		  opac.TauAbsGeo[0][rfield.ipEnergyBremsThin-1]);
		fac = MAX2(fac,0.);
		rfield.EnergyBremsThin = (float)(rfield.anu[rfield.ipEnergyBremsThin-1] + rfield.widflx[rfield.ipEnergyBremsThin-1]*
		  fac);
	}
	else
	{
		rfield.EnergyBremsThin = 0.;
	}

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

	if( rfield.ipPlasma > 0 )
	{
		/* >>chng 02 jul 25, increase index for plasma frequency until within proper cell */
		while( rfield.plsfrq > rfield.anu[rfield.ipPlasma]+rfield.widflx[rfield.ipPlasma]/2. )
			++rfield.ipPlasma;

		/* >>chng 02 jul 25, decrease index for plasma frequency until within proper cell */
		while( rfield.ipPlasma>2 && rfield.plsfrq < rfield.anu[rfield.ipPlasma]-rfield.widflx[rfield.ipPlasma]/2. )
			--rfield.ipPlasma;
	}

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

	/* is plasma frequency within energy grid? */
	if( rfield.plsfrq > rfield.anu[0] )
	{
		rfield.lgPlasNu = TRUE;
	}

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

	/* now increment ipEnergyBremsThin still further, until above plasma frequency */
	while( rfield.ipEnergyBremsThin < rfield.nflux && 
		rfield.anu[rfield.ipEnergyBremsThin] <= rfield.EnergyBremsThin )
	{
		++rfield.ipEnergyBremsThin;
	}

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

/*velset set thermal velocities for all particles in gas */
static void velset(void)
{
	long int nelem;
	double turb2;

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

	/* usually TurbVel =0, reset with turbulence command
	 * cm/s here, but was entered in km/s with command */
	turb2 = DoppVel.TurbVel*DoppVel.TurbVel;

	/* this is option to dissipate the turbulence.  DispScale is entered with
	 * dissipate keyword on turbulence command.  The velocity is reduced here,
	 * by an assumed exponential scale, and also adds heat */
	if( DoppVel.DispScale > 0. )
	{
		/* square of exp depth dependence */
		turb2 *= sexp( 2.*radius.depth / DoppVel.DispScale );
	}

	/* in case of D-Critical flow include inital velocity as
	 * a component of turbulence */
	if( wind.windv0 < 0. )
	{
		turb2 += POW2(wind.windv0);
	}

	/* computes one doppler width, in cm/sec,
	 * for each element with atomic number the array index*/
	for( nelem=0; nelem < LIMELM; nelem++ )
	{
		DoppVel.doppler[nelem] = 
			/*(float)sqrt(1.651e8*phycon.te/dense.AtomicWeight[nelem]+*/
			/* >>chng 00 may 02 to physical constants */
			(float)sqrt(2.*BOLTZMANN/ATOMIC_MASS_UNIT*phycon.te/dense.AtomicWeight[nelem]+
		  turb2);
		/* this is average (NOT rms) particle speed for Maxwell distribution, Mihalas 70, 9-70 */
		DoppVel.AveVel[nelem] = sqrt(8.*BOLTZMANN/PI/ATOMIC_MASS_UNIT*phycon.te/dense.AtomicWeight[nelem]);
	}

	/* DoppVel.doppler[LIMELM] is CO, vector is dim LIMELM+1 */
	/* C12O16 */
	DoppVel.doppler[LIMELM] = 
		(float)sqrt(2.*BOLTZMANN/ATOMIC_MASS_UNIT*phycon.te/
		(dense.AtomicWeight[5]+dense.AtomicWeight[7]) + turb2);
	DoppVel.AveVel[LIMELM] = 
		(float)sqrt(8.*BOLTZMANN/ATOMIC_MASS_UNIT*phycon.te/
		(dense.AtomicWeight[5]+dense.AtomicWeight[7]) + turb2);

	/* C13O16 */
	DoppVel.doppler[LIMELM+1] = 
		(float)sqrt(2.*BOLTZMANN/ATOMIC_MASS_UNIT*phycon.te/
		(dense.AtomicWeight[5]*13./12.+dense.AtomicWeight[7]) + turb2);
	DoppVel.AveVel[LIMELM+1] = 
		(float)sqrt(8.*BOLTZMANN/ATOMIC_MASS_UNIT*phycon.te/
		(dense.AtomicWeight[5]*13./12.+dense.AtomicWeight[7]) + turb2);

	/* H2 */
	DoppVel.doppler[LIMELM+2] = 
		(float)sqrt(2.*BOLTZMANN/ATOMIC_MASS_UNIT*phycon.te/
		(2.*dense.AtomicWeight[0]) + turb2);
	DoppVel.AveVel[LIMELM+2] = 
		(float)sqrt(8.*BOLTZMANN/ATOMIC_MASS_UNIT*phycon.te/
		(2.*dense.AtomicWeight[0]) + turb2);

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

static void FillGFF( void )
{

	long i,i1,i2,i3,j,charge,GffMAGIC = 20613;
	double Temp, photon;
	int lgEOL;

#	define chLine_LENGTH 1000
	char chLine[chLine_LENGTH] , 
		/* this must be longer than chDataPath, set in path.h */
		chFilename[FILENAME_PATH_LENGTH_2];

	FILE *ioDATA;

	for( i = 0; i < N_TE_GFF; i++ )
	{
		TeGFF[i] = 0.5f*i;
	}

	for( i = 0; i< N_PHOTON_GFF; i++ )
	{
		PhoGFF[i] = 0.125f*i - 8.f;
	}

	if( (GauntFF = (float***)MALLOC((size_t)((LIMELM+1)*sizeof(float**)) ) ) == NULL )
		BadMalloc();
	for( i = 1; i <= LIMELM; i++ )
	{
		if( (GauntFF[i] = (float**)MALLOC((size_t)(N_PHOTON_GFF*sizeof(float*)) ) ) == NULL )
			BadMalloc();
		for( j = 0; j < N_PHOTON_GFF; j++ )
		{
			if( (GauntFF[i][j] = (float*)MALLOC((size_t)(N_TE_GFF*sizeof(float)) ) ) == NULL )
				BadMalloc();
		}
	}

	if( (GauntFF_T = (float**)MALLOC((size_t)((LIMELM+1)*sizeof(float*)) ) ) == NULL )
		BadMalloc();
	for( i = 1; i <= LIMELM; i++ )
	{
		if( (GauntFF_T[i] = (float*)MALLOC((size_t)(N_PHOTON_GFF*sizeof(float)) ) ) == NULL )
			BadMalloc();
	}
	
	if( !rfield.lgCompileGauntFF )
	{
		if( lgDataPathSet == TRUE )
		{
			/*path set, so look only there */
			strcpy( chFilename , chDataPath );
			strcat( chFilename , "gauntff.dat" );
		}
		else
		{
			/* path not set, check local space only */
			strcpy( chFilename , "gauntff.dat" );
		}

		if( trace.lgTrace )
			fprintf( ioQQQ," FillGFF opening gauntff.dat:");

		if( ( ioDATA = fopen( chFilename , "r" ) ) == NULL )
		{
			fprintf( ioQQQ, " FillGFF could not open gauntff.dat\n" );
			if( lgDataPathSet == TRUE )
				fprintf( ioQQQ, " even tried path\n" );

			if( lgDataPathSet == TRUE )
			{
				fprintf( ioQQQ, " FillGFF could not open gauntff.dat\n");
				fprintf( ioQQQ, " path is ==%s==\n",chDataPath );
				fprintf( ioQQQ, " final path is ==%s==\n",chFilename );
				fprintf( ioQQQ, " Defaulting to on-the-fly computation, ");
				fprintf( ioQQQ, "but the code runs much faster if you compile gauntff.dat!\n");
			}
			
			/* Do on the fly computation of Gff's instead.	*/
			for( charge=1; charge<=LIMELM; charge++ )
			{
				for( i=0; i<N_PHOTON_GFF; i++ )
				{
					photon = pow(10.,PhoGFF[i]);
					for (j=0; j<N_TE_GFF; j++)
					{
						Temp = pow(10.,TeGFF[j]);
						GauntFF[charge][i][j] = (float)CalcThermAveGaunt( Temp, (double)charge, photon );
					}
				}
			}
		}
		else 
		{
			/* check that magic number is ok */
			if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
			{
				fprintf( ioQQQ, " FillGFF could not read first line of gauntff.dat.\n");
				puts( "[Stop in FillGFF]" );
				cdEXIT(EXIT_FAILURE);
			}
			i = 1;
			i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
			i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
			i3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
			
			if( i1 !=GffMAGIC || i2 !=N_PHOTON_GFF || i3 !=N_TE_GFF )
			{
				fprintf( ioQQQ, 
					" FillGFF: the version of gauntff.dat is not the current version.\n" );
				fprintf( ioQQQ, 
					" FillGFF: I expected to find the numbers  %li %li %li and got %li %li %li instead.\n" ,
					GffMAGIC ,
					N_PHOTON_GFF,
					N_TE_GFF,
					i1 , i2 , i3 );
				fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
				fprintf( ioQQQ, 
					" FillGFF: please recompile the data file with the COMPILE GAUNT command.\n" );
				puts( "[Stop in FillGFF]" );
				cdEXIT(EXIT_FAILURE);
			}

			/* now read in the data */
			for( charge = 1; charge <= LIMELM; charge++ )
			{
				for( i = 0; i<N_PHOTON_GFF; i++ )
				{
					/* get next line image */
					if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
					{
						fprintf( ioQQQ, " FillGFF could not read first line of gauntff.dat.\n");
						puts( "[Stop in FillGFF]" );
						cdEXIT(EXIT_FAILURE);
					}
					/* each line starts with charge and energy level ( index in rfield ) */
					i3 = 1;
					i1 = (long)FFmtRead(chLine,&i3,INPUT_LINE_LENGTH,&lgEOL);
					i2 = (long)FFmtRead(chLine,&i3,INPUT_LINE_LENGTH,&lgEOL);
					/* check that these numbers are correct */
					if( i1!=charge || i2!=i )
					{
						fprintf( ioQQQ, " FillGFF detected insanity in gauntff.dat.\n");
						fprintf( ioQQQ, 
							" FillGFF: please recompile the data file with the COMPILE GAUNT command.\n" );
						puts( "[Stop in FillGFF]" );
						cdEXIT(EXIT_FAILURE);
					}

					/* loop over temperatures to produce array of free free gaunt factors	*/
					for (j = 0; j < N_TE_GFF; j++)
					{
						GauntFF[charge][i][j] = (float)FFmtRead(chLine,&i3,chLine_LENGTH,&lgEOL);

						if( lgEOL )
						{
							fprintf( ioQQQ, " FillGFF detected insanity in gauntff.dat.\n");
							fprintf( ioQQQ, 
								" FillGFF: please recompile the data file with the COMPILE GAUNT command.\n" );
							puts( "[Stop in FillGFF]" );
							cdEXIT(EXIT_FAILURE);
						}
					}
				}

			}

			/* check that magic number is ok */
			if( fgets( chLine , (int)sizeof(chLine) , ioDATA ) == NULL )
			{
				fprintf( ioQQQ, " FillGFF could not read first line of gauntff.dat.\n");
				puts( "[Stop in FillGFF]" );
				cdEXIT(EXIT_FAILURE);
			}
			i = 1;
			i1 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
			i2 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
			i3 = (long)FFmtRead(chLine,&i,INPUT_LINE_LENGTH,&lgEOL);
			
			if( i1 !=GffMAGIC || i2 !=N_PHOTON_GFF || i3 !=N_TE_GFF )
			{
				fprintf( ioQQQ, 
					" FillGFF: the version of gauntff.dat is not the current version.\n" );
				fprintf( ioQQQ, 
					" FillGFF: I expected to find the numbers  %li %li %li and got %li %li %li instead.\n" ,
					GffMAGIC ,
					N_PHOTON_GFF,
					N_TE_GFF,
					i1 , i2 , i3 );
				fprintf( ioQQQ, "Here is the line image:\n==%s==\n", chLine );
				fprintf( ioQQQ, 
					" FillGFF: please recompile the data file with the COMPILE GAUNT command.\n" );
				puts( "[Stop in FillGFF]" );
				cdEXIT(EXIT_FAILURE);
			}

			/* close the data file */
			fclose( ioDATA );
		}
		if( trace.lgTrace )
			fprintf( ioQQQ,"  - opened and read ok.\n");
	}
	else
	{
		/* option to create table of gaunt factors,
		 * executed with the compile gaunt command */
		FILE *ioGFF;
		
		/*GffMAGIC the magic number for the table of recombination coefficients */
		if( ( ioGFF = fopen( "gauntff.dat" , "w" ) ) == NULL )
		{
			fprintf( ioQQQ, " FillGFF could not open gauntff.dat for writing.\n" );
			puts( "[Stop in FillGFF]" );
			cdEXIT(EXIT_FAILURE);
		}
		fprintf(ioGFF,"%li\t%li\t%li\tGFF isoelectronic sequence recomb data, created by COMPILE GAUNT command, with %li energy levels and %li temperatures.\n",
			GffMAGIC ,
			N_PHOTON_GFF,
			N_TE_GFF,
			N_PHOTON_GFF,
			N_TE_GFF );

		for( charge = 1; charge <= LIMELM; charge++ )
		{
			for( i=0; i < N_PHOTON_GFF; i++ )
			{
				fprintf(ioGFF, "%li\t%li", charge, i );
				/* loop over temperatures to produce array of gaunt factors	*/
				for (j = 0; j < N_TE_GFF; j++)
				{
					/* Store gaunt factor in N_TE_GFF half dec steps */
					Temp = pow(10.,TeGFF[j]);
					photon = pow(10.,PhoGFF[i]);
					GauntFF[charge][i][j] = (float)CalcThermAveGaunt( Temp, (double)charge, photon );
					fprintf(ioGFF, "\t%f", GauntFF[charge][i][j] );
				}
				fprintf(ioGFF, "\n" );
			}
		}
			
		/* end the file with the same information */
		fprintf(ioGFF,"%li\t%li\t%li\tGFF isoelectronic sequence recomb data, created by COMPILE GAUNT command, with %li energy levels and %li temperatures.\n",
			GffMAGIC ,
			N_PHOTON_GFF,
			N_TE_GFF,
			N_PHOTON_GFF,
			N_TE_GFF );

		fclose( ioGFF );

		fprintf( ioQQQ, "FillGFF: compilation complete, gauntff.dat created.\n" );
	}

	lgGffNotFilled = FALSE;

	{
		/* We have already checked the validity of CalcThermAveGaunt in sanitycheck.c.
		 * Now we check to see if the InterpolateGff routine also works correctly.	*/
		/*@-redef@*/
		enum {DEBUG_LOC=FALSE};
		/* if set true there will be two problems at very low temps */
		/*@+redef@*/
		if( DEBUG_LOC )
		{
			double gaunt, error;
			float tempsave = phycon.te;
			long logu, loggamma2;

			for( logu=-4; logu<=4; logu++)
			{
				/* Uncommenting each of the three print statements in this bit
				 * will produce a nice table comparable to Table 2 of
				 * >>refer	free-free	gaunts	Sutherland, R.S., 1998, MNRAS, 300, 321-330 */
				/* fprintf(ioQQQ,"%li\t", logu);*/
				for(loggamma2=-4; loggamma2<=4; loggamma2++)
				{ 
					double SutherlandGff[9][9]=
					{	{5.5243, 5.5213, 5.4983, 5.3780, 5.0090, 4.4354, 3.8317, 3.2472, 2.7008},
						{4.2581, 4.2577, 4.2403, 4.1307, 3.7816, 3.2436, 2.7008, 2.2126, 1.8041},
						{3.0048, 3.0125, 3.0152, 2.9434, 2.6560, 2.2131, 1.8071, 1.4933, 1.2771},
						{1.8153, 1.8367, 1.8880, 1.9243, 1.7825, 1.5088, 1.2886, 1.1507, 1.0747},
						{0.8531, 0.8815, 0.9698, 1.1699, 1.2939, 1.1988, 1.1033, 1.0501, 1.0237},
						{0.3101, 0.3283, 0.3900, 0.5894, 0.9725, 1.1284, 1.0825, 1.0419, 1.0202},
						{0.1007, 0.1080, 0.1335, 0.2281, 0.5171, 0.9561, 1.1065, 1.0693, 1.0355},
						{0.0320, 0.0344, 0.0432, 0.0772, 0.1997, 0.5146, 0.9548, 1.1042, 1.0680},
						{0.0101, 0.0109, 0.0138, 0.0249, 0.0675, 0.1987, 0.5146, 0.9547, 1.1040}};

					phycon.te = (float)(TE1RYD/pow(10.,loggamma2));
					phycon.alogte = (float)log10((double)phycon.te);
					gaunt = InterpolateGff( 1, pow(10.,(double)(logu-loggamma2)) );
					error = fabs( gaunt - SutherlandGff[logu+4][loggamma2+4] ) /gaunt;
					/*fprintf(ioQQQ,"%1.3f\t", gaunt);*/
					if( error>0.05 )
					{
						fprintf(ioQQQ," Tfidle found insane gff. log(u) %li, log(gamma2) %li, error %.3e\n",
							logu, loggamma2, error );
					}
				}
				/*fprintf(ioQQQ,"\n");*/
			}
			phycon.te = tempsave;
			phycon.alogte = (float)log10(phycon.te);
		}
	}
	
	return;
}

/* Interpolate Gff at some temperature */
static float InterpolateGff( long charge , double ERyd )
{
	double GauntAtLowerPho, GauntAtUpperPho;
	static long int ipTe=-1, ipPho=-1;
	double gaunt = 0.;
	double xmin , xmax;
	long i;

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

	if( ipTe < 0 )
	{
		/* te totally unknown */
		if( phycon.alogte < TeGFF[0] || phycon.alogte > TeGFF[N_TE_GFF-1] )
		{
			fprintf(ioQQQ," InterpolateGff called with te out of bounds \n");
			puts( "[Stop in InterpolateGff]" );
			cdEXIT(EXIT_FAILURE);
		}
		/* now search for temperature */
		for( i=0; i<N_TE_GFF-1; ++i )
		{
			if( phycon.alogte > TeGFF[i] && phycon.alogte <= TeGFF[i+1] )
			{
				/* found the temperature - use it */
				ipTe = i;
				break;
			}
		}
		ASSERT( (ipTe >=0) && (ipTe < N_TE_GFF-1)  );
			
	}
	else if( phycon.alogte < TeGFF[ipTe] )
	{
		/* temp is too low, must also lower ipTe */
		ASSERT( phycon.alogte > TeGFF[0] );
		/* decrement ipTe until it is correct */
		while( phycon.alogte < TeGFF[ipTe] && ipTe > 0)
			--ipTe;
	}
	else if( phycon.alogte > TeGFF[ipTe+1] )
	{
		/* temp is too high */
		ASSERT( phycon.alogte <= TeGFF[N_TE_GFF-1] );
		/* increment ipTe until it is correct */
		while( phycon.alogte > TeGFF[ipTe+1] && ipTe < N_TE_GFF-1)
			++ipTe;
	}

	ASSERT( (ipTe >=0) && (ipTe < N_TE_GFF-1)  );

	/* ipTe should now be valid */
	ASSERT( phycon.alogte >= TeGFF[ipTe] && phycon.alogte <= TeGFF[ipTe+1] && ipTe < N_TE_GFF-1 );

	/***************/
	/* This bit is completely analogous to the above, but for the photon vector instead of temp.	*/
	if( ipPho < 0 )
	{
		if( log10(ERyd) < PhoGFF[0] || log10(ERyd) > PhoGFF[N_PHOTON_GFF-1] )
		{
			fprintf(ioQQQ," InterpolateGff called with photon energy out of bounds \n");
			puts( "[Stop in InterpolateGFF]" );
			cdEXIT(EXIT_FAILURE);
		}
		for( i=0; i<N_PHOTON_GFF-1; ++i )
		{
			if( log10(ERyd) > PhoGFF[i] && log10(ERyd) <= PhoGFF[i+1] )
			{
				ipPho = i;
				break;
			}
		}
		ASSERT( (ipPho >=0) && (ipPho < N_PHOTON_GFF-1)  );
			
	}
	else if( log10(ERyd) < PhoGFF[ipPho] )
	{
		ASSERT( log10(ERyd) >= PhoGFF[0] );
		while( log10(ERyd) < PhoGFF[ipPho] && ipPho > 0)
			--ipPho;
	}
	else if( log10(ERyd) > PhoGFF[ipPho+1] )
	{
		ASSERT( log10(ERyd) <= PhoGFF[N_PHOTON_GFF-1] );
		while( log10(ERyd) > PhoGFF[ipPho+1] && ipPho < N_PHOTON_GFF-1)
			++ipPho;
	}
	ASSERT( (ipPho >=0) && (ipPho < N_PHOTON_GFF-1)  );
	ASSERT( log10(ERyd)>=PhoGFF[ipPho] 
		&& log10(ERyd)<=PhoGFF[ipPho+1] && ipPho<N_PHOTON_GFF-1 );

	/* Calculate the answer...must interpolate on two variables.
	 * First interpolate on T, at both the lower and upper photon energies.
	 * Then interpolate between these results for the right photon energy.	*/

	GauntAtLowerPho = (phycon.alogte - TeGFF[ipTe]) / (TeGFF[ipTe+1] - TeGFF[ipTe]) *
		(GauntFF[charge][ipPho][ipTe+1] - GauntFF[charge][ipPho][ipTe]) + GauntFF[charge][ipPho][ipTe];

	GauntAtUpperPho = (phycon.alogte - TeGFF[ipTe]) / (TeGFF[ipTe+1] - TeGFF[ipTe]) *
		(GauntFF[charge][ipPho+1][ipTe+1] - GauntFF[charge][ipPho+1][ipTe]) + GauntFF[charge][ipPho+1][ipTe];

	gaunt = (log10(ERyd) - PhoGFF[ipPho]) / (PhoGFF[ipPho+1] - PhoGFF[ipPho]) * 
		(GauntAtUpperPho - GauntAtLowerPho) + GauntAtLowerPho;

	xmax = MAX4( GauntFF[charge][ipPho][ipTe+1], GauntFF[charge][ipPho+1][ipTe+1],
		GauntFF[charge][ipPho][ipTe], GauntFF[charge][ipPho+1][ipTe] );
	ASSERT( gaunt <= xmax );

	xmin = MIN4( GauntFF[charge][ipPho][ipTe+1], GauntFF[charge][ipPho+1][ipTe+1],
		GauntFF[charge][ipPho][ipTe], GauntFF[charge][ipPho+1][ipTe] );
	ASSERT( gaunt >= xmin );

	ASSERT( gaunt > 0. );

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

	return (float)gaunt;
}

/* Interpolate in table t[lta][ltb], with physical values for the
	 second index given by v[ltb], for values x, and put results in
	 a[lta]; store the index found if that's useful; assumes v[] is
	 sorted */
static int LinterpTable(float **t, float *v, long int lta, long int ltb, float x, float *a, long int *pipx)
{
	long int ipx=-1;
	float frac;
	long i;

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

	if (pipx != NULL)
		ipx = *pipx;
	
	fhunt (v,ltb,x,&ipx); 		/* search for index */
	if (pipx != NULL)
		*pipx = ipx;
	
	if( ipx == -1 || ipx == ltb )
	{
		return -1;
	}
	
	ASSERT( (ipx >=0) && (ipx < ltb-1)  );
	ASSERT( x >= v[ipx] && x <= v[ipx+1]);
	
	frac = (x - v[ipx]) / (v[ipx+1] - v[ipx]);
	for( i=0; i<lta; i++ )
	{
		a[i] = frac*t[i][ipx+1]+(1.f-frac)*t[i][ipx];
	}
	
#	ifdef DEBUG_FUN
	fputs( " <->LinterpTable()\n", debug_fp );
#	endif
	
	return 0 ;
}

/* Interpolate in table t[lta][ltb], with physical values for the second index given by v[ltb],
	 for values yy[ny], and put results in a[lta][ly] */
static int LinterpVector(float **t, float *v, long lta , long ltb, float *yy, long ly, float **a)
{
	float yl, frac;
	long i, j, n;

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

	if( yy[0] < v[0] || yy[ly-1] > v[ltb-1] )
	{
		return -1;
	}

	n = 0;
	yl = yy[n];
	for (j = 1; j < ltb && n < ly; j++) {
		while (yl < v[j] && n < ly) {
			frac = (yl-v[j-1])/(v[j]-v[j-1]);
			for (i = 0; i < lta; i++)
				a[i][n] = frac*t[i][j]+(1.f-frac)*t[i][j-1];
			n++;
			if (n == ly)
				break;
			assert (yy[n] > yy[n-1]);
			yl = yy[n];
		}
	}

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

	return 0;
}
void fhunt(float *xx, long int n, float x, long int *j)
{
	/*lint -e731 boolean argument to equal / not equal */
	long int jl, jm, jh, in;
	int up;

	jl = *j;
	up = (xx[n-1] > xx[0]);
	if (jl < 0 || jl >= n) 
	{
		jl = -1;
		jh = n;
	} 
	else 
	{
		in = 1;
		if ((x >= xx[jl]) == up) 
		{
			if (jl == n-1) 
			{
				*j = jl;
				return;
			}
			jh = jl + 1;
			while ((x >= xx[jh]) == up)
			{
				jl = jh;
				in += in;
				jh += in;
				if (jh >= n)
				{
					jh = n;
					break;
				}
			}
		}
		else
		{
			if (jl == 0)
			{
				*j = -1;
				return;
			}
			jh = jl--;
			while ((x < xx[jl]) == up)
			{
				jh = jl;
				in += in;
				jl -= in;
				if (jl <= 0)
				{
					jl = 0;
					break;
				}
			}
		}
	}
	while (jh-jl != 1)
	{
		jm = (jh+jl)/2;
		if ((x > xx[jm]) == up)
			jl = jm;
		else
			jh = jm;
	}
	*j = jl;
	return;
}
	/*lint +e731 boolean argument to equal / not equal */
