/* 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 */
#include "cddefines.h"
#include "cddrive.h"
/* #include "optimize.h" */
#include "grid.h"
#include "punch.h"
#include "rfield.h"
#include "prt.h"
#include "input.h"
#include "version.h"

#define RECORDSIZE	2880
#define LINESIZE	80

#if defined(_BIG_ENDIAN) 
	/* the value of A will not be manipulated */
#	define HtoNL(A) (A)	
/*
#	define HtoNS(A) (A)
#	define NtoHS(A) (A)
#	define NtoHL(A) (A)
*/
#else
/* defined(_LITTLE_ENDIAN) */
/* the value of A will be byte swapped */
#	define HtoNL(A) ((((A) & 0xff000000) >> 24) | \
		(((A) & 0x00ff0000) >> 8) | \
		(((A) & 0x0000ff00) << 8) | \
		(((A) & 0x000000ff) << 24))
/*
#	define HtoNS(A) ((((A) & 0xff00) >> 8) | (((A) & 0x00ff) << 8))
#	define NtoHS HtoNS
#	define NtoHL HtoNL
*/
/*#else
error One of BIG_ENDIAN or LITTLE_ENDIAN must be #defined.*/
#endif

#define ByteSwap5(x) ByteSwap((unsigned char *) &x,sizeof(x))

void ByteSwap(unsigned char * b, int n)
{
	register int i = 0;
	register int j = n-1;
	while (i<j)
	{
		char temp = b[i];
		b[i] = b[j];
		b[j] = temp;
		/* std::swap(b[i], b[j]); */
		i++, j--;
	}
	return;
}

static FILE *fits_output;
static long bytesAdded = 0;
static long bitpix = 8;
static long pcount = 0;
static long gcount = 1;

void punchFITS_PrimaryHeader( int lgAddModel );
void punchFITS_ParamHeader( long numParamValues, long nintparm, long naddparm );
void punchFITS_ParamData( char **paramNames,
						 long *paramMethods,
						 float **paramRange,
						 float **paramData,
						 long nintparm,
						 long naddparm,
						 long numParamValues );
void punchFITS_EnergyHeader( long numEnergies );
void punchFITS_EnergyData( float *Energies, long numEnergies );
void punchFITS_SpectraHeader( long nintparm, long naddparm, long totNumModels, long numEnergies );
void punchFITS_SpectraData( float **interpParameters,
						   float **theSpectrum,
						   long totNumModels,
						   long numEnergies,
						   long nintparm,
						   long naddparm  );
void punchFITS_GenericHeader( long numEnergies );
void punchFITS_GenericData( long numEnergies );
void writeCloudyDetails( void );
long addComment( char *CommentToAdd );
long addKeyword_txt( char *theKeyword, const void *theValue, char *theComment, long Str_Or_Log );
long addKeyword_num( char *theKeyword, long theValue, char *theComment);
float swap_float(float f);

void punchFITSfile( FILE* io, char *type_of_file)
{

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

	fits_output = io;

	if( lgMatch("XSPA",type_of_file) )
	{
		punchFITS_PrimaryHeader( TRUE );
		punchFITS_ParamHeader( grid.numParamValues, grid.nintparm, grid.naddparm );
		punchFITS_ParamData( grid.paramNames, grid.paramMethods, grid.paramRange, grid.paramData,
			grid.nintparm, grid.naddparm, grid.numParamValues );
		punchFITS_EnergyHeader( grid.numEnergies );
		punchFITS_EnergyData( grid.Energies, grid.numEnergies );
		punchFITS_SpectraHeader( grid.nintparm, grid.naddparm, grid.totNumModels, grid.numEnergies);
		punchFITS_SpectraData( grid.interpParameters, grid.Spectra,
			grid.totNumModels, grid.numEnergies, grid.nintparm, grid.naddparm );
	}
	else if( lgMatch("XSPM",type_of_file) )
	{
		punchFITS_PrimaryHeader( FALSE );
		punchFITS_ParamHeader( grid.numParamValues, grid.nintparm, grid.naddparm );
		punchFITS_ParamData( grid.paramNames, grid.paramMethods, grid.paramRange, grid.paramData,
			grid.nintparm, grid.naddparm, grid.numParamValues );
		punchFITS_EnergyHeader( grid.numEnergies );
		punchFITS_EnergyData( grid.Energies, grid.numEnergies );
		punchFITS_SpectraHeader( grid.nintparm, grid.naddparm, grid.totNumModels, grid.numEnergies);
		punchFITS_SpectraData( grid.interpParameters, grid.ExpMinusTau,
			grid.totNumModels, grid.numEnergies, grid.nintparm, grid.naddparm );
	}
	/* in this case, we just have a regular run and are punching the predicted spectrum in
	 * the form of a FITS file. */
	else if( lgMatch("FITS",type_of_file) )
	{
		punchFITS_PrimaryHeader( FALSE );
		punchFITS_GenericHeader( rfield.nflux - 1 );
		punchFITS_GenericData( rfield.nflux -1 );
	}
	else 
	{
		fprintf( ioQQQ, "PROBLEM:  This is not an option. \n" );
		cdEXIT( EXIT_FAILURE );
	}

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

	return;
}

