/home66/gary/public_html/cloudy/c08_branch/source/punch_fits.cpp

Go to the documentation of this file.

/* This file is part of Cloudy and is copyright (C)1978-2008 by Gary J. Ferland and
 * others.  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"
#include "physconst.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))

#if !defined(_BIG_ENDIAN) 
STATIC 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;
}
#endif

static FILE *ioFITS_OUTPUT;
static long bytesAdded = 0;
static long bitpix = 8;
static long pcount = 0;
static long gcount = 1;
static long maxParamValues = 0;
const char ModelUnits[2][17] = {"'dimensionless '", "'photons/cm^2/s'" };

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

void punchFITSfile( FILE* ioPUN, int option )
{

        long i;

        DEBUG_ENTRY( "punchFITSfile()" );

        if( !grid.lgGridDone && option != NUM_OUTPUT_TYPES )
        {
                /* don't do any of this until flag set true. */
                return;
        }

        ioFITS_OUTPUT = ioPUN;

#if     0
        {

                long i,j;
                FILE *asciiDump;

                if( (asciiDump = fopen(  "gridspectra.con","w" ) ) == NULL )
                {
                        fprintf( ioQQQ, "punchFITSfile could not open punch file for writing.\nSorry.\n" );
                        cdEXIT(EXIT_FAILURE);
                }

                for( i=0; i<grid.numEnergies-1; i++ )
                {
                        fprintf( asciiDump, "%.12e\t", grid.Energies[i] );
                        for( j=0; j<grid.totNumModels; j++ )
                        {
                                fprintf( asciiDump, "%.12e\t", grid.Spectra[10][j][i] );
                        }
                        fprintf( asciiDump, "\n" );
                }
                fclose( asciiDump );
        }
#endif

        /* This is generic FITS option */
        if( option == NUM_OUTPUT_TYPES )
        {
                punchFITS_PrimaryHeader( false );
                punchFITS_GenericHeader( rfield.nflux - 1 );
                punchFITS_GenericData( rfield.nflux -1 );
        }
        /* These are specially designed XSPEC outputs. */
        else if( option >= 1 && option < NUM_OUTPUT_TYPES )
        {
                bool lgAdditiveModel;

                /* option 10 is exp(-tau). */
                if( option == 10 )
                {
                        /* false says not an additive model */
                        lgAdditiveModel = false;
                }
                else
                {
                        lgAdditiveModel = true;                 
                }

                punchFITS_PrimaryHeader( lgAdditiveModel );

                for( i=0; i<grid.nintparm+grid.naddparm; i++ )
                {
                        maxParamValues = MAX2( maxParamValues, grid.numParamValues[i] );
                }

                ASSERT( maxParamValues >= 2 );

                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( lgAdditiveModel, grid.nintparm, grid.naddparm, grid.totNumModels, grid.numEnergies);
                punchFITS_SpectraData( grid.interpParameters, grid.Spectra[option],
                        grid.totNumModels, grid.numEnergies, grid.nintparm, grid.naddparm );
        }
        else 
        {
                fprintf( ioQQQ, "PROBLEM - undefined option encountered in punchFITSfile. \n" );
                cdEXIT( EXIT_FAILURE );
        }
        return;
}

