/*DoCoStar read in and interpolate on costar grid of windy O atmospheres */
/*CompileCoStar rebin costar stellar atmospheres to match cloudy energy grid, 
 *called by the compile stars command */
/*RebinCostar - called by CompileCoStar to do much of the work */
/*lint -e713 loss of sign in promotion */
/*lint -e737 loss of sign in promotion */
#include <string.h>
#include <float.h>
#include "cddefines.h"
#include "physconst.h"
#include "rfield.h"
#include "called.h"
#include "incidspec.h"
#include "path.h"
#include "varypar.h"
#include "con0.h"
#include "costar.h"

/* internal routines */
void SortCoStar(float[],float[],float[],float[],char[],int[],float[],int[],long int,long int*);
void FindHCoStar(long int,double,int,float[],int[],float[],int[],long int,long int);
void FindVCoStar(double,int,float[],float[],int[],long int);
void ReadCoStar(FILE*,long int,long int,float[],float[],float[],char[],int[]);
void InterpolateCoStar(float[],float,float[],float);
void RebinCostar(long,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


/*DoCoStar read in and interpolate on costar grid of windy O atmospheres */
void DoCoStar(
	/* 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 )
{
	/* will be used for reading in werner.mod */
	FILE *ioIN ;

	char chLine[81];
	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;

	long int j, 
	  ModelBase,
	  nModels,
	  nTracks,
	  res;

	float *Age,
	  *Gravity,
	  *Mass,
	  *Teff,
	  wl[NCELL],
	  *Work,
	  *Work2,
	  Work3[5];

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

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

	/* structure of costar.mod 
	 *
	 * the first 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, " DoCoStar called with insane value for intmode: %d.\n",intmode );
		puts( "[Stop in DoCoStar]" );
		exit(1);
	}
	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);

	/* 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,  "costar.mod" );
	}
	else
	{
		/* path not set, look here */
		strcpy( chLine , "costar.mod" );
	}

	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 DoCoStar]" );
		exit(1);
	}

	/* read in the saved cloudy energy scale so we can confirm this is a good image */
	if( abs( fread( wl, 1, sizeof(wl),ioIN ) - sizeof(wl) ) !=0 )
	{
		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 DoCoStar]" );
		exit(1);
	}

	/* now get number of models */
	if( abs( fread( &nModels, 1, sizeof(nModels),ioIN ) - sizeof(nModels) ) !=0 )
	{
		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 DoCoStar]" );
		exit(1);
	}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

	/* 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 */
	res = fseek(ioIN,0,SEEK_END);
	if( res == 0 )
	{
		res = ftell( ioIN );
		if( (res - ModelBase) != nModels*(long)sizeof(rfield.flux) )
		{
			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*sizeof(rfield.flux),res );
			fprintf( ioQQQ, " please re-compile costar.mod\n" );
			puts( "[Stop in DoCoStar]" );
			exit(1);
		}
	}

	/* 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 DoCoStar]" );
			exit(1);
		}
	}			

	/* 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 DoCoStar\n" );
		puts( "[Stop in DoCoStar]" );
		exit(1);
	}

	if( (iWork2 = (int *)malloc( sizeof(int)*nTracks*2 )) == NULL )
	{ 
		fprintf( ioQQQ, " not enough memory to allocate iWork2 in DoCoStar\n" );
		puts( "[Stop in DoCoStar]" );
		exit(1);
	}

	/* 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",
					       '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 DoCoStar is called in non-optimizing mode,
	 * when optimizing DoCoStar should report back to 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, " requested CoStar model is out of range.\n" );
		puts( "[Stop in DoCoStar]" );
		exit(1);
	}
	assert( iWork3[0] >= 0 && iWork3[0] < nTracks );
	assert( iWork3[1] >= 0 && iWork3[1] < nTracks );
	assert( IWORK2(0,iWork3[0]) >= 0 && IWORK2(0,iWork3[0]) < nModels );
	assert( IWORK2(1,iWork3[0]) >= 0 && IWORK2(1,iWork3[0]) < nModels );
	assert( IWORK2(0,iWork3[1]) >= 0 && IWORK2(0,iWork3[1]) < nModels );
	assert( IWORK2(1,iWork3[1]) >= 0 && IWORK2(1,iWork3[1]) < 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( VaryPar.lgVarOn )
	{
		i2 = (intmode >= 3 ? 3 : 1);
		if( intmode != 4 )
		{
			VaryPar.varang[VaryPar.nparm][0] = +FLT_MAX;
			VaryPar.varang[VaryPar.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));
					VaryPar.varang[VaryPar.nparm][0] =
						MIN2(VaryPar.varang[VaryPar.nparm][0],temp);
					VaryPar.varang[VaryPar.nparm][1] =
						MAX2(VaryPar.varang[VaryPar.nparm][1],temp);
				}
			}
		}
		else
		{
			int ptr0,ptr1;
			ptr0 = IWORK(0,IWORK(1,iWork3[0]));
			ptr1 = IWORK(0,IWORK(1,iWork3[1]));
			VaryPar.varang[VaryPar.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);
			VaryPar.varang[VaryPar.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",VaryPar.varang[VaryPar.nparm][0],VaryPar.varang[VaryPar.nparm][1]);
#		endif

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

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

		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(rfield.flux,frac1,flux1,frac2);
	}
	if( iWork3[1] != iWork3[0] )
	{
		float flux2[NCELL],frac3,frac4;

		ReadCoStar(ioIN,IWORK2(0,iWork3[1]),ModelBase,flux2,Teff,Gravity,chGrid,modid);
		if( IWORK2(1,iWork3[1]) != IWORK2(0,iWork3[1]) )
		{
			float flux3[NCELL],frac5,frac6;

			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(rfield.flux,frac3,flux2,frac4);
	}

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

	/* now stuff what we read into the arrays that will read them later with interpolate */
	for( j=0; j < NCELL ; j++ )
	{
		IncidSpec.tslop[j][IncidSpec.nspec-1] = (float)rfield.flux[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 */
		IncidSpec.tnu[j][IncidSpec.nspec-1] = wl[j];
	}

	*nstar = NCELL;
	fclose( ioIN );
	free( Teff);
	free( modid);
	free( Gravity);
	free( Mass);
	free( Age);
	free( Work);
	free( Work2);
	free( iWork);
	free( iWork2);
	free( chGrid);

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

/* sort CoStar models according to track and index number, i.e. A1, A2, ..., An, B1, ..., Bm, etc... */
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(0,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 */
	long 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) = 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) = 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 */
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 */
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[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[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] = j;
			iWork3[1] = 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] = j;
				iWork3[1] = 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 */
	con0.lgCoStarInterpolationCaution = (iWork3[1]-iWork3[0] > 1);

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

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

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

	/* position file */
	if( fseek(ioIN,ModelBase+mod_num*NCELL*sizeof(float), SEEK_SET )!= 0 )
	{
		fprintf( ioQQQ, " Error seeking costar atmosphere%4ld\n", 
		  mod_num );
		puts( "[Stop in ReadCoStar]" );
		exit(1);
	}

	/* read model */
	if( abs( fread( arr, 1 , NCELL*sizeof(float), ioIN ) - NCELL*sizeof(float)) != 0 )
	{
		fprintf( ioQQQ, " Error trying to read costar atmosphere%4ld\n", 
		  mod_num );
		puts( "[Stop in ReadCoStar]" );
		exit(1);
	}

	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 */