void punchFITS_PrimaryHeader( int lgAddModel )
{
	char *ModelName = "'CLOUDY'";

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

	bytesAdded = 0;
	
	bytesAdded += addKeyword_txt( "SIMPLE"	, "T",					"file does conform to FITS standard", 1 );
	bytesAdded += addKeyword_num( "BITPIX"	, bitpix,				"number of bits per data pixel" );
	bytesAdded += addKeyword_num( "NAXIS"	, 0,					"number of data axes" );
	bytesAdded += addKeyword_txt( "EXTEND"	, "T",					"FITS dataset may contain extensions", 1 );
	bytesAdded += addKeyword_txt( "CONTENT" , "'MODEL   '",			"spectrum file contains time intervals and event", 0 );
	bytesAdded += addKeyword_txt( "MODLNAME", ModelName,			"Model name", 0 );
	bytesAdded += addKeyword_txt( "MODLUNIT", "'photons/cm^2/s'",	"Model units", 0 );
	bytesAdded += addKeyword_txt( "REDSHIFT", "T",				"If true then redshift will be included as a par", 1 );
	if( lgAddModel == TRUE )
	{
		bytesAdded += addKeyword_txt( "ADDMODEL", "T",				"If true then this is an additive table model", 1 );
	}
	else
	{
		bytesAdded += addKeyword_txt( "ADDMODEL", "F",				"If true then this is an additive table model", 1 );
	}

	/* bytes are added here as well */
	writeCloudyDetails();

	bytesAdded += addKeyword_txt( "HDUCLASS", "'OGIP    '",			"Format conforms to OGIP/GSFC conventions", 0 );
	bytesAdded += addKeyword_txt( "HDUCLAS1", "'XSPEC TABLE MODEL'","Extension contains an image", 0 );
	bytesAdded += addKeyword_txt( "HDUVERS"	, "'1.0.0   '",			"Version of format (OGIP memo OGIP-92-001)", 0 );
		/* After everything else */
	bytesAdded += fprintf(fits_output, "%-80s", "END" );
	
	ASSERT( bytesAdded%LINESIZE == 0 );

	/* Now add blanks */
	while( bytesAdded%RECORDSIZE > 0 )
	{
		bytesAdded += fprintf(fits_output, "%-1s", " " );
	}
#	ifdef DEBUG_FUN
	fputs( "<->punchFITS_PrimaryHeader()\n", debug_fp );
#	endif

	return;
}