STATIC void punchFITS_PrimaryHeader( bool lgAddModel )
{
        const char *ModelName = "'CLOUDY'";

        DEBUG_ENTRY( "punchFITS_PrimaryHeader()" );

        bytesAdded = 0;

        fixit(); /* bitpix is wrong when realnum is double? */

        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", ModelUnits[lgAddModel],       "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(ioFITS_OUTPUT, "%-80s", "END" );

        ASSERT( bytesAdded%LINESIZE == 0 );

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

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

        DEBUG_ENTRY( "punchFITS_ParamHeader()" );

        ASSERT( nintparm+naddparm <= LIMPAR );

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

        naxis = 2;
        /* >>chng 06 aug 23, change to maximum number of parameter values. */
        naxis1 = 44+4*maxParamValues;
        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  );

        /* >>chng 06 aug 23, use maxParamValues instead of numParamValues */
        /* The size of this array is dynamic, set to size of the maximum of the numParamValues vector */
        sprintf( theValue_temp,         "%ld%s", maxParamValues, "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(ioFITS_OUTPUT, "%-80s", "END" );

        ASSERT( bytesAdded%LINESIZE == 0 );

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

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

        DEBUG_ENTRY( "punchFITS_ParamData()" );

        ASSERT( nintparm+naddparm <= LIMPAR );

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

#define LOG2LINEAR 0

                paramMethods[i] = HtoNL(paramMethods[i]);
                /* >>chng 06 aug 23, numParamValues is now an array.  */
                numTemp = HtoNL(numParamValues[i]);

#if LOG2LINEAR
                /* change to linear */
                paramRange[i][0] = (realnum)pow( 10., (double)paramRange[i][0] );
                paramRange[i][1] = (realnum)pow( 10., (double)paramRange[i][1] );
                paramRange[i][2] = (realnum)pow( 10., (double)paramRange[i][2] );
                paramRange[i][3] = (realnum)pow( 10., (double)paramRange[i][3] );
                paramRange[i][4] = (realnum)pow( 10., (double)paramRange[i][4] );
                paramRange[i][5] = (realnum)pow( 10., (double)paramRange[i][5] );
#endif

#if !defined(_BIG_ENDIAN) 
                ByteSwap5( paramRange[i][0] );
                ByteSwap5( paramRange[i][1] );
                ByteSwap5( paramRange[i][2] );
                ByteSwap5( paramRange[i][3] );
                ByteSwap5( paramRange[i][4] );
                ByteSwap5( paramRange[i][5] );
#endif

                /* >>chng 06 aug 23, numParamValues is now an array.  */
                for( j=0; j<numParamValues[i]; j++ )
                {

#if LOG2LINEAR
                        paramData[i][j] = (realnum)pow( 10., (double)paramData[i][j] );
#endif

#if !defined(_BIG_ENDIAN) 
                        ByteSwap5( paramData[i][j] );
#endif
                }

                bytesAdded += fprintf(ioFITS_OUTPUT, "%-12s", paramNames[i] );
                bytesAdded += (long)fwrite( &paramMethods[i],   1,                                sizeof(int32),   ioFITS_OUTPUT );
                bytesAdded += (long)fwrite( paramRange[i],              1,                              6*sizeof(realnum),   ioFITS_OUTPUT );
                bytesAdded += (long)fwrite( &numTemp,                   1,                                sizeof(int32),   ioFITS_OUTPUT );
                /* >>chng 06 aug 23, numParamValues is now an array.  */
                bytesAdded += (long)fwrite( paramData[i],               1, (unsigned)numParamValues[i]*sizeof(realnum),   ioFITS_OUTPUT );

                for( j=numParamValues[i]+1; j<=maxParamValues; j++ )
                {
                        realnum filler = -10.f;
                        bytesAdded += (long)fwrite( &filler,            1, sizeof(realnum),   ioFITS_OUTPUT );
                }
        }

        /* Switch the endianness again */
        for( i=0; i<nintparm+naddparm; i++ )
        {
                paramMethods[i] = HtoNL(paramMethods[i]);

#if !defined(_BIG_ENDIAN) 
                ByteSwap5( paramRange[i][0] );
                ByteSwap5( paramRange[i][1] );
                ByteSwap5( paramRange[i][2] );
                ByteSwap5( paramRange[i][3] );
                ByteSwap5( paramRange[i][4] );
                ByteSwap5( paramRange[i][5] );
#endif

                /* >>chng 06 aug 23, numParamValues is now an array.  */
                for( j=0; j<numParamValues[i]; j++ )
                {
#if !defined(_BIG_ENDIAN) 
                        ByteSwap5( paramData[i][j] );
#endif
                }
        }

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

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

        DEBUG_ENTRY( "punchFITS_EnergyHeader()" );

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

        naxis = 2;
        naxis1 = 2*sizeof(realnum);
        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(ioFITS_OUTPUT, "%-80s", "END" );

        ASSERT( bytesAdded%LINESIZE == 0 );

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

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

        DEBUG_ENTRY( "punchFITS_EnergyData()" );

        /* Now add the energies data */
        for( i=0; i<numEnergies; i++ )
        {
                realnum EnergyLow, EnergyHi;
                /* Convert all of these to kev */
                EnergyLow = 0.001f*(realnum)EVRYD*(Energies[i]-rfield.widflx[i]/2.f);

                if( i == numEnergies-1 )
                {
                        EnergyHi = 0.001f*(realnum)EVRYD*(Energies[i] + rfield.widflx[i]/2.f);
                }
                else
                {
                        EnergyHi = 0.001f*(realnum)EVRYD*(Energies[i+1] - rfield.widflx[i+1]/2.f);
                }

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

                bytesAdded += (long)fwrite( &EnergyLow  , 1, sizeof(realnum), ioFITS_OUTPUT );
                bytesAdded += (long)fwrite( &EnergyHi   , 1, sizeof(realnum), ioFITS_OUTPUT );
        }

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

