/* This file is part of Cloudy and is copyright (C)1978-2006 by Gary J. Ferland
 * For conditions of distribution and use see copyright notice in license.txt */
/*CoStarInterpolate read in and interpolate on costar grid of windy O atmospheres */
/*CoStarCompile rebin costar stellar atmospheres to match cloudy energy grid, 
 *called by the compile stars command */
#include "cddefines.h"
#include "physconst.h"
#include "rfield.h"
#include "called.h"
#include "path.h"
#include "optimize.h"
#include "continuum.h"
#include "stars.h"

/*lint -e713 loss of sign in promotion */
/*lint -e737 loss of sign in promotion */

/* will use used to remember amount of header information */
static size_t nBlocksize;

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

/* internal routines */
static void SortCoStar(float[],float[],float[],float[],char[],int[],float[],int[],long int,int*);
static void FindHCoStar(long int,double,int,float[],int[],float[],int[],long int,long int);
static void FindVCoStar(double,int,float[],float[],int[],long int);
static void ReadCoStar(FILE*,long int,long int,float[],float[],float[],char[],int[]);
static void InterpolateCoStar(float[],float,float[],float);

#define WORK(I,J) (Work[(I)*nModels+(J)])
#define IWORK(I,J) (iWork[(I)*nModels+(J)])
#define WORK2(I,J) (Work2[(I)*nTracks+(J)])
#define IWORK2(I,J) (iWork2[(I)*nTracks+(J)])

/* this is to turn on debug print statements */
#define DEBUGPRT 0