void punchFITS_ParamHeader( long numParamValues, long nintparm, long naddparm )
{
	long numFields = 10;
	long naxis, naxis1, naxis2;
	char theValue[20];
	char theValue_temp[20];

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

	/* Make sure the previous blocks are the right size */
	ASSERT( bytesAdded%RECORDSIZE == 0 );

	naxis = 2;
	naxis1 = 44+4*numParamValues;
	naxis2 = nintparm+naddparm;

	bytesAdded += addKeyword_txt( "XTENSION", "'BINTABLE'",			"binary table extension", 0  );
	bytesAdded += addKeyword_num( "BITPIX"	, bitpix,					"8-bit bytes" );
	bytesAdded += addKeyword_num( "NAXIS"	, naxis,					"2-dimensional binary table" );
	bytesAdded += addKeyword_num( "NAXIS1"	, naxis1,				"width of table in bytes" );
	bytesAdded += addKeyword_num( "NAXIS2"	, naxis2,					"number of rows in table" );
	bytesAdded += addKeyword_num( "PCOUNT"	, pcount,					"size of special data area" );
	bytesAdded += addKeyword_num( "GCOUNT"	, gcount,					"one data group (required keyword)" );
	bytesAdded += addKeyword_num( "TFIELDS"	, numFields,			"number of fields in each row" );
	bytesAdded += addKeyword_txt( "TTYPE1"	, "'NAME    '",			"label for field   1", 0  );
	bytesAdded += addKeyword_txt( "TFORM1"	, "'12A     '",			"data format of the field: ASCII Character", 0  );
	bytesAdded += addKeyword_txt( "TTYPE2"	, "'METHOD  '",			"label for field   2", 0  );
	bytesAdded += addKeyword_txt( "TFORM2"	, "'J       '",			"data format of the field: 4-byte INTEGER", 0  );
	bytesAdded += addKeyword_txt( "TTYPE3"	, "'INITIAL '",			"label for field   3", 0  );
	bytesAdded += addKeyword_txt( "TFORM3"	, "'E       '",			"data format of the field: 4-byte REAL", 0  );
	bytesAdded += addKeyword_txt( "TTYPE4"	, "'DELTA   '",			"label for field   4", 0  );
	bytesAdded += addKeyword_txt( "TFORM4"	, "'E       '",			"data format of the field: 4-byte REAL", 0  );
	bytesAdded += addKeyword_txt( "TTYPE5"	, "'MINIMUM '",			"label for field   5", 0  );
	bytesAdded += addKeyword_txt( "TFORM5"	, "'E       '",			"data format of the field: 4-byte REAL", 0  );
	bytesAdded += addKeyword_txt( "TTYPE6"	, "'BOTTOM  '",			"label for field   6", 0  );
	bytesAdded += addKeyword_txt( "TFORM6"	, "'E       '",			"data format of the field: 4-byte REAL", 0  );
	bytesAdded += addKeyword_txt( "TTYPE7"	, "'TOP     '",			"label for field   7", 0  );
	bytesAdded += addKeyword_txt( "TFORM7"	, "'E       '",			"data format of the field: 4-byte REAL", 0  );
	bytesAdded += addKeyword_txt( "TTYPE8"	, "'MAXIMUM '",			"label for field   8", 0  );
	bytesAdded += addKeyword_txt( "TFORM8"	, "'E       '",			"data format of the field: 4-byte REAL", 0  );
	bytesAdded += addKeyword_txt( "TTYPE9"	, "'NUMBVALS'",			"label for field   9", 0  );
	bytesAdded += addKeyword_txt( "TFORM9"	, "'J       '",			"data format of the field: 4-byte INTEGER", 0  );
	bytesAdded += addKeyword_txt( "TTYPE10"	, "'VALUE   '",			"label for field  10", 0  );
	
	/* The size of this array is dynamic, set to size of numParamValues */
	sprintf( theValue_temp,		"%ld%s", numParamValues, "E" );
	sprintf( theValue,		"%s%-8s%s", "'",theValue_temp,"'" );
	bytesAdded += addKeyword_txt( "TFORM10"	, theValue,			"data format of the field: 4-byte REAL", 0  );

	bytesAdded += addKeyword_txt( "EXTNAME"	, "'PARAMETERS'",		"name of this binary table extension", 0  );
	bytesAdded += addKeyword_txt( "HDUCLASS", "'OGIP    '",			"Format conforms to OGIP/GSFC conventions", 0  );
	bytesAdded += addKeyword_txt( "HDUCLAS1", "'XSPEC TABLE MODEL'","model spectra for XSPEC", 0  );
	bytesAdded += addKeyword_txt( "HDUCLAS2", "'PARAMETERS'",		"Extension containing paramter info", 0  );
	bytesAdded += addKeyword_txt( "HDUVERS"	, "'1.0.0   '",			"Version of format (OGIP memo OGIP-92-001)", 0  );
	bytesAdded += addKeyword_num( "NINTPARM", nintparm,				"Number of interpolation parameters" );
	bytesAdded += addKeyword_num( "NADDPARM", naddparm,				"Number of additional parameters" );
	/* After everything else */
	bytesAdded += fprintf(fits_output, "%-80s", "END" );
	
	ASSERT( bytesAdded%LINESIZE == 0 );

	/* Now add blanks */
	while( bytesAdded%RECORDSIZE > 0 )
	{
		bytesAdded += fprintf(fits_output, "%-1s", " " );
	}

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

	return;
}

