/* 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 */
/*RauchInitialize02 K Volk's routine to compact the 66 Rauch atmospheres into a large ascii file.*/
/*RauchInterpolate02 get one of the Rauch model atmospheres, original version by K. Volk */
/*RauchCompile02 rebin Rauch stellar models to match energy grid of code*/
#include "cddefines.h"
#include "physconst.h"
#include "path.h"
#include "called.h"
#include "rfield.h"
#include "stars.h"
/* this is the number of points in each Rauch stellar continuum */
#define	NRAUCH	19951
/* number of atmospheres */
#ifdef NMODS
#undef NMODS
#endif
#define	NMODS	51
/* The number of models in the initial release of the new H-Fe set, June/July 2002 */


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

/* offset due to storage of content id */
static size_t nOffset;

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

/*========================================================================*/
void RauchInitialize02(void)
{
	char chFName[21], 
		chLine[INPUT_LINE_LENGTH],		/* used for getting input lines */
		chFileName[FILENAME_PATH_LENGTH_2] ,	/* file name for star input atmos */
	   chSuff[22];			/* the suffix for the star file */

	/* As the wavelengths are always the same only the fluxes
	 * will be written to the output file.
	 *
	 * Before running this program issue the following command where the Rauch
	 * model atmosphere files are kept (050000_50_solar_iron.bin.0.1 and so on)
	 *
	 *   ls *solar*.1 > rauchmods.list
	 *
	 * and check to see that there are 51 lines in the file.
	 */

	/* concatenated output will go to this file */
	FILE *ioOut,	/* pointer to output file we came here to create*/
		*ioIn ;		/* pointer to input files we will reaed */

	long int i, 
	  j, 
	  k;

	float wl;
	/* this is optional string that will come before star name, so that
	 * this can be executed in the debugger*/
	/*char chPathName[]="C:\\projects\\Cloudy\\temp_star\\solar\\";*/
	char chPathName[]="";

	static char chFile[NMODS][12]=
{"0050000_50","0050000_60","0050000_70","0050000_80",
"0060000_50","0060000_60","0060000_70","0060000_80",
"0070000_50","0070000_60","0070000_70","0070000_80",
"0080000_50","0080000_60","0080000_70","0080000_80",
"0090000_50","0090000_60","0090000_70","0090000_80",
"0100000_50","0100000_60","0100000_70","0100000_80",
"0110000_60","0110000_70","0110000_80",
"0120000_60","0120000_70","0120000_80",
"0130000_60","0130000_70","0130000_80",
"0140000_60","0140000_70","0140000_80",
"0150000_60","0150000_70","0150000_80",
"0160000_60","0160000_70","0160000_80",
"0170000_60","0170000_70","0170000_80",
"0180000_60","0180000_70","0180000_80",
"0190000_60","0190000_70","0190000_80"};

	double *fluxes;

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

	fprintf( ioQQQ, "RauchInitialize02 on the job.  This will take a while.\n" );

	/*grab some space for the fluxes */
	if( (fluxes = (double *)MALLOC( sizeof(double)*NRAUCH)) == NULL )
		BadMalloc();

	/* loop over two types of stellar atmospheres, solar and halo */
	for( k=0; k < 2; k++ )
	{
		if( k == 0 )
		{
			/* name of the output file */
			strcpy( chFName, "rauch02_solar.ascii" );
			strcpy( chSuff, "_solar_iron.bin_0.1" );
			fprintf( ioQQQ, "Doing solar abundance set....\n" );
		}
		else
		{
			/* name of the output file */
			strcpy( chFName, "rauch02_halo.ascii" );
			strcpy( chSuff, "_halo__iron.bin_0.1" );
			fprintf( ioQQQ, "Doing halo abundance set....\n" );
		}
		ioOut = fopen( chFName , "w" );
		/* check that we got it ok */
		if( ioOut == NULL )
		{
			fprintf(ioQQQ,"error opening output file for Rauch atmospheres.\n");
			fprintf(ioQQQ,"This is routine RauchInitialize02.\n");
			puts( "[Stop in RauchInitialize02]" );
			cdEXIT(EXIT_FAILURE);
		}

		for( i=0; i < NMODS; i++ )
		{
			/* must create name of next stellar atmosphere */
			strcpy( chFileName , chPathName );
			strcat( chFileName , chFile[i] );
			strcat( chFileName , chSuff );
			/* now open next stellar atmosphere for reading*/
			ioIn = fopen( chFileName , "r" );
			if( ioIn == NULL )
			{
				fprintf(ioQQQ,"error opening input file for Rauch atmospheres.\n");
				fprintf(ioQQQ,"name was %s\n",chFileName);
				fprintf(ioQQQ,"This is routine RauchInitialize02.\n");
				puts( "[Stop in RauchInitialize02]" );
				cdEXIT(EXIT_FAILURE);
			}
			
			/* print to screen to show progress */
			fprintf( ioQQQ, " star%3li %s\n", i, chFile[i]);

			/* throw away the header information, all lines that start with '*' */
#			if 0
			for( j=0; j < 34; j++ )
			{
				if(fgets( chLine , (int)sizeof(chLine) , ioIn )==NULL )
				{
					fprintf( ioQQQ, " RauchInitialize02 error in header atmosphile file %4ld%4ld\n", 
					  i, j );
					puts( "[Stop in RauchInitialize02]" );
					cdEXIT(EXIT_FAILURE);
				}
			}
#			endif
			/* get first line */
			j = 0;
			if(fgets( chLine , (int)sizeof(chLine) , ioIn )==NULL )
			{
				fprintf( ioQQQ, " RauchInitialize error in header atmosphile file %4ld%4ld\n", 
					i, j );
				puts( "[Stop in RauchInitialize]" );
				cdEXIT(EXIT_FAILURE);
			}
			/* >>chng 02 nov 20, now keep reading them until don't hit the *
			 * since number of comments may change */
			while( chLine[0] == '*' )
			{
				if(fgets( chLine , (int)sizeof(chLine) , ioIn )==NULL )
				{
					fprintf( ioQQQ, " RauchInitialize error in header atmosphile file %4ld%4ld\n", 
						i, j );
					puts( "[Stop in RauchInitialize]" );
					cdEXIT(EXIT_FAILURE);
				}
				++j;
			}

			for( j=0; j < NRAUCH; j++ )
			{
				float ttemp;
				/* get the input line */
				/* >>chng 02 nov 20, don't reread very first line image since we got it above */
				if( j > 0 )
				{
					if(fgets( chLine , (int)sizeof(chLine) , ioIn )==NULL )
					{
						fprintf( ioQQQ, " RauchInitialize02 error in header atmosphere file %4ld%4ld\n", 
						i, j );
						puts( "[Stop in RauchInitialize02]" );
						cdEXIT(EXIT_FAILURE);
					}
				}
				/* scan off wavelength (don't save) and flux)*/
				sscanf( chLine , "%6f%10e", &wl, &ttemp ); 
				fluxes[j] = ttemp; 
			}

			/* finished - close the unit */
			fclose(ioIn);

			/* now write to output file */
			PrintE93( ioOut, fluxes[0] );
			for( j=1; j < NRAUCH; j++ )
			{
				PrintE93( ioOut, fluxes[j] );
				/* want to have 9 numbers per line */
				if(!((j+1)%9))
					fprintf(ioOut,"\n");
			}
			/* need to throw a newline if we did not end on an exact line */
			if(((j+1)%9))
				fprintf(ioOut,"\n");
		}

		fclose(ioOut);
	}
	/* free the space grabbed above */
	free(fluxes);

	fprintf( ioQQQ, " RauchInitialize02 init went OK - now compile the ascii files.\n");
	puts( "[Stop in RauchInitialize02]" );
	cdEXIT(EXIT_FAILURE); 
}