void InterpolateCoStar(
	float flux1[NCELL],
	float fac1,
	float flux2[NCELL],
	float fac2)
{
	long int j;
	double sum;

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

	for( j=0; j < NCELL ; j++ )
	{
		if( flux1[j] > SMALLFLOAT && flux2[j] > SMALLFLOAT )
		{
			sum = log10(flux1[j])*fac1 + log10(flux2[j])*fac2;
			flux1[j] = (float)pow(10., sum );
		}
		else
		{
			flux1[j] = 0.f;
		}
	}

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

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

	/* these will be malloced into large work arrays*/
	float *StarEner , *StarFlux , *CloudyFlux , *Teff , *Gravity , *Mass, *Age;
	/* 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 ;

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

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

	/* malloc some workspace */
	if( (StarEner = (float *)malloc( sizeof(float)*NCELL )) == NULL )
	{ 
		fprintf( ioQQQ, " not enough memory to allocate StarEner in CompileCoStar\n" );
		puts( "[Stop in CompileCoStar]" );
		exit(1);
	}

	if( (StarFlux = (float *)malloc( sizeof(float)*NCELL )) == NULL )
	{ 
		fprintf( ioQQQ, " not enough memory to allocate StarFlux in CompileCoStar\n" );
		puts( "[Stop in CompileCoStar]" );
		exit(1);
	}

	if( (CloudyFlux = (float *)malloc( sizeof(float)*NCELL )) == NULL )
	{ 
		fprintf( ioQQQ, " not enough memory to allocate CloudyFlux in CompileCoStar\n" );
		puts( "[Stop in CompileCoStar]" );
		exit(1);
	}

	/* 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( (ioIN = fopen( "costar.ascii", "r" ) ) == NULL )
	{
		fprintf( ioQQQ, " CompileCoStar fails opening costar.ascii\n" );
		fprintf( ioQQQ, " Where is this file?\n" );
		puts( "[Stop in CompileCoStar]" );
		exit(1);
	}

	/* this will be the file we create, that will be read to compute models, 
	 * open to write binary */
	if( (ioOUT = fopen( "costar.mod", "wb" ) ) == NULL )
	{
		fprintf( ioQQQ, " CompileCoStar fails creating costar.mod\n" );
		fprintf( ioQQQ, " This is impossible??\n" );
		puts( "[Stop in CompileCoStar]" );
		exit(1);
	}
	fprintf( ioQQQ, " CompileCoStar creating costar.mod\n" );

	/* first dump the Cloudy energy scale in the binary file.   This will
	 * be used for sanity checks when the file is read*/
	if( abs( fwrite( rfield.AnuOrg,		1,sizeof(rfield.AnuOrg),ioOUT ) - 
		sizeof(rfield.AnuOrg) ) !=0 )
	{
		fprintf( ioQQQ, " problem trying to write costar.mod\n" );
		fprintf( ioQQQ, " I expected to write %li words, but fwrite was short\n",
			(long)sizeof(rfield.AnuOrg) );
		puts( "[Stop in CompileCoStar]" );
		exit(1);
	}

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

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

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

	/* now get number of models and number of wavelengths */
	if( fgets( chLine , sizeof(chLine) , ioIN ) == NULL )
	{
		fprintf( ioQQQ, " CompileCoStar fails reading energy grid1.\n" );
		puts( "[Stop in CompileCoStar]" );
		exit(-1);
	}
	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, " CompileCoStar scanned off impossible values of nModels or nWL=%li %li\n",
			nModels, nWL );
		puts( "[Stop in CompileCoStar]" );
		exit(1);
	}

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

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

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

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

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

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

	/* get all Teff's for the stars */
	for( i=0; i<nModels; ++i )
	{
		if( fgets( chLine , sizeof(chLine) , ioIN ) == NULL )
		{
			fprintf( ioQQQ, " CompileCoStar fails reading energy grid1.\n" );
			puts( "[Stop in CompileCoStar]" );
			exit(-1);
		}
		/* 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 );
	}

	/* now save the number of models in the binary file */
	if( abs( fwrite( &nModels, 1, sizeof(long),ioOUT ) - 
		sizeof(long) ) !=0 )
	{
		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 CompileCoStar]" );
		exit(1);
	}

	/* now save the grid of temperatures in the binary file */
	if( abs( fwrite( Teff, 1, nModels*sizeof(float),ioOUT ) - 
		nModels*sizeof(float) ) !=0 )
	{
		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 CompileCoStar]" );
		exit(1);
	}

	/* now save the grid of gravities in the binary file */
	if( abs( fwrite( Gravity, 1, nModels*sizeof(float),ioOUT ) - 
		nModels*sizeof(float) ) !=0 )
	{
		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 CompileCoStar]" );
		exit(1);
	}

	/* now save the grid of gravities in the binary file */
	if( abs( fwrite( Mass, 1, nModels*sizeof(float),ioOUT ) - 
		nModels*sizeof(float) ) !=0 )
	{
		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 CompileCoStar]" );
		exit(1);
	}

	/* now save the grid of model ages in the binary file */
	if( abs( fwrite( Age, 1, nModels*sizeof(float),ioOUT ) - 
		nModels*sizeof(float) ) !=0 )
	{
		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 CompileCoStar]" );
		exit(1);
	}

	/* now save the grid of model char letters in the binary file */
	if( abs( fwrite( chGrid, 1, nModels*sizeof(char),ioOUT ) - 
		nModels*sizeof(char) ) !=0 )
	{
		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 CompileCoStar]" );
		exit(1);
	}

	/* now save the grid of model ids in the binary file */
	if( abs( fwrite( modid, 1, nModels*sizeof(int),ioOUT ) - 
		nModels*sizeof(int) ) !=0 )
	{
		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 CompileCoStar]" );
		exit(1);
	}

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

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

		/* 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 , sizeof(chLine) , ioIN ) == NULL )
			{
				fprintf( ioQQQ, " CompileCoStar fails reading the star data.\n" );
				puts( "[Stop in CompileCoStar]" );
				exit(-1);
			}
			sscanf( chLine , "%g %g" , 
				&StarEner[j],&StarFlux[j] );

			/* continuum flux was log, convert to linear, */
			StarFlux[j] = (float)pow(10.,StarFlux[j]);
			/* cannot pass a zero to Volk's routine */
			StarFlux[j] = MAX2(1e-30f,	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] );*/
		}

		/* 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 */
		RebinCostar(nWL-1, StarEner+1, StarFlux+1, CloudyFlux );

		/* this is a debugging option */
		{
			enum {DEBUG = FALSE};
			if( DEBUG )
			{
				FILE *ioBUG;
				if( (ioBUG = fopen("test.txt","w")) !=NULL )
				{
					for( i=0; i<NCELL; ++i)
					{
						fprintf(ioBUG,"%e %e\n", rfield.anu[i], CloudyFlux[i]);
					}
					fclose( ioBUG );
				}
				exit(1);
			}
		}
		/* write the continuum out as a binary file */
		if( abs( fwrite( CloudyFlux, 1, NCELL*sizeof(float),ioOUT ) - 
			NCELL*sizeof(float) ) !=0 )
		{
			fprintf( ioQQQ, " problem2 trying to write costar.mod\n" );
			fprintf( ioQQQ, " I expected to write %li words, but fwrite was short\n",
				NCELL*sizeof(float) );
			puts( "[Stop in CompileCoStar]" );
			exit(1);
		}

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

	/* up to here */
	/* The Klaus Werner model wavelength grid was converted to Rydberg units 
	 * when I originally created the direct access file, but here they are 
	 * converted back to Hz. */

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

	fclose( ioIN );
	fclose( ioOUT );

	fprintf( ioQQQ, " CompileCoStar completed ok, costar.ascii can be deleted.\n" );
	free( StarEner );
	free( StarFlux );
	free( CloudyFlux );
	free( Teff );
	free( modid);
	free( Gravity );
	free( Mass );
	free( Age );
	free( chGrid);

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