void punchFITS_ParamData( char **paramNames,
					   long *paramMethods,
					   float **paramRange,
					   float **paramData,
					   long nintparm,
					   long naddparm,
					   long numParamValues )
{
	long i, j;

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

	/* Now add the parameters data */
	for( i=0; i<nintparm+naddparm; i++ )
	{
		int32 numTemp;

		paramMethods[i] = HtoNL(paramMethods[i]);
		numTemp = HtoNL(numParamValues);
		paramRange[i][0] = swap_float( paramRange[i][0] );
		paramRange[i][1] = swap_float( paramRange[i][1] );
		paramRange[i][2] = swap_float( paramRange[i][2] );
		paramRange[i][3] = swap_float( paramRange[i][3] );
		paramRange[i][4] = swap_float( paramRange[i][4] );
		paramRange[i][5] = swap_float( paramRange[i][5] );

		for( j=0; j<numParamValues; j++ )
		{
			paramData[i][j] = swap_float( paramData[i][j] );
		}
		
		bytesAdded += fprintf(fits_output, "%-12s", paramNames[i] );
		bytesAdded += (long)fwrite( &paramMethods[i],	1,				  sizeof(int32),   fits_output );
		bytesAdded += (long)fwrite( paramRange[i],		1,				6*sizeof(float),   fits_output );
		bytesAdded += (long)fwrite( &numTemp,			1,				  sizeof(int32),   fits_output );
		bytesAdded += (long)fwrite( paramData[i],		1, (unsigned)numParamValues*sizeof(float),   fits_output );
	}

	/* Switch the endianness again */
	for( i=0; i<nintparm+naddparm; i++ )
	{
		paramMethods[i] = HtoNL(paramMethods[i]);
		paramRange[i][0] = swap_float( paramRange[i][0] );
		paramRange[i][1] = swap_float( paramRange[i][1] );
		paramRange[i][2] = swap_float( paramRange[i][2] );
		paramRange[i][3] = swap_float( paramRange[i][3] );
		paramRange[i][4] = swap_float( paramRange[i][4] );
		paramRange[i][5] = swap_float( paramRange[i][5] );

		for( j=0; j<numParamValues; j++ )
		{
			paramData[i][j] = swap_float( paramData[i][j] );
		}
	}

	while( bytesAdded%RECORDSIZE > 0 )
	{
		int	tempInt = 0;
		bytesAdded += (long)fwrite( &tempInt, 1, 1,   fits_output );
	}

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

	return;
}

void punchFITS_EnergyHeader( long numEnergies )
{
	long numFields = 2;
	long naxis, naxis1, naxis2;

#	ifdef DEBUG_FUN
	fputs( "<+>punchFITS_EnergyHeader()\n", debug_fp );
#	endif
	
	/* Make sure the previous blocks are the right size */
	ASSERT( bytesAdded%RECORDSIZE == 0 );

	naxis = 2;
	naxis1 = 2*sizeof(float);
	naxis2 = numEnergies;

	bytesAdded += addKeyword_txt( "XTENSION", "'BINTABLE'",			"binary table extension", 0 );
	bytesAdded += addKeyword_num( "BITPIX"	, bitpix,					"8-bit bytes" );
	bytesAdded += addKeyword_num( "NAXIS"	, naxis,					"2-dimensional binary table" );
	bytesAdded += addKeyword_num( "NAXIS1"	, naxis1,					"width of table in bytes" );
	bytesAdded += addKeyword_num( "NAXIS2"	, naxis2,					"number of rows in table" );
	bytesAdded += addKeyword_num( "PCOUNT"	, pcount,					"size of special data area" );
	bytesAdded += addKeyword_num( "GCOUNT"	, gcount,					"one data group (required keyword)" );
	bytesAdded += addKeyword_num( "TFIELDS"	, numFields,			"number of fields in each row" );
	bytesAdded += addKeyword_txt( "TTYPE1"	, "'ENERG_LO'",			"label for field   1", 0  );
	bytesAdded += addKeyword_txt( "TFORM1"	, "'E       '",			"data format of the field: 4-byte REAL", 0  );
	bytesAdded += addKeyword_txt( "TTYPE2"	, "'ENERG_HI'",			"label for field   2", 0  );
	bytesAdded += addKeyword_txt( "TFORM2"	, "'E       '",			"data format of the field: 4-byte REAL", 0  );
	bytesAdded += addKeyword_txt( "EXTNAME"	, "'ENERGIES'",			"name of this binary table extension", 0  );
	bytesAdded += addKeyword_txt( "HDUCLASS", "'OGIP    '",			"Format conforms to OGIP/GSFC conventions", 0  );
	bytesAdded += addKeyword_txt( "HDUCLAS1", "'XSPEC TABLE MODEL'","model spectra for XSPEC", 0  );
	bytesAdded += addKeyword_txt( "HDUCLAS2", "'ENERGIES'",			"Extension containing energy bin info", 0  );
	bytesAdded += addKeyword_txt( "HDUVERS"	, "'1.0.0   '",			"Version of format (OGIP memo OGIP-92-001)", 0  );
	/* After everything else */
	bytesAdded += fprintf(fits_output, "%-80s", "END" );
	
	ASSERT( bytesAdded%LINESIZE == 0 );

	while( bytesAdded%RECORDSIZE > 0 )
	{
		bytesAdded += fprintf(fits_output, "%-1s", " " );
	}

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

	return;
}