/* RauchCompile02 rebin Rauch stellar models to match energy grid of code
 * limit to number of cells in frequency array
 * NB is this every is changed, also need to change formats in
 * Kevin Volk's stellar atmospheres codes */
/* helper routine located below */

/*========================================================================*/
/* >>chng 01 feb 12, added return value to indicate success (0) or failure (1) */
int RauchCompile02(void)
{
	FILE 
		/* used for input */
		* ioIN ,
		/* used for output */
		* ioOUT ;

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

	/* used for converting F_lambda -> F_nu */
	double convert;

	/* used to save line image */
	char chLine[INPUT_LINE_LENGTH];

	long int i, 
	  imod, 
	  j, 
	  k;

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

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

	nBlocksize = (size_t)rfield.nupper*sizeof(float );
	ASSERT( nBlocksize > 0 );

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

	/* following has extra 8 since we must read last line, which can be partly empty*/
	scratch = (float *)MALLOC( sizeof(float)*(NRAUCH+8) );
	if( scratch == NULL )
	{ 
		printf( " not enough memory to allocate scratch in RauchCompile02\n" );
		puts( "[Stop in RauchCompile02]" );
		cdEXIT(EXIT_FAILURE);
	}

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

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

	/* This subroutine re-bins the Rauch atmosphere fluxes to match the Cloudy
	 * grid.  The wavelengths of the Rauch models are 5.0 to 2000.0 angstroms 
	 * with a step size of 0.1 angstroms, 19951 points each.  Thus the wavelength
	 * grid is not saved in the file of ascii values. */

	/* With some many points I need to average over the Cloudy grid rather than 
	 * interpolating as is the case for the Kurucz or Werner models.  At short 
	 * wavelengths the two grids are of similar spacing. */

	/* nvals is the number of values per Rauch model.  nmods is the total
	 * number of models in the set.*/
	fprintf( ioQQQ, " RauchCompile02 on the job.\n" );
	fprintf( ioQQQ, " There are a total of %i models, and I will need to do this twice.\n\n",NMODS );

	/* 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;

	/* StarFlux is the Rauch stellar atmosphere wavelength grid (used as temp_star storage) */
	for( i=0; i < NRAUCH; i++ )
	{
		/* notice that this is in decreasing wavelength order */
		StarFlux[i] = (float)(2000. - 0.1*(float)(i));
	}

	/* convert continuum over to increasing energy in Ryd */
	for( i=0; i < NRAUCH; i++ )
	{
		/* >>chng 00 oct 16, change to infinite mass Rydberg, by PvH
		   StarEner[i] = (float)(1e10/(StarFlux[i]*10967758.30)); */
		StarEner[i] = (float)(RYDLAM/StarFlux[i]);
	}

	/* this is used to convert the Rauch "astrophysical" fluxes to cgs F_nu units
	 * The values in the file are F_lambda values in erg/s/cm^2/cm. They have to
	 * be multiplied by pi to get the "real" surface flux. This constant is needed
	 * in the conversion to F_nu in erg/s/cm^2/Hz */
	convert = SPEEDLIGHT/POW2(FR1RYD)*PI;

	for( k=1; k <= 2; k++ )
	{
		if( k == 1 )
		{
			/* this is the solar rauch grid */
			/*if( (ioIN = fopen( "rauch02_solar.ascii", "r" ) ) == NULL )*/
			if( (ioIN = fopen( "rauch02_solar.ascii", "r" ) ) == NULL )
			{
				fprintf( ioQQQ, " NewRauchCompile could not find rauch02_solar.ascii.\n" );
				free( CloudyFlux ), 
				free( StarFlux ); 
				free( StarEner );
				free( scratch );
				return(1);
			}
			fprintf( ioQQQ, " NewRauchCompile got rauch02_solar.ascii.\n" );

			 
			/* the frequencies are in increasing order; when the fluxes are 
			 * read in they will have to be put in reverse order since they go from 
			 * 5 to 2000 angstroms. */
			if( (ioOUT = fopen( "rauch02_solar.mod", "wb" ) ) == NULL )
			{
				fprintf( ioQQQ, " NewRauchCompile failed creating rauch02_solar.mod.\n" );
				free( CloudyFlux ), 
				free( StarFlux ); 
				free( StarEner );
				free( scratch );
				return(1);
			}

		}
		else
		{
			/* this is the halo rauch grid */
			if( (ioIN = fopen( "rauch02_halo.ascii", "r" ) ) == NULL )
			{
				fprintf( ioQQQ, " NewRauchCompile could not find rauch02_halo.ascii.\n" );
				free( CloudyFlux ), 
				free( StarFlux ); 
				free( StarEner );
				free( scratch );
				return(1);
			}
			fprintf( ioQQQ, " NewRauchCompile got rauch02_halo.ascii.\n" );

			if( (ioOUT = fopen( "rauch02_halo.mod", "wb" ) ) == NULL )
			{
				free( CloudyFlux ), 
				free( StarFlux ); 
				free( StarEner );
				free( scratch );
				return(1);
			}

		}


		/* >>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);
		}

		/* write out the cloudy energy grid for later sanity checks */

		if( fwrite( rfield.AnuOrg, 1, nBlocksize,ioOUT ) - (size_t)nBlocksize )
		{
			fprintf( ioQQQ, " NewRauchCompile failed writing anu array.\n" );
			free( CloudyFlux ), 
			free( StarFlux ); 
			free( StarEner );
			free( scratch );
			return(1);
		}

		for( imod=0; imod < NMODS; imod++ )
		{
			/* now read the stellar fluxes */
			i = 0;
			while( i < NRAUCH )
			{
				if( fgets( chLine , (int)sizeof(chLine) , ioIN ) == NULL )
				{
					fprintf( ioQQQ, " NewRauchCompile failed reading star flux.\n" );
					puts( "[Stop in newrauchcompile]" );
					cdEXIT(EXIT_FAILURE);
				}
				sscanf( chLine , "%9f %9f %9f %9f %9f %9f %9f %9f %9f" , 
					&scratch[i  ],
					&scratch[i+1],
					&scratch[i+2],
					&scratch[i+3],
					&scratch[i+4],
					&scratch[i+5],
					&scratch[i+6],
					&scratch[i+7],
					&scratch[i+8]);
				/* increment counter, nine numbers per line*/
				i += 9;
			}

			for( j=0; j < NRAUCH; j++ )
			{
				/* The Rauch models go down to 5 Angstroms (about 182 Rydb) but I choose to 
				 * terminate the array at 150 Rydb since the fluxes are extremely small at
				 * the shortest wavelengths -- 17 orders of magnitude below peak for the 
				 * 500000 K model.  This differs a bit from the handling of the Werner models
				 * even though the temperature range is similar.*/

				/* >>chng 00 oct 26, use full range since reason for above not clear, PvH */

				/* flip values around (wavelength to frequency) and convert to f_nu */
				StarFlux[j] = scratch[NRAUCH-1-j]*(float)convert/POW2(StarEner[j]);
			}

			/* the re-binned values are returned in the "CloudyFlux" array */
			RebinAtmosphere(NRAUCH,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.AnuOrg[i], CloudyFlux[i]);
				}
				cdEXIT(EXIT_FAILURE);
			}*/

			/* write the continuum out as a binary file */
			if( fwrite( CloudyFlux,		1,nBlocksize,ioOUT ) - nBlocksize )
			{
				fprintf( ioQQQ, " RauchCompile02 failed writing star flux.\n" );
				free( CloudyFlux ), 
				free( StarFlux ); 
				free( StarEner );
				free( scratch );
				return(1);
			}

			/* write every 10th model so we are clearly still alive */
			if( imod%10 == 0 )
			{
				fprintf( ioQQQ, " NewRauchCompile wrote stellar model number%4ld\n", 
				  imod );
			}
		}

		fclose(ioIN);
		fclose(ioOUT);

		if( k == 1 )
		{
			fprintf( ioQQQ, " NewRauchCompile finished set 1, rauch02_solar.ascii file may be deleted.\n\n" );
		}
		else
		{
			fprintf( ioQQQ, " NewRauchCompile completed ok, rauch02_halo.ascii file may be deleted.\n\n" );
		}
	}

	/* now free up the memory we claimed*/
	free( CloudyFlux ), 
	free( StarFlux ); 
	free( StarEner );
	free( scratch );

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

