/** @file  tst_psFits.c
*
*  @brief Contains the tests for psFits.[ch]
*
*
*  @author Robert DeSonia, MHPCC
*
*  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
*  @date $Date: 2005-01-18 19:35:24 $
*
*  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
*/

#include "psTest.h"
#include "pslib.h"

#include <unistd.h>
#include <string.h>
#include <sys/stat.h>
#include <sys/types.h>


#define GENIMAGE(img,c,r,TYP, valueFcn) \
img = psImageAlloc(c,r,PS_TYPE_##TYP); \
for (psU32 row=0;row<r;row++) { \
    ps##TYP* imgRow = img->data.TYP[row]; \
    for (psU32 col=0;col<c;col++) { \
        imgRow[col] = (ps##TYP)(valueFcn); \
    } \
}

static bool makeMulti(void);  // implicitly tests psFitsSetExtName
static bool makeTable(void);
const char* multiFilename = "multi.fits";
const char* tableFilename = "table.fits";
const int tableNumRows = 10;


// N.B., the tests to Image read/write was liberally taken from the now
// deprecated psImageReadSection/psImageWriteSection function tests.
static psS32 testImageRead(void);
static psS32 testImageWrite(void);

static psS32 tst_psFitsAlloc( void );
static psS32 tst_psFitsMoveExtName( void ); // also tests psFitsGetExtName
static psS32 tst_psFitsMoveExtNum( void );  // also tests psFitsGetExtNum, psFitsGetSize
static psS32 tst_psFitsReadHeader( void );
static psS32 tst_psFitsReadHeaderSet( void );
static psS32 tst_psFitsReadTable( void );

testDescription tests[] = {
                              {tst_psFitsAlloc, 801, "psFitsAlloc", 0, false},
                              {tst_psFitsMoveExtName, 802, "psFitsMoveExtName", 0, false},
                              {tst_psFitsMoveExtName, 802, "psFitsGetExtName", 0, true},
                              {tst_psFitsMoveExtNum, 803, "psFitsMoveExtNum", 0, false},
                              {tst_psFitsMoveExtNum, 803, "psFitsGetExtNum", 0, true},
                              {tst_psFitsMoveExtNum, 803, "psFitsGetSize", 0, true},
                              {tst_psFitsReadHeader, 804, "psFitsReadHeader", 0, false},
                              {tst_psFitsReadHeaderSet,805, "psFitsReadHeaderSet", 0, false},
                              {tst_psFitsReadTable,809, "psFitsReadTable", 0, false},
                              {testImageRead,567, "psFitsReadImage", 0, false},
                              {testImageWrite,569, "psFitsWriteImage", 0, false},
                              {NULL}
                          };

psS32 main( psS32 argc, char* argv[] )
{
    psLogSetLevel( PS_LOG_INFO );

    return ( ! runTestSuite( stderr, "psImage", tests, argc, argv ) );
}

bool makeMulti(void)
{

    remove
        (multiFilename);

    psFits* fitsFile = psFitsAlloc(multiFilename);

    if (fitsFile == NULL) {
        psError(PS_ERR_UNKNOWN, false,
                "Could not create 'multi' FITS file.");
        return false;
    }

    psImage* image = psImageAlloc(16,16,PS_TYPE_F32);

    char extname[80];
    for (int lcv = 0; lcv < 8; lcv++) {
        snprintf(extname,80,"ext-%d", lcv);

        psMetadata* header = psMetadataAlloc();

        psMetadataAdd(header,PS_LIST_TAIL, "MYINT",
                      PS_TYPE_S32, PS_META_PRIMITIVE,
                      "psS32 Item", (psS32)lcv);

        psMetadataAdd(header,PS_LIST_TAIL, "MYFLT",
                      PS_TYPE_F32, PS_META_PRIMITIVE,
                      "psF32 Item", (float)(1.0f/(float)(1+lcv)));

        psMetadataAdd(header,PS_LIST_TAIL, "MYDBL",
                      PS_TYPE_F64, PS_META_PRIMITIVE,
                      "psF64 Item", (double)(1.0/(double)(1+lcv)));

        psMetadataAdd(header,PS_LIST_TAIL, "MYBOOL",
                      PS_TYPE_BOOL, PS_META_PRIMITIVE,
                      "psBool Item",
                      (lcv%2 == 0));

        psMetadataAdd(header,PS_LIST_TAIL, "MYSTR",
                      PS_TYPE_PTR, PS_META_STR,
                      "String Item",
                      extname);

        // set the pixels in the image
        psBinaryOp(image,image,"=",psScalarAlloc(lcv,PS_TYPE_F32));
        if (! psFitsWriteImage(fitsFile,header,image,1, extname) ) {
            psError(PS_ERR_UNKNOWN, false,
                    "Could not write image.");
            return false;
        }

        psFree(header);
    }
    psFree(image);
    psFree(fitsFile);

    return true;
}

bool makeTable(void)
{

    remove
        (tableFilename);

    psFits* fitsFile = psFitsAlloc(tableFilename);

    if (fitsFile == NULL) {
        psError(PS_ERR_UNKNOWN, false,
                "Could not create 'table' FITS file.");
        return false;
    }

    // make the PHU an image (per FITS standard, it must be)
    psImage* image = psImageAlloc(16,16,PS_TYPE_F32);

    if (! psFitsWriteImage(fitsFile,NULL,image,1,"primary") ) {
        psError(PS_ERR_UNKNOWN, false,
                "Could not write PHU image.");
        return false;
    }

    psFree(image);

    // build a table structure
    psArray* table = psArrayAlloc(tableNumRows);
    table->n = 0;
    for (int row = 0; row < tableNumRows; row++) {
        psMetadata* header = psMetadataAlloc();

        psMetadataAdd(header,PS_LIST_TAIL, "MYINT",
                      PS_TYPE_S32, PS_META_PRIMITIVE,
                      "psS32 Item", (psS32)row);

        psMetadataAdd(header,PS_LIST_TAIL, "MYFLT",
                      PS_TYPE_F32, PS_META_PRIMITIVE,
                      "psF32 Item", (float)(1.0f/(float)(1+row)));

        psMetadataAdd(header,PS_LIST_TAIL, "MYDBL",
                      PS_TYPE_F64, PS_META_PRIMITIVE,
                      "psF64 Item", (double)(1.0/(double)(1+row)));

        psMetadataAdd(header,PS_LIST_TAIL, "MYBOOL",
                      PS_TYPE_BOOL, PS_META_PRIMITIVE,
                      "psBool Item",
                      (row%2 == 0));

        psArrayAdd(table,0,header);
    }

    psFitsWriteTable(fitsFile, NULL, table, "table-1");

    psFree(table);
    psFree(fitsFile);

    return true;
}

psS32 tst_psFitsAlloc( void )
{

    if (! makeMulti() ) {
        return 1;
    }

    psFits* fitsFile = psFitsAlloc(multiFilename);

    if (fitsFile == NULL) {
        psError(PS_ERR_UNKNOWN, false,
                "psFitsAlloc returned NULL on existing file.");
        return 1;
    }

    int extNum = psFitsGetExtNum(fitsFile);
    if (extNum != 0) {
        psError(PS_ERR_UNKNOWN, false,
                "psFitsAlloc was not queued to the PHU, but to extension #%d.",
                extNum);
        return 2;
    }

    psFree(fitsFile);

    // make sure the file doesn't already exist.
    if (access("new.fits", F_OK) == 0) {
        if (remove
                ("new.fits") != 0) {
            psError(PS_ERR_UNKNOWN, false,
                    "Couldn't delete the new.fits file.");
            return 3;
        }
    }

    fitsFile = psFitsAlloc("new.fits");

    if (fitsFile == NULL) {
        psError(PS_ERR_UNKNOWN, false,
                "psFitsAlloc returned NULL on existing file.");
        return 4;
    }

    psFree(fitsFile);

    // now, if psFitsAlloc actually created the file, I shouldn't error in removing it.
    if (remove
            ("new.fits") != 0) {
        psError(PS_ERR_UNKNOWN, false,
                "psFitsAlloc seemed to not have created a new file.");
        return 5;
    }


    return 0;
}

psS32 tst_psFitsMoveExtName( void )
{

    if (! makeMulti() ) {
        return 1;
    }

    psFits* fits = psFitsAlloc(multiFilename);

    if (fits == NULL) {
        psError(PS_ERR_UNKNOWN, false,
                "psFitsAlloc returned NULL on existing file.");
        return 1;
    }

    int numHDUs = psFitsGetSize(fits);

    if (numHDUs < 2) {
        psError(PS_ERR_UNKNOWN,true,
                "The 'multi' FITS file does not have multiple HDUs.");
        return 2;
    }

    char extName[80];
    psRegion region = {0,0,0,0};

    for (int lcv = 0; lcv < numHDUs; lcv++) {
        snprintf(extName,80,"ext-%d",lcv);
        // try to move to the named extension.
        if (! psFitsMoveExtName(fits, extName) ) {
            psError(PS_ERR_UNKNOWN, false,
                    "Failed to move to ext-%d.",
                    lcv);
            return 3;
        }

        // check to see if I can retrieve the name back from the psFits object.
        char* nameFromFile = psFitsGetExtName(fits);
        if (strcmp(nameFromFile,extName) != 0) { // hey, it didn't move?
            psError(PS_ERR_UNKNOWN, false,
                    "Failed to retrieve the extension name back ('%s' vs '%s'",
                    nameFromFile, extName);
            return 3;
        }
        psFree(nameFromFile);

        // check that the image is associated to the extension moved, i.e.,
        // did we really move to the proper extension?
        psImage* image = psFitsReadImage(NULL, fits,region,0);

        if (image == NULL || abs(image->data.F32[0][0] - (float)lcv) > FLT_EPSILON) {
            psError(PS_ERR_UNKNOWN, true,
                    "The image pixel 0,0 of ext-%d was %g, expected %d.",
                    lcv,image->data.F32[0][0],lcv);
            return 4;
        }
        psFree(image);
    }

    for (int lcv = numHDUs-1; lcv >= 0; lcv--) {
        snprintf(extName,80,"ext-%d",lcv);
        // try to move to the named extension.
        if (! psFitsMoveExtName(fits, extName) ) {
            psError(PS_ERR_UNKNOWN, false,
                    "Failed to move to ext-%d.",
                    lcv);
            return 5;
        }

        // check to see if I can retrieve the name back from the psFits object.
        char* nameFromFile = psFitsGetExtName(fits);
        if (strcmp(nameFromFile,extName) != 0) { // hey, it didn't move?
            psError(PS_ERR_UNKNOWN, false,
                    "Failed to retrieve the extension name back ('%s' vs '%s'",
                    nameFromFile, extName);
            return 5;
        }
        psFree(nameFromFile);

        // check that the image is associated to the extension moved, i.e.,
        // did we really move to the proper extension?
        psImage* image = psFitsReadImage(NULL, fits,region,0);

        if (abs(image->data.F32[0][0] - (float)lcv) > FLT_EPSILON) {
            psError(PS_ERR_UNKNOWN, true,
                    "The image pixel 0,0 of ext-%d was %g, expected %d.",
                    lcv,image->data.F32[0][0],lcv);
            return 6;
        }
        psFree(image);
    }

    // check to see if given a bogus extension name, it errors.
    psLogMsg(__func__,PS_LOG_INFO, "Following should be an error.");
    if (psFitsMoveExtName(fits, "bogus") || psErrorGetStackSize() != 1) {
        psError(PS_ERR_UNKNOWN, false,
                "Moving to non-existant HDU didn't fail.");
        return 7;
    }

    // check to see if given a NULL psFits, it errors.
    psLogMsg(__func__,PS_LOG_INFO, "Following should be an error.");
    if (psFitsMoveExtName(NULL, "bogus") || psErrorGetStackSize() != 1) {
        psError(PS_ERR_UNKNOWN, false,
                "Operation of NULL psFits didn't fail.");
        return 8;
    }

    // check to see if given a NULL extname, it errors.
    psLogMsg(__func__,PS_LOG_INFO, "Following should be an error.");
    if (psFitsMoveExtName(fits, NULL) || psErrorGetStackSize() != 1) {
        psError(PS_ERR_UNKNOWN, false,
                "Operation of NULL extname didn't fail.");
        return 9;
    }

    psFree(fits);

    return 0;
}

psS32 tst_psFitsMoveExtNum( void )
{

    if (! makeMulti() ) {
        return 1;
    }

    psFits* fits = psFitsAlloc(multiFilename);

    if (fits == NULL) {
        psError(PS_ERR_UNKNOWN, false,
                "psFitsAlloc returned NULL on existing file.");
        return 1;
    }

    int numHDUs = psFitsGetSize(fits);

    // as a side test, let's make sure psFitsGetSize can handle NULL.
    psLogMsg(__func__,PS_LOG_INFO,
             "Following should be an error.");
    psErrorClear();
    if (psFitsGetSize(NULL) != 0 || psErrorGetStackSize() != 1) {
        psError(PS_ERR_UNKNOWN,true,
                "The 'multi' FITS file does not have multiple HDUs.");
        return 2;
    }

    if (numHDUs != 8) {
        psError(PS_ERR_UNKNOWN,true,
                "The 'multi' FITS file does not have multiple HDUs.");
        return 2;
    }

    psRegion region = {0,0,0,0};

    // test absolute positioning
    for (int lcv = 0; lcv < numHDUs; lcv++) {
        // try to move to the extension
        if (! psFitsMoveExtNum(fits, lcv, false) ) {
            psError(PS_ERR_UNKNOWN, false,
                    "Failed to move to extension %d.",
                    lcv);
            return 3;
        }

        // check to see if I can retrieve the number back from the psFits object.
        if (psFitsGetExtNum(fits) != lcv) { // hey, it didn't move?
            psError(PS_ERR_UNKNOWN, false,
                    "Failed to retrieve the extension number back (%d vs %d)",
                    psFitsGetExtNum(fits), lcv);
            return 5;
        }

        // check that the image is associated to the extension moved, i.e.,
        // did we really move to the proper extension?
        psImage* image = psFitsReadImage(NULL, fits,region,0);

        if (image == NULL || abs(image->data.F32[0][0] - (float)lcv) > FLT_EPSILON) {
            psError(PS_ERR_UNKNOWN, true,
                    "The image pixel 0,0 of ext-%d was %g, expected %d.",
                    lcv,image->data.F32[0][0],lcv);
            return 4;
        }
        psFree(image);
    }

    for (int lcv = numHDUs-1; lcv >= 0; lcv--) {
        // try to move to the extension
        if (! psFitsMoveExtNum(fits, lcv, false) ) {
            psError(PS_ERR_UNKNOWN, false,
                    "Failed to move to extension %d.",
                    lcv);
            return 5;
        }

        // check that the image is associated to the extension moved, i.e.,
        // did we really move to the proper extension?
        psImage* image = psFitsReadImage(NULL, fits,region,0);

        if (abs(image->data.F32[0][0] - (float)lcv) > FLT_EPSILON) {
            psError(PS_ERR_UNKNOWN, true,
                    "The image pixel 0,0 of ext-%d was %g, expected %d.",
                    lcv,image->data.F32[0][0],lcv);
            return 6;
        }
        psFree(image);
    }

    // test relative positioning
    psFitsMoveExtNum(fits,0,false);
    for (int lcv = 1; lcv < numHDUs; lcv++) {
        // try to move to the extension
        if (! psFitsMoveExtNum(fits, 1, true) ) {
            psError(PS_ERR_UNKNOWN, false,
                    "Failed to move to extension %d.",
                    lcv);
            return 13;
        }

        // check to see if I can retrieve the number back from the psFits object.
        if (psFitsGetExtNum(fits) != lcv) { // hey, it didn't move?
            psError(PS_ERR_UNKNOWN, false,
                    "Failed to retrieve the extension number back (%d vs %d)",
                    psFitsGetExtNum(fits), lcv);
            return 13;
        }

        // check that the image is associated to the extension moved, i.e.,
        // did we really move to the proper extension?
        psImage* image = psFitsReadImage(NULL, fits,region,0);

        if (image == NULL || abs(image->data.F32[0][0] - (float)lcv) > FLT_EPSILON) {
            psError(PS_ERR_UNKNOWN, true,
                    "The image pixel 0,0 of ext-%d was %g, expected %d.",
                    lcv,image->data.F32[0][0],lcv);
            return 14;
        }
        psFree(image);
    }

    for (int lcv = numHDUs-2; lcv >= 0; lcv--) {
        // try to move to the extension
        if (! psFitsMoveExtNum(fits, -1, true) ) {
            psError(PS_ERR_UNKNOWN, false,
                    "Failed to move to extension %d.",
                    lcv);
            return 15;
        }

        // check to see if I can retrieve the number back from the psFits object.
        if (psFitsGetExtNum(fits) != lcv) { // hey, it didn't move?
            psError(PS_ERR_UNKNOWN, false,
                    "Failed to retrieve the extension number back (%d vs %d)",
                    psFitsGetExtNum(fits), lcv);
            return 15;
        }

        // check that the image is associated to the extension moved, i.e.,
        // did we really move to the proper extension?
        psImage* image = psFitsReadImage(NULL, fits,region,0);

        if (abs(image->data.F32[0][0] - (float)lcv) > FLT_EPSILON) {
            psError(PS_ERR_UNKNOWN, true,
                    "The image pixel 0,0 of ext-%d was %g, expected %d.",
                    lcv,image->data.F32[0][0],lcv);
            return 16;
        }
        psFree(image);
    }

    // check to see if given a negative extension number, it errors.
    psLogMsg(__func__,PS_LOG_INFO, "Following should be an error.");
    psErrorClear();
    if (psFitsMoveExtNum(fits, -1, false) || psErrorGetStackSize() != 1) {
        psError(PS_ERR_UNKNOWN, false,
                "Moving to negative HDU didn't fail.");
        return 21;
    }

    // check to see if relative positioning beyond PHU, it errors.
    psFitsMoveExtNum(fits,0,false);
    psLogMsg(__func__,PS_LOG_INFO, "Following should be an error.");
    psErrorClear();
    if (psFitsMoveExtNum(fits, -1, true) || psErrorGetStackSize() != 1) {
        psError(PS_ERR_UNKNOWN, false,
                "Moving to negative HDU didn't fail.");
        return 22;
    }


    // check to see if given a extension greater than the total #HDUs, it errors.
    psLogMsg(__func__,PS_LOG_INFO, "Following should be an error.");
    psErrorClear();
    if (psFitsMoveExtNum(fits, numHDUs, false) || psErrorGetStackSize() != 1) {
        psError(PS_ERR_UNKNOWN, false,
                "Moving to a HDU beyond the file's contents didn't fail.");
        return 31;
    }

    // check to see if relative positioning beyond PHU, it errors.
    psFitsMoveExtNum(fits,numHDUs-1,false);
    psLogMsg(__func__,PS_LOG_INFO, "Following should be an error.");
    psErrorClear();
    if (psFitsMoveExtNum(fits, 1, true) || psErrorGetStackSize() != 1) {
        psError(PS_ERR_UNKNOWN, false,
                "Moving to negative HDU didn't fail.");
        return 32;
    }

    // check to see if given a NULL psFits, it errors.
    psLogMsg(__func__,PS_LOG_INFO, "Following should be an error.");
    psErrorClear();
    if (psFitsMoveExtNum(NULL, 0, false) || psErrorGetStackSize() != 1) {
        psError(PS_ERR_UNKNOWN, false,
                "Operation of NULL psFits didn't fail.");
        return 40;
    }

    psFree(fits);

    return 0;
}



static psS32 tst_psFitsReadHeader( void )
{
    if (! makeMulti() ) {
        return 1;
    }

    psFits* fits = psFitsAlloc(multiFilename);

    if (fits == NULL) {
        psError(PS_ERR_UNKNOWN, false,
                "psFitsAlloc returned NULL on existing file.");
        return 1;
    }

    int numHDUs = psFitsGetSize(fits);

    if (numHDUs < 8) {
        psError(PS_ERR_UNKNOWN,true,
                "The 'multi' FITS file does not have multiple HDUs.");
        return 2;
    }

    char extname[80];
    for (int hdunum = 0; hdunum < numHDUs; hdunum++) {
        snprintf(extname,80,"ext-%d",hdunum);

        psFitsMoveExtNum(fits,hdunum,false);

        psMetadata* header = psFitsReadHeader(NULL,fits);
        if (header == NULL) {
            psError(PS_ERR_UNKNOWN, false,
                    "Failed to read header");
            return 3;
        }

        psMetadata* header2 = psMetadataAlloc();
        header2 = psFitsReadHeader(header2,fits);
        if (header2 == NULL) {
            psError(PS_ERR_UNKNOWN, false,
                    "Failed to read header");
            return 11;
        }

        if (header->list->size < 1 || header->list->size != header2->list->size) {
            psError(PS_ERR_UNKNOWN, true,
                    "Reading the header given a NULL input psMetadata differed "
                    "from giving an existing psMetadata.");
            return 12;
        }

        // check for the extra metadata items
        psS32 intItem = psMetadataLookupS32(header, "MYINT",NULL);
        psF32 fltItem = psMetadataLookupF32(header, "MYFLT",NULL);
        psF64 dblItem = psMetadataLookupF64(header, "MYDBL",NULL);
        psMetadataItem* boolItem = psMetadataLookup(header, "MYBOOL");
        char* strItem = psMetadataLookupPtr(header, "MYSTR",NULL);

        if (intItem != hdunum) {
            psError(PS_ERR_UNKNOWN, true,
                    "Failed to retrieve psS32 metadata item from file.");
            return 20;
        }

        if (fabsf(fltItem - 1.0f/(float)(1+hdunum)) > FLT_EPSILON) {
            psError(PS_ERR_UNKNOWN, true,
                    "Failed to retrieve psF32 metadata item from file.  Got %f vs %f",
                    fltItem,1.0f/(float)(1+hdunum));
            return 21;
        }

        if (abs(dblItem - 1.0/(double)(1+hdunum)) > DBL_EPSILON) {
            psError(PS_ERR_UNKNOWN, true,
                    "Failed to retrieve psF64 metadata item from file.  Got %g vs %g",
                    dblItem, 1.0/(double)(1+hdunum));
            return 22;
        }

        if (boolItem == NULL ||
                boolItem->type != PS_META_PRIMITIVE ||
                boolItem->pType != PS_TYPE_BOOL) {
            psError(PS_ERR_UNKNOWN, true,
                    "Failed to retrieve psBool metadata item from file.");
            return 23;
        }

        if (strItem == NULL || strncmp(strItem,extname,strlen(extname)) != 0) {
            psError(PS_ERR_UNKNOWN, true,
                    "Failed to retrieve string metadata item from file.  Got '%s' vs '%s' (%d)",
                    strItem,extname,strlen(extname));
            return 24;
        }

        psFree(header);
        psFree(header2);
    }

    psLogMsg(__func__,PS_LOG_INFO,"following should be an error (input psFits = NULL)");
    psMetadata* header = psFitsReadHeader(NULL,NULL);

    if (header != NULL || psErrorGetStackSize() != 1) {
        psError(PS_ERR_UNKNOWN, true,
                "psFitsReadHeader didn't error on a NULL psFits.");
        return 30;
    }


    psFree(fits);

    return 0;
}

static psS32 tst_psFitsReadHeaderSet( void )
{
    if (! makeMulti() ) {
        return 1;
    }

    psFits* fits = psFitsAlloc(multiFilename);

    if (fits == NULL) {
        psError(PS_ERR_UNKNOWN, false,
                "psFitsAlloc returned NULL on existing file.");
        return 1;
    }

    int numHDUs = psFitsGetSize(fits);

    if (numHDUs < 8) {
        psError(PS_ERR_UNKNOWN,true,
                "The 'multi' FITS file does not have multiple HDUs.");
        return 2;
    }

    // move to the middle
    psFitsMoveExtNum(fits,numHDUs/2, false);

    psHash* headerSet = psFitsReadHeaderSet(NULL, fits);

    if (headerSet == NULL) {
        psError(PS_ERR_UNKNOWN, false,
                "psFitsReadHeaderSet returned NULL unexpectedly.");
        return 3;
    }

    if (psFitsGetExtNum(fits) != numHDUs/2) {
        psError(PS_ERR_UNKNOWN, false,
                "psFitsReadHeaderSet changed the CHU.");
        return 4;
    }

    char extname[80];
    for (int i = 0; i < numHDUs; i++) {
        if (i == 0) {
            snprintf(extname, 80, "PHU");
        } else {
            snprintf(extname, 80, "ext-%d", i);
        }

        psMetadata* header = psHashLookup(headerSet, extname);

        if (header == NULL) {
            psError(PS_ERR_UNKNOWN, false,
                    "psFitsReadHeader returned NULL unexpectedly for HDU#%d.",
                    i);
            return 5;
        }

        psS32 intItem = psMetadataLookupS32(header, "MYINT",NULL);

        if (intItem != i) {
            psError(PS_ERR_UNKNOWN, false,
                    "psFitsReadHeader for HDU#%d had a MYINT of %d, expected %d.",
                    intItem, i);
            return 6;
        }

    }

    psHash* set2 = psHashAlloc(20);

    if (psFitsReadHeaderSet(set2,NULL) != NULL) {
        psError(PS_ERR_UNKNOWN, false,
                "psFitsReadHeaderSet returned non-NULL given a NULL psFits.");
        return 10;
    }

    psHash* set3 = psHashAlloc(20);
    psFitsReadHeaderSet(set3,fits);

    if (set3 == NULL) {
        psError(PS_ERR_UNKNOWN, false,
                "psFitsReadHeaderSet returned NULL unexpectedly.");
        return 11;
    }

    for (int i = 0; i < numHDUs; i++) {
        if (i == 0) {
            snprintf(extname, 80, "PHU");
        } else {
            snprintf(extname, 80, "ext-%d", i);
        }

        psMetadata* header = psHashLookup(set3, extname);

        if (header == NULL) {
            psError(PS_ERR_UNKNOWN, false,
                    "psFitsReadHeader returned NULL unexpectedly for HDU#%d.",
                    i);
            return 5;
        }

        psS32 intItem = psMetadataLookupS32(header, "MYINT",NULL);

        if (intItem != i) {
            psError(PS_ERR_UNKNOWN, false,
                    "psFitsReadHeader for HDU#%d had a MYINT of %d, expected %d.",
                    intItem, i);
            return 6;
        }

    }



    psFree(headerSet);
    psFree(set3);
    // set2 should have been freed by psFitsReadHeaderSet upon an error.

    psFree(fits);

    return 0;
}

static psS32 tst_psFitsReadTable( void )
{


    if (! makeTable() ) {
        return 1;
    }

    psFits* fits = psFitsAlloc(tableFilename);

    if (fits == NULL) {
        psError(PS_ERR_UNKNOWN, false,
                "psFitsAlloc returned NULL on existing file.");
        return 1;
    }

    psFitsMoveExtNum(fits,1,false);

    psArray* table = psFitsReadTable(fits);

    if (table == NULL) {
        psError(PS_ERR_UNKNOWN, false,
                "psFitsReadTable returned NULL unexpectedly.");
        return 2;
    }

    if (table->n != tableNumRows) {
        psError(PS_ERR_UNKNOWN, false,
                "Expected %d rows, but read %d.",
                tableNumRows, table->n);
        return 3;
    }


    for (int row = 0; row < table->n; row++) {
        psMetadata* rowData = table->data[row];

        psS32 intItem = psMetadataLookupS32(rowData, "MYINT",NULL);
        psF32 fltItem = psMetadataLookupF32(rowData, "MYFLT",NULL);
        psF64 dblItem = psMetadataLookupF64(rowData, "MYDBL",NULL);
        psBool boolItem = psMetadataLookupBool(rowData, "MYBOOL",NULL);

        if (intItem != row) {
            psError(PS_ERR_UNKNOWN, true,
                    "Failed to retrieve psS32 metadata item from file (row=%d).  Got %d vs %d",
                    row, intItem, row);
            return 20;
        }

        if (fabsf(fltItem - 1.0f/(float)(1+row)) > FLT_EPSILON) {
            psError(PS_ERR_UNKNOWN, true,
                    "Failed to retrieve psF32 metadata item from file (row=%d).  Got %f vs %f",
                    row, fltItem,1.0f/(float)(1+row));
            return 21;
        }

        if (abs(dblItem - 1.0/(double)(1+row)) > DBL_EPSILON) {
            psError(PS_ERR_UNKNOWN, true,
                    "Failed to retrieve psF64 metadata item from file (row=%d).  Got %g vs %g",
                    row, dblItem, 1.0/(double)(1+row));
            return 22;
        }

        if ( boolItem != ((row&0x01) == 0)) {
            psError(PS_ERR_UNKNOWN, true,
                    "Failed to retrieve psBool metadata item from file (row=%d). Got %d vs %d",
                    row, boolItem, ((row&0x01) == 0));
            return 23;
        }

    }

    psFree(table);
    psFree(fits);

    psArray* nullTest = psFitsReadTable(NULL);

    if (nullTest != NULL || psErrorGetStackSize() != 1) {
        psError(PS_ERR_UNKNOWN, true,
                "psFitsReadTable returned non-NULL when given NULL.");
        return 30;
    }

    return 0;
}

psS32 testImageRead(void)
{
    psS32 N = 256;
    psS32 M = 128;

    /*
        This function shall open the specified FITS file, read the specified data
        and place the data into a psImage structure. This function shall generate
        an error message and return NULL if any of the input parameters are out of
        range, image file doesn't exist or image is zero or one dimensional.

        Verify the returned psImage structure contains expected values, if the input
        parameter filename specifies an available FITS file with known 2dimensional
        data, input parameters col, row, ncol, nrow, z specify data range with the
        FITS file. Cases should include 1x1, Nx1, 1xN, NxN and MxN sub images and
        total FITS file image. (done in macro)

        Verify the returned psImage structure is equal to the input parameter
        'output', if specified. (done in macro)

        */

    /* generate FITS file to read */

    #define testReadTypeSize(m, n, readM0, readN0, readM, readN, TYP, filename) \
    { \
        psImage* img = NULL; \
        psImage* img2 = NULL; \
        psImage* img3 = NULL; \
        psImage* img4 = NULL; \
        /*        psImagimge* img_ref = NULL; */ \
        \
        GENIMAGE(img,m,n,TYP,row+2*col); \
        img2 = psImageCopy(img2,img,PS_TYPE_##TYP); \
        GENIMAGE(img3,m,n,TYP,row+2*col); \
        psImageClip(img3,32.0,32.0,120.0,120.0); \
        img4 = psImageCopy(img4,img3,PS_TYPE_##TYP); \
        remove(filename); \
        psFits* fits = psFitsAlloc(filename); \
        psRegion region = {0,0,0,0}; \
        if (! psFitsWriteImage(fits, NULL, img, 2, "primary")) { \
            psError(PS_ERR_UNKNOWN, true,"Failed to write test image %s",filename); \
            return 1; \
        } \
        if (! psFitsUpdateImage(fits,img3, region, 1)) { \
            psError(PS_ERR_UNKNOWN, true,"Failed to write test image %s",filename); \
            return 2; \
        } \
        if (! psFitsWriteImage(fits,NULL, img3, 2, "extension")) { \
            psError(PS_ERR_UNKNOWN, true,"Failed to write test image %s",filename); \
            return 3; \
        } \
        if (! psFitsUpdateImage(fits,img,region, 1)) { \
            psError(PS_ERR_UNKNOWN, true,"Failed to write test image %s",filename); \
            return 4; \
        } \
        psFree(img); \
        psFree(fits); \
        img = NULL; \
        psFree(img3); \
        img3 = NULL; \
        fits = psFitsAlloc(filename); \
        psRegion reg = {readM0, readM, readN0, readN}; \
        img = psFitsReadImage(img, fits, reg, 0); \
        img3 = psFitsReadImage(img3, fits, reg, 1); \
        if (img3 == NULL) { \
            psError(PS_ERR_UNKNOWN, true,"Failed to read test image %s",filename); \
            return 6; \
        } \
        for (psU32 row = readN0; row < readN; row++) { \
            ps##TYP* imgRow = img->data.TYP[row-readN0]; \
            ps##TYP* img2Row = img2->data.TYP[row]; \
            ps##TYP* img3Row = img3->data.TYP[row-readN0]; \
            ps##TYP* img4Row = img4->data.TYP[row]; \
            for (psU32 col = readM0; col < readM; col++) { \
                if (fabsf(imgRow[col-readM0]-img2Row[col]) > FLT_EPSILON) { \
                    psError(PS_ERR_UNKNOWN, true,"Image changed in I/O operation at %d,%d,0 (%.2f vs %.2f) for %s", \
                            col,row,(psF32)imgRow[col-readM0],(psF32)img2Row[col],filename); \
                    return 7; \
                } \
                if (fabsf(img3Row[col-readM0]-img4Row[col]) > FLT_EPSILON) { \
                    psError(PS_ERR_UNKNOWN, true,"Image changed in I/O operation at %d,%d,1 (%.2f vs %.2f) for %s", \
                            col,row,(psF32)img3Row[col-readM0],(psF32)img4Row[col],filename); \
                    return 8; \
                } \
            } \
        } \
        psFree(img); \
        img = NULL; \
        psFree(img3); \
        img3 = NULL; \
        psFitsMoveExtNum(fits,1, false); \
        img3 = psFitsReadImage(img3, fits, reg, 0); \
        img = psFitsReadImage(img, fits, reg, 1); \
        if (img == NULL) { \
            psError(PS_ERR_UNKNOWN, true,"Failed to read test image %s",filename); \
            return 9; \
        } \
        for (psU32 row = readN0; row < readN; row++) { \
            ps##TYP* imgRow = img->data.TYP[row-readN0]; \
            ps##TYP* img2Row = img2->data.TYP[row]; \
            ps##TYP* img3Row = img3->data.TYP[row-readN0]; \
            ps##TYP* img4Row = img4->data.TYP[row]; \
            for (psU32 col = readM0; col < readM; col++) { \
                if (fabsf(imgRow[col-readM0]-img2Row[col]) > FLT_EPSILON) { \
                    psError(PS_ERR_UNKNOWN, true,"Image changed in I/O operation at %d,%d,0 (%.2f vs %.2f) for %s", \
                            col,row,(psF32)imgRow[col-readM0],(psF32)img2Row[col],filename); \
                    return 10; \
                } \
                if (fabsf(img3Row[col-readM0]-img4Row[col]) > FLT_EPSILON) { \
                    psError(PS_ERR_UNKNOWN, true,"Image changed in I/O operation at %d,%d,1 (%.2f vs %.2f) for %s", \
                            col,row,(psF32)img3Row[col-readM0],(psF32)img4Row[col],filename); \
                    return 11; \
                } \
            } \
        } \
        psFree(img); \
        psFree(img2); \
        psFree(img3); \
        psFree(img4); \
        psFree(fits); \
    }

    #define testReadType(TYP,filename) \
    testReadTypeSize(1,1,0,0,0,0,TYP,"tmpImages/1x1_" filename); \
    testReadTypeSize(M,1,M/4,0,M*3/4,0,TYP,"tmpImages/Mx1_" filename); \
    testReadTypeSize(1,N,0,N/4,0,N*3/4,TYP,"tmpImages/1xN_" filename); \
    testReadTypeSize(M,N,M/4,N/4,M*3/4,N*3/4,TYP,"tmpImages/MxN_" filename);

    mkdir("tmpImages",0777);

    testReadTypeSize(1,1,0,0,0,0,U8,"tmpImages/1x1_" "U8.fits");
    testReadTypeSize(M,1,M/4,0,M*3/4,0,U8,"tmpImages/Mx1_" "U8.fits");
    testReadTypeSize(1,N,0,N/4,0,N*3/4,U8,"tmpImages/1xN_" "U8.fits");
    testReadTypeSize(M,N,M/4,N/4,M*3/4,N*3/4,U8,"tmpImages/MxN_" "U8.fits");


    testReadType(U8,"U8.fits");
    testReadType(S8,"S8.fits");   // Not a requirement
    testReadType(S16,"S16.fits");
    testReadType(U16,"U16.fits"); // Not a requirement
    testReadType(S32,"S32.fits");
    testReadType(U32,"U32.fits"); // Not a requirement
    testReadType(F32,"F32.fits");
    testReadType(F64,"F64.fits");

    return 0;
}

psS32 testImageWrite(void)
{
    psImage* img = NULL;
    psImage* img2 = NULL;
    psS32 m = 64;
    psS32 n = 96;

    /*
    This function shall write the specified section within a psImage structure
    to a FITS file. If the specifiedfile exists, then data should overwrite the
    section to write. If the specified file doesn't exist, it shall be created.
    If an extenstion is specified, then a basic primary header data unit shall
    be created.
    */

    /*
    Verify a FITS file named filename is generated and contains expected
    values, if the input parameter input contains known data values, input
    parameters col, row, ncol, nrow specify a valid data region within psImage
    structure.

    Verify a FITS file named filename is generated and contains a primary
    header data unit with extension with expected values, if the input
    parameter input contains known data values, input parameters col, row,
    ncol, nrow specify a valid data region within psImage structure and
    extname and/or extnum specify an extenstion to write.

    N.B. : these are done in testImageRead tests, see above.
    */

    /*
    Verify a FITS file named filename is overwritten and contains
    expected values, if the input parameter input contains known data values,
    input parameters col, row, ncol, nrow specify a valid data region within
    psImage structure.
    */

    GENIMAGE(img,m,n,F32,0);
    GENIMAGE(img2,m,n,F32,row+2*col);
    mkdir("tmpImages",0777);
    remove
        ("tmpImages/writeTest.fits");
    psRegion region = {
                          0,0,0,0
                      };
    psFits* fits = psFitsAlloc("tmpImages/writeTest.fits");

    if (! psFitsWriteImage(fits, NULL, img,1,NULL)) {
        psError(PS_ERR_UNKNOWN, true,"Couldn't write writeTest.fits.");
        return 14;
    }
    if (! psFitsUpdateImage(fits, img2, region, 0)) {
        psError(PS_ERR_UNKNOWN, true,"Couldn't update writeTest.fits.");
        return 15;
    }
    psFree(img);
    psFree(img2);

    // Did it really overwrite the pixel values?  Let's read it in and see.
    psFree(fits);
    fits = psFitsAlloc("tmpImages/writeTest.fits");
    img = psFitsReadImage(NULL, fits, region, 0);
    if (img == NULL) {
        psError(PS_ERR_UNKNOWN, true,"Could not read in writeTest.fits.");
        return 16;
    }
    for (psU32 row=0;row<n;row++) {
        psF32* imgRow = img->data.F32[row];
        for (psU32 col=0;col<m;col++) {
            if (fabsf(imgRow[col] - (row+2*col)) > FLT_EPSILON) {
                psError(PS_ERR_UNKNOWN, true,"The image values were not overwritten at %d,%d (%.2f vs %.2f)",
                        col,row,imgRow[col],(row+2*col));
                return 17;
            }
        }
    }

    psFree(img);

    /*
    Verify false is returned and program execution is not stopped, if the input image
    is null.
    */
    psLogMsg(__func__,PS_LOG_INFO,"Following should generate an error message because input image is null.");
    if ( psFitsWriteImage(fits,NULL,NULL, 1, NULL) ) {
        psError(PS_ERR_UNKNOWN, true,"psImageWriteSection did not return false when input image is NULL.");
        return 20;
    }

    psFree(fits);

    return 0;
}