void punchFITS_EnergyData( float *Energies, long numEnergies )
{
	long i;

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

	/* Now add the energies data */
	for( i=0; i<numEnergies; i++ )
	{
		float EnergyLow, EnergyHi;
		EnergyLow = Energies[i]-rfield.widflx[i]/2.f;
		
		if( i == numEnergies-1 )
		{
			EnergyHi = Energies[i] + rfield.widflx[i]/2.f;
		}
		else
		{
			EnergyHi = Energies[i+1]-rfield.widflx[i+1]/2.f;
		}

#if !defined(_BIG_ENDIAN) 
		ByteSwap5(EnergyLow);
		ByteSwap5(EnergyHi);
#endif

		bytesAdded += (long)fwrite( &EnergyLow	, 1, sizeof(float), fits_output );
		bytesAdded += (long)fwrite( &EnergyHi	, 1, sizeof(float), fits_output );
	}

	while( bytesAdded%RECORDSIZE > 0 )
	{
		int	tempInt = 0;
		bytesAdded += (long)fwrite( &tempInt, 1, 1,   fits_output );
	}

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

	return;
}

void punchFITS_SpectraHeader( long nintparm, long naddparm, long totNumModels, long numEnergies )
{
	long i, numFields = 2+naddparm;
	long naxis, naxis1, naxis2;
	char theKeyword1[8];
	char theKeyword2[8];
	char theKeyword3[8];
	char theValue1[20];
	char theValue2[20];
	char theValue2temp[20];
	char theValue[20];
	char theValue_temp[20];
	char theComment1[47];

#	ifdef DEBUG_FUN
	fputs( "<+>punchFITS_SpectraHeader()\n", debug_fp );
#	endif
	
	/* Make sure the previous blocks are the right size */
	ASSERT( bytesAdded%RECORDSIZE == 0 );

	naxis = 2;
	naxis1 = ( numEnergies*(naddparm+1) + nintparm ) * (long)sizeof(float);
	naxis2 = totNumModels; 

	bytesAdded += addKeyword_txt( "XTENSION", "'BINTABLE'",			"binary table extension", 0  );
	bytesAdded += addKeyword_num( "BITPIX"	, bitpix,				"8-bit bytes" );
	bytesAdded += addKeyword_num( "NAXIS"	, naxis,				"2-dimensional binary table" );
	bytesAdded += addKeyword_num( "NAXIS1"	, naxis1,				"width of table in bytes" );
	bytesAdded += addKeyword_num( "NAXIS2"	, naxis2,				"number of rows in table" );
	bytesAdded += addKeyword_num( "PCOUNT"	, pcount,				"size of special data area" );
	bytesAdded += addKeyword_num( "GCOUNT"	, gcount,				"one data group (required keyword)" );
	bytesAdded += addKeyword_num( "TFIELDS"	, numFields,			"number of fields in each row" );
	
	/******************************************/
	/* These are the interpolation parameters */
	/******************************************/
	bytesAdded += addKeyword_txt( "TTYPE1"	, "'PARAMVAL'",			"label for field   1", 0 );
	/* The size of this array is dynamic, set to size of nintparm */
	sprintf( theValue2temp,		"%ld%s", nintparm, "E" );
	sprintf( theValue2,		"%s%-8s%s", "'",theValue2temp,"'" );
	bytesAdded += addKeyword_txt( "TFORM1"	, theValue2,			"data format of the field: 4-byte REAL", 0  );

	/******************************************/
	/* This is the interpolated spectrum      */	
	/******************************************/
	bytesAdded += addKeyword_txt( "TTYPE2"	, "'INTPSPEC'",	"label for field 2", 0  );
	/* The size of this array is dynamic, set to size of numEnergies */
	sprintf( theValue_temp,		"%ld%s", numEnergies, "E" );
	sprintf( theValue,		"%s%-8s%s", "'",theValue_temp,"'" );
	bytesAdded += addKeyword_txt( "TFORM2"	, theValue,	"data format of the field: 4-byte REAL", 0  );
	bytesAdded += addKeyword_txt( "TUNIT2"	, "'photons/cm^2/s'",	"physical unit of field", 0  );

	/******************************************/
	/* These are the additional parameters    */
	/******************************************/
	for( i=1; i<=naddparm; i++ )
	{
		sprintf( theKeyword1,	"%s%ld", "TTYPE", i+2 );
		sprintf( theKeyword2,	"%s%ld", "TFORM", i+2 );
		sprintf( theKeyword3,	"%s%ld", "TUNIT", i+2 );
		
		sprintf( theValue1,		"%s%2.2ld%s", "'ADDSP", i, "'" );
		/* The size of this array is dynamic, set to size of numEnergies */
		sprintf( theValue2temp,		"%ld%s", numEnergies, "E" );
		sprintf( theValue2,		"%s%-8s%s", "'",theValue2temp,"'" );

		sprintf( theComment1,	"%s%ld", "label for field ", i+2 );

		bytesAdded += addKeyword_txt( theKeyword1	, theValue1,		theComment1, 0  );
		bytesAdded += addKeyword_txt( theKeyword2	, theValue2,		"data format of the field: 4-byte REAL", 0  );
		bytesAdded += addKeyword_txt( theKeyword3	, "'photons/cm^2/s'",	"physical unit of field", 0  );
	}
	
	bytesAdded += addKeyword_txt( "EXTNAME"	, "'SPECTRA '",			"name of this binary table extension", 0  );
	bytesAdded += addKeyword_txt( "HDUCLASS", "'OGIP    '",			"Format conforms to OGIP/GSFC conventions", 0  );
	bytesAdded += addKeyword_txt( "HDUCLAS1", "'XSPEC TABLE MODEL'","model spectra for XSPEC", 0  );
	bytesAdded += addKeyword_txt( "HDUCLAS2", "'MODEL SPECTRA'",	"Extension containing model spectra", 0  );
	bytesAdded += addKeyword_txt( "HDUVERS"	, "'1.0.0   '",			"Version of format (OGIP memo OGIP-92-001)", 0  );
	/* After everything else */
	bytesAdded += fprintf(fits_output, "%-80s", "END" );

	ASSERT( bytesAdded%LINESIZE == 0 );
	
	while( bytesAdded%RECORDSIZE > 0 )
	{
		bytesAdded += fprintf(fits_output, "%-1s", " " );
	}

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

	return;
}