/*RauchInterpolate02 get one of the Rauch model atmospheres, coded by K. Volk */
/*lint -e713 unsigned to signed */
/*lint -e737 signed to unsigned */

/*========================================================================*/
void RauchInterpolate02(long int *nstar, 
  double temp_star, 
  double alogg, 
  long int iabund)
{
	/* pointer to the rebinned binary grids */
	FILE *ioIN ;
	char chLine[FILENAME_PATH_LENGTH_2];

	long int i, 
	  ipGrav, 
	  ipGravUp,
	  j, 
	  ipTeff ,
	  ipTeffUp,
	  k , 
	  version;

	/* these are used to get the binary files, 
	 * they must not be changed to double without change those files too */
	float
	  *flux1/*[rfield.nupper]*/, 
	  *flux2/*[rfield.nupper]*/,
	  *FluxInterp;

	double fact1, 
	  fact2, 
	  fr1, 
	  fr2, 
	  fr3, 
	  fr4, 
	  lumi,
	  chk,
	  xval;

	static double teff[NMODS]=
{50000.,50000.,50000.,50000.,60000.,60000.,60000.,60000.,
70000.,70000.,70000.,70000.,80000.,80000.,80000.,80000.,
90000.,90000.,90000.,90000.,100000.,100000.,100000.,100000.,
110000.,110000.,110000.,120000.,120000.,120000.,130000.,130000.,
130000.,140000.,140000.,140000.,150000.,150000.,150000.,160000.,
160000.,160000.,170000.,170000.,170000.,180000.,180000.,180000.,
190000.,190000.,190000};

	static double xlogg[NMODS]=
{5.0,6.0,7.0,8.0,5.0,6.0,7.0,8.0,
5.0,6.0,7.0,8.0,5.0,6.0,7.0,8.0,
5.0,6.0,7.0,8.0,5.0,6.0,7.0,8.0,
6.0,7.0,8.0,6.0,7.0,8.0,6.0,7.0,
8.0,6.0,7.0,8.0,6.0,7.0,8.0,6.0,
7.0,8.0,6.0,7.0,8.0,6.0,7.0,8.0,
6.0,7.0,8.0};

#	define NTVAL	15
	static double tval[NTVAL]={50000.,60000.,70000.,80000.,90000.,100000.,
	  110000.,120000.,130000.,140000.,150000.,160000.,170000.,180000.,
	  190000.};

	/* The 0th offset block is the wavelength scale, 
	 * the first model is number 2 in jval, so you must subtract
	 * 2 when using jval as index into teff or xlogg */
	static long jval[NTVAL][5]={
		{2,3,4,5,5},
		{6,7,8,9,9},
		{10,11,12,13,13},
		{14,15,16,17,17},
		{18,19,20,21,21},
		{22,23,24,25,25},
		{26,26,27,28,28},
		{29,29,30,31,31},
		{32,32,33,34,34},
		{35,35,36,37,37},
		{38,38,39,40,40},
		{41,41,42,43,43},
		{44,44,45,46,46},
		{47,47,48,49,49},
		{50,50,51,52,52}
	};

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

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

	/* get one of the Rauch model atmospheres, coded by K. Volk */

	if( lgDataPathSet )
	{
		/* path is set, generate full path name with file */
		strcpy( chLine , chDataPath );
		if( iabund == 0 )
		{
			/* solar abundances, the default */
			strcat( chLine,  "rauch02_solar.mod" );
		}
		else
		{
			/* halo abundances, set with halo keyword */
			strcat( chLine,  "rauch02_halo.mod" );
		}
	}
	else
	{
		/* path not set, look here */
		if( iabund == 0 )
		{
			strcpy( chLine,  "rauch02_solar.mod" );
		}
		else
		{
			strcpy( chLine,  "rauch02_halo.mod" );
		}
	}

	if( (ioIN = fopen( chLine , "rb" )) == NULL )
	{
		/* something went wrong */
		fprintf( ioQQQ, "Error: Rauch stellar atmosphere file rauch02_solar.mod or rauch02_halo.mod not found.\n" );
		fprintf( ioQQQ, " The path I tried: ==%s==\n", chLine );
		fprintf( ioQQQ, " And here comes its hexadecimal representation:\n" );
		for( i=0; i < FILENAME_PATH_LENGTH_2; i++ )
		{
			fprintf( ioQQQ, " '%c'=%#02x", chLine[i], (unsigned int)chLine[i] );
			if( chLine[i] == '\0' ) 
			{
				break;
			}
		}
		fprintf( ioQQQ, "\n" );
		puts( "[Stop in NewRauchInterpolate]" );
		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, " Rauch02Interpolate failed reading nBlocksize.\n" );
		puts( "[Stop in Rauch02Interpolate]" );
		cdEXIT(EXIT_FAILURE);
	}
	ASSERT( nBlocksize > 0 );

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

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

	/* create some space */
	if( (flux1 = (float*)MALLOC(nBlocksize ) ) == NULL )
		BadMalloc();
	if( (flux2 = (float*)MALLOC(nBlocksize ) ) == NULL )
		BadMalloc();
	if( (FluxInterp = (float*)MALLOC(nBlocksize ) ) == NULL )
		BadMalloc();

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

	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++ )
	{
		/* save energy scale for check against code's in conorm
		 * (scale not defined when this routine called */
		rfield.tNuRyd[rfield.nspec][j] = flux1[j];
	}

	/* check whether the models have the correct effective temperature, for debugging only */
	if( FALSE ) 
	{
		for( j=0; j < NMODS; j++ ) 
		{
			if( fread( FluxInterp, 1, nBlocksize, ioIN ) - nBlocksize  ) 
			{
				fprintf( ioQQQ, " problem trying to read rauch model %li\n", j+1 );
				fprintf( ioQQQ, " I expected to read %li words, but fread was short\n",
					 (long)nBlocksize );
				puts( "[Stop in NewRauchInterpolate]" );
				cdEXIT(EXIT_FAILURE);
			}
			lumi = 0.;
			/* rebinned models are in cgs F_nu units */
			for( k=1; k < rfield.nupper; k++ ) {
				lumi += (rfield.tNuRyd[rfield.nspec][k] - rfield.tNuRyd[rfield.nspec][k-1])*
					(FluxInterp[k] + FluxInterp[k-1])/2.;
			}
			/* now convert luminosity to effective temperature */
			chk = pow(lumi*FR1RYD/STEFAN_BOLTZ,0.25);
			/* allow a tolerance of 2% */
			if( fabs((teff[j]-chk)/teff[j]) > 0.02 ) {
				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 );
			}
		}
	}
		
	/* determine which models are to be interpolated */
	ipTeff = -1;
	for( j=0; j < NTVAL-1; j++ )
	{
		if( tval[j] <= temp_star && tval[j+1] > temp_star )
		{
			/* this points to model at or below desired temperature */
			ipTeff = j;
			break;
		}
	}

	/* this happens if above loop never got a hit */
	if( ipTeff<0 )
	{
		/* above loop could not catch highest temperature */
		/*lint -e777 test floats for equal */
		if( tval[NTVAL-1] == temp_star )
		/*lint +e777 test floats for equal */
		{
			/* >>chng 06 feb 08, val of 23 had been copied over from old Rauch file, not
			 * correct for this set.  Bug found by Marcin Hajduk */
			/*ipTeff = 23;*/
			ipTeff = 14;
		}
		else
		{
			fprintf( ioQQQ, 
				" Requested temperature of %.2e is not within the range"
				"%.0f to %.0f\n", 
			  temp_star, tval[0], tval[NTVAL-1] );
			puts( "[Stop in NewRauchInterpolate]" );
			cdEXIT(EXIT_FAILURE);
		}
	}
	/* we now have fully validated pointers to lower and upper limits to
	 * stellar temperature scale */
	ipTeffUp = MIN2( ipTeff+1 , NTVAL-1 );

	/* now generate same types of pointers within gravity scale, 
	 * which has 1 dex steps, log 5 = 5, 6, 7, and 8 for cooler stars,
	 * hottest stars only 8 */
	/* xval = (alogg - 4.0); */
	/* >>>chng 99 apr 28, from above to below, apparently left on Fort scale,
	 * move onto c scale for addressing array */
	xval = (alogg - 5.0);

	ipGrav = (long)(xval);
	ipGrav = MAX2(0, ipGrav);
	ipGrav = MIN2( ipGrav , 4 );
	ipGravUp = MIN2( ipGrav+1, 4);

	/* Interpolation is carried out first in Teff and then in log(g).
	 * Two Rauch models are read in for the lower log(g) value first, which
	 * I denote models A and B with B the higher T model [ipTeff+1 rather than ipTeff].
	 * unless alogg is 8.0 a second pair of Rauch models are read in for the
	 * higher log(g) value, denoted models C and D.  Then the interpolation is
	 *
	 *  RESULT = fr2*(fr3*B+fr4*A)+fr1*(fr3*D+fr4*C)
	 *
	 * when one has 4 models and
	 *
	 * RESULT = (fr3*B+fr4*A)
	 *
	 * when one has just 2 models. */

	/* this sets up interpolation fractions in log g */
	fr1 = xval - (float)(ipGrav);
	fr2 = (1. - fr1);

	/* now set up interplation fractions in Teff */
	if( ipTeffUp == ipTeff )
	{
		fr3 = 0.;
	}
	else
	{
		fr3 = ((temp_star - tval[ipTeff])/(tval[ipTeffUp] - tval[ipTeff]));
	}
	fr4 = 1. - fr3;

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

	if( fr2>1. || fr2 < 0. )
	{
		fprintf( ioQQQ, 
			"The specified surface gravity is outside the range of the Rauch stellar atmospheres.\n" );
		fprintf( ioQQQ, 
			"Please consult Hazy Vol 1 to find the stellar parameters that are included.\nSorry.\n" );
		puts( "[Stop in NewRauchInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

	/* skip over jval-1 stars */
	/* >>chng 01 oct 18, add nOffset */
	if( fseek(ioIN,  (long)((jval[ipTeff][ipGrav]-1)*nBlocksize+nOffset), SEEK_SET )!= 0 )
	{
		fprintf( ioQQQ, " Error1 seeking Rauch atmosphere%4ld\n", 
		  jval[ipTeff][ipGrav] );
		puts( "[Stop in NewRauchInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

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

	if( called.lgTalk )
	{
		/* say we read in the first atmosphere */
		fprintf( ioQQQ, 
			"                       * <<< new Rauch model %3ld read in"
			".  T_eff = %10.2f log(g) = %8.5f  >>>     *\n", 
			/* -1 below to be same as C90 */
		  jval[ipTeff][ipGrav]-1, 
		   teff[jval[ipTeff][ipGrav]-2], 
		  xlogg[jval[ipTeff][ipGrav]-2] );
	}

	/* >>chng 01 oct 18, add nOffset */
	if( fseek(ioIN,  (long)((jval[ipTeffUp][ipGrav]-1)*nBlocksize+nOffset), SEEK_SET )!= 0 )
	{
		fprintf( ioQQQ, " Error2 seeking new Rauch atmosphere%4ld\n", 
		  jval[ipTeffUp][ipGrav] );
		puts( "[Stop in NewRauchInterpolate]" );
		cdEXIT(EXIT_FAILURE);
	}

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

	/* convert to logs since we will interpolate in log flux */
	for( i=0; i<rfield.nupper; ++i )
	{
		flux1[i] = (float)log10( MAX2(1e-37, flux1[i] ) );
		flux2[i] = (float)log10( MAX2(1e-37, flux2[i] ) );
	}

	if( called.lgTalk )
	{
		fprintf( ioQQQ, "                       * <<< new Rauch model %3ld read in.  T_eff = %10.2f log(g) = %8.5f  >>>     *\n", 
		  jval[ipTeffUp][ipGrav]-1, 
		   teff[jval[ipTeffUp][ipGrav]-2], 
		  xlogg[jval[ipTeffUp][ipGrav]-2] );
	}

	/* The following is an approximate scaling to use for the range of 
	 * temperatures above 200000 K where the temperature steps are large 
	 * and thus the interpolations are over large ranges.  For the lower 
	 * temperatures I assume that there is no need for this since the 
	 * step of 10000 K is not so large.
	 *
	 * It should be remembered that this interpolation is not exact, and 
	 * the possible error at high temperatures might be large enough to 
	 * matter.
	 */
	if( temp_star > 200000. )
	{
		fact1 = log10(temp_star/tval[ipTeff])*4.;
		fact2 = log10(temp_star/tval[ipTeffUp])*4.;
	}
	else
	{
		fact1 = 0.;
		fact2 = 0.;
	}

	/* save interpolated continuum here for now, will use again below */
	for( j=0; j < rfield.nupper; j++ )
	{
		FluxInterp[j] = 
			(float)(fr4*(flux1[j] + fact1) + fr3*(flux2[j] + fact2));
	}

	/* ipGrav is offset from lowest gravity, so ranges from 0 through 3 */
	if( ipGrav < 4 )
	{
		/* for log(g) < 8.0 repeat the interpolation at the larger log(g) then
		 * carry out the interpolation in log(g) */
		/* >>chng 01 oct 18, add nOffset */
		if( fseek(ioIN,  (long)((jval[ipTeff][ipGravUp]-1)*nBlocksize+nOffset), SEEK_SET )!= 0 )
		{
			fprintf( ioQQQ, " Error3 seeking new Rauch atmosphere%4ld\n", 
			  jval[ipTeff][ipGravUp] );
			puts( "[Stop in NewRauchInterpolate]" );
			cdEXIT(EXIT_FAILURE);
		}

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

		/* convert to logs since we will interpolate in log flux */
		for( i=0; i<rfield.nupper; ++i )
		{
			flux1[i] = (float)log10( MAX2(1e-37, flux1[i] ) );
		}

		if( called.lgTalk )
		{
			fprintf( ioQQQ, 
				"                       * <<< new Rauch model %3ld"
				" read in.  T_eff = %10.2f log(g) = %8.5f  >>>     *\n", 
			  jval[ipTeff][ipGravUp]-1, 
			   teff[jval[ipTeff][ipGravUp]-2], 
			  xlogg[jval[ipTeff][ipGravUp]-2] );
		}

		/* >>chng 01 oct 18, add nOffset */
		if( fseek(ioIN,  (long)((jval[ipTeffUp][ipGravUp]-1)*nBlocksize+nOffset), SEEK_SET )!= 0 )
		{
			fprintf( ioQQQ, " Error2 seeking new Rauch atmosphere%4ld\n", 
			  jval[ipTeffUp][ipGrav] );
			puts( "[Stop in NewRauchInterpolate]" );
			cdEXIT(EXIT_FAILURE);
		}

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

		/* convert to logs since we will interpolate in log flux */
		for( i=0; i<rfield.nupper; ++i )
		{
			flux2[i] = (float)log10( MAX2(1e-37, flux2[i] ) );
		}

		if( called.lgTalk )
		{
			fprintf( ioQQQ, 
				"                       * <<< new Rauch model %3ld"
				" read in.  T_eff = %10.2f log(g) = %8.5f  >>>     *\n", 
			  jval[ipTeffUp][ipGravUp]-1, 
			   teff[jval[ipTeffUp][ipGravUp]-2], 
			  xlogg[jval[ipTeffUp][ipGravUp]-2] );
		}

		/* at this stage things are still log flux, so interpolate */
		for( j=0; j < rfield.nupper; j++ )
		{
			flux1[j] = (float)((flux1[j] + fact1)*fr4 + (flux2[j] + fact2)*fr3);
		}

		/* now convert back to linear, set to zero if very small */
		for( j=0; j < rfield.nupper; j++ )
		{
			FluxInterp[j] = 
				(float)pow(10.,fr2*FluxInterp[j] + fr1*flux1[j]);
			if( FluxInterp[j] < 1e-37 )
				FluxInterp[j] = 0.;
		}
	}
	else
	{
		/* One gets here if log(g) is 8.0 and there is no higher log(g) pair of
		 * models to interpolate.*/
		for( j=0; j < rfield.nupper; j++ )
		{
			FluxInterp[j] = (float)pow(10.f,FluxInterp[j]);
			if( FluxInterp[j] < 1e-37 )
				FluxInterp[j] = 0.;
		}
	}

	/* 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 += (rfield.tNuRyd[rfield.nspec][k] - rfield.tNuRyd[rfield.nspec][k-1])*
			(FluxInterp[k] + FluxInterp[k-1])/2.;
	}
	/* now convert luminosity to effective temperature */
	chk = pow(lumi*FR1RYD/STEFAN_BOLTZ,0.25);
	/* allow a tolerance of 5% */
	/* one of the Rauch models has a known error of 7%, increase the
	 * threshold for declaring insanity so that this model passes */
	if( fabs((temp_star-chk)/temp_star) > 0.10 && called.lgTalk )
	{
		/* this means that the claimed effective temperature is not equal to the
		 * actual effective temperature.  For this model the difference would be
		 * a logical error */
		fprintf( ioQQQ,
			 "*** PROBLEM , T_eff discrepancy for this model, expected teff %.2f, "
			 "log(g) %.2f, integration yielded teff %.2f, delta %.2f%%\n",
			 temp_star, alogg, chk, (chk-temp_star)/temp_star*100. );
		TotalInsanity();
	}

	/* Note on the interpolation: 26 October 2000 (Peter van Hoof)
	 *
	 * I computed the effective temperature for a random sample of interpolated
	 * atmospheres by integrating the flux as shown above and compared the results
	 * with the expected effective temperature using DELTA = (COMP-EXPEC)/EXPEC.
	 *
	 * I found that the average discrepacy was:
	 *
	 *     DELTA = -0.66% +/- 1.00% (SOLAR models, sample size 5000)
	 *     DELTA = -0.63% +/- 0.92% (HALO models, sample size 5000)
	 *
	 * The most extreme discrepancies for the SOLAR models were
	 *     -3.80% <= DELTA <= 1.68%
	 *
	 * The most negative discrepancies were for T_eff = 260 kK, log(g) = 5 - 7
	 * The most positive discrepancies were for T_eff =  70 kK, log(g) = 6
	 *
	 * The most extreme discrepancies for the HALO models were
	 *     -3.56% <= DELTA <= 6.46%
	 *
	 * The most negative discrepancies were for T_eff = 255 kK, log(g) = 5 - 6.5
	 * The most positive discrepancies were for T_eff = 200 kK, log(g) = 9
	 *
	 * NB - There is a clear problem with the HALO T_eff = 200 kK, log(g) = 9 model
	 *      it may not have converged properly.
	 *
	 * 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 */

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

	*nstar = rfield.nupper ;
	fclose( ioIN );

	free(flux1);
	free(flux2);
	free(FluxInterp);

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

/*lint +e713 unsigned to signed */
/*lint +e737 signed to unsigned */