/*====================================================================== */
/*RebinCostar rebin Costar atmospheres onto Cloudy grid */
void RebinCostar(
  /* the number of points in the incident continuum*/
  long nCont , 
  /* the energy grid for the stellar continuum in Ryd*/
  float *StarEner , 
  /* the predicted continuum flux on the stellar continuum scale */
  float *StarFlux , 
  /* what we want, the continuum on the cloudy grid */
  float *CloudyFlux )
{
	long int i, 
	  ipHi, 
	  ipLo, 
	  j;
	double EdgeHigh, 
	  EdgeLow, 
	  power, 
	  slope2, 
	  sum, 
	  v1, 
	  val, 
	  x1, 
	  x2;

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

	/* this is the Rayleigh-Jeans slope in F_nu, for low frequency extrapolation */
	slope2 = 2.0e00;

	for( j=0; j < NCELL; j++ )
	{
		/* EdgeHigh is upper bound of this continnum cell */
		EdgeHigh = rfield.anu[j] + 0.5e00*rfield.widflx[j];

		/* EdgeLow is lower bound of this continnum cell */
		EdgeLow = rfield.anu[j] - 0.5e00*rfield.widflx[j];

		if( EdgeHigh < StarEner[0] )
		{
			/* this is case where cloudy's continuum is below stellar continuum,
			 * so we do Rayleigh Jeans extrapolation */
			CloudyFlux[j] = (float)(StarFlux[0]*(pow(rfield.anu[j]/
			  StarEner[0],slope2)));
		}
		else
		{
			if( EdgeLow > StarEner[nCont-1] )
			{
				/* case where cloudy continuum above highest stellar point */
				CloudyFlux[j] = 0.0e00;
			}
			else
			{
				/* now go through stellar continuum to find bins corresponding to
				 * this cloudy cell, stellar continuum defined through nCont cells */
				ipLo = 0;
				ipHi = 0;
				for( i=0; i < nCont; i++ )
				{
					if( StarEner[i+1] > EdgeLow && StarEner[i] <= EdgeLow  )
					{
						ipLo = i;
						break;
					}
				}
				for( i=MAX2(0,ipLo-1); i < nCont; i++ )
				{
					if( StarEner[i+1] > EdgeHigh && StarEner[i] <= EdgeHigh )
					{
						ipHi = i;
						break;
					}
				}

				/* Do the case where the grid point and its edges are between two of the
				 * stellar model points:  add segments with power-law interpolation up to
				 * do the averaging.*/

				/* When EdgeLow < anu3(1) ipLo is equal to -1 and ipHi is equal to 1, so do it 
				 * as if ipHi = ipLo to avoid going outside array bounds. */
				if( ipHi == ipLo || ipLo>ipHi)
				{
					power = log(StarFlux[ipLo+1]/StarFlux[ipLo])/log(StarEner[ipLo+1]/
					  StarEner[ipLo]);

					val = power*log(rfield.anu[j]/StarEner[ipLo]);

					CloudyFlux[j] = (float)(StarFlux[ipLo]*exp(val));
				}
				else
				{
					sum = 0.0e00;

					/* ipHi points to stellar point at high end of cloudy continuum cell,
					 * ipLo points to low end */
					for( i=ipLo; i <= ipHi; i++ )
					{
						power = log(StarFlux[i+1]/StarFlux[i])/
						  log(StarEner[i+1]/StarEner[i]);

						if( i == ipLo )
						{
							x1 = EdgeLow;
							x2 = StarEner[i+1];
							v1 = StarFlux[i]*exp(power*log(x1/StarEner[i]));
							/*v2 = StarFlux[i+1];*/
						}

						else if( i == ipHi )
						{
							x2 = EdgeHigh;
							x1 = StarEner[i];
							/*v2 = StarFlux[i]*exp(power*log(x1/StarEner[i]));*/
							v1 = StarFlux[i];
						}

						/*if( i > ipLo && i < ipHi )*/
						else
						{
							x1 = StarEner[i];
							x2 = StarEner[i+1];
							v1 = StarFlux[i];
							/*v2 = StarFlux[i+1];*/
						}

						if( fabs(power+1.0) < 0.001 )
						{
							val = x1*v1*log(x2/x1);
						}
						else
						{
							val = pow(x2/x1,power + 1.e00) - 1.e00;
							val = val*x1*v1/(power + 1.e00);
						}
						sum += val;
					}

					CloudyFlux[j] = (float)(sum/rfield.widflx[j]);
				}
			}
		}
	}

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