void punchFITS_SpectraData( float **interpParameters, float **theSpectrum,
						   long totNumModels, long numEnergies, long nintparm, long naddparm )
{
	long i, j;
	long naxis2 = totNumModels;

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

	/* Now add the spectra data */
	for( i=0; i<naxis2; i++ )
	{

#if !defined(_BIG_ENDIAN)
		for( j = 0; j<numEnergies; j++ )
		{
			ByteSwap5( theSpectrum[i][j] );
			/* theSpectrum[i][j] = swap_float(  theSpectrum[i][j] ); */
		}

		for( j = 0; j<nintparm; j++ )
		{
			interpParameters[i][j] = swap_float( interpParameters[i][j] );
		}
#endif

		/* The interpolation parameters vector */
		bytesAdded += (long)fwrite( interpParameters[i],	1,	 (unsigned)nintparm*sizeof(float), fits_output );
		/* The interpolated spectrum */
		bytesAdded += (long)fwrite( theSpectrum[i],			1, (unsigned)numEnergies*sizeof(float), fits_output );
		/* The additional parameters */
		for( j=1; j<=naddparm; j++ )
		{
			bytesAdded += (long)fwrite( theSpectrum[i],		1, (unsigned)numEnergies*sizeof(float), fits_output );
		}

#if !defined(_BIG_ENDIAN)
		/* Switch the endianness back to native. */
		for( j = 0; j<numEnergies; j++ )
		{
			ByteSwap5( theSpectrum[i][j] );
			/* theSpectrum[i][j] = swap_float(  theSpectrum[i][j] ); */
		}

		for( j = 0; j<nintparm; j++ )
		{
			interpParameters[i][j] = swap_float( interpParameters[i][j] );
		}
#endif
	}

	while( bytesAdded%RECORDSIZE > 0 )
	{
		int	tempInt = 0;
		bytesAdded += (long)fwrite( &tempInt, 1, 1,   fits_output );
	}

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

	return;
}

