/* 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 */
/*AtlasCompile rebin Kurucz stellar models to match energy grid of code */
/*AtlasGetOne get one of the Atlas model atmospheres, coded by K. Volk */
/*AtlasInterpolate interpolate on atlas model atmospheres, by K Volk */
#include "cddefines.h"
#include "physconst.h"
#include "rfield.h"
#include "path.h"
#include "called.h"
#include "stars.h"
/*lint -e613 possible use of null pointer */
/* NVALS is the number of values per Kurucz model. NATLAS is the total 
 * number of models in the set.  The last two are the Solar and Vega models.*/
#define	NATLAS	412
#define	NVALS	1221

static size_t nOffset;

/* the version number for this way to write out the atmospheres */
static long int VERSION=011017 ;

/*AtlasGetOne get one of the Atlas model atmospheres, 
 * original version by K. Volk */
static void AtlasGetOne(double temp, 
  double alogg, 
  long int jval[][11], 
  float tval[], 
  long int *ierr, 
  float teff[], 
  float xlogg[],
  FILE *ioIN ,
  float *fluxarray ,
  size_t nBlocksize );

void AtlasCompile(void)
{
	char chLine[INPUT_LINE_LENGTH], chC7[8], chC12[13];

	long int i, 
	  imod, 
	  j;

	/* these will be malloced into large work arrays*/
	float *StarEner , *StarFlux , *CloudyFlux , *scratch;
	/* these contain frequencies for the major absorption edges */
	float Edges[3];
	long nmalloc;

	/* we will write the binary results to this file*/
	FILE *ioOUT ,
		/* will get the ascii from this file */
		* ioIN ,
		/* this is short list of models in main file */
		* ioLIST;
	size_t nBlocksize = (size_t)rfield.nupper*sizeof(float );

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

	/* This is a program to re-bin the Kurucz stellar models spectrum to match the 
	 * CLOUDY grid.  For wavelengths shorter than supplied in the Kurucz files,
	 * the flux will be set to zero.  At long wavelengths a Rayleigh-Jeans
	 * extrapolation will be used. */

	/* This version uses power-law interpolation between the points of the stellar
	 * model.*/

	fprintf( ioQQQ, " AtlasCompile on the job.\n" );

	/* make sure continuum arrays have been allocated */
	ASSERT( rfield.nupper > 0 );

	/* define the major absorption edges that require special attention during rebinning
	 *
	 * NB the frequencies should be chosen here such that they are somewhere inbetween
	 * the two frequency points that straddle the edge in the atmosphere model, the
	 * software in RebinAtmosphere will seek out the exact values of those two points
	 * e.g.: in the CoStar models the H I edge is straddled by wavelength points at
	 * 911.67 and 911.85 A, so Edges[0] should be chosen somewhere inbetween (e.g. at 911.76A).
	 *
	 * NB beware not to choose edges too close to one another (i.e. on the order of the
	 * resolution of the Cloudy frequency grid). E.g. the He II Balmer edge nearly coincides
	 * with the H I Ly edge, they should be treated as one edge. Trying to separate them will
	 * almost certainly lead to erroneous behaviour in RebinAtmosphere */
	Edges[0] = (float)(RYDLAM/911.76);
	Edges[1] = (float)(RYDLAM/504.26);
	Edges[2] = (float)(RYDLAM/227.84);

	/* MALLOC some workspace */
	StarEner = (float *)MALLOC( sizeof(float)*NVALS );
	if( StarEner == NULL )
	{ 
		printf( " not enough memory to allocate StarEner in AtlasCompile\n" );
		puts( "[Stop in AtlasCompile]" );
		cdEXIT(EXIT_FAILURE);
	}
	ASSERT( StarEner !=NULL );

	StarFlux = (float *)MALLOC( sizeof(float)*NVALS );
	if( StarFlux == NULL )
	{ 
		printf( " not enough memory to allocate StarFlux in AtlasCompile\n" );
		puts( "[Stop in AtlasCompile]" );
		cdEXIT(EXIT_FAILURE);
	}
	ASSERT( StarEner != StarFlux );

	CloudyFlux = (float *)MALLOC( nBlocksize );
	if( CloudyFlux == NULL )
	{ 
		printf( " not enough memory to allocate CloudyFlux in AtlasCompile\n" );
		puts( "[Stop in AtlasCompile]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* use NVALS+4 because data are read in in multiples of 5 */
	nmalloc = MAX3(rfield.nupper,NATLAS,NVALS+4);
	scratch = (float *)MALLOC( sizeof(float)* (size_t)nmalloc);
	if( scratch == NULL )
	{ 
		printf( " not enough memory to allocate scratch in AtlasCompile\n" );
		puts( "[Stop in AtlasCompile]" );
		cdEXIT(EXIT_FAILURE);
	}
	ASSERT( scratch!=NULL);
	for( i=0; i<nmalloc; ++i )
	{
		scratch[i] = -1.f;
	}

	/* open main file of stellar atmospheres for reading */
	if( (ioIN = fopen( "kurucz.ascii", "r" ) ) == NULL )
	{
		char chFilename[FILENAME_PATH_LENGTH_2] = "junk";
		/* did not find it - try on the path */
		if( lgDataPathSet == TRUE )
		{
			/*path set, so look only there */
			strcpy( chFilename , chDataPath );
			strcat( chFilename , "kurucz.ascii" );
		}
		if( (ioIN = fopen( chFilename, "r" ) ) == NULL )
		{
			fprintf( ioQQQ, " AtlasCompile could not find kurucz.ascii.\n" );
			puts( "[Stop in getatlas]" );
			cdEXIT(EXIT_FAILURE);
		}
	}
	fprintf( ioQQQ, " AtlasCompile got kurucz.ascii.\n" );

	/* open much smaller list of stellar atmospheres */
	if( (ioLIST = fopen( "kurucz.list", "r" ) ) == NULL )
	{
		char chFilename[FILENAME_PATH_LENGTH_2] = "junk";
		/* did not find it - try on the path */
		if( lgDataPathSet == TRUE )
		{
			/*path set, so look only there */
			strcpy( chFilename , chDataPath );
			strcat( chFilename , "kurucz.list" );
		}
		if( (ioLIST = fopen( chFilename, "r" ) ) == NULL )
		{
			fprintf( ioQQQ, " AtlasCompile could not find kurucz.list.\n" );
			puts( "[Stop in getatlas]" );
			cdEXIT(EXIT_FAILURE);
		}
	}
	fprintf( ioQQQ, " AtlasCompile got kurucz.list.\n" );

	/* create the binary output file */
	if( (ioOUT = fopen( "atlas.mod", "wb" ) ) == NULL )
	{
		fprintf( ioQQQ, " AtlasCompile failed creating kurucz.mod.\n" );
		puts( "[Stop in getatlas]" );
		cdEXIT(EXIT_FAILURE);
	}
	fprintf( ioQQQ, " AtlasCompile got atlas.mod.\n" );

	/* >>chng 01 oct 17, add version and size to this array */
	/* first write out a version number for version checks when reading in */

	/* >>chng 01 oct 17, add version and size to this array */
	/* now write out nBlocksize, the size of the continuum array */
	if( fwrite( &nBlocksize,	1, sizeof(nBlocksize) ,ioOUT ) - sizeof(nBlocksize) )
	{
		fprintf( ioQQQ, " AtlasCompile failed writng nBlocksize.\n" );
		puts( "[Stop in getatlas]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* now write out rfield.nupper, the size of the continuum array */
	if( fwrite( &rfield.nupper,	1, sizeof(rfield.nupper) ,ioOUT ) - sizeof(rfield.nupper) )
	{
		fprintf( ioQQQ, " AtlasCompile failed writng rfield.nupper.\n" );
		puts( "[Stop in getatlas]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( fwrite( &VERSION,	1, sizeof(VERSION) ,ioOUT ) - sizeof(VERSION) )
	{
		fprintf( ioQQQ, " AtlasCompile failed writng VERSION.\n" );
		puts( "[Stop in getatlas]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* write out the cloudy energy grid for later sanity checks */
	if( fwrite( rfield.AnuOrg,	1, nBlocksize,ioOUT ) - nBlocksize  )
	{
		fprintf( ioQQQ, " AtlasCompile failed writng anu array.\n" );
		puts( "[Stop in getatlas]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* first read wavelenghth grid from start of large data file*/
	i = 0;
	while( i < NVALS )
	{
		if( fgets( chLine , (int)sizeof(chLine) , ioIN ) == NULL )
		{
			fprintf( ioQQQ, " AtlasCompile could not read kurucz wavelengths.\n" );
			fprintf( ioQQQ, " AtlasCompile: there is a problem in the file kurucz.ascii.\n" );
			puts( "[Stop in AtlasCompile]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* this is the wavelength grid */
		ASSERT( i+4 < nmalloc-1);
		sscanf( chLine , "%f %f %f %f %f" , 
			&scratch[i],&scratch[i+1],&scratch[i+2],&scratch[i+3],&scratch[i+4] );

		/* increment counter, five numbers per line*/
		i += 5;
	}

	fprintf( ioQQQ, 
		" AtlasCompile got kurucz wavelengths, low to high is %11.2e to %11.2e\n", 
		 scratch[0], scratch[NVALS-1] );

	/* the wavelength values scratch are in nm units, convert to freq (Hz)
	 * the energy range is 5.7e-4 to 10 Ryd */
	/* 	c = 299792458.0f; */
	/* 	rydb_h = 10967758.30f; */
	/* These are the speed of light and H-atom Rydberg constant in MKS units, for 
	 * conversion of frequency to rydberg units.  The infinite mass
	 * R value is */ 
	/*      rydb_inf = 10973731.534; */
	/* and the proton to electron mass ratio is 1836.152701 for this calculation. */
	/* 	con1 = 1.f/(c*rydb_inf); */

	for( i=0; i < NVALS; i++ )
	{
		/* convert wavelength (in nm) grid to inf. mass rydberg units */
		/* scratch[i] = (float)(299792458.e09/scratch[i]*con1); */
		ASSERT( scratch[i] > 0. );
		scratch[i] = (float)(0.1*RYDLAM/scratch[i]);
	}

	/* flip them over */
	for( j=0; j < NVALS; j++ )
	{
		/* StarEner now has correct energy grid, in Rydbergs, in increasing order */
		StarEner[j] = scratch[NVALS-j-1];
		ASSERT( StarEner[j] > 0. );

		/* sanity check */
		if( j > 0 )
		{
			ASSERT( StarEner[j] > StarEner[j-1] );
		}
	}

	/* now get grid of effective temperatures and surface gravities from smaller file */
	/* first line is header */
	if( fgets( chLine , (int)sizeof(chLine) , ioLIST ) == NULL )
	{
		fprintf( ioQQQ, " AtlasCompile could not read empty header.\n" );
		puts( "[Stop in AtlasCompile]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* sanity check - NATLAS is set 412 above*/
	ASSERT( rfield.nupper >= NATLAS );

	/* now read in NATLAS pairs of effective temperatures and surface gravities */
	for( j=0; j < NATLAS; j++ )
	{
		if( fgets( chLine , (int)sizeof(chLine) , ioLIST ) == NULL )
		{
			fprintf( ioQQQ, " AtlasCompile failed readin teff, logg, for model%4ld\n", 
			  j );
			puts( "[Stop in AtlasCompile]" );
			cdEXIT(EXIT_FAILURE);
		}
		
		/* parse the temperature and gravity off the line */
		sscanf( chLine , "%6li%s%8f %s %f " , 
			&i,
			chC12, 
			/* the temperature */
			&scratch[j], 
			chC7, 
			/* the gravity */
			&CloudyFlux[j] );
		/* print temp and grav */
		/*fprintf(ioQQQ," %.2f\t%.2f\n", scratch[j],CloudyFlux[j] );*/

	}
	fprintf( ioQQQ, " Got grid of temperatures and gravities.\n" );


	/* dump the temperatures into the binary file */
	if( fwrite( scratch , 1, nBlocksize, ioOUT ) - nBlocksize  )
	{
		fprintf( ioQQQ, " AtlasCompile failed writing teff.\n" );
		fprintf( ioQQQ, " I expected to write %li words, but fwrite was short\n",
			(rfield.nupper*(long)sizeof(float )) );
		puts( "[Stop in AtlasCompile]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* dump the gravities into the binary file */
	if( fwrite( CloudyFlux , 1, nBlocksize, ioOUT ) - nBlocksize )
	{
		fprintf( ioQQQ, " AtlasCompile failed writing grav.\n" );
		fprintf( ioQQQ, " I expected to write %li words, but fwrite was short\n",
			(long)nBlocksize );
		puts( "[Stop in AtlasCompile]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* >>chng 00 oct 24, removed redundant second copy of AnuOrg array, PvH */

	fprintf( ioQQQ, "\n\n\n\n" );
	fprintf( ioQQQ, " Now would be a good time to get on the Cloudy mailing list, to\n" );
	fprintf( ioQQQ, " find out about corrections and improvements to the code.\n" );
	fprintf( ioQQQ, " Go to http://www.nublado.org to get on the list.\n" );
	fprintf( ioQQQ, "\n\n\n\n" );
	fprintf( ioQQQ, " I will print every 10th model so you know I am still alive.\n" );
	fprintf( ioQQQ, " (there are %d models)\n", NATLAS );
	
	for( imod=0; imod < NATLAS; imod++ )
	{
		i = 0;
		while( i < NVALS )
		{
			if( fgets( chLine , (int)sizeof(chLine) , ioIN ) == NULL )
			{
				fprintf( ioQQQ, " AtlasCompile failed reading star flux.\n" );
				fprintf( ioQQQ, " AtlasCompile: there is a problem in the file kurucz.ascii.\n" );
				puts( "[Stop in AtlasCompile]" );
				cdEXIT(EXIT_FAILURE);
			}

			/* actually scan off the fluxes, 5 per line */
			sscanf( chLine , "%f %f %f %f %f" , 
				&scratch[i],&scratch[i+1],&scratch[i+2],&scratch[i+3],&scratch[i+4] );

			/* increment the counter */
			i += 5;
		}

		for( j=0; j < NVALS; j++ )
		{
			/* original models are in cgs H_nu units, so multiply with 4pi to get F_nu */
			StarFlux[j] = (float)(PI4*scratch[NVALS-1-j]);
		}

		/* actually do the rebinning. */
		RebinAtmosphere(NVALS, StarEner, StarFlux, CloudyFlux, 3L, Edges );

		{
			/* following should be set true to print contributors */
			/*@-redef@*/
			enum {DEBUG_LOC=FALSE};
			/*@+redef@*/
			if( DEBUG_LOC )
			{
				FILE *ioBUG;
				ioBUG = fopen("test.txt","w");
				for( i=0; i<rfield.nupper; ++i)
				{
					fprintf(ioBUG,"%e %e\n", rfield.AnuOrg[i], CloudyFlux[i]);
				}
				cdEXIT(EXIT_FAILURE);
			}
		}

		/* write out the cloudy flux */
		if( fwrite( CloudyFlux , 1, nBlocksize, ioOUT ) - nBlocksize  )
		{
			fprintf( ioQQQ, " AtlasCompile failed writing CloudyFlux star flux.\n" );
			fprintf( ioQQQ, " I expected to write %li words, but fwrite was short\n",
				(long)nBlocksize );
			puts( "[Stop in AtlasCompile]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* let the world know we are still alive, every 10th model */
		if( imod%10 == 0 )
		{
			fprintf( ioQQQ, " AtlasCompile wrote%4ld\n", imod );
		}
	}

	fclose( ioOUT );
	fclose( ioIN );
	fclose( ioLIST );

	free( StarEner );
	free( StarFlux );
	free( scratch );
	free( CloudyFlux );

	fprintf( ioQQQ, " AtlasCompile completed OK, kurucz.ascii and kurucz.list files may be deleted.\n" );

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


/*AtlasInterpolate interpolate on atlas model atmospheres, 
 * original version by K Volk */
void AtlasInterpolate(long int *nstar, 
  double temp, 
  double alogg, 
  int lgTrace)
{
	/* will point to binary file */
	FILE *ioIN;
	char chLine[FILENAME_PATH_LENGTH];
	long int i, 
	  ierr, 
	  j, 
	  jj, 
	  jmod, 
	  jval[61][11], 
	  k;
	double chk,
	  dev1, 
	  dev2,
	  xLuminosity;
	float *ReadArray;

	/* >>chng 01 oct 17, will get this from file */
	size_t nBlocksize /*= (size_t)rfield.nupper*sizeof(float )*/;

	long int version_atlas;
	float 
	  teff[NATLAS], 
	  xlogg[NATLAS];
	static float tval[61]={3500.,3750.,4000.,4250.,4500.,4750.,5000.,
	  5250.,5500.,5750.,6000.,6250.,6500.,6750.,7000.,7250.,7500.,
	  7750.,8000.,8250.,8500.,8750.,9000.,9250.,9500.,9750.,10000.,
	  10500.,11000.,11500.,12000.,12500.,13000.,14000.,15000.,16000.,
	  17000.,18000.,19000.,20000.,21000.,22000.,23000.,24000.,25000.,
	  26000.,27000.,28000.,29000.,30000.,31000.,32000.,33000.,34000.,
	  35000.,37500.,40000.,42500.,45000.,47500.,50000.};

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

	/* NATLAS is the number of accessable models in the new model set
	 * record 0 is the AnuOrg grid, record 1 holds the T_eff values,
	 * and record 2 holds the log(g) values, so there are NATLAS+3 records
	 * in the main file */

	if( lgDataPathSet )
	{
		/* path is set, generate full path name with file */
		strcpy( chLine , chDataPath );
		strcat( chLine,  "atlas.mod" );
	}
	else
	{
		/* path not set, look here */
		strcpy( chLine , "atlas.mod" );
	}

	/* put junk in indices so we know if set */
	for( j=0; j < 11; j++ )
	{
		for( k=0; k < 61; k++ )
		{
			jval[k][j] = -1;
		}
	}

	/* gjf open(unit=92,file='atlas.mod', */
	if( lgTrace )
	{
		fprintf( ioQQQ, " About to open file=%s\n", chLine );
	}

	if( (ioIN = fopen( chLine , "rb" )) == NULL )
	{
		fprintf( ioQQQ, "Error: Kurucz stellar atmosphere file atlas.mod not found.\n" );
		fprintf( ioQQQ, " The path I tried: ==%s==\n", chLine );
		fprintf( ioQQQ, " And here comes its hexadecimal representation:\n" );
		/* write(qq,'(1X,A1,O4)') (chLine(i:i),chLine(i:i),i=1,80) */
		for( i=0; i < FILENAME_PATH_LENGTH; i++ )
		{
			fprintf( ioQQQ, " '%c'=%#02x", chLine[i], (unsigned int)chLine[i] );
			if( chLine[i] == '\0' ) {
				break;
			}
		}
		fprintf( ioQQQ, "\n" );
		puts( "[Stop in AtlasInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* >>chng 01 oct 17, add version and size to this array */

	/* >>chng 01 oct 17, add version and size to this array */
	/* now write out nBlocksize, the size of the continuum array */
	if( fread( &nBlocksize,	1, sizeof(nBlocksize) ,ioIN ) - sizeof(nBlocksize) )
	{
		fprintf( ioQQQ, " AtlasInterpolate failed reading nBlocksize.\n" );
		puts( "[Stop in AtlasInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}
	if( fread( &rfield.nupper,	1, sizeof(rfield.nupper) ,ioIN ) - sizeof(rfield.nupper) )
	{
		fprintf( ioQQQ, " AtlasInterpolate failed reading nBlocksize.\n" );
		puts( "[Stop in AtlasInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}
	/* first write out a version number for version checks when reading in */
	if( fread( &version_atlas,	1, sizeof(VERSION) ,ioIN ) - sizeof(VERSION) )
	{
		fprintf( ioQQQ, " AtlasInterpolate failed reading VERSION.\n" );
		puts( "[Stop in AtlasInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}
	/* this had better be zero */
	if( VERSION-version_atlas )
	{
		fprintf( ioQQQ, " AtlasInterpolate: there is a version mismatch between the compiled atlas I expected and the one I found.\n" );
		fprintf( ioQQQ, " AtlasInterpolate: Please recompile the atlas stellar atmospheres with the COMPILE STARS ONLY ATLAS command.\n" );
		puts( "[Stop in AtlasInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* this is th total offset, used for the seeks */
	nOffset= sizeof(nBlocksize) + sizeof(rfield.nupper) + sizeof(VERSION);

	if( lgTrace )
	{
		fprintf( ioQQQ, " About to get wavelengths\n" );
	}

	/* MALLOC scratch space */
	if( (ReadArray = (float*)MALLOC(nBlocksize) )==NULL )
		BadMalloc();

	/* read in the saved cloudy energy scale so we can confirm this is a good image */
	if( fread( ReadArray, 1, nBlocksize,ioIN ) - nBlocksize  )
	{
		fprintf( ioQQQ, " problem trying to read atlas.mod wavelengths \n" );
		fprintf( ioQQQ, " I expected to read %li words, but fread was short\n",
			(long)nBlocksize );
		puts( "[Stop in AtlasInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* tnu array will be used to hold freq grid (should be the same as AnuOrg) -
	 * this will be verified AFTER continuum energy mesh is known (it is not
	 * known when this routine is executed) */
	for( j=0; j < rfield.nupper; j++ )
	{
		rfield.tNuRyd[j][rfield.nspec] = ReadArray[j];
	}

	if( lgTrace )
	{
		fprintf( ioQQQ, " About to get temps\n" );
	}

	/* next we read the the teff and log(g) values for the models in the
	 * set; if the parameters are in range, interpolation is done first in
	 * temperature and then in log(g)
	 *
	 * start by reading in the effective temperatures */

	if( fread( ReadArray, 1, nBlocksize,ioIN ) -  nBlocksize )
	{
		fprintf( ioQQQ, " problem trying to read atlas.mod temperatures \n" );
		fprintf( ioQQQ, " I expected to read %li words, but fread was short\n",
			(long)nBlocksize );
		puts( "[Stop in AtlasInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* now copy over to actual teff array */
	for( j=0; j < NATLAS; j++ )
	{
		teff[j] = ReadArray[j];
	}

	if( lgTrace )
	{
		fprintf( ioQQQ, " About to get gravs\n" );
	}

	/* now read the log(g) values */
	if( fread( ReadArray, 1, nBlocksize,ioIN ) - nBlocksize )
	{
		fprintf( ioQQQ, " problem trying to read atlas.mod gravities \n" );
		fprintf( ioQQQ, " I expected to read %li words, but fread was short\n",
			(long)nBlocksize );
		puts( "[Stop in AtlasInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* and copy them over as well */
	for( j=0; j < NATLAS; j++ )
	{
		xlogg[j] = ReadArray[j];
	}

	/* check whether the models have the correct effective temperature, for debugging only */
	if( FALSE ) 
	{
		for( j=0; j < NATLAS; j++ ) 
		{
			if( fread( ReadArray, 1, nBlocksize, ioIN ) - nBlocksize ) 
			{
				fprintf( ioQQQ, " problem trying to read atlas model %li\n", j+1 );
				fprintf( ioQQQ, " I expected to read %li words, but fread was short\n",
					 (long)nBlocksize );
				puts( "[Stop in AtlasInterpolate]" );
				cdEXIT(EXIT_FAILURE);
			}
			xLuminosity = 0.;
			/* rebinned models are in cgs F_nu units */
			for( k=1; k < rfield.nupper; k++ ) 
			{
				xLuminosity += (rfield.tNuRyd[k][rfield.nspec] - rfield.tNuRyd[k-1][rfield.nspec])*
					(ReadArray[k] + ReadArray[k-1])/2.;
			}
			/* now convert luminosity to effective temperature */
			chk = pow(xLuminosity*FR1RYD/STEFAN_BOLTZ,0.25);
			/* allow a tolerance of 0.5% */
			if( fabs((teff[j]-chk)/teff[j]) > 0.005 ) 
			{
				printf("*** WARNING, T_eff discrepancy for model %li, expected teff %.2f, "
				       "log(g) %.2f, integration yielded teff %.2f\n",
				       j+1, teff[j], xlogg[j], chk );
			}
		}
	}
		
	for( j=0; j < NATLAS; j++ )
	{
		jj = (long)(2.*xlogg[j]+0.001) + 1;

		if( jj < 1 || jj > 11 )
		{
			fprintf( ioQQQ, " AtlasInterpolate finds insane jj:%5ld j, xlogg(j),=%5ld%10.2e\n", 
			  jj, j, xlogg[j] );
			puts( "[Stop in AtlasInterpolate]" );
			cdEXIT(EXIT_FAILURE);
		}
		for( k=0; k < 61; k++ )
		{
			/* >>chng 00 oct 24, use C scale, PvH */
			/* this confirms that indices are within bounds - pclint was not sure of this */
			ASSERT( (k < 61) && (jj-1)<11 );
			/*lint -e777 test float equal */
			if( teff[j] == tval[k] )
			/*lint +e777 test float equal */
				jval[k][jj-1] = j;
		}
	}

	/*  determine which model or models are needed.*/

	/* TEMP and ALOGG are desired temp and log g */

	jmod = -1;
	for( j=0; j < NATLAS; ++j )
	{
		dev1 = fabs(temp/teff[j]-1.);
		/* this is because log(g) can be zero for the new set of models */
		if( xlogg[j] > 0.4 )
		{
			dev2 = fabs(alogg/xlogg[j]-1.);
		}
		else
		{
			dev2 = fabs(alogg-xlogg[j]);
		}

		/* this is test for exact model, result is jmod still -1 if not hit */
		if( dev1 <= 10.*FLT_EPSILON && dev2 <= 10.*FLT_EPSILON )
		{
			jmod = j;
			break;
		}
	}

	if( alogg < 0. || alogg > 5. )
	{
		/* this rejects temps outside the range of models */
		fprintf( ioQQQ, " Gravity outside range of 0.0 to 5.0\n" );
		puts( "[Stop in AtlasInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* following blocks either a) interpolate or b) get exact model */
	if( jmod == -1 )
	{
		if( temp < 3500. || temp > 50000. )
		{
			/*  this rejects values outside the range of model temperatures */
			fprintf( ioQQQ, " Temperature outside range of 3500 to 50,000K\n" );
			puts( "[Stop in AtlasInterpolate]" );
			cdEXIT(EXIT_FAILURE);
		}
		if( lgTrace )
		{
			fprintf( ioQQQ, " About to call GETATL.\n" );
		}

		/*  AtlasGetOne carries out the interpolation in t_eff and log(g) */
		AtlasGetOne(temp,alogg,jval,tval,&ierr,teff,xlogg , ioIN , ReadArray , nBlocksize);

		if( ierr == 1 )
		{
			fprintf( ioQQQ, " GETATL returned ierr=1, stop.\n" );
			puts( "[Stop in AtlasInterpolate]" );
			cdEXIT(EXIT_FAILURE);
		}
	}
	else
	{
		if( lgTrace )
		{
			fprintf( ioQQQ, " Exact model about to be read in.\n" );
		}

		/* this branch, we hit a proper model, use exactly that one */
		/* +3 below is to skip first three records, which are cloudy energy
		 * scale, teff, and log g */
		if( fseek(ioIN,  (long)((jmod+3)*nBlocksize+nOffset), SEEK_SET ) )
		{
			fprintf( ioQQQ, " Error seeking exact Atlas atmosphere%4ld\n", 
			  jmod );
			puts( "[Stop in AtlasInterpolate]" );
			cdEXIT(EXIT_FAILURE);
		}

		if( fread( ReadArray, 1 , nBlocksize, ioIN ) - nBlocksize )
		{
			fprintf( ioQQQ, " Error trying to read exact Atlas atmosphere%4ld\n", 
			  jmod );
			puts( "[Stop in AtlasInterpolate]" );
			cdEXIT(EXIT_FAILURE);
		}

		if( called.lgTalk )
		{
			fprintf( ioQQQ, 
				"                       * <<< K. Volk s Kurucz model %3ld"
				" read in. T_eff = %8.2f log(g) = %8.5f  >>> *\n", 
			  jmod+1, teff[jmod], xlogg[jmod] );
		}

	}

	/* sanity check: see whether this model has the correct effective temperature */
	xLuminosity = 0.;
	/* rebinned models are in cgs F_nu units */
	for( k=1; k < rfield.nupper; k++ ) 
	{
		xLuminosity += (rfield.tNuRyd[k][rfield.nspec] - rfield.tNuRyd[k-1][rfield.nspec])*
			(ReadArray[k] + ReadArray[k-1])/2.;
	}
	/* now convert luminosity to effective temperature */
	chk = pow(xLuminosity*FR1RYD/STEFAN_BOLTZ,0.25);

	/* allow a tolerance of 1% */
	if( fabs((temp-chk)/temp) > 0.01 && called.lgTalk ) 
	{
		/* this means that the claimed effective temperature is not equal to the
		 * actual effective temperature.  For this model the difference would be
		 * a logical error */
		fprintf( ioQQQ,
			 "*** PROBLEM, T_eff discrepancy for this model, expected teff %.2f, "
			 "log(g) %.2f, integration yielded teff %.2f, delta %.2f%%\n",
			 temp, 
			 alogg, 
			 chk, 
			 (chk-temp)/temp*100. );
		TotalInsanity();
	}

	/* Note on the interpolation: 26 October 2000 (Peter van Hoof)
	 *
	 * I computed the effective temperature for a random sample of interpolated
	 * atmospheres by integrating the flux as shown above and compared the results
	 * with the expected effective temperature using DELTA = (COMP-EXPEC)/EXPEC.
	 *
	 * I found that the average discrepacy was:
	 *
	 *     DELTA = -0.10% +/- 0.06% (sample size 5000)
	 *
	 * The most extreme discrepancies were
	 *     -0.30% <= DELTA <= 0.21%
	 *
	 * The most negative discrepancies were for T_eff =  36 -  39 kK, log(g) = 4.5 - 5
	 * The most positive discrepancies were for T_eff = 3.5 - 4.0 kK, log(g) = 0 - 1
	 *
	 * The interpolation in the ATLAS grid is clearly very accurate */

	for( j=0; j < rfield.nupper; j++ )
	{
		/* TSLOP is f_nu for each TNU point */
		rfield.tslop[j][rfield.nspec] = ReadArray[j];
	}

	*nstar = rfield.nupper;

	free( ReadArray );

	fclose(ioIN);

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

/*AtlasGetOne get one of the Atlas model atmospheres, coded by K. Volk */
static void AtlasGetOne(double temp, 
  double alogg, 
  long int jval[][11], 
  float tval[], 
  long int *ierr, 
  float teff[], 
  float xlogg[],
  FILE *ioIN ,
  float *fluxarray ,
  size_t nBlocksize )
{
	long int i, 
	  ipGrav, 
	  j, 
	  ipTeff,
	  ipGravUp ,
	  ipTeffUp ;

	/* following used in binary read with fread, do not change to
	 * double unless binary file changed too */
	float *flux1/*[rfield.nupper]*/, 
	  *flux2/*[rfield.nupper]*/;

	 double  fr1, 
	  fr2, 
	  fr3, 
	  fr4, 
	  xval;
	 /* this is offset in call to fseek to skip initial records */
	 const long iOff = 3;

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

	 ASSERT( fluxarray != NULL );

	if( (flux1 = (float*)MALLOC(nBlocksize ) ) == NULL )
		BadMalloc();
	if( (flux2 = (float*)MALLOC(nBlocksize ) ) == NULL )
		BadMalloc();

	/* parameter (NVALS=1221,NATLAS=412) */
	for( j=0; j < 61; j++ )
	{
		for( i=9; i >= 0; i-- )
		{
			if( jval[j][i] == -1 )
				jval[j][i] = jval[j][i+1];
		}
	}

	/* determine which models are to be interpolated */
	ipTeff = -1;
	for( j=0; j < 60; j++ )
	{
		if( tval[j] <= (float)temp && tval[j+1] > (float)temp )
		{
			ipTeff = j;
			break;
		}
	}

	/* nb - ipTeff is pointer within c array */
	if( ipTeff < 0 )
	{
		if( tval[60] == (float)temp )
		{
			ipTeff = 60;
		}
		else
		{
			fprintf( ioQQQ, " Requested temperature of%10.2e is not within the range%10.2e%10.2e\n", 
			  temp, tval[0], tval[60] );
			puts( "[Stop in AtlasGetOne]" );
			cdEXIT(EXIT_FAILURE);
		}
	}
	ipTeffUp = MIN2( ipTeff+1 , 60 );

	xval = 2.*alogg ;
	ipGrav = (long)(xval);
	ipGrav = MIN2( ipGrav , 10 );
	ipGravUp = MIN2( ipGrav+1 , 10 );

	/* Interpolation is carried out first in Teff and then in LOG(g).
	 * Two Kurucz models are read in for the lower LOG(g) value first, which
	 * I denote models A and B with B the higher T model [ipTeff+1 rather than ipTeff].
	 * unless alogg is 5.0 a second pair of Kurucz models are read in for the
	 * higher LOG(g) value, denoted models C and D.  Then the interpolation is
	 *
	 *  RESULT = fr2*(fr3*B+fr4*A)+fr1*(fr3*D+fr4*C)
	 *
	 * when one has 4 models and 
	 *
	 * RESULT = (fr3*B+fr4*A) 
	 *
	 * when one has just 2 models. */
	fr1 = xval - (float)(ipGrav);
	fr2 = 1. - fr1;
	if( ipTeffUp==ipTeff )
	{
		fr3 = 0.;
	}
	else
	{
		fr3 = ((float)temp - tval[ipTeff])/(tval[ipTeffUp] - tval[ipTeff]);
	}
	fr4 = 1. - fr3;

	/* this is check that we are within rational bounds */
	if( fr3 < 0. || fr3 > 1. )
	{
		fprintf( ioQQQ, "fr3 insanity in AtlasGetOne\n" );
		ShowMe();
		puts( "[Stop in AtlasGetOne]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( fr2 < 0. || fr2 > 1. )
	{
		fprintf( ioQQQ, "fr2 insanity in AtlasGetOne\n" );
		ShowMe();
		puts( "[Stop in AtlasGetOne]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* at this point we should have valid gravity (ipGrav) and
	 * teff (ipTeff) pointers */
	if( fseek(ioIN,  (long)((jval[ipTeff][ipGrav]+iOff)*(nBlocksize)+nOffset), SEEK_SET ) )
	{
		fprintf( ioQQQ, " Error1 seeking Atlas atmosphere%4ld\n", 
		  jval[ipTeff][ipGrav]+1 );
		puts( "[Stop in AtlasGetOne]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( fread( flux1, 1 , nBlocksize, ioIN ) - nBlocksize )
	{
		fprintf( ioQQQ, " Error1 trying to read Atlas atmosphere%4ld\n", 
		  jval[ipTeff][ipGrav]+1 );
		puts( "[Stop in AtlasGetOne]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( called.lgTalk )
	{
		fprintf( ioQQQ, 
			"                       * <<< K. Volk s Kurucz model %3ld"
			" read in. T_eff = %8.2f log(g) = %8.5f  >>> *\n", 
		  jval[ipTeff][ipGrav]+1, 
		  /* >>chng 00 oct 24, corrected index for teff (repeated 3 more times below), PvH */
		  teff[jval[ipTeff][ipGrav]], 
		  xlogg[jval[ipTeff][ipGrav]] );
	}

	if( fseek(ioIN,  (long)((jval[ipTeffUp][ipGrav]+iOff)*nBlocksize+nOffset), SEEK_SET ) )
	{
		fprintf( ioQQQ, " Error1 seeking Atlas atmosphere%4ld\n", 
		  jval[ipTeff][ipGrav]+1 );
		puts( "[Stop in AtlasGetOne]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( fread( flux2, 1 , nBlocksize, ioIN ) - nBlocksize )
	{
		fprintf( ioQQQ, " Error2 trying to read Atlas atmosphere%4ld\n", 
		  jval[ipTeff][ipGrav]+1 );
		puts( "[Stop in AtlasGetOne]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( called.lgTalk )
	{
		fprintf( ioQQQ, 
			"                       * <<< K. Volk s Kurucz model %3ld"
			" read in. T_eff = %8.2f log(g) = %8.5f  >>> *\n", 
		   jval[ipTeffUp][ipGrav]+1, 
		   teff[jval[ipTeffUp][ipGrav]], 
		  xlogg[jval[ipTeffUp][ipGrav]] );
	}


	for( j=0; j < rfield.nupper; j++ )
	{
		fluxarray[j] = (float)(
			fr4*log10(MAX2(1e-37,flux1[j]) )  + 
			fr3*log10(MAX2(1e-37,flux2[j]) ) );
	}

	if( ipGrav < 10 )
	{
		/* for LOG(g) < 5.0 repeat the interpolation at the larger LOG(g) then
		 * carry out the interpolation in log(g)
		 * >>chng 97 jul 24, get rid of go to */
		if( fseek(ioIN,  (long)((jval[ipTeff][ipGravUp]+iOff)*nBlocksize+nOffset), SEEK_SET ) )
		{
			fprintf( ioQQQ, " Error3 seeking Atlas atmosphere%4ld\n", 
			  jval[ipTeff][ipGrav]+1 );
			puts( "[Stop in AtlasGetOne]" );
			cdEXIT(EXIT_FAILURE);
		}

		if( fread( flux1, 1 , nBlocksize, ioIN ) - nBlocksize )
		{
			fprintf( ioQQQ, " Error3 trying to read Atlas atmosphere%4ld\n", 
			  jval[ipTeff][ipGravUp]+1 );
			puts( "[Stop in AtlasGetOne]" );
			cdEXIT(EXIT_FAILURE);
		}

		if( called.lgTalk )
		{
			fprintf( ioQQQ, 
				"                       * <<< K. Volk s Kurucz model %3ld"
				" read in. T_eff = %8.2f log(g) = %8.5f  >>> *\n", 
			        jval[ipTeff][ipGravUp]+1, 
			   teff[jval[ipTeff][ipGravUp]], 
			  xlogg[jval[ipTeff][ipGravUp]] );
		}

		if( fseek(ioIN,  (long)((jval[ipTeffUp][ipGravUp]+iOff)*nBlocksize+nOffset), SEEK_SET ) )
		{
			fprintf( ioQQQ, " Error4 seeking Atlas atmosphere%4ld\n", 
			  jval[ipTeff][ipGrav]+1 );
			puts( "[Stop in AtlasGetOne]" );
			cdEXIT(EXIT_FAILURE);
		}

		if( fread( flux2, 1 , nBlocksize, ioIN ) - nBlocksize )
		{
			fprintf( ioQQQ, " Error4 trying to read Atlas atmosphere%4ld\n", 
			  jval[ipTeff][ipGrav]+1 );
			puts( "[Stop in AtlasGetOne]" );
			cdEXIT(EXIT_FAILURE);
		}

		if( called.lgTalk )
		{
			fprintf( ioQQQ, 
				"                       * <<< K. Volk s Kurucz model %3ld"
				" read in. T_eff = %8.2f log(g) = %8.5f  >>> *\n", 
			       jval[ipTeffUp][ipGravUp]+1, 
			   teff[jval[ipTeffUp][ipGravUp]], 
			  xlogg[jval[ipTeffUp][ipGravUp]] );
		}

		for( j=0; j < rfield.nupper; j++ )
		{
			flux1[j] = (float)(
				fr4*log10(MAX2(1e-37,flux1[j] ) ) + 
				fr3*log10(MAX2(1e-37,flux2[j] ) ) );
		}

		for( j=0; j < rfield.nupper; j++ )
		{
			/* >>chng 96 dec 27, as per K Volk change, 
			 * to interpolation in log of space*/
			if( fluxarray[j]>-36. && flux1[j]>-36. )
			{
				fluxarray[j] = 
					(float)pow(10.,fr2*fluxarray[j] + fr1*flux1[j]);
			}
			else
			{
				fluxarray[j] = 0.;
			}
		}
	}
	else
	{
		/* One gets here if log(g) is 5.0 and there is no higher log(g) pair of 
		 * models to interpolate.
		 * */
		for( j=0; j < rfield.nupper; j++ )
		{
			fluxarray[j] = (float)pow(10.f,fluxarray[j]);
			if( fluxarray[j]<1e-36 )
				fluxarray[j] = 0.;
		}
	}

	*ierr = 0;
	free( flux1);
	free( flux2);

#	ifdef DEBUG_FUN
	fputs( " <->AtlasGetOne()\n", debug_fp );
#	endif
	return;
}
/*lint -e613 possible use of null pointer */