/*CoStarInterpolate read in and interpolate on costar grid of windy O atmospheres */
void CoStarInterpolate(
	/* number of points in the output continuum, will be ncell */
	long int *nstar, 
	/* which interpolation mode is requested
	 * intmode = 1: use Teff and nmodid
	 * intmode = 2: use Teff and log(g)
	 * intmode = 3: use M_ZAMS and age
	 * intmode = 4: use age and M_ZAMS */
	int intmode,
	/* Teff for intmode = 1,2; M_ZAMS for intmode = 3, age for intmode = 4 */
	double par1, 
	/* not valid for intmode = 1; log(g) for intmode = 2; age for intmode = 3, M_ZAMS for intmode = 4 */
	double par2,
	/* the index number for the model sequence, only valid for intmode = 1 */
	int nmodid ,
	/* flag indicating whether solar (==0) or halo (==1) abundances */
	int iabund )
{
	/* will be used for reading in werner.mod */
	FILE *ioIN ;

	long int version;

	char chLine[FILENAME_PATH_LENGTH_2],
		chInFile[FILENAME_PATH_LENGTH_2];
	char *chGrid;

	int *iWork,	  
	  *iWork2,
	  iWork3[2],
	  i2,
	/* this will be an integer between 1 and 7 saying which age sequence model is */
	  *modid,
	  ptr,
	  nTracks,
	  res;

	long int End,
	  j,
	  k,
	  ModelBase,
	  /* >>chng 00 dec 06, nModels promoted to long; this MUST match definition in CoStarCompile ! PvH */
	  nModels;

	float *Age,
	  *Gravity,
	  *Mass,
	  *Teff,
	  *scratch,
	  *wl/*[rfield.nupper]*/,
	  *Work,
	  *Work2,
	  Work3[5]={0.,0.,0.,0.,0.};

	double chk,
	  lumi;

	/*size_t nBlocksize = (size_t)rfield.nupper*sizeof(float );*/

	const float SECURE = (1.f + 20.f*FLT_EPSILON);

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

	/* structure of costar.mod 
	 *
	 * first comes three numbers, the blocksize, nupper, and version number,
	 * the next record contaings an image of the cloudy energy
	 * grid, to use for sanity checks
	 * next is the number of stellar continua, nModels
	 * followed by a list of temperatures, log(g)'s, ZAMS
	 * masses, ages and model ids for these models
	 * the remaining records are rebinned stellar continua
	 * there are a total of nModels of these
	 * ModelBase is pointing exactly at the beginning of
	 * the rebinned stellar continua in costar.mod
	 */

	if( intmode <= 0 || intmode > 4 )
	{
		fprintf( ioQQQ, " CoStarInterpolate called with insane value for intmode: %d.\n",intmode );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}
	else if( intmode == 4 )
	{
		/* swap parameters, hence mimic intmode = 3 */
		double temp;
		temp = par1;
		par1 = par2;
		par2 = temp;
	}
	/* use log10(M_ZAMS) internally */
	if( intmode == 3 || intmode == 4 )
		par1 = log10(par1);

	if( iabund == 0 )
	{
		/* the solar abundance star set */
		strcpy( chInFile , "costar_sol.mod");
	}
	else if( iabund ==1 )
	{
		/* the halo abundance star set */
		strcpy( chInFile , "costar_halo.mod");
	}
	else
	{
		strcpy( chInFile , "trash4lint");
		fprintf( ioQQQ, " CoStarInterpolate insane iabund.\n" );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* true is special path set, false if file in same dir */
	if( lgDataPathSet )
	{
		/* path is set, generate full path name with file */
		strcpy( chLine , chDataPath );
		strcat( chLine,  chInFile );
	}
	else
	{
		/* path not set, look here */
		strcpy( chLine , chInFile );
	}

	if( (ioIN = fopen( chLine , "rb" )) == NULL )
	{
		/* something went wrong */
		fprintf( ioQQQ, "ERROR: The CoStar stellar atmosphere file was not found.\nSorry.\n" );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* >>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, " CoStarInterpolate failed reading nBlocksize.\n" );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}
	if( fread( &rfield.nupper,	1, sizeof(rfield.nupper) ,ioIN ) - sizeof(rfield.nupper) )
	{
		fprintf( ioQQQ, " CoStarInterpolate failed reading nBlocksize.\n" );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}
	/* first write out a version number for version checks when reading in */
	if( fread( &version,	1, sizeof(VERSION) ,ioIN ) - sizeof(VERSION) )
	{
		fprintf( ioQQQ, " CoStarInterpolate failed reading VERSION.\n" );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}
	/* this had better be zero */
	if( VERSION-version )
	{
		fprintf( ioQQQ, " CoStarInterpolate: there is a version mismatch between the compiled costar I expected and the one I found.\n" );
		fprintf( ioQQQ, " CoStarInterpolate: Please recompile the costar stellar atmospheres with the command COMPILE STARS ONLY COSTAR.\n" );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( (wl = (float*)MALLOC(nBlocksize ) ) == NULL )
		BadMalloc();
	if( (scratch = (float*)MALLOC(nBlocksize ) ) == NULL )
		BadMalloc();

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

	/* now get number of models */
	/* >>chng 00 dec 06, nModels changed from int to long above, it caused problems here
	 * problem caught on Tru64 alpha, PvH */
	if( fread( &nModels, 1, sizeof(nModels),ioIN ) - (int)sizeof(nModels) )
	{
		fprintf( ioQQQ, " problem trying to read costar.mod wavelengths \n" );
		fprintf( ioQQQ, " I expected to read %li words, but fread was short\n",
			(long)sizeof(nModels) );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* allocate space for the stellar temperatures */
	if( (Teff = (float *)MALLOC( sizeof(float)*nModels )) == NULL )
	{ 
		fprintf( ioQQQ, " not enough memory to allocate Teff in CoStarInterpolate\n" );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* allocate space for the stellar gravity */
	if( (Gravity = (float *)MALLOC( sizeof(float)*nModels )) == NULL )
	{ 
		fprintf( ioQQQ, " not enough memory to allocate Gravity in CoStarInterpolate\n" );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* allocate space for ZAMS masses */
	if( (Mass = (float *)MALLOC( sizeof(float)*nModels )) == NULL )
	{ 
		fprintf( ioQQQ, " not enough memory to allocate Mass in CoStarInterpolate\n" );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* allocate space for model ages */
	if( (Age = (float *)MALLOC( sizeof(float)*nModels )) == NULL )
	{ 
		fprintf( ioQQQ, " not enough memory to allocate Age in CoStarInterpolate\n" );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* allocate temporary space */
	if( (Work = (float *)MALLOC( sizeof(float)*nModels*4 )) == NULL )
	{ 
		fprintf( ioQQQ, " not enough memory to allocate Work in CoStarInterpolate\n" );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* allocate space for the model integer ids */
	if( (modid = (int *)MALLOC( sizeof(int)*nModels )) == NULL )
	{ 
		fprintf( ioQQQ, " not enough memory to allocate modid in CoStarInterpolate nModels=%li\n",nModels );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* allocate temporary space */
	if( (iWork = (int *)MALLOC( sizeof(int)*nModels*3 )) == NULL )
	{ 
		fprintf( ioQQQ, " not enough memory to allocate iWork in CoStarInterpolate nModels=%li\n",nModels );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* allocate space for the model char ids */
	if( (chGrid = (char *)MALLOC( sizeof(char)*nModels )) == NULL )
	{ 
		fprintf( ioQQQ, " not enough memory to allocate chGrid in CoStarInterpolate nModels=%li\n",nModels );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* now read the grid of temperatures in the binary file */
	if( fread( Teff, 1, nModels*sizeof(float),ioIN ) - nModels*(int)sizeof(float)  )
	{
		fprintf( ioQQQ, " problem trying to read temperature grid for costar models\n" );
		fprintf( ioQQQ, " I expected to read %li words, but fread was short\n",
			nModels*sizeof(float) );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* now read the grid of gravities in the binary file */
	if( fread( Gravity, 1, nModels*sizeof(float),ioIN ) - nModels*(int)sizeof(float) )
	{
		fprintf( ioQQQ, " problem trying to read gravity for costar models\n" );
		fprintf( ioQQQ, " I expected to read %li words, but fread was short\n",
			nModels*sizeof(float) );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* now read the grid of ZAMS masses in the binary file */
	if( fread( Mass, 1, nModels*sizeof(float),ioIN ) - nModels*(int)sizeof(float) )
	{
		fprintf( ioQQQ, " problem trying to read ZAMS masses for costar models\n" );
		fprintf( ioQQQ, " I expected to read %li words, but fread was short\n",
			nModels*sizeof(float) );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* now read the grid of model ages in the binary file */
	if( fread( Age, 1, nModels*sizeof(float),ioIN ) - nModels*(int)sizeof(float) )
	{
		fprintf( ioQQQ, " problem trying to read model ages for costar models\n" );
		fprintf( ioQQQ, " I expected to read %li words, but fread was short\n",
			nModels*sizeof(float) );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* now read the grid of model char letters in the binary file */
	if( fread( chGrid, 1, nModels*sizeof(char),ioIN ) - nModels*(int)sizeof(char) )
	{
		fprintf( ioQQQ, " problem trying to read chGrid grid for costar models\n" );
		fprintf( ioQQQ, " I expected to read %li words, but fread was short\n",
			nModels*sizeof(char) );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* now read the grid of model ids in the binary file */
	if( fread( modid, 1, nModels*sizeof(int),ioIN ) - nModels*(int)sizeof(int) )
	{
		fprintf( ioQQQ, " problem trying to read model id grid for costar models\n" );
		fprintf( ioQQQ, " I expected to read %li words, but fread was short\n",
			nModels*sizeof(float) );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* after this point the model atmospheres start, so remember it... */
	ModelBase = ftell( ioIN );

	/* sanity check: does the file have the correct length ? */
	/* NOTE: this operation is not necessarily supported by all operating systems */
#ifdef SEEK_END
	res = fseek(ioIN,0,SEEK_END);
	if( res == 0 )
	{
		End = ftell( ioIN );
		if( (End - ModelBase) != (long)nModels*(long)nBlocksize )
		{
			fprintf( ioQQQ, " problem performing sanity check for size of costar.mod\n" );
			fprintf( ioQQQ, " I expected to find %ld words, but actually found %ld words\n",
				 ModelBase+nModels*nBlocksize,End );
			fprintf( ioQQQ, " please re-compile costar.mod\n" );
			puts( "[Stop in CoStarInterpolate]" );
			cdEXIT(EXIT_FAILURE);
		}
	}
#endif

	/* check whether the models have the correct effective temperature, for debugging only */
	if( FALSE ) 
	{
		fseek( ioIN, ModelBase, SEEK_SET );
		for( j=0; j < nModels; j++ ) {
			if( fread( scratch, 1, nBlocksize, ioIN ) - (int)nBlocksize  ) 
			{
				fprintf( ioQQQ, " problem trying to read CoStar model %li\n", j+1 );
				fprintf( ioQQQ, " I expected to read %li words, but fread was short\n",
					 (long)nBlocksize );
				puts( "[Stop in CoStarInterpolate]" );
				cdEXIT(EXIT_FAILURE);
			}
			lumi = 0.;
			/* rebinned models are in cgs F_nu units */
			for( k=1; k < rfield.nupper; k++ ) {
				lumi += (wl[k] - wl[k-1])*(scratch[k] + scratch[k-1])/2.;
			}
			/* now convert luminosity to effective temperature */
			chk = pow(lumi*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], Gravity[j], chk );
			}
		}
	}
		
	/* convert ZAMS masses to logarithms */
	for( j=0; j < nModels; j++ )
	{
		if( Mass[j] > 0.f )
		{
			Mass[j] = (float)log10(Mass[j]);
		}
		else
		{
			fprintf( ioQQQ, " I found a non-positive ZAMS mass in costar.mod (Mass=%f)\n", Mass[j] );
			fprintf( ioQQQ, " The file costar.mod is probably corrupt, please re-compile.\n" );
			puts( "[Stop in CoStarInterpolate]" );
			cdEXIT(EXIT_FAILURE);
		}
	}			

	/* sort the models according to track */
	SortCoStar(Teff,Gravity,Mass,Age,chGrid,modid,Work,iWork,nModels,&nTracks);

	/* now allocate some more temp workspace */
	if( (Work2 = (float *)MALLOC( sizeof(float)*nTracks*5 )) == NULL )
	{ 
		fprintf( ioQQQ, " not enough memory to allocate Work2 in CoStarInterpolate\n" );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( (iWork2 = (int *)MALLOC( sizeof(int)*nTracks*2 )) == NULL )
	{ 
		fprintf( ioQQQ, " not enough memory to allocate iWork2 in CoStarInterpolate\n" );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* first do horizontal search, i.e. search along individual tracks */
	for( j=0; j < nTracks; j++ )
	{
		if( intmode == 1 )
		{
			ptr = IWORK(1,j);
			if( IWORK(2,j) >= nmodid ) {
				ptr = ptr + nmodid - 1;
				IWORK2(0,j) = IWORK(0,ptr);
				IWORK2(1,j) = IWORK(0,ptr);
				WORK2(0,j) = 1.f;
				WORK2(1,j) = WORK(0,ptr);
				WORK2(2,j) = WORK(1,ptr);
				WORK2(3,j) = WORK(2,ptr);
				WORK2(4,j) = WORK(3,ptr);
			}
			else
			{
				IWORK2(0,j) = -1;
				WORK2(0,j) = -1.f;
			}
		}
		else
		{
			/*void FindHCoStar(long int,double,int,float[],int[],float[],int[],
			 * long int,long int);*/
			FindHCoStar(j,par2,intmode,Work,iWork,Work2,iWork2,nModels,nTracks);
		}
	}

#	if DEBUGPRT
	for( j=0; j<nTracks; j++ ) 
	{
		if( WORK2(0,j) >= 0.f ) 
			printf("track %c: models %c%d, %c%d, frac %f, vals %g %g %g %g\n",
				(char)('A'+j),
				chGrid[IWORK2(0,j)],
				modid[IWORK2(0,j)],
				chGrid[IWORK2(1,j)],
				modid[IWORK2(1,j)],
				WORK2(0,j),
				WORK2(1,j),
				WORK2(2,j),
				WORK2(3,j),
				WORK2(4,j));
	}
#	endif

	/* now do vertical search, i.e. interpolate between tracks */

	/*void FindVCoStar(double,int,float[],float[],int[],long int);*/
	FindVCoStar(par1,intmode,Work2,Work3,iWork3,nTracks);

	/* This should only happen when CoStarInterpolate is called in non-optimizing mode,
	 * when optimizing CoStarInterpolate should report back to optimize_func()...
	 * The fact that FindVCoStar allows interpolation between non-adjoining tracks
	 * should guarantee that this will not happen. */
	if( Work3[0] < 0.f )
	{
		fprintf( ioQQQ, " The parameters for the requested CoStar model are out of the range of tables we have.\n" );
		fprintf( ioQQQ, " A table in Hazy 1 gives the range of temperature, gravity, and age that we have.\n" );
		puts( "[Stop in CoStarInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}
	ASSERT( iWork3[0] >= 0 && iWork3[0] < nTracks );
	ASSERT( iWork3[1] >= 0 && iWork3[1] < nTracks );
	ASSERT( IWORK2(0,iWork3[0]) >= 0 && IWORK2(0,iWork3[0]) < (int)nModels );
	ASSERT( IWORK2(1,iWork3[0]) >= 0 && IWORK2(1,iWork3[0]) < (int)nModels );
	ASSERT( IWORK2(0,iWork3[1]) >= 0 && IWORK2(0,iWork3[1]) < (int)nModels );
	ASSERT( IWORK2(1,iWork3[1]) >= 0 && IWORK2(1,iWork3[1]) < (int)nModels );
#	if DEBUGPRT
	printf("interpolate between tracks %c and %c, frac %f, vals %g %g %g %g\n",'A'+iWork3[0],'A'+iWork3[1],
	       Work3[0],Work3[1],Work3[2],Work3[3],Work3[4]);
#	endif

	/* set limits for optimizer */
	if( optimize.lgVarOn )
	{
		i2 = (intmode >= 3 ? 3 : 1);
		if( intmode != 4 )
		{
			optimize.varang[optimize.nparm][0] = +FLT_MAX;
			optimize.varang[optimize.nparm][1] = -FLT_MAX;
			for( j=0; j < nTracks; j++ )
			{
				if( WORK2(0,j) >= 0.f )
				{
					float temp;
					if( i2 == 3 )
						temp = WORK2(i2,j);
					else
						temp = (float)log10(WORK2(i2,j));
					optimize.varang[optimize.nparm][0] =
						MIN2(optimize.varang[optimize.nparm][0],temp);
					optimize.varang[optimize.nparm][1] =
						MAX2(optimize.varang[optimize.nparm][1],temp);
				}
			}
		}
		else
		{
			int ptr0,ptr1;
			ptr0 = IWORK(0,IWORK(1,iWork3[0]));
			ptr1 = IWORK(0,IWORK(1,iWork3[1]));
			optimize.varang[optimize.nparm][0] = (float)log10(MAX2(Age[ptr0],Age[ptr1]));
#			if DEBUGPRT
			printf("set limit 0: (models %d, %d) %f %f\n",ptr0,ptr1,Age[ptr0],Age[ptr1]);
#			endif
			ptr0 = IWORK(0,IWORK(1,iWork3[0])+IWORK(2,iWork3[0])-1);
			ptr1 = IWORK(0,IWORK(1,iWork3[1])+IWORK(2,iWork3[1])-1);
			optimize.varang[optimize.nparm][1] = (float)log10(MIN2(Age[ptr0],Age[ptr1]));
#			if DEBUGPRT
			printf("set limit 1: (models %d, %d) %f %f\n",ptr0,ptr1,Age[ptr0],Age[ptr1]);
#			endif
		}
#		if DEBUGPRT
		printf("set limits: %f %f\n",optimize.varang[optimize.nparm][0],optimize.varang[optimize.nparm][1]);
#		endif

		/* check sanity of optimization limits */
		if( optimize.varang[optimize.nparm][1] <= optimize.varang[optimize.nparm][0] )
		{
			fprintf( ioQQQ, " CoStar: no room to optimize parameter: lower limit %.4f, upper limit %.4f.\n",
				 optimize.varang[optimize.nparm][0],optimize.varang[optimize.nparm][1] );
			puts( "[Stop in CoStarInterpolate]" );
			cdEXIT(EXIT_FAILURE);
		}
		else
		{
			/* make a bit of room for round-off errors */
			optimize.varang[optimize.nparm][0] *= SECURE;
			optimize.varang[optimize.nparm][1] /= SECURE;
		}
	}

	/* now read the actual models and interpolate if necessary */
	ReadCoStar(ioIN,IWORK2(0,iWork3[0]),ModelBase,scratch,Teff,Gravity,chGrid,modid);
	if( IWORK2(1,iWork3[0]) != IWORK2(0,iWork3[0]) )
	{
		float *flux1/*[rfield.nupper]*/,frac1,frac2;
		if( (flux1 = (float*)MALLOC(nBlocksize ) ) == NULL )
			BadMalloc();

		ReadCoStar(ioIN,IWORK2(1,iWork3[0]),ModelBase,flux1,Teff,Gravity,chGrid,modid);
		frac1 = WORK2(0,iWork3[0]);
		frac2 = 1.f - frac1;
#		if DEBUGPRT
		printf("Interpolate first track %g\n",frac1);
#		endif
		InterpolateCoStar(scratch,frac1,flux1,frac2);
		free(flux1);
	}
	if( iWork3[1] != iWork3[0] )
	{
		float *flux2/*[rfield.nupper]*/,frac3,frac4;
		if( (flux2 = (float*)MALLOC(nBlocksize ) ) == NULL )
			BadMalloc();

		ReadCoStar(ioIN,IWORK2(0,iWork3[1]),ModelBase,flux2,Teff,Gravity,chGrid,modid);
		if( IWORK2(1,iWork3[1]) != IWORK2(0,iWork3[1]) )
		{
			float *flux3/*[rfield.nupper]*/,frac5,frac6;
			if( (flux3 = (float*)MALLOC((size_t)(nBlocksize) ) ) == NULL )
				BadMalloc();

			ReadCoStar(ioIN,IWORK2(1,iWork3[1]),ModelBase,flux3,Teff,Gravity,chGrid,modid);
			frac5 = WORK2(0,iWork3[1]);
			frac6 = 1.f - frac5;
#			if DEBUGPRT
			printf("Interpolate second track %g\n",frac5);
#			endif
			InterpolateCoStar(flux2,frac5,flux3,frac6);
		}
		frac3 = Work3[0];
		frac4 = 1.f - frac3;
#		if DEBUGPRT
		printf("Interpolate inbetween tracks %g\n",frac3);
#		endif
		InterpolateCoStar(scratch,frac3,flux2,frac4);
		free( flux2);
	}

	/* now write some final info */
	if( called.lgTalk )
	{
		fprintf( ioQQQ, 
			"                       * <<< costar T_eff = %7.1f, log(g) = %4.2f, M(ZAMS) = %5.1f, age = ", 
			Work3[1],Work3[2],pow(10.f,Work3[3]));
		fprintf( ioQQQ, PrintEfmt("%8.2e", Work3[4]) );
		fprintf( ioQQQ, "  >>> *\n" );
	}

	/* sanity check: see whether this model has the correct effective temperature */
	lumi = 0.;
	/* rebinned models are in cgs F_nu units */
	for( k=1; k < rfield.nupper; k++ ) 
	{
		lumi += (wl[k] - wl[k-1])*(scratch[k] + scratch[k-1])/2.;
	}
	/* now convert luminosity to effective temperature */
	chk = pow(lumi*FR1RYD/STEFAN_BOLTZ,0.25);
	/* allow a tolerance of 5% */
	if( fabs((Work3[1]-chk)/Work3[1]) > 0.05 ) 
	{
		/* 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",
			 Work3[1], Work3[2], chk, (chk-Work3[1])/Work3[1]*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 = -1.16% +/- 0.69% (SOLAR models, sample size 4590)
	 *     DELTA = -1.17% +/- 0.70% (HALO models, sample size 4828)
	 *
	 * The most extreme discrepancies for the SOLAR models were
	 *     -3.18% <= DELTA <= -0.16%
	 *
	 * The most negative discrepancies were for  T_eff = 35 kK, log(g) = 3.5
	 * The least negative discrepancies were for T_eff = 50 kK, log(g) = 4.1
	 *
	 * The most extreme discrepancies for the HALO models were
	 *     -2.90% <= DELTA <= -0.13%
	 *
	 * The most negative discrepancies were for  T_eff = 35 kK, log(g) = 3.5
	 * The least negative discrepancies were for T_eff = 50 kK, log(g) = 4.1
	 *
	 * Since Cloudy checks the scaling elsewhere there is no need to re-scale 
	 * things here, but this inaccuracy should be kept in mind since it could
	 * indicate problems with the flux distribution */

	/* now stuff what we read into the arrays that will read them later with interpolate */
	/* create arrays and zero unused elements */
	rfield.lgContMalloc[rfield.nspec] = TRUE;
	if( (rfield.tNuRyd[rfield.nspec] = (float*)MALLOC((size_t)(NCELL*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.tslop[rfield.nspec] = (float*)MALLOC((size_t)(NCELL*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.tFluxLog[rfield.nspec] = (float*)MALLOC((size_t)(NCELL*sizeof(float)) ) ) == NULL )
		BadMalloc();

	for( j=0; j < rfield.nupper; j++ )
	{
		rfield.tslop[rfield.nspec][j] = scratch[j];
		/* we save the wavelength scale so that it can be cross checked against
		 * one inside cloudy - this checks whether the grid is for the current
		 * version of the wavelength scale */
		rfield.tNuRyd[rfield.nspec][j] = wl[j];
	}

	*nstar = rfield.nupper;
	fclose( ioIN );
	free( Teff);
	free( modid);
	free( Gravity);
	free( Mass);
	free( Age);
	free( Work);
	free( Work2);
	free( iWork);
	free( iWork2);
	free( chGrid);
	free( wl);
	free( scratch);

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

/* sort CoStar models according to track and index number, i.e. A1, A2, ..., An, B1, ..., Bm, etc... */
static void SortCoStar(
	/* first parameter: Teff[nModels] */
	float par1[],
	/* second parameter: Gravity[nModels] */
	float par2[],
	/* third parameter: Mass[nModels] */
	float par3[],
	/* fourth parameter: Age[nModels] */
	float par4[],
	/* chGrid[nModels]: identifier for evolutionary track */
	char chGrid[],
	/* modid[nModels]: step number in evolutionary track */
	int modid[],
	/* Work[4*nModels]: float workspace
	 * WORK(0,*): Teff of sorted models
	 * WORK(1,*): log(g) of sorted models
	 * WORK(2,*): ZAMS mass of sorted models
	 * WORK(3,*): Age of sorted models */
	float Work[],
	/* iWork[3*nModels]: int workspace
	 * IWORK(0,*): index number of sorted models (i.e., pointers into Teff, Gravity etc.)
	 * IWORK(1,j): start of track no. j (i.e., pointer into IWORK(0,*)), not defined for j >= nTracks
	 * IWORK(2,j): no. of models in track no. j, not defined for j >= nTracks */
	int iWork[],
	/* total number of atmosphere models */
	long int nModels,
	/* total number of tracks found */
	int *nTracks)
{
	int index;
	long int found,
	  i,
	  ptr,
	  tlen;
	char track;

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

	ptr = 0;
	*nTracks = 0;
	do
	{
		track = (char)('A'+*nTracks);
		index = 1;
		tlen = 0;
		IWORK(1,*nTracks) = (int)ptr;
		do
		{
			found = 0;
			for( i=0; i<nModels; i++ )
			{
				if( chGrid[i] == track && modid[i] == index )
				{
					WORK(0,ptr) = par1[i];
					WORK(1,ptr) = par2[i];
					WORK(2,ptr) = par3[i];
					WORK(3,ptr) = par4[i];
					IWORK(0,ptr) = (int)i;
					ptr++;
					found = 1;
					break;
				}
			}
			index++;
			tlen += found;
		} while( found > 0 );
		IWORK(2,*nTracks) = (int)tlen;
		(*nTracks)++;
	} while( tlen > 0 );
	(*nTracks)--;

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

/* find which models to use for interpolation along a given evolutionary track */
static void FindHCoStar(
	long int track,
	/* requested logg or age */
	double par2,
	/* interpolation mode, should be 2 (using Teff, logg) or 3 (using mass, age) */
	int intmode,
	/* internal workspace */
	float Work[],
	int iWork[],
	/* Work2 and iWork2 contain information needed for interpolation along tracks:
	 * the interpolated model will be calculated as:
	 *    WORK2(0,j)*MODEL(IWORK2(0,j)) + (1.-WORK2(0,j))*MODEL(IWORK2(1,j))
	 *
	 * Work2[5*nTracks]: float workspace
	 * WORK2(0,j): fractional coefficient for interpolation along track j
	 * WORK2(1,j): Teff for interpolated model along track j
	 * WORK2(2,j): log(g) for interpolated model along track j
	 * WORK2(3,j): Mass for interpolated model along track j
	 * WORK2(4,j): Age for interpolated model along track j */
	float Work2[],
	/* iWork2[2*nTracks]: int workspace
	 * IWORK2(0,j): model number for first model used in interpolation (i.e., index into Teff, etc.)
	 * IWORK2(1,j): model number for second model used in interpolation */
	int iWork2[],
	/* total number of atmosphere models */
	long int nModels,
	/* total number of evolutionary tracks */
	long int nTracks)
{
	int ptr;
	long int i2,
	  j;

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

	i2 = (intmode >= 3 ? 3 : 1);
	IWORK2(0,track) = -1;
	WORK2(0,track) = -1.f;
	ptr = IWORK(1,track);
	for( j=0; j < IWORK(2,track); j++ )
	{
		/* do we have an exact match ? */
		/* Caution: this implicitly assumes that none of the model
		 * parameters are zero. As of June 1999 this is the case. */
		if( fabs(par2-(double)WORK(i2,ptr+j))/(double)WORK(i2,ptr+j) <= 10.*FLT_EPSILON )
		{
			IWORK2(0,track) = IWORK(0,ptr+j);
			IWORK2(1,track) = IWORK(0,ptr+j);
			WORK2(0,track) = 1.f;
			WORK2(1,track) = WORK(0,ptr+j);
			WORK2(2,track) = WORK(1,ptr+j);
			WORK2(3,track) = WORK(2,ptr+j);
			WORK2(4,track) = WORK(3,ptr+j);
			break;
		}
	}

	if( IWORK2(0,track) >= 0 )
	{
#	        ifdef DEBUG_FUN
		fputs( " <->FindHCoStar()\n", debug_fp );
#	        endif
		return;
	}

	for( j=0; j < IWORK(2,track)-1; j++ )
	{
		/* do we interpolate ? */
		if( ((float)par2 - WORK(i2,ptr+j))*((float)par2 - WORK(i2,ptr+j+1)) < 0.f )
		{
			float frac;

			IWORK2(0,track) = IWORK(0,ptr+j);
			IWORK2(1,track) = IWORK(0,ptr+j+1);
			frac = ((float)par2 - WORK(i2,ptr+j+1))/(WORK(i2,ptr+j) - WORK(i2,ptr+j+1));
			WORK2(0,track) = frac;
			WORK2(1,track) = frac*WORK(0,ptr+j) + (1.f-frac)*WORK(0,ptr+j+1);
			WORK2(2,track) = frac*WORK(1,ptr+j) + (1.f-frac)*WORK(1,ptr+j+1);
			WORK2(3,track) = frac*WORK(2,ptr+j) + (1.f-frac)*WORK(2,ptr+j+1);
			WORK2(4,track) = frac*WORK(3,ptr+j) + (1.f-frac)*WORK(3,ptr+j+1);
			break;
		}
	}

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

/* find which tracks to use for interpolation inbetween tracks */
static void FindVCoStar(
	/* requested Teff or ZAMS mass */
	double par1,
	/* interpolation mode, should be 2 (using Teff, logg) or 3 (using mass, age) */
	int intmode,
	/* internal workspace */
	float Work2[],
	/* Work3 and iWork3 contain information needed for interpolation inbetween tracks:
	 * the interpolated model will be calculated as:
	 *    Work3[0]*MODEL_ALONG_TRACK(iWork3[0]) + (1.-Work3[0])*MODEL_ALONG_TRACK(iWork3[1])
	 *
	 * Work3[0]: fractional coefficient for interpolation inbetween tracks
	 * Work3[1]: Teff for interpolated model
	 * Work3[2]: log(g) for interpolated model
	 * Work3[3]: Mass for interpolated model
	 * Work3[4]: Age for interpolated model */
	float Work3[], /* Work3[5] */
	/* iWork3[0]: track number for first track used in interpolation (i.e., 0 is 'A', etc.)
	 * iWork3[1]: track number for second track used in interpolation
	 *
	 * NOTE: FindVCoStar raises a caution when interpolating between non-adjoining tracks,
	 *       i.e. when (iWork3[1]-iWork3[0]) > 1 */
	int iWork3[], /*iWork3[2] */
	/* total number of evolutionary tracks */
	long int nTracks)
{
	long int i1,
	  j;

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

	i1 = (intmode >= 3 ? 3 : 1);
	iWork3[0] = -1;
	Work3[0] = -1.f;
	for( j=0; j < nTracks; j++ )
	{
		/* do we have an exact match ? */
		if( WORK2(0,j) >= 0.f && fabs(par1-(double)WORK2(i1,j))/(double)WORK2(i1,j) <= 10.*FLT_EPSILON )
		{
			iWork3[0] = (int)j;
			iWork3[1] = (int)j;
			Work3[0] = 1.f;
			Work3[1] = WORK2(1,j);
			Work3[2] = WORK2(2,j);
			Work3[3] = WORK2(3,j);
			Work3[4] = WORK2(4,j);
			break;
		}
	}

	if( iWork3[0] >= 0 )
	{
#	        ifdef DEBUG_FUN
		fputs( " <->FindVCoStar()\n", debug_fp );
#	        endif
		return;
	}

	for( j=0; j < nTracks-1; j++ )
	{
		if( WORK2(0,j) >= 0.f )
		{
			long int i,j2;

			/* find next valid track */
			j2 = 0;
			for( i = j+1; i < nTracks; i++ )
			{
				if( WORK2(0,i) >= 0.f )
				{
					j2 = i;
					break;
				}
			}

			/* do we interpolate ? */
			if( j2 > 0 && ((float)par1-WORK2(i1,j))*((float)par1-WORK2(i1,j2)) < 0.f )
			{
				float frac;

				iWork3[0] = (int)j;
				iWork3[1] = (int)j2;
				frac = ((float)par1 - WORK2(i1,j2))/(WORK2(i1,j) - WORK2(i1,j2));
				Work3[0] = frac;
				Work3[1] = frac*WORK2(1,j) + (1.f-frac)*WORK2(1,j2);
				Work3[2] = frac*WORK2(2,j) + (1.f-frac)*WORK2(2,j2);
				Work3[3] = frac*WORK2(3,j) + (1.f-frac)*WORK2(3,j2);
				Work3[4] = frac*WORK2(4,j) + (1.f-frac)*WORK2(4,j2);
				break;
			}
		}
	}

	/* raise caution when we interpolate between non-adjoining tracks */
	continuum.lgCoStarInterpolationCaution = (iWork3[1]-iWork3[0] > 1);

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

static void ReadCoStar(
	FILE *ioIN,
	long int mod_num,
	long int ModelBase,
	float arr[], /* arr[NC_ELL] */
	float Teff[],
	float Gravity[],
	char chGrid[],
	int modid[])
{

	/*size_t nBlocksize = (size_t)rfield.nupper*sizeof(float );*/

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

	ASSERT( nBlocksize > 0 );

	/* position file */
	if( fseek(ioIN,(long)(ModelBase+mod_num*nBlocksize), SEEK_SET ) )
	{
		fprintf( ioQQQ, " Error seeking costar atmosphere%4ld\n", 
		  mod_num );
		puts( "[Stop in ReadCoStar]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* read model */
	if( fread( arr, 1 , nBlocksize, ioIN ) - nBlocksize )
	{
		fprintf( ioQQQ, " Error trying to read costar atmosphere%4ld\n", 
		  mod_num );
		puts( "[Stop in ReadCoStar]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( called.lgTalk )
	{
		fprintf( ioQQQ, 
			"                       * <<< costar     model %c%1d (%2ld) read in.  T_eff = %7.1f, log(g) = %4.2f       >>> *\n", 
			 chGrid[mod_num],modid[mod_num],mod_num+1,Teff[mod_num],Gravity[mod_num] );
	}

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

/* Interpolate between two CoStar models */
static void InterpolateCoStar(float flux1[], /* flux1[NC_ELL] */
			      float fac1,
			      float flux2[], /* flux2[NC_ELL] */
			      float fac2)
{
	long int j;
	double sum;

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

	for( j=0; j < rfield.nupper; j++ )
	{
#		if DEBUGPRT
		if( j >= 840 && j <= 844 ) {
			printf("Interpolation test %ld, flux1 %g, flux2 %g\n",j,flux1[j],flux2[j]);
		}
#		endif
		if( flux1[j] > SMALLFLOAT && flux2[j] > SMALLFLOAT )
		{
			sum = log(flux1[j])*fac1 + log(flux2[j])*fac2;
			flux1[j] = (float)exp(sum);
		}
		else
		{
			flux1[j] = 0.f;
		}
	}

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

/*======================================================================*/
/*CoStarCompile rebin costar stellar atmospheres to match cloudy energy grid, 
 *called by the compile stars command */
void CoStarCompile(void)
{
	char chLine[FILENAME_PATH_LENGTH_2];
	char chInFile[FILENAME_PATH_LENGTH_2];
	char chOutFile[FILENAME_PATH_LENGTH_2];
	long int i, 
		j ,
	  nskip , 
	  nModels , 
	  nWL ,
	  nAbund ;

	/* these will be malloced into large work arrays*/
	float *StarEner , *StarFlux , *CloudyFlux , *Teff , *Gravity , *Mass, *Age;
	/* these contain frequencies for the major absorption edges */
	float Edges[3];
	/* this will be an integer between 1 and 7 saying which age sequence model is */
	int *modid;
	char *chGrid;

	/* we will write the binary results to this file*/
	FILE *ioOUT ,
		/* will get the ascii from this file */
		* ioIN ;

	/*size_t nBlocksize = (size_t)rfield.nupper*sizeof(float );*/

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

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

	nBlocksize = (size_t)rfield.nupper*sizeof(float );
	ASSERT( nBlocksize>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);

	for( nAbund=0; nAbund<2; ++nAbund )
	{
		/* MALLOC some workspace */
		if( (CloudyFlux = (float *)MALLOC( nBlocksize )) == NULL )
		{ 
			fprintf( ioQQQ, " not enough memory to allocate CloudyFlux in CoStarCompile\n" );
			puts( "[Stop in CoStarCompile]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* This is a program to re-bin the costar stellar model spectra to match the 
		 * Cloudy grid.  For short wavelengths I will use a power law extrapolation
		 * of the model values (which should be falling rapidly) if needed.  At long
		 * wavelengths I will assume Rayleigh-Jeans from the last stellar model point
		 * to extrapolate to 1 cm wavelength. */

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

		/* costar.ascii is the original data file obtained off the web, 
		 * open as read only */
		if( nAbund == 0 )
		{
			/* the solar abundance star set */
			strcpy( chInFile , "Sc1_costar_z020_lb.fluxes");
			strcpy( chOutFile , "costar_sol.mod");
		}
		else if( nAbund ==1 )
		{
			/* the halo abundance star set */
			strcpy( chInFile , "Sc1_costar_z004_lb.fluxes");
			strcpy( chOutFile , "costar_halo.mod");
		}
		else
		{
			strcpy( chOutFile , "trash4lint");
			strcpy( chInFile , "trash4lint");
			fprintf( ioQQQ, " CoStarCompile insane nAbund.\n" );
			puts( "[Stop in CoStarCompile]" );
			cdEXIT(EXIT_FAILURE);
		}
		if( (ioIN = fopen( chInFile, "r" ) ) == NULL )
		{
			fprintf( ioQQQ, " CoStarCompile fails opening %s\n", chInFile );
			fprintf( ioQQQ, " Where is this file?\n" );
			puts( "[Stop in CoStarCompile]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* this will be the file we create, that will be read to compute models, 
		 * open to write binary */
		if( (ioOUT = fopen( chOutFile, "wb" ) ) == NULL )
		{
			fprintf( ioQQQ, " CoStarCompile fails creating %s\n",chOutFile );
			fprintf( ioQQQ, " This is impossible??\n" );
			puts( "[Stop in CoStarCompile]" );
			cdEXIT(EXIT_FAILURE);
		}
		fprintf( ioQQQ, " CoStarCompile reading %s and creating %s\n",chInFile ,chOutFile );


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

		ASSERT( rfield.nupper > 0 );
		/* 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);
		}

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

		/* next dump the Cloudy energy scale in the binary file.   This will
		 * be used for sanity checks when the file is read*/
		if( fwrite( rfield.AnuOrg, 1,nBlocksize,ioOUT ) -  nBlocksize )
		{
			fprintf( ioQQQ, " problem trying to write %s\n",chOutFile );
			fprintf( ioQQQ, " I expected to write %li words, but fwrite was short\n",
				(long)nBlocksize );
			puts( "[Stop in CoStarCompile]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* ======================================================================== */

		/* get first line and see how many more to skip */
		if( fgets( chLine , (int)sizeof(chLine) , ioIN ) == NULL )
		{
			fprintf( ioQQQ, " CoStarCompile fails reading energy grid1.\n" );
			puts( "[Stop in CoStarCompile]" );
			cdEXIT(EXIT_FAILURE);
		}
		sscanf( chLine , "%li" , &nskip );

		/* now skip the header information */
		for( i=0; i<nskip; ++i )
		{
			if( fgets( chLine , (int)sizeof(chLine) , ioIN ) == NULL )
			{
				fprintf( ioQQQ, " CoStarCompile fails skipping block1.\n" );
				puts( "[Stop in CoStarCompile]" );
				cdEXIT(EXIT_FAILURE);
			}
		}

		/* now get number of models and number of wavelengths */
		if( fgets( chLine , (int)sizeof(chLine) , ioIN ) == NULL )
		{
			fprintf( ioQQQ, " CoStarCompile fails reading energy grid1.\n" );
			puts( "[Stop in CoStarCompile]" );
			cdEXIT(EXIT_FAILURE);
		}
		sscanf( chLine , "%li%li" , &nModels , &nWL );

		fprintf(ioQQQ," there are %li models and %li wavelengths.\n", nModels , nWL );

		if( nModels <= 0 || nWL <= 0 )
		{
			fprintf( ioQQQ, " CoStarCompile scanned off impossible values of nModels or nWL=%li %li\n",
				nModels, nWL );
			puts( "[Stop in CoStarCompile]" );
			cdEXIT(EXIT_FAILURE);
		}

		if( (StarEner = (float *)MALLOC( sizeof(float)*nWL )) == NULL )
		{ 
			fprintf( ioQQQ, " not enough memory to allocate StarEner in CoStarCompile\n" );
			puts( "[Stop in CoStarCompile]" );
			cdEXIT(EXIT_FAILURE);
		}

		if( (StarFlux = (float *)MALLOC( sizeof(float)*nWL )) == NULL )
		{ 
			fprintf( ioQQQ, " not enough memory to allocate StarFlux in CoStarCompile\n" );
			puts( "[Stop in CoStarCompile]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* allocate space for the stellar temperatures */
		if( (Teff = (float *)MALLOC( sizeof(float)*nModels )) == NULL )
		{ 
			fprintf( ioQQQ, " not enough memory to allocate Teff in CoStarCompile nModels=%li\n",nModels );
			puts( "[Stop in CoStarCompile]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* allocate space for the stellar gravities */
		if( (Gravity = (float *)MALLOC( sizeof(float)*nModels )) == NULL )
		{ 
			fprintf( ioQQQ, " not enough memory to allocate Gravity in CoStarCompile nModels=%li\n",nModels );
			puts( "[Stop in CoStarCompile]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* allocate space for the ZAMS masses */
		if( (Mass = (float *)MALLOC( sizeof(float)*nModels )) == NULL )
		{ 
			fprintf( ioQQQ, " not enough memory to allocate Mass in CoStarCompile nModels=%li\n",nModels );
			puts( "[Stop in CoStarCompile]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* allocate space for the model ages */
		if( (Age = (float *)MALLOC( sizeof(float)*nModels )) == NULL )
		{ 
			fprintf( ioQQQ, " not enough memory to allocate Age in CoStarCompile nModels=%li\n",nModels );
			puts( "[Stop in CoStarCompile]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* allocate space for the model ids */
		if( (modid = (int *)MALLOC( sizeof(int)*nModels )) == NULL )
		{ 
			fprintf( ioQQQ, " not enough memory to allocate modid in CoStarCompile nModels=%li\n",nModels );
			puts( "[Stop in CoStarCompile]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* allocate space for the model letters */
		if( (chGrid = (char *)MALLOC( sizeof(char)*nModels )) == NULL )
		{ 
			fprintf( ioQQQ, " not enough memory to allocate chGrid in CoStarCompile nModels=%li\n",nModels );
			puts( "[Stop in CoStarCompile]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* get all Teff's for the stars */
		for( i=0; i<nModels; ++i )
		{
			if( fgets( chLine , (int)sizeof(chLine) , ioIN ) == NULL )
			{
				fprintf( ioQQQ, " CoStarCompile fails reading energy grid1.\n" );
				puts( "[Stop in CoStarCompile]" );
				cdEXIT(EXIT_FAILURE);
			}
			/* first letter on line is indicator of grid */
			chGrid[i] = chLine[0];
			/* get the model id number */
			sscanf( chLine+1 , "%i" , modid+i );
			/* get the temperature */
			sscanf( chLine+23 , "%g" , Teff+i );
			/* get the surface gravity */
			sscanf( chLine+31 , "%g" , Gravity+i );
			/* get the ZAMS mass */
			sscanf( chLine+7 , "%g" , Mass+i );
			/* get the model age */
			sscanf( chLine+15 , "%g" , Age+i );
		}

		ASSERT( nModels > 0 );

		/* now save the number of models in the binary file */
		if( fwrite( &nModels, 1, sizeof(long),ioOUT ) - (int)sizeof(long) )
		{
			fprintf( ioQQQ, " problem trying to write number of costar models\n" );
			fprintf( ioQQQ, " I expected to write %li words, but fwrite was short\n",
				(long)sizeof(long) );
			puts( "[Stop in CoStarCompile]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* now save the grid of temperatures in the binary file */
		if( fwrite( Teff, 1, nModels*sizeof(float),ioOUT ) - nModels*(int)sizeof(float) )
		{
			fprintf( ioQQQ, " problem trying to write temperature grid for costar models\n" );
			fprintf( ioQQQ, " I expected to write %ld words, but fwrite was short\n",
				nModels*sizeof(float) );
			puts( "[Stop in CoStarCompile]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* now save the grid of gravities in the binary file */
		if( fwrite( Gravity, 1, nModels*sizeof(float),ioOUT ) - nModels*(int)sizeof(float) )
		{
			fprintf( ioQQQ, " problem trying to write gravity for costar models\n" );
			fprintf( ioQQQ, " I expected to write %ld words, but fwrite was short\n",
				nModels*sizeof(float) );
			puts( "[Stop in CoStarCompile]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* now save the grid of gravities in the binary file */
		if( fwrite( Mass, 1, nModels*sizeof(float),ioOUT ) - nModels*(int)sizeof(float) )
		{
			fprintf( ioQQQ, " problem trying to write ZAMS masses for costar models\n" );
			fprintf( ioQQQ, " I expected to write %ld words, but fwrite was short\n",
				nModels*sizeof(float) );
			puts( "[Stop in CoStarCompile]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* now save the grid of model ages in the binary file */
		if( fwrite( Age, 1, nModels*sizeof(float),ioOUT ) - nModels*(int)sizeof(float) )
		{
			fprintf( ioQQQ, " problem trying to write ages for costar models\n" );
			fprintf( ioQQQ, " I expected to write %ld words, but fwrite was short\n",
				nModels*sizeof(float) );
			puts( "[Stop in CoStarCompile]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* now save the grid of model char letters in the binary file */
		if( fwrite( chGrid, 1, nModels*sizeof(char),ioOUT ) - nModels*(int)sizeof(char) )
		{
			fprintf( ioQQQ, " problem trying to write chGrid grid for costar models\n" );
			fprintf( ioQQQ, " I expected to write %ld words, but fwrite was short\n",
				nModels*sizeof(char) );
			puts( "[Stop in CoStarCompile]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* now save the grid of model ids in the binary file */
		if( fwrite( modid, 1, nModels*sizeof(int),ioOUT ) - nModels*(int)sizeof(int) )
		{
			fprintf( ioQQQ, " problem trying to write model id grid for costar models\n" );
			fprintf( ioQQQ, " I expected to write %ld words, but fwrite was short\n",
				nModels*sizeof(float) );
			puts( "[Stop in CoStarCompile]" );
			cdEXIT(EXIT_FAILURE);
		}

		/* get the star data */
		for( i=0; i<nModels; ++i )
		{
			/* get number to skip */
			if( fgets( chLine , (int)sizeof(chLine) , ioIN ) == NULL )
			{
				fprintf( ioQQQ, " CoStarCompile fails reading the skip in star data.\n" );
				puts( "[Stop in CoStarCompile]" );
				cdEXIT(EXIT_FAILURE);
			}
			sscanf( chLine , "%li" , &nskip );

			for( j=0; j<nskip; ++j )
			{
				if( fgets( chLine , (int)sizeof(chLine) , ioIN ) == NULL )
				{
					fprintf( ioQQQ, " CoStarCompile fails doing the skip.\n" );
					puts( "[Stop in CoStarCompile]" );
					cdEXIT(EXIT_FAILURE);
				}
			}

			/* now read in the wavelength and flux for this star, read in 
			 * backwards since we want to be in increasing energy order rather
			 * than wavelength */
			for( j=nWL-1; j>=0; --j )
			{
				if( fgets( chLine , (int)sizeof(chLine) , ioIN ) == NULL )
				{
					fprintf( ioQQQ, " CoStarCompile fails reading the star data.\n" );
					puts( "[Stop in CoStarCompile]" );
					cdEXIT(EXIT_FAILURE);
				}
				sscanf( chLine , "%g %g" , 
					&StarEner[j],&StarFlux[j] );

				/* continuum flux was log, convert to linear, also do
				 * conversion from "astrophysical" flux to F_nu in cgs units */
				StarFlux[j] = (float)(PI*pow(10.f,StarFlux[j]));
				/* StarEner was in Angstroms, convert to Ryd */
				StarEner[j] = (float)(RYDLAM/StarEner[j]);
				/*fprintf(ioQQQ,"%g\t%g\n", StarEner[j], StarFlux[j] );*/

				/* sanity check */
				if( j < nWL-1 )
				{
					ASSERT( StarEner[j] < StarEner[j+1] );
				}
			}

			/* this will do the heavy lifting, and define arrays used below,
			 * NB the lowest energy point in these grids appears to be bogus.
			 * tell rebin about nWL-1 */
			RebinAtmosphere(nWL-1, StarEner+1, StarFlux+1, CloudyFlux, 3L, Edges );

			/* this is a debugging option */
			{
				/*@-redef@*/
				enum {DEBUG_LOC=FALSE};
				/*@+redef@*/
				if( DEBUG_LOC )
				{
					FILE *ioBUG;
					if( (ioBUG = fopen("test.txt","w")) !=NULL )
					{
						for( i=0; i < rfield.nupper; ++i )
						{
							fprintf(ioBUG,"%e %e\n", rfield.anu[i], CloudyFlux[i]);
						}
						fclose( ioBUG );
					}
					cdEXIT(EXIT_FAILURE);
				}
			}
			/* write the continuum out as a binary file */
			if( fwrite( CloudyFlux, 1, nBlocksize, ioOUT ) - nBlocksize )
			{
				fprintf( ioQQQ, " problem2 trying to write %s\n",chOutFile );
				fprintf( ioQQQ, " I expected to write %li words, but fwrite was short\n",
					(long)nBlocksize );
				puts( "[Stop in CoStarCompile]" );
				cdEXIT(EXIT_FAILURE);
			}

			fprintf( ioQQQ, " CoStarCompile wrote stellar atmosphere number%3ld\n", i );
		}

		/* up to here */
		/* ======================================================================== */

		fclose( ioIN );
		fclose( ioOUT );

		fprintf( ioQQQ, "\n CoStarCompile completed ok, %s can be deleted.\n\n", chInFile );
		free( StarEner );
		free( StarFlux );
		free( CloudyFlux );
		free( Teff );
		free( modid );
		free( Gravity );
		free( Mass );
		free( Age );
		free( chGrid );
	}

#	ifdef DEBUG_FUN
	fputs( " <->CoStarCompile()\n", debug_fp );
#	endif
	return;
}
/*lint +e713 loss of sign in promotion */
/*lint +e737 loss of sign in promotion */