void punchFITS_GenericHeader( long numEnergies )
{
	long numFields = 2;
	long naxis, naxis1, naxis2;

#	ifdef DEBUG_FUN
	fputs( "<+>punchFITS_GenericHeader()\n", debug_fp );
#	endif
	
	/* Make sure the previous blocks are the right size */
	ASSERT( bytesAdded%RECORDSIZE == 0 );

	naxis = 2;
	naxis1 = numFields*(long)sizeof(float);
	naxis2 = numEnergies;

	bytesAdded += addKeyword_txt( "XTENSION", "'BINTABLE'",			"binary table extension", 0 );
	bytesAdded += addKeyword_num( "BITPIX"	, bitpix,					"8-bit bytes" );
	bytesAdded += addKeyword_num( "NAXIS"	, naxis,					"2-dimensional binary table" );
	bytesAdded += addKeyword_num( "NAXIS1"	, naxis1,					"width of table in bytes" );
	bytesAdded += addKeyword_num( "NAXIS2"	, naxis2,					"number of rows in table" );
	bytesAdded += addKeyword_num( "PCOUNT"	, pcount,					"size of special data area" );
	bytesAdded += addKeyword_num( "GCOUNT"	, gcount,					"one data group (required keyword)" );
	bytesAdded += addKeyword_num( "TFIELDS"	, numFields,			"number of fields in each row" );
	bytesAdded += addKeyword_txt( "TTYPE1"	, "'ENERGY  '",			"label for field   1", 0  );
	bytesAdded += addKeyword_txt( "TFORM1"	, "'E       '",			"data format of the field: 4-byte REAL", 0  );
	bytesAdded += addKeyword_txt( "TTYPE2"	, "'TRN_SPEC'",			"label for field   2", 0  );
	bytesAdded += addKeyword_txt( "TFORM2"	, "'E       '",			"data format of the field: 4-byte REAL", 0  );
	bytesAdded += addKeyword_txt( "EXTNAME"	, "'SPECTRA '",			"name of this binary table extension", 0  );
	bytesAdded += addKeyword_txt( "HDUCLASS", "'OGIP    '",			"Format conforms to OGIP/GSFC conventions", 0  );
	bytesAdded += addKeyword_txt( "HDUCLAS1", "'XSPEC TABLE MODEL'","model spectra for XSPEC", 0  );
	bytesAdded += addKeyword_txt( "HDUCLAS2", "'ENERGIES'",			"Extension containing energy bin info", 0  );
	bytesAdded += addKeyword_txt( "HDUVERS"	, "'1.0.0   '",			"Version of format (OGIP memo OGIP-92-001)", 0  );
	/* After everything else */
	bytesAdded += fprintf(fits_output, "%-80s", "END" );

	ASSERT( bytesAdded%LINESIZE == 0 );

	while( bytesAdded%RECORDSIZE > 0 )
	{
		bytesAdded += fprintf(fits_output, "%-1s", " " );
	}

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

	return;
}

void punchFITS_GenericData( long numEnergies )
{
	long i;

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

	float *TransmittedSpectrum;

	if( (TransmittedSpectrum = (float*)MALLOC(sizeof(float)*(unsigned)(numEnergies) ) )==NULL )
		BadMalloc();

	cdSPEC2( 2, numEnergies, TransmittedSpectrum );
	
	/* Now add the energies data */
	for( i=0; i<numEnergies; i++ )
	{
		float Energy;
		Energy = rfield.AnuOrg[i];

#if !defined(_BIG_ENDIAN) 
		ByteSwap5(Energy);
		ByteSwap5(TransmittedSpectrum[i]);
#endif

		bytesAdded += (long)fwrite( &Energy	,				1, sizeof(float), fits_output );
		bytesAdded += (long)fwrite( &TransmittedSpectrum[i],1, sizeof(float), fits_output );
	}

	while( bytesAdded%RECORDSIZE > 0 )
	{
		int	tempInt = 0;
		bytesAdded += (long)fwrite( &tempInt, 1, 1,   fits_output );
	}

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

	return;
}

