/* 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 */
/*WernerCompile rebin Werner stellar atmospheres to match cloudy energy grid */
/*WernerInterpolate_one get a single Werner PN atmosphere, by K Volk */
/*WernerInterpolate read in and interpolate on Werner grid of PN atmospheres, by K Volk */
#include "cddefines.h"
#include "rfield.h"
#include "path.h"
#include "called.h"
#include "physconst.h"
#include "stars.h"
/* number of models in the stream */
#	ifdef NMODS
#undef NMODS
#	endif
#define	NMODS	20
/* number of frequency points in Werner model,
 * will allocate 4 more than this is MALLOC */
#ifdef NVALS
#	undef NVALS
#endif
#define	NVALS 514

static size_t nOffset;
static size_t nBlocksize;

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

/*WernerInterpolate_one get a single Werner PN atmosphere, by K Volk */
static void WernerInterpolate_one(double temp_star, 
  double alogg, 
  long int *ierr, 
  FILE * ioIN );

/*WernerCompile rebin Werner stellar atmospheres to match cloudy energy grid */
void WernerCompile(void)
{
	char chLine[INPUT_LINE_LENGTH];
	long int i, 
	  imod,
	  line,
	  nAllocated;

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

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

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

	fprintf( ioQQQ, " WernerCompile 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] = 0.99946789f;
	Edges[1] = 1.8071406f;
	Edges[2] = 3.9996377f;

	/* this is size of vectors and will be checked in asserts for array overwrite */
	nAllocated = NVALS + 4;
	/* MALLOC some workspace */
	StarEner = (float *)MALLOC( sizeof(float)*(unsigned)nAllocated );
	if( StarEner == NULL )
	{ 
		printf( " not enough memory to allocate StarEner in WernerCompile\n" );
		puts( "[Stop in WernerCompile]" );
		cdEXIT(EXIT_FAILURE);
	}

	StarFlux = (float *)MALLOC( sizeof(float)*(unsigned)nAllocated );
	if( StarFlux == NULL )
	{ 
		printf( " not enough memory to allocate StarFlux in WernerCompile\n" );
		puts( "[Stop in WernerCompile]" );
		cdEXIT(EXIT_FAILURE);
	}

	CloudyFlux = (float *)MALLOC( sizeof(float)*(unsigned)(rfield.nupper+1) );
	if( CloudyFlux == NULL )
	{ 
		printf( " not enough memory to allocate CloudyFlux in WernerCompile\n" );
		puts( "[Stop in WernerCompile]" );
		cdEXIT(EXIT_FAILURE);
	}


	/* This is a program to re-bin the Werner stellar model spectra to match the 
	 * CLOUDY grid.  For short wavelengths I will assume zero flux. At long
	 * wavelengths I will assume Rayleigh-Jeans from the last stellar model point.
	 * This version uses power-law interpolation between the points of the stellar
	 * model. */

	/* This "werner.ascii" file is a straight ascii dump of the Klaus Werner 
	 * stellar model files which he gave to me in 1992.  The first set of values 
	 * is the frequency grid (in Ryd) followed by the atmosphere models in order
	 * of increasing temperature and log(g). */

	if( (ioIN = fopen( "werner.ascii", "r" ) ) == NULL )
	{
		fprintf( ioQQQ, " WernerCompile fails opening werner.ascii\n" );
		fprintf( ioQQQ, " Where is this file?\n" );
		puts( "[Stop in getwerner]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( (ioOUT = fopen( "werner.mod", "wb" ) ) == NULL )
	{
		fprintf( ioQQQ, " WernerCompile fails creating werner.mod\n" );
		fprintf( ioQQQ, " This is impossible??\n" );
		puts( "[Stop in getwerner]" );
		cdEXIT(EXIT_FAILURE);
	}

	fprintf( ioQQQ, " WernerCompile reading werner.ascii and creating werner.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);
	}

	/* first 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, " problem1 trying to write werner.mod\n" );
		fprintf( ioQQQ, " I expected to write %li words, but fwrite was short\n",
			(long)nBlocksize );
		puts( "[Stop in WernerCompile]" );
		cdEXIT(EXIT_FAILURE);
	}

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

	/* get energy grid, */
	i = 0;
	while( i < NVALS )
	{
		if( fgets( chLine , (int)sizeof(chLine) , ioIN ) == NULL )
		{
			fprintf( ioQQQ, " WernerCompile fails reading energy grid1.\n" );
			puts( "[Stop in GetData]" );
			cdEXIT(EXIT_FAILURE);
		}

		ASSERT( (i+4) < (nAllocated-1) );
		sscanf( chLine , "%f %f %f %f %f" , 
			&StarEner[i],&StarEner[i+1],&StarEner[i+2],&StarEner[i+3],&StarEner[i+4] );

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

	/* >>chng 00 oct 20, StarEner[172] is out of sequence. As per the Klaus Werner comment,
	 * it will be omitted. The energy grid is very dense in this region and was most likely
	 * intended to sample an absorption line which was not included in this particular grid.
	 * StarFlux[172] is therefore always equal to the flux in neighbouring points (at least
	 * those with slightly smaller energies). It is therefore safe to ignore this point. PvH */
	for( i=172; i < NVALS-1; i++ )
	{
		StarEner[i] = StarEner[i+1];
	}

	/* read in the store stellar continuum energy array */
	for( i=1; i < NVALS-1; ++i )
	{
		/* following is due to there being fluxes above and below certain 
		 * wavelengths where  the opacity changes 
		 * (i.e. the Lyman and Balmer limits for example) which are 
		 * assigned the same wavelength in the Klaus Werner files. */
		if( StarEner[i] == StarEner[i-1] )
		{
			float val1,val2;

			/* >>chng 00 oct 18, take care that StarEner remains monotonic, PvH */
			val1 = 0.99999f*StarEner[i-1];
			val2 = (i > 1) ? StarEner[i-1] - 0.1f*(StarEner[i-1]-StarEner[i-2]) : 0.f;
			StarEner[i-1] = MAX2(val1,val2);
			val1 = 1.00001f*StarEner[i];
			val2 = (i < NVALS-1) ? StarEner[i] + 0.1f*(StarEner[i+1]-StarEner[i]) : FLT_MAX;
			StarEner[i] = MIN2(val1,val2);
		}

		/* sanity check */
		ASSERT( StarEner[i] > StarEner[i-1] );
	}

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

	line = 0;
	for( imod=0; imod < NMODS; imod++ )
	{
		/* now get the models */
		i = 0;
		while( i < NVALS )
		{
			if( fgets( chLine , (int)sizeof(chLine) , ioIN ) == NULL )
			{
				fprintf( ioQQQ, " WernerCompile fails reading star file1.\n" );
				puts( "[Stop in GetData]" );
				cdEXIT(EXIT_FAILURE);
			}
			++line;
			{
				/* following should be set true to print input lines */
				/*@-redef@*/
				enum {DEBUG_LOC=FALSE};
				/*@+redef@*/
				if( DEBUG_LOC )
				{
					fprintf(ioQQQ,"%li %li %s", imod, line,chLine );
					fflush(ioQQQ);
				}
			}

			ASSERT( (i+4) < (nAllocated-1) );
			sscanf( chLine , "%f %f %f %f %f" , 
				&StarFlux[i],&StarFlux[i+1],&StarFlux[i+2],
				&StarFlux[i+3],&StarFlux[i+4] );
			/* increment counter, five numbers per line*/
			i += 5;
		}

		/* >>chng 00 oct 20, As per the Klaus Werner comment,
		 * StarFlux[172] is deleted. See also the comment above. PvH */
		for( i=172; i < NVALS-1; i++ )
		{
			StarFlux[i] = StarFlux[i+1];
		}

		/* continuum flux was log, convert to linear, 
		 * and make sure no doubled energy points */
		for( i=0; i < NVALS-1; i++ )
		{
			StarFlux[i] = (float)pow(10.f,StarFlux[i]);
		}

		/* this will do the heavy lifting, and define arrays used below
		 * NVALS is the number of energy points, minus 1 because StarEner[172] was deleted */
		RebinAtmosphere(NVALS-1, StarEner, StarFlux, CloudyFlux, 3L, Edges );

		/*{
			FILE *ioBUG;
			ioBUG = fopen("test.txt","w");
			for( i=0; i<rfield.nupper; ++i)
			{
				fprintf(ioBUG,"%e %e\n", rfield.anu[i], CloudyFlux[i]);
			}
			cdEXIT(EXIT_FAILURE);
		}*/
		/* write the continuum out as a binary file */
		if( fwrite( CloudyFlux, 1, nBlocksize, ioOUT ) - nBlocksize )
		{
			fprintf( ioQQQ, " problem2 trying to write werner.mod\n" );
			fprintf( ioQQQ, " I expected to write %li words, but fwrite was short\n",
				(long)nBlocksize );
			puts( "[Stop in WernerCompile]" );
			cdEXIT(EXIT_FAILURE);
		}

		fprintf( ioQQQ, " WernerCompile wrote stellar atmosphere number%3ld\n", 
		  imod+1 );
	}

	fclose( ioIN );
	fclose( ioOUT );

	fprintf( ioQQQ, "\n WernerCompile completed ok, werner.ascii can be deleted.\n\n" );
	free( StarEner );
	free( StarFlux );
	free( CloudyFlux );

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

static float *scratch;

/*WernerInterpolate read in and interpolate on Werner grid of PN atmospheres, by K Volk */
void WernerInterpolate(long int *nstar, 
  double temp_star, 
  double alogg)
{
	/* will be used for reading in werner.mod */
	FILE *ioIN ;

	char chLine[FILENAME_PATH_LENGTH_2];

	int lgAgain;

	long int ierr, 
	  j, 
	  jmod,
	  k , 
	  version;

	double dev1, 
	  dev2,
	  lumi,
	  chk;

	float *wl/*[rfield.nupper]*/;

	static float teff[NMODS]={80000.,80000.,80000.,80000.,100000.,
	  100000.,100000.,100000.,120000.,120000.,120000.,140000.,140000.,
	  140000.,160000.,160000.,180000.,180000.,200000.,200000.};

	static float xlogg[NMODS]={5.,6.,7.,8.,5.,6.,7.,8.,6.,7.,8.,6.,
	  7.,8.,7.,8.,7.,8.,7.,8.};

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

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

	/* This subroutine was added (28 dec 1992) to read from the set of
	 * hot white dwarf model atmospheres from Klaus Werner at Kiel. The 
	 * values are read in (energy in rydberg units, f_nu in cgs units)
	 * for any of NMODS=20 models.  T_eff values of 80, 100, 120, 140, 160,
	 * 180 and 200 thousand K and log(g) values of 5.0, 6.0, 7.0, and 8.0 are
	 * available.  This routine does the interpolation.  If a model has a T_eff
	 * and a log(g) value within 10 epsilon of the values specified in the
	 * parameters, that model is chosen for use.  Each model had NVALS=514
	 * points before rebinning. The Rayleigh-Jeans tail was extrapolated.
	 *
	 * file used: werner.mod (contains the wavelength array, then the models)
	 *
	 * NMODS is the number of accessable models in the werner model set.
	 * Data set 0 is the wavelength grid, so there are NMODS+1 data sets
	 * in the main file.  Each set has rfield.nupper re-binned values, increasing in 
	 * frequency.  The original models had NVALS=514 values.
	 * The file is an unformatted direct access file. */

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

	if( (ioIN = fopen( chLine , "rb" )) == NULL )
	{
		/* something went wrong */
		fprintf( ioQQQ, "ERROR: The Klaus Werner stellar atmosphere file was not found.\nSorry.\n" );
		puts( "[Stop in WernerInterpolate]" );
		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, " WernerInterpolate failed reading nBlocksize.\n" );
		puts( "[Stop in WernerInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}
	if( fread( &rfield.nupper,	1, sizeof(rfield.nupper) ,ioIN ) - sizeof(rfield.nupper) )
	{
		fprintf( ioQQQ, " WernerInterpolate failed reading nBlocksize.\n" );
		puts( "[Stop in WernerInterpolate]" );
		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, " WernerInterpolate failed reading VERSION.\n" );
		puts( "[Stop in WernerInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}
	/* this had better be zero */
	if( VERSION-version )
	{
		fprintf( ioQQQ, " WernerInterpolate: there is a version mismatch between the compiled werner atmospheres I expected and the one I found.\n" );
		fprintf( ioQQQ, " WernerInterpolate: Please recompile the werner stellar atmospheres with the command COMPILE STARS ONLY WERNER.\n" );
		puts( "[Stop in WernerInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* true is special path set, false if file in same dir */
	if( (wl = (float*)MALLOC(nBlocksize ) ) == NULL )
		BadMalloc();
	if( (scratch = (float*)MALLOC(nBlocksize ) ) == NULL )
		BadMalloc();

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

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

	/* check whether the models have the correct effective temperature, for debugging only */
	if( FALSE ) 
	{
		for( j=0; j < NMODS; j++ ) 
		{
			if( fread( scratch, 1, nBlocksize, ioIN ) - nBlocksize ) 
			{
				fprintf( ioQQQ, " problem trying to read werner model %li\n", j+1 );
				fprintf( ioQQQ, " I expected to read %li words, but fread was short\n",
					 (long)nBlocksize );
				puts( "[Stop in WernerInterpolate]" );
				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)/chk) > 0.005 ) 
			{
				printf( "*** WARNING, T_eff discrepancy for model %li, expected teff %.2f, ",
					j+1, teff[j]);
				printf("log(g) %.2f, integration yielded teff %.2f\n",
					xlogg[j], chk );
			}
		}
	}
		
	/* TEMP and ALOGG are desired temp_star and log g */

	/* this may (or may not) catch bad input parameters */
	if( temp_star < 80000. || temp_star > 200000. )
	{
		fprintf( ioQQQ, " Entered temperature is outside range of grid.\n" );
		puts( "[Stop in WernerInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( alogg < 5. || alogg > 8. )
	{
		fprintf( ioQQQ, " Entered gravity is outside range of grid.\n" );
		puts( "[Stop in WernerInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	lgAgain = TRUE;
	j = 0;
	/* -1 will be indication that interpolation is needed */
	jmod = -1;

	while( lgAgain && j < NMODS )
	{
		dev1 = fabs(temp_star/teff[j]-1.);
		dev2 = fabs(alogg-xlogg[j]);

		/* this is a test for how close we are to the tabulated models */
		if( dev1 <= 10.*FLT_EPSILON && dev2 <= 10.*FLT_EPSILON )
		{
			/* jmod is C-scale pointer into teff and xlogg arrays */
			jmod = j;
			lgAgain = FALSE;
		}
		++j;
	}

	/* if jmod still -1 then we need to interpolate */
	if( jmod == -1 )
	{
		/* was called getmod but name changed to following */
		WernerInterpolate_one( temp_star, alogg, &ierr, ioIN );
		if( ierr > 0 )
		{
			fprintf( ioQQQ, " Error occurred in WernerInterpolate.\n" );
			puts( "[Stop in WernerInterpolate]" );
			cdEXIT(EXIT_FAILURE);
		}
	}
	else
	{
		/* this branch, we hit a proper model, use exactly that one */
		/* the jmod+1 below is to offset the fseek by one; this is to skip
		 * the set of energies, which are stored in the first record but not
		 * included in the list of models */
		/* >>chng 01 oct 18, added nOffset */
		if( fseek( ioIN, (long)((jmod+1)*nBlocksize+nOffset) , SEEK_SET )!= 0 )
		{
			fprintf( ioQQQ, " Error seeking Werner atmosphere%4ld\n", 
			  jmod );
			puts( "[Stop in WernerInterpolate]" );
			cdEXIT(EXIT_FAILURE);
		}

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

		if( called.lgTalk )
		{
			fprintf( ioQQQ, 
				 "                       * <<< Klaus Werner stellar model %2ld read in. T_eff = %8.1f log(g) = %6.3f >>> *\n", 
			  jmod+1, teff[jmod], xlogg[jmod] );
		}
	}

	/* 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((temp_star-chk)/temp_star) > 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",
			 temp_star, alogg, chk, (chk-temp_star)/temp_star*100. );
		/* TotalInsanity exits after announcing a problem */
		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.71% +/- 0.71% (sample size 5000)
	 *
	 * The most extreme discrepancies were
	 *     -4.37% <= DELTA <= 0.24%
	 *
	 * The most negative discrepancies were for T_eff =  95 kK, log(g) = 5
	 * The most positive discrepancies were for T_eff = 160 kK, log(g) = 8
	 *
	 * 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 */
	rfield.lgContMalloc[rfield.nspec] = TRUE;
	if( (rfield.tNuRyd[rfield.nspec] = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.tslop[rfield.nspec] = (float*)MALLOC((size_t)(rfield.nupper*sizeof(float)) ) ) == NULL )
		BadMalloc();
	if( (rfield.tFluxLog[rfield.nspec] = (float*)MALLOC((size_t)(rfield.nupper*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
		 * the 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( wl );
	free( scratch );

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

/*WernerInterpolate_one interpolate T and G on the Werner PN atmosphere, by K Volk */
static void WernerInterpolate_one(
						  /* desired stellar temperature */
						  double temp_star, 
						  /* desired surface gravity */
						  double alogg, 
						  /* this will be zero if ended ok */
						  long int *ierr, 
						  /* the file with the stellar continuum */
						  FILE * ioIN )
{
	long int 
		ipGrav, 
		ipGravUp ,
		j, 
		ipTeff,
		ipTeffUp;

	/* following two used in binary read, must remain float 
	 * or binary file changed with it */
	float *flux1/*[rfield.nupper]*/, 
	  *flux2/*[rfield.nupper]*/; 

	double fr1,
	  fr2, 
	  fr3, 
	  fr4;

	/* offset for fseek to first star file */
	const long iOff=-1;

	/* the computed gravities */
	static float gval[4]={5.,6.,7.,8.};

	/* the computed temperatures*/
	static float tval[7]=
	{80000.,100000.,120000.,140000.,160000.,180000., 200000.};

	/* these are array indices for specific models, organized by grav and teff */
	static long modnum[4][7]=
	{ 
	{2,6,10,13,16,18,20},
	{3,7,10,13,16,18,20},
	{4,8,11,14,16,18,20},
	{5,9,12,15,17,19,21}
	};

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

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

	/* this is error indicator, error condition detected if >0 upon exit */
	*ierr = 1;

	/* this subroutine gets stellar atmospheres from the Werner file
	 * and interpolates on them, returning the flux in FLUX as F_nu in cgs units */

	/* pointer to temp_star will range from 0 to 6, there are seven temps stored */
	ipTeff = (long)((temp_star-80000.)/20000.);
	ipTeff = MAX2(0, ipTeff );
	ipTeff = MIN2( ipTeff , 6);

	/* upper model, make sure grid not overrun */
	ipTeffUp = MIN2( 6, ipTeff + 1 );

	/* check whether we can deal with the surface gravity */
	/* tabulated alogg is 5 to 8, 
	 * pointer will range from 0 to 3 */
	ipGrav = (long)(alogg) - 5;
	ipGrav = MAX2( 0 , ipGrav );
	ipGrav = MIN2( 3 , ipGrav );

	/* upper model, make sure grid not overrun */
	ipGravUp = MIN2( 3, ipGrav + 1 );

	/* at this point we should have valid gravity (ipGrav) and
	 * teff (ipTeff) pointers */

	/* Read in first model ==================== */

	/* >>chng 01 oct 18, added nOffset */
	if( fseek(ioIN,  (long)((modnum[ipGrav][ipTeff]+iOff)*nBlocksize+nOffset), SEEK_SET ) )
	{
		fprintf( ioQQQ, " Error seeking Werner atmosphere%4ld\n", 
		  modnum[ipGrav][ipTeff] );
		puts( "[Stop in WernerInterpolate_one]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( fread( flux1, 1 , nBlocksize, ioIN ) -	nBlocksize )
	{
		fprintf( ioQQQ, " Error trying to read Werner atmosphere%4ld\n", 
		  modnum[ipGrav][ipTeff] );
		puts( "[Stop in WernerInterpolate_one]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* say what we got */
	if( called.lgTalk )
	{
		fprintf( ioQQQ, 
		  "                       * <<< Klaus Werner stellar model %2ld read in. T_eff = %8.1f log(g) = %6.3f >>> *\n", 
		  modnum[ipGrav][ipTeff]-1, tval[ipTeff], gval[ipGrav] );
	}

	/* >>chng 01 oct 18, added nOffset */
	if( fseek(ioIN,  (long)((modnum[ipGravUp][ipTeff]+iOff)*nBlocksize+nOffset), SEEK_SET ) )
	{
		fprintf( ioQQQ, " Error seeking Werner atmosphere%4ld\n", 
		  modnum[ipGravUp][ipTeff] );
		puts( "[Stop in WernerInterpolate_one]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* Read in second model ==================== */

	if( fread( flux2, 1 , nBlocksize, ioIN ) - nBlocksize )
	{
		fprintf( ioQQQ, " Error2 trying to read Werner atmosphere%4ld\n", 
		  modnum[ipGravUp][ipTeff] );
		puts( "[Stop in WernerInterpolate_one]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( called.lgTalk )
	{
		fprintf( ioQQQ,  
			"                       * <<< Klaus Werner stellar model %2ld read in."
			" T_eff = %8.1f log(g) = %6.3f >>> *\n", 
		  modnum[ipGravUp][ipTeff]-1 , tval[ipTeff], gval[ipGravUp] );
	}

	fr1 = alogg - gval[ipGrav];
	fr2 = 1. - fr1;
	if( ipTeffUp == ipTeff )
	{
		fr3 = 1.;
	}
	else
	{
		fr3 = (temp_star - tval[ipTeff])/(tval[ipTeffUp] - tval[ipTeff]);
	}
	fr4 = 1. - fr3;

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

	if( fr2 < 0. || fr2 > 1. )
	{
		fprintf( ioQQQ, "fr2 insanity in WernerInterpolate_one\n" );
		ShowMe();
		puts( "[Stop in WernerInterpolate_one]" );
		cdEXIT(EXIT_FAILURE);
	}
	/* save inerpolated continuum in flux array */
	for( j=0; j < rfield.nupper; j++ )
	{
		/* this is the interpolation between the two read in so far */
		scratch[j] = (float)(
			fr2*log10(MAX2(1e-37,flux1[j])) + 
			fr1*log10(MAX2(1e-37,flux2[j])));
	}
	

	/* Read in third model ==================== */

	/* >>chng 01 oct 18, added nOffset */
	if( fseek(ioIN,  (long)((modnum[ipGrav][ipTeffUp]+iOff)*nBlocksize+nOffset), SEEK_SET ) )
	{
		fprintf( ioQQQ, " Error seeking Werner atmosphere%4ld\n", 
		  modnum[ipGrav][ipTeffUp] );
		puts( "[Stop in WernerInterpolate_one]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( ( fread( flux1, 1 , nBlocksize, ioIN ) -
		nBlocksize) != 0 )
	{
		fprintf( ioQQQ, " Error trying to read Werner atmosphere%4ld\n", 
		  modnum[ipGrav][ipTeffUp] );
		puts( "[Stop in WernerInterpolate_one]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( called.lgTalk )
	{
		fprintf( ioQQQ, 
		  "                       * <<< Klaus Werner stellar model %2ld read in. T_eff = %8.1f log(g) = %6.3f >>> *\n", 
		  modnum[ipGrav][ipTeffUp]-1 , tval[ipTeffUp], gval[ipGrav] );
	}

	/* Read in fourth model ==================== */

	/* >>chng 01 oct 18, added nOffset */
	if( fseek(ioIN,  (long)((modnum[ipGravUp][ipTeffUp]+iOff)*nBlocksize+nOffset), SEEK_SET ) )
	{
		fprintf( ioQQQ, " Error seeking Werner atmosphere%4ld\n", 
		  modnum[ipGravUp][ipTeffUp] );
		puts( "[Stop in WernerInterpolate_one]" );
		cdEXIT(EXIT_FAILURE);
	}

	if( ( fread( flux2, 1 , nBlocksize, ioIN ) -
		nBlocksize) != 0 )
	{
		fprintf( ioQQQ, " Error3 trying to read Werner atmosphere%4ld\n", 
		  modnum[ipGravUp][ipTeffUp] );
		puts( "[Stop in WernerInterpolate_one]" );
		cdEXIT(EXIT_FAILURE);
	}


	if( called.lgTalk )
	{
		fprintf( ioQQQ, 
		  "                       * <<< Klaus Werner stellar model %2ld read in. T_eff = %8.1f log(g) = %6.3f >>> *\n", 
		  modnum[ipGravUp][ipTeffUp]-1, tval[ipTeffUp], gval[ipGravUp] );
	}

	/* now do the actual interpolation */
	for( j=0; j < rfield.nupper; j++ )
	{
		flux1[j] = (float)(
			fr2*log10(MAX2(1e-37,flux1[j])) + 
			fr1*log10(MAX2(1e-37,flux2[j])));
	}

	for( j=0; j < rfield.nupper; j++ )
	{
		scratch[j] = (float)pow(10.,(fr4*scratch[j] + fr3*flux1[j]));
		if( scratch[j]<1e-36 )
			scratch[j] = 0.;
	}

	/* say that things turned out ok */
	*ierr = 0;

	free( flux1 );
	free( flux2);

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