STATIC void punchFITS_SpectraHeader( bool lgAddModel, 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];

        DEBUG_ENTRY( "punchFITS_SpectraHeader()" );

        ASSERT( nintparm + naddparm <= LIMPAR );

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

        naxis = 2;
        naxis1 = ( numEnergies*(naddparm+1) + nintparm ) * (long)sizeof(realnum);
        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"  , ModelUnits[lgAddModel],       "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       , ModelUnits[lgAddModel],       "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(ioFITS_OUTPUT, "%-80s", "END" );

        ASSERT( bytesAdded%LINESIZE == 0 );

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

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

        DEBUG_ENTRY( "punchFITS_SpectraData()" );

        ASSERT( nintparm + naddparm <= LIMPAR );

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

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

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

                /* The interpolation parameters vector */
                bytesAdded += (long)fwrite( interpParameters[i],        1,       (unsigned)nintparm*sizeof(realnum), ioFITS_OUTPUT );
                /* The interpolated spectrum */
                bytesAdded += (long)fwrite( theSpectrum[i],                     1, (unsigned)numEnergies*sizeof(realnum), ioFITS_OUTPUT );

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

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

                /* >>chng 06 aug 23, disable additional parameters for now */
                if( naddparm > 0 )
                {
                        /* The additional parameters */

                        /* bytesAdded += (long)fwrite( theSpectrum[i],          1, (unsigned)numEnergies*sizeof(realnum), ioFITS_OUTPUT ); */
                        /* \todo 2      must create another array if we are to punch additional parameter information. */
                        fprintf( ioQQQ, " Additional parameters not currently supported.\n" );
                        cdEXIT( EXIT_FAILURE );
                }
        }

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

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

        DEBUG_ENTRY( "punchFITS_GenericHeader()" );

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

        naxis = 2;
        naxis1 = numFields*(long)sizeof(realnum);
        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(ioFITS_OUTPUT, "%-80s", "END" );

        ASSERT( bytesAdded%LINESIZE == 0 );

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

STATIC void punchFITS_GenericData( long numEnergies )
{
        long i;

        DEBUG_ENTRY( "punchFITS_GenericData()" );

        realnum *TransmittedSpectrum;

        TransmittedSpectrum = (realnum*)MALLOC(sizeof(realnum)*(unsigned)(numEnergies) );

        cdSPEC2( 8, numEnergies, TransmittedSpectrum );

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

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

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

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

        free( TransmittedSpectrum );
        return;
}

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

        /* 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 " );
        // strncat guarantees that terminating 0 byte will always be written...
        strncat( tempString, t_version::Inst().chVersion, sizeof(tempString)-strlen(tempString) );
        bytesAdded += addComment( tempString );
        bytesAdded += addComment( t_version::Inst().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< MIN2(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;
}

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

        DEBUG_ENTRY( "addKeyword_txt()" );

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

        ASSERT( numberOfBytesWritten%LINESIZE == 0 );
        return numberOfBytesWritten;
}

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

        DEBUG_ENTRY( "addKeyword_num()" );

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

        ASSERT( numberOfBytesWritten%LINESIZE == 0 );
        return numberOfBytesWritten;
}

long addComment( const char *CommentToAdd )
{
        long i, numberOfBytesWritten = 0;
        char tempString[70] = "                                                                     ";

        DEBUG_ENTRY( "addComment()" );

        strncpy( &tempString[0], CommentToAdd, 69 );
        ASSERT( (int)strlen( tempString ) <= 70 );      

        /* tabs violate FITS standard, replace them with spaces. */
        for( i=0; i<69; i++ )
        {
                if( tempString[i] == '\t' )
                {
                        tempString[i] = ' ';
                }
        }

        numberOfBytesWritten = fprintf(ioFITS_OUTPUT, "COMMENT   %-70s", tempString );

        ASSERT( numberOfBytesWritten%LINESIZE == 0 );
        return numberOfBytesWritten;
}

---