void writeCloudyDetails( void )
{
	char chVer[10];
	char timeString[30]="";
	char tempString[70];
	time_t now;
	long i, j, k;

	cdVersion(chVer ) ;

	/* usually print date and time info - do not if "no times" command entered, 
	 * which set this flag false */
	now = time(NULL);
	if( prt.lgPrintTime ) 
	{
		/* now add date of this run */
		/* now print this time at the end of the string.  the system put cr at the end of the string */
		strcpy( timeString , ctime(&now) );
	}
	/* ctime puts a carriage return at the end, but we can't have that in a fits file.
	 * remove the carriage return here. */
	for( i=0; i<30; i++ )
	{
		if( timeString[i] == '\n' )
		{
			timeString[i] = ' ';
		}
	}

	strcpy( tempString, "Generated by Cloudy version " );
	strcat( tempString, chVer );
	bytesAdded += addComment( tempString );
	bytesAdded += addComment( version.chInfo );
	strcpy( tempString, "--- " );
	strcat( tempString, timeString );
	bytesAdded += addComment( tempString );
	bytesAdded += addComment( "Input string was as follows: " );
	/* >>chng 05 nov 24, from <nSave to <=nSave bug caught by PvH */
	for( i=0; i<=input.nSave; i++ )
	{
		char firstLine[70], extraLine[65];

		for( j=0; j<INPUT_LINE_LENGTH; j++ )
		{
			if( input.chCardSav[i][j] =='\0' )
				break;
		}
		
		ASSERT( j < 200 );
		for( k=0; k< MIN(69, j); k++ )
		{
			firstLine[k] = input.chCardSav[i][k];
		}
		firstLine[k] = '\0';
		bytesAdded += addComment( firstLine );
		if( j >= 69 )
		{
			for( k=69; k< 133; k++ )
			{
				extraLine[k-69] = input.chCardSav[i][k];
			}
			/* >> chng 06 jan 05, this was exceeding array bounds. */
			extraLine[64] = '\0'; 
			strcpy( tempString, "more " );
			strcat( tempString, extraLine );
			bytesAdded += addComment( tempString );
			if( j >= 133 )
			{
				for( k=133; k< 197; k++ )
				{
					extraLine[k-133] = input.chCardSav[i][k];
				}
				extraLine[64] = '\0';
				strcpy( tempString, "more " );
				strcat( tempString, extraLine );
				bytesAdded += addComment( tempString );
			}
		}
	}

	return;
}

long addKeyword_txt( char *theKeyword, const void *theValue, char *theComment, long Str_Or_Log )
{
	long numberOfBytesWritten = 0;

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

	/* False means string, true means logical */
	if( Str_Or_Log == 0 )
	{
		numberOfBytesWritten = fprintf(fits_output, "%-8s%-2s%-20s%3s%-47s",
			theKeyword,
			"= ",
			(char *)theValue,
			" / ",
			theComment );
	}
	else
	{
		ASSERT( Str_Or_Log == 1 );
		numberOfBytesWritten = fprintf(fits_output, "%-8s%-2s%20s%3s%-47s",
			theKeyword,
			"= ",
			(char *)theValue,
			" / ",
			theComment );
	}

	ASSERT( numberOfBytesWritten%LINESIZE == 0 );

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

	return numberOfBytesWritten;
}

long addKeyword_num( char *theKeyword, long theValue, char *theComment)
{
	long numberOfBytesWritten = 0;

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

	numberOfBytesWritten = fprintf(fits_output, "%-8s%-2s%20ld%3s%-47s",
		theKeyword,
		"= ",
		theValue,
		" / ",
		theComment );

	ASSERT( numberOfBytesWritten%LINESIZE == 0 );

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

	return numberOfBytesWritten;
}

long addComment( char *CommentToAdd )
{
	long numberOfBytesWritten = 0;

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

	numberOfBytesWritten = fprintf(fits_output, "COMMENT   %-70s", CommentToAdd );

	ASSERT( numberOfBytesWritten%LINESIZE == 0 );

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

	return numberOfBytesWritten;
}

float swap_float(float f)
{
	union
	{
		unsigned int integer;
		float real;
	} integer_real;

	ASSERT(sizeof(float) == 4);
	ASSERT(sizeof(unsigned int) == 4);

	integer_real.real = f;
	integer_real.integer = HtoNL(integer_real.integer);

	return integer_real.real;
}
