/** @file  psFits.c
 *
 *  @brief Contains Fits I/O routines
 *
 *  @ingroup FileIO
 *
 *  @author Robert DeSonia, MHPCC
 *
 *  @version $Revision: 1.50 $ $Name: not supported by cvs2svn $
 *  @date $Date: 2005-09-28 20:02:46 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
 */

#include <unistd.h>

#include "psFits.h"
#include "string.h"
#include "psError.h"
#include "psErrorText.h"
#include "psImageStructManip.h"
#include "psMemory.h"
#include "psString.h"
#include "psLogMsg.h"
#include "psTrace.h"
#include "psVector.h"

#define MAX_STRING_LENGTH 256  // maximum length string for FITS routines

// list of FITS header keys to ignore.
static char* standardFitsKeys[] = {
                                      NULL
                                  };

static psElemType convertFitsToPsType(int datatype)
{
    switch (datatype) {
    case TBYTE:
        return PS_TYPE_U8;
    case TSBYTE:
        return PS_TYPE_S8;
    case TSHORT:
        return PS_TYPE_S16;
    case TUSHORT:
        return PS_TYPE_U16;
    case TLONG:
        if (sizeof(long) == 8) {
            return PS_TYPE_S64;
        }
        // no break
    case TINT:
        return PS_TYPE_S32;
    case TULONG:
        if (sizeof(unsigned long) == 8) {
            return PS_TYPE_U64;
        }
        // no break
    case TUINT:
        return PS_TYPE_U32;
    case TLONGLONG:
        return PS_TYPE_S64;
    case TFLOAT:
        return PS_TYPE_F32;
    case TDOUBLE:
        return PS_TYPE_F64;
    case TCOMPLEX:
        return PS_TYPE_C32;
    case TDBLCOMPLEX:
        return PS_TYPE_C64;
    case TLOGICAL:
        return PS_TYPE_BOOL;
    default:
        psError(PS_ERR_IO, true,
                "Unknown FITS datatype, %d.",
                datatype);
        return 0;
    }
}

static bool convertPsTypeToFits(int type, int* bitPix, double* bZero, int* dataType)
{

    int bitpix;
    int datatype;
    double bzero = 0.0;

    switch (type) {

    case PS_TYPE_U8:
        bitpix = BYTE_IMG;
        datatype = TBYTE;
        break;

    case PS_TYPE_S8:
        bitpix = BYTE_IMG;
        bzero = INT8_MIN;
        datatype = TSBYTE;
        break;

    case PS_TYPE_U16:
        bitpix = SHORT_IMG;
        bzero = -1.0 * INT16_MIN;
        datatype = TUSHORT;
        break;

    case PS_TYPE_S16:
        bitpix = SHORT_IMG;
        datatype = TSHORT;
        break;

    case PS_TYPE_U32:
        bitpix = LONG_IMG;
        bzero = -1.0 * INT32_MIN;
        datatype = TUINT;
        break;

    case PS_TYPE_S32:
        bitpix = LONG_IMG;
        datatype = TINT;
        break;

    case PS_TYPE_F32:
        bitpix = FLOAT_IMG;
        datatype = TFLOAT;
        break;

    case PS_TYPE_F64:
        bitpix = DOUBLE_IMG;
        datatype = TDOUBLE;
        break;

    case PS_DATA_STRING:
        bitpix = BYTE_IMG;
        datatype = TSTRING;
        break;

    case PS_DATA_BOOL:
        bitpix = BYTE_IMG;
        datatype = TLOGICAL;
        break;

    default: {
            char* typeStr;
            PS_TYPE_NAME(typeStr,type);
            psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                    PS_ERRORTEXT_psFits_TYPE_UNSUPPORTED,
                    typeStr);
            return false;
        }
    }

    // pass the requested parameters  (NULL parameters are not set, of course).
    if (bitPix != NULL) {
        *bitPix = bitpix;
    }

    if (dataType != NULL) {
        *dataType = datatype;
    }

    if (bZero != NULL) {
        *bZero = bzero;
    }

    return true;
}

static char getTForm(int type)
{
    switch (type) {
    case PS_TYPE_U8:
    case PS_TYPE_S8:
        return 'B';
    case PS_TYPE_S16:
        return 'I';
    case PS_TYPE_S32:
        return 'J';
    case PS_TYPE_U64:
    case PS_TYPE_S64:
        return 'K';
    case PS_TYPE_U16:
        return 'U';
    case PS_TYPE_U32:
        return 'V';
    case PS_TYPE_F32:
        return 'E';
    case PS_TYPE_F64:
        return 'D';
    case PS_TYPE_C32:
        return 'C';
    case PS_TYPE_C64:
        return 'M';
    case PS_TYPE_BOOL:
        return 'L';
    case PS_DATA_STRING:
        return 'A';
    default:
        return '?';
    }
}


static bool getMetadataTForm(psMetadataItem* item, char** fitsType, int repeat)
{
    if (item == NULL) {
        return false;
    }

    *fitsType = NULL;
    char tform = getTForm(item->type);

    switch (item->type) {
    case PS_TYPE_U8:
    case PS_TYPE_S8:
    case PS_TYPE_S16:
    case PS_TYPE_S32:
    case PS_TYPE_U64:
    case PS_TYPE_S64:
    case PS_TYPE_U16:
    case PS_TYPE_U32:
    case PS_TYPE_F32:
    case PS_TYPE_F64:
    case PS_TYPE_C32:
    case PS_TYPE_C64:
    case PS_TYPE_BOOL:
        psStringAppend(fitsType, "%d%c",repeat,tform);
        break;

    case PS_DATA_STRING:
        psStringAppend(fitsType,"%dA",
                       strlen(item->data.V)*repeat);
        break;

    case PS_DATA_VECTOR: {
            psVector* vec = item->data.V;
            tform = getTForm(vec->type.type);
            psStringAppend(fitsType,"%d%c",repeat*vec->n,tform);
        }
        break;

    default:
        return false;
    }

    return true;
}

static bool isHDUEmpty(const psFits* fits)
{
    /* check for keys - no keys means this is really an empty HDU */
    int keysexist = -1;
    int morekeys;
    int status = 0;

    fits_get_hdrspace(fits->fd, &keysexist, &morekeys, &status);

    // if no keys exist and not primary HDU, this really is an empty HDU
    if (keysexist == 0) {
        return true;
    }

    return false;

}

static void fitsFree(psFits* fits)
{
    int status = 0;

    if (fits != NULL) {
        (void)fits_close_file(fits->fd, &status);
        psFree(fits->filename);
    }
}

psFits* psFitsAlloc(const char* name)
{
    int status = 0;
    fitsfile *fptr = NULL;      /* Pointer to the FITS file */

    if (name == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_FILENAME_NULL);
        return NULL;
    }

    /* Open/Create the FITS file */
    if (access(name, F_OK) == 0) {     // file exists
        (void)fits_open_file(&fptr, name, READWRITE, &status);
        if (fptr == NULL) { // if failed, try openning as just read-only
            status = 0;
            (void)fits_open_file(&fptr, name, READONLY, &status);
        }
        if (fptr == NULL || status != 0) {
            char fitsErr[MAX_STRING_LENGTH];
            fits_get_errstatus(status, fitsErr);
            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                    PS_ERRORTEXT_psFits_FILENAME_INVALID,
                    name, fitsErr);
            return NULL;
        }
    } else {  // file does not exist, so create.
        (void)fits_create_file(&fptr, name, &status);
        if (fptr == NULL || status != 0) {
            char fitsErr[MAX_STRING_LENGTH];
            fits_get_errstatus(status, fitsErr);
            psError(PS_ERR_IO, true,
                    PS_ERRORTEXT_psFits_FILENAME_CREATE_FAILED,
                    name, fitsErr);
            return NULL;
        }
    }

    psFits* fits = psAlloc(sizeof(psFits));
    fits->filename = psAlloc(strlen(name)+1);
    fits->fd = fptr;
    strcpy((char*)fits->filename,name);
    psMemSetDeallocator(fits,(psFreeFunc)fitsFree);

    return fits;
}


bool psMemCheckFits(psPtr ptr)
{
    return ( psMemGetDeallocator(ptr) == (psFreeFunc)fitsFree );
}


bool psFitsMoveExtName(const psFits* fits,
                       const char* extname)
{
    int status = 0;

    if (fits == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_NULL);
        return false;
    }

    if (extname == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_EXTNAME_NULL);
        return false;
    }


    if (fits_movnam_hdu(fits->fd, ANY_HDU, (char*)extname, 0, &status) != 0) {
        char fitsErr[MAX_STRING_LENGTH];
        fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_LOCATION_INVALID, true,
                PS_ERRORTEXT_psFits_EXTNAME_INVALID,
                extname, fits->filename, fitsErr);
        return false;
    }

    return true;
}

bool psFitsMoveExtNum(const psFits* fits,
                      int extnum,
                      bool relative)
{
    if (fits == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_NULL);
        return false;
    }

    int status = 0;
    int hdutype = 0;

    if (relative) {
        fits_movrel_hdu(fits->fd, extnum, &hdutype, &status);
        if (status != 0) {
            char fitsErr[MAX_STRING_LENGTH];
            fits_get_errstatus(status, fitsErr);
            psError(PS_ERR_LOCATION_INVALID, true,
                    PS_ERRORTEXT_psFits_EXTNUM_REL_MOVE_FAILED,
                    extnum, fits->filename, fitsErr);
            return false;
        }
    } else {
        fits_movabs_hdu(fits->fd, extnum+1, &hdutype, &status);
        if (status != 0) {
            char fitsErr[MAX_STRING_LENGTH];
            fits_get_errstatus(status, fitsErr);
            psError(PS_ERR_LOCATION_INVALID, true,
                    PS_ERRORTEXT_psFits_EXTNUM_ABS_MOVE_FAILED,
                    extnum, fits->filename, fitsErr);
            return false;
        }
    }

    return true;
}

int psFitsGetExtNum(const psFits* fits)
{
    int hdunum;

    if (fits == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_NULL);
        return PS_FITS_TYPE_NONE;
    }


    return fits_get_hdu_num(fits->fd,&hdunum) - 1;
}

psString psFitsGetExtName(const psFits* fits)
{
    if (fits == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_NULL);
        return NULL;
    }

    int status = 0;
    char name[MAX_STRING_LENGTH];

    if (fits_read_key_str(fits->fd, "EXTNAME", name, NULL, &status) != 0) {
        status = 0;
        if (fits_read_key_str(fits->fd, "HDUNAME", name, NULL, &status) != 0) {
            int num = psFitsGetExtNum(fits);
            snprintf(name, MAX_STRING_LENGTH, "EXT-%3d",num);
        }
    }
    return psStringCopy(name);
}

bool psFitsSetExtName(psFits* fits, const char* name)
{
    if (fits == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_NULL);
        return false;
    }

    if (name == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_EXTNAME_NULL);
        return false;
    }

    int status = 0;

    if (fits_update_key_str(fits->fd, "EXTNAME", (char*)name, NULL, &status) != 0) {
        char fitsErr[MAX_STRING_LENGTH];
        (void)fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_WRITE_FAILED,
                fits->filename, fitsErr);
        return false;
    }

    return true;
}

int psFitsGetSize(const psFits* fits)
{
    if (fits == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_NULL);
        return 0;
    }

    int num = 0;
    int status = 0;

    if (fits_get_num_hdus(fits->fd, &num, &status) != 0) {
        char fitsErr[MAX_STRING_LENGTH];
        fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_LOCATION_INVALID, true,
                PS_ERRORTEXT_psFits_GETNUMHDUS_FAILED,
                fits->filename, fitsErr);
        return 0;
    }

    return num;
}

psFitsType psFitsGetExtType(const psFits* fits)
{
    if (fits == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_NULL);
        return PS_FITS_TYPE_NONE;
    }

    int status = 0;
    int hdutype = PS_FITS_TYPE_NONE;

    if (fits_get_hdu_type(fits->fd, &hdutype, &status) != 0) {
        char fitsErr[MAX_STRING_LENGTH];
        fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_LOCATION_INVALID, true,
                PS_ERRORTEXT_psFits_GETHDUTYPE_FAILED,
                fits->filename, fitsErr);
        return PS_FITS_TYPE_NONE;
    }

    if (hdutype == PS_FITS_TYPE_IMAGE &&
            psFitsGetExtNum(fits) > 0 &&
            isHDUEmpty(fits)) {
        return PS_FITS_TYPE_ANY;
    }

    return hdutype;
}

psMetadata* psFitsReadHeader(psMetadata* out,
                             const psFits* fits)
{
    if (fits == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_NULL);
        return NULL;
    }

    if (out == NULL) {
        out = psMetadataAlloc();
        if (out == NULL) {
            psError(PS_ERR_UNKNOWN, false,
                    "Failed to allocate a new psMetadata container.");
            return NULL;
        }
    }

    // Get number of key names
    int numKeys = 0;
    int keyNum = 0;
    int status = 0;
    fits_get_hdrpos(fits->fd, &numKeys, &keyNum, &status);

    // Get each key name. Keywords start at one.
    char keyType;
    char keyName[MAX_STRING_LENGTH];
    char keyValue[MAX_STRING_LENGTH];
    char keyComment[MAX_STRING_LENGTH];
    psBool tempBool;
    psBool success;
    psBool stdKey;
    for (int i = 1; i <= numKeys; i++) {

        fits_read_keyn(fits->fd, i, keyName, keyValue, keyComment, &status);

        stdKey = false;

        int stdKeyIdx = 0;
        while (standardFitsKeys[stdKeyIdx] != NULL && ! stdKey) {
            if (strcmp(keyName,standardFitsKeys[stdKeyIdx++]) == 0) {
                stdKey = true;
            }
        }

        if (keyValue[0] != 0) { // blank values are not handled by fits_get_keytype
            fits_get_keytype(keyValue, &keyType, &status);
        } else {
            keyType = 'C';
        }
        if (status != 0) {
            break;
        }

        if (! stdKey) {
            switch (keyType) {
            case 'X': // bit
            case 'I': // short int.
            case 'J': // int.
            case 'B': // byte
                success = psMetadataAdd(out,
                                        PS_LIST_TAIL,
                                        keyName,
                                        PS_DATA_S32 | PS_META_DUPLICATE_OK,
                                        keyComment,
                                        atoi(keyValue));
                break;
            case 'U': // unsigned int. may not fit in a psS32
            case 'K': // long int. can't all fit in a psS32
            case 'F':
                success = psMetadataAdd(out,
                                        PS_LIST_TAIL,
                                        keyName,
                                        PS_DATA_F64 | PS_META_DUPLICATE_OK,
                                        keyComment,
                                        atof(keyValue));
                break;
            case 'C':
                // remove the single-quotes at front/end
                if (keyValue[0] == '\'' && keyValue[strlen(keyValue)-1] == '\'') {
                    keyValue[strlen(keyValue)-1] = '\0';
                    success = psMetadataAdd(out,
                                            PS_LIST_TAIL,
                                            keyName,
                                            PS_DATA_STRING | PS_META_DUPLICATE_OK,
                                            keyComment,
                                            keyValue+1);
                } else {
                    success = psMetadataAdd(out,
                                            PS_LIST_TAIL,
                                            keyName,
                                            PS_DATA_STRING | PS_META_DUPLICATE_OK,
                                            keyComment,
                                            keyValue);
                }
                break;
            case 'L':
                tempBool = (keyValue[0] == 'T') ? 1 : 0;
                success = psMetadataAdd(out,
                                        PS_LIST_TAIL,
                                        keyName,
                                        PS_DATA_BOOL | PS_META_DUPLICATE_OK,
                                        keyComment,
                                        tempBool);
                break;
            default:
                psError(PS_ERR_IO, true,
                        PS_ERRORTEXT_psFits_METATYPE_INVALID,
                        keyType);
                return out;
            }

            if (!success) {
                psError(PS_ERR_UNKNOWN, false,
                        PS_ERRORTEXT_psFits_METADATA_ADD_FAILED,
                        keyName);
                return out;
            }
        }

    }

    if ( status != 0) {
        char fitsErr[MAX_STRING_LENGTH];
        (void)fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_METADATA_ADD_FAILED,
                fitsErr);
        return false;
    }

    return out;
}

psMetadata* psFitsReadHeaderSet(psMetadata* out,
                                const psFits* fits)
{
    if (fits == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_NULL);
        psFree(out);
        return NULL;
    }

    if (out == NULL) {
        out = psMetadataAlloc();
        if (out == NULL) {
            psError(PS_ERR_UNKNOWN, false,
                    "Failed to allocate a new psMetadata container.");
            return NULL;
        }
    }

    int size = psFitsGetSize(fits);

    int origPosition = psFitsGetExtNum(fits);

    for (int lcv=0; lcv < size; lcv++) {
        psFitsMoveExtNum(fits, lcv, false);

        char* name = NULL;
        if (lcv == 0) {
            name = psStringCopy("PHU");
        } else {
            name = psFitsGetExtName(fits);
        }

        psMetadata* header = psFitsReadHeader(NULL, fits);
        if (name != NULL && header != NULL) {
            psMetadataAddMetadata(out, PS_LIST_HEAD, name, 0, "FITS Header",
                                  header);
        } else {
            psWarning("Failed to read HDU#%d header data.",
                      lcv);
        }

        psFree(name);
        psFree(header);
    }

    // reposition to the original position
    psFitsMoveExtNum(fits, origPosition, false);

    return out;
}

psImage* psFitsReadImage(psImage* output, // a psImage to recycle.
                         const psFits* fits,    // the psFits object
                         psRegion region, // the region in the FITS image to read
                         int z)           // the z-plane in the FITS image cube to read
{
    psS32 status = 0;           /* CFITSIO file vars */
    psS32 nAxis = 0;
    psS32 anynull = 0;
    psS32 bitPix = 0;           /* Pixel type */
    long nAxes[3];
    long firstPixel[3];         /* lower-left corner of image subset */
    long lastPixel[3];          /* upper-right corner of image subset */
    long increment[3];          /* increment for image subset */
    char fitsErr[80] = "";      /* CFITSIO error message string */
    psS32 fitsDatatype = 0;
    psS32 datatype = 0;

    if (fits == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_NULL);
        psFree(output);
        return NULL;
    }

    // check to see if we even are positioned on an image HDU
    int hdutype;
    if ( fits_get_hdu_type(fits->fd,&hdutype, &status) != 0) {
        char fitsErr[MAX_STRING_LENGTH];
        (void)fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_GET_HDU_TYPE_FAILED,
                fitsErr);
        return NULL;
    }
    if (hdutype != IMAGE_HDU) {
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_NOT_IMAGE_TYPE);
        return NULL;
    }

    /* Get the data type 'bitPix' from the FITS image */
    if (fits_get_img_equivtype(fits->fd, &bitPix, &status) != 0) {
        fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_DATATYPE_UNKNOWN,
                fitsErr);
        psFree(output);
        return NULL;
    }

    /* Get the dimensions 'nAxis' from the FITS image */
    if (fits_get_img_dim(fits->fd, &nAxis, &status) != 0) {
        (void)fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_IMAGE_DIM_UNKNOWN,
                fitsErr);
        psFree(output);
        return NULL;
    }

    /* Validate the number of axis */
    if ((nAxis < 2) || (nAxis > 3)) {
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_IMAGE_DIMENSION_UNSUPPORTED,
                nAxis);
        psFree(output);
        return NULL;
    }

    /* Get the Image size from the FITS file */
    if (fits_get_img_size(fits->fd, nAxis, nAxes, &status) != 0) {
        (void)fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_IMAGE_SIZE_UNKNOWN,
                fitsErr);
        psFree(output);
        return NULL;
    }

    firstPixel[0] = region.x0 + 1;
    firstPixel[1] = region.y0 + 1;
    firstPixel[2] = z + 1;

    if (region.x1 > 0) {
        lastPixel[0] = region.x1;
    } else {
        lastPixel[0] = nAxes[0] + region.x1; // n.b., region.x1 < 0
    }
    if (region.y1 > 0) {
        lastPixel[1] = region.y1;
    } else {
        lastPixel[1] = nAxes[1] + region.y1; // n.b., region.y1 < 0
    }
    lastPixel[2] = z + 1;

    increment[0] = 1;
    increment[1] = 1;
    increment[2] = 1;

    switch (bitPix) {
    case BYTE_IMG:
        datatype = PS_TYPE_U8;
        fitsDatatype = TBYTE;
        break;
    case SBYTE_IMG:
        datatype = PS_TYPE_S8;
        fitsDatatype = TSBYTE;
        break;
    case USHORT_IMG:
        datatype = PS_TYPE_U16;
        fitsDatatype = TUSHORT;
        break;
    case SHORT_IMG:
        datatype = PS_TYPE_S16;
        fitsDatatype = TSHORT;
        break;
    case ULONG_IMG:
        datatype = PS_TYPE_U32;
        fitsDatatype = TUINT;
        break;
    case LONG_IMG:
        datatype = PS_TYPE_S32;
        fitsDatatype = TINT;
        break;
    case LONGLONG_IMG:
        datatype = PS_TYPE_S64;
        fitsDatatype = TLONGLONG;
        break;
    case FLOAT_IMG:
        datatype = PS_TYPE_F32;
        fitsDatatype = TFLOAT;
        break;
    case DOUBLE_IMG:
        datatype = PS_TYPE_F64;
        fitsDatatype = TDOUBLE;
        break;
    default:
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_FITS_TYPE_UNSUPPORTED,
                bitPix);
        psFree(output);
        return NULL;
    }

    output = psImageRecycle(output,
                            lastPixel[0]-firstPixel[0]+1,
                            lastPixel[1]-firstPixel[1]+1,
                            datatype);

    if (output == NULL) {
        psError(PS_ERR_UNKNOWN, false,
                "Failed to allocate a properly sized image.");
        return false;
    }

    // n.b., this assumes contiguous image buffer
    if (fits_read_subset(fits->fd, fitsDatatype, firstPixel, lastPixel, increment,
                         NULL, output->data.V[0], &anynull, &status) != 0) {
        psFree(output);
        (void)fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_READ_FAILED,
                fitsErr);
        return NULL;
    }

    return output;

}

bool psFitsWriteImage(psFits* fits,
                      psMetadata* header,
                      const psImage* input,
                      int numZPlanes)
{

    if (fits == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_NULL);
        return false;
    }

    if (input == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_IMAGE_NULL);
        return false;
    }
    int numCols = input->numCols;
    int numRows = input->numRows;

    int status = 0;

    // determine the FITS-equivalent parameters
    int bitPix;
    double bZero;
    int dataType;
    if (! convertPsTypeToFits(input->type.type, &bitPix, &bZero, &dataType) ) {
        return false;
    }

    int naxis = 3;
    long naxes[3];

    naxes[0] = numCols;
    naxes[1] = numRows;
    naxes[2] = numZPlanes;

    if (numZPlanes < 2) {
        naxis = 2;
    }

    fits_create_img(fits->fd, bitPix, naxis, naxes, &status);

    if (bZero != 0) {        // set the bscale/bzero
        fits_write_key_dbl(fits->fd, "BZERO", bZero, 12, "Pixel Value Offset", &status);
        fits_write_key_dbl(fits->fd, "BSCALE", 1.0, 12, "Pixel Value Scale", &status);
        fits_set_bscale(fits->fd, 1.0, bZero, &status);
    }

    // write the header, if any.
    if (header != NULL) {
        psFitsWriteHeader(header, (psPtr)fits);
    }

    if (input->parent == NULL) { // if no parent, assume that the image data is contiguous
        fits_write_img(fits->fd,
                       dataType,              // datatype
                       1,                     // writing to the first z-plane
                       numCols*numRows,       // number of elements to write, i.e., the whole image
                       input->data.V[0],      // the data
                       &status);
    } else { // image data may not be contiguous; write one row at a time
        int firstPixel = 1;
        for (int row = 0; row < numRows; row++) {
            fits_write_img(fits->fd,
                           dataType,          // datatype
                           firstPixel,
                           numCols,           // number of elements to write, i.e., one row's worth
                           input->data.V[row],// the raw row data
                           &status);
            firstPixel += numCols;  // move to next row
        }
    }

    if ( status != 0) {
        char fitsErr[MAX_STRING_LENGTH];
        (void)fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_WRITE_FAILED,
                fits->filename, fitsErr);
        return false;
    }

    return true;

}

bool psFitsUpdateImage(psFits* fits,
                       const psImage* input,
                       psRegion region,
                       int z)
{
    int status = 0;

    if (fits == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_NULL);
        return false;
    }

    if (input == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_IMAGE_NULL);
        return false;
    }

    // check to see if we are positioned on an image HDU
    int hdutype;
    if ( fits_get_hdu_type(fits->fd,&hdutype, &status) != 0) {
        char fitsErr[MAX_STRING_LENGTH];
        (void)fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_GET_HDU_TYPE_FAILED,
                fitsErr);
        return NULL;
    }
    if (hdutype != IMAGE_HDU) {
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_NOT_IMAGE_TYPE);
        return NULL;
    }

    int numCols = input->numCols;
    int numRows = input->numRows;

    // determine the FITS-equivalent parameters
    int bitPix;
    double bZero;
    int dataType;
    if (! convertPsTypeToFits(input->type.type, &bitPix, &bZero, &dataType) ) {
        return false;
    }

    //check to see if the HDU has the same datatype
    int fileBitpix;
    int naxis;
    long nAxes[3];
    nAxes[2] = 1;
    fits_get_img_param(fits->fd, 3, &fileBitpix, &naxis, nAxes, &status);

    //check to see if the HDU has the same datatype
    if (bitPix != fileBitpix) {
        char* fitsTypeStr;
        char* imageTypeStr;
        PS_TYPE_NAME(fitsTypeStr,fileBitpix);
        PS_TYPE_NAME(imageTypeStr,input->type.type);
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_IMAGE_UPDATE_TYPE_MISMATCH,
                fitsTypeStr, imageTypeStr);
        return false;
    }

    //check if the HDU has the z-plane requested
    if (z >= nAxes[2]) {
        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
                PS_ERRORTEXT_psFits_FITS_Z_SMALL,
                nAxes[2],z);
        return false;
    }

    // determine the region in the FITS file domain
    long firstPixel[3];
    long lastPixel[3];

    firstPixel[0] = region.x0 + 1;
    firstPixel[1] = region.y0 + 1;
    firstPixel[2] = z + 1;

    if (region.x1 > 0) {
        lastPixel[0] = region.x1;
    } else {
        lastPixel[0] = nAxes[0] + region.x1; // n.b., region.x1 < 0
    }
    if (region.y1 > 0) {
        lastPixel[1] = region.y1;
    } else {
        lastPixel[1] = nAxes[1] + region.y1; // n.b., region.y1 < 0
    }
    lastPixel[2] = z + 1;

    if (firstPixel[0] < 1 || firstPixel[0] > nAxes[0] ||
            firstPixel[1] < 1 || firstPixel[1] > nAxes[1] ||
            lastPixel[0] < 1 || lastPixel[0] > nAxes[0] ||
            lastPixel[1] < 1 || lastPixel[1] > nAxes[1]) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                "Specified region [%d:%d,%d:%d], is not valid given the %dx%d FITS image.",
                region.y0,region.y1-1,region.x0,region.x1-1);
        return false;

    }

    int dx = lastPixel[0] - firstPixel[0];
    int dy = lastPixel[1] - firstPixel[1];
    if (dx > numCols ||
            dy > numRows) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                "The region [%d:%d,%d:%d], is not valid given the input %dx%d image.",
                firstPixel[1]-1,lastPixel[1]-1,
                firstPixel[0]-1,lastPixel[0]-1,
                numCols, numRows);
        return false;
    }

    psImage* subset;
    if (dx != numCols || dy != numRows) {
        // the input image needs to be subsetted
        subset = psImageSubset((psImage*)input, psRegionSet(0,dx+1,0,dy+1));
    } else {
        subset = psMemIncrRefCounter((psImage*)input);
    }

    fits_write_subset(fits->fd, dataType, firstPixel, lastPixel, subset->data.V[0], &status);

    if ( status != 0) {
        char fitsErr[MAX_STRING_LENGTH];
        (void)fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_WRITE_FAILED,
                fits->filename, fitsErr);
        return false;
    }

    return true;
}

bool psFitsWriteHeader(const psMetadata* output,
                       psFits* fits)
{

    if (fits == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_NULL);
        return false;
    }

    if (output == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_METADATA_NULL);
        return false;
    }

    int status = 0;

    //transverse the metadata list and add each key.

    psListIterator* iter = psListIteratorAlloc(output->list,PS_LIST_HEAD,true);
    psMetadataItem* item;
    while ( (item=psListGetAndIncrement(iter)) != NULL ) {
        switch (item->type) {
        case PS_DATA_BOOL: {
                int value = item->data.B;
                fits_update_key(fits->fd,
                                TLOGICAL,
                                item->name,
                                &value,
                                item->comment,
                                &status);
                break;
            }
        case PS_DATA_S32:
            fits_update_key(fits->fd,
                            TINT,
                            item->name,
                            &item->data.S32,
                            item->comment,
                            &status);
            break;
        case PS_DATA_F32:
            fits_update_key(fits->fd,
                            TFLOAT,
                            item->name,
                            &item->data.F32,
                            item->comment,
                            &status);
            break;
        case PS_DATA_F64:
            fits_update_key(fits->fd,
                            TDOUBLE,
                            item->name,
                            &item->data.F64,
                            item->comment,
                            &status);
            break;
        case PS_DATA_STRING:
            fits_update_key(fits->fd,
                            TSTRING,
                            item->name,
                            item->data.V,
                            item->comment,
                            &status);
            break;
        default:  // all other META types are ignored
            break;
        }

        if ( status != 0) {
            char fitsErr[MAX_STRING_LENGTH];
            (void)fits_get_errstatus(status, fitsErr);
            psError(PS_ERR_IO, true,
                    PS_ERRORTEXT_psFits_WRITE_FAILED,
                    fits->filename, fitsErr);
            return false;
        }
    }

    return true;
}

psMetadata* psFitsReadTableRow(const psFits* fits,
                               int row)
{
    long numRows;
    int numCols;
    int status = 0;

    if (fits == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_NULL);
        return NULL;
    }

    // check to see if we even are positioned on a table HDU
    int hdutype;
    fits_get_hdu_type(fits->fd,&hdutype, &status);
    if ( status != 0) {
        char fitsErr[MAX_STRING_LENGTH];
        (void)fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_GET_HDU_TYPE_FAILED,
                fitsErr);
        return NULL;
    }
    if (hdutype != ASCII_TBL && hdutype != BINARY_TBL) {
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_NOT_TABLE_TYPE);
        return NULL;
    }

    // get the size of the FITS table
    fits_get_num_rows(fits->fd, &numRows, &status);
    fits_get_num_cols(fits->fd, &numCols, &status);
    if ( status != 0) {
        char fitsErr[MAX_STRING_LENGTH];
        (void)fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_GET_TABLE_SIZE_FAILED,
                fitsErr);
        return NULL;
    }

    psTrace(".psFits.psFitsReadTableRow",5,"Table size is %ix%i\n",numCols, numRows);
    // the row parameter in the proper range?
    if (row < 0 || row >= numRows) {
        psError(PS_ERR_BAD_PARAMETER_VALUE, true,
                PS_ERRORTEXT_psFits_ROW_INVALID,
                row,numRows);
        return NULL;
    }

    psMetadata* data = psMetadataAlloc();

    int typecode;
    long repeat;
    long width;
    char name[60];
    for (int col = 1; col <= numCols; col++) {
        // get the column name
        if (hdutype == BINARY_TBL) {
            fits_get_bcolparms(fits->fd, col, name,
                               NULL, NULL, NULL, NULL, NULL, NULL, NULL, &status);
        } else {
            fits_get_acolparms(fits->fd, col, name,
                               NULL, NULL, NULL, NULL, NULL, NULL, NULL, &status);
        }
        // get the column type
        fits_get_coltype(fits->fd, col, &typecode, &repeat, &width, &status);

        if (status == 0) {

            #define READ_TABLE_ROW_CASE(FITSTYPE, NATIVETYPE, TYPE) \
        case FITSTYPE: { \
                if (repeat == 1) { \
                    NATIVETYPE value; \
                    int anynul = 0; \
                    fits_read_col(fits->fd, FITSTYPE, col,row+1, \
                                  1, 1, NULL, &value, &anynul, &status); \
                    psTrace(".psFits.psFitsReadTableRow",5,"Column #%i, '%s', is type %i, repeat %i, Value = %g\n", \
                            col, name, typecode, repeat, (double)value); \
                    psMetadataAdd(data,PS_LIST_TAIL, name, \
                                  PS_DATA_##TYPE, \
                                  "", (ps##TYPE)value); \
                } else { \
                    psVector* value = psVectorAlloc(repeat,PS_DATA_##TYPE); \
                    int anynul = 0; \
                    fits_read_col(fits->fd, FITSTYPE, col,row+1, \
                                  1, repeat, NULL, value->data.U8, &anynul, &status); \
                    psMetadataAdd(data,PS_LIST_TAIL, name, \
                                  PS_DATA_VECTOR, \
                                  "", value); \
                    psFree(value); \
                } \
                break; \
            }

            switch (typecode) {
            case TBYTE:
            case TSHORT:
            case TLONGLONG:
                READ_TABLE_ROW_CASE(TLONG, long, S32)
                READ_TABLE_ROW_CASE(TFLOAT, float, F32)
                READ_TABLE_ROW_CASE(TDOUBLE, double, F64)
                READ_TABLE_ROW_CASE(TLOGICAL, bool, BOOL);
            case TSTRING: {
                    char* value = psAlloc(repeat+1);
                    int anynul = 0;
                    fits_read_col(fits->fd, TSTRING, col,row+1,
                                  1, 1, NULL, &value, &anynul, &status);
                    psTrace(".psFits.psFitsReadTableRow",5,"Column #%i, '%s', is type %i, repeat %i, value = %s\n",
                            col, name, typecode, repeat, value);
                    if (anynul == 0) {
                        psMetadataAdd(data,PS_LIST_TAIL, name,
                                      PS_DATA_STRING,
                                      "", value);
                    }
                    psFree(value);
                    break;
                }
            default:
                psWarning("Data type (%d) not supportted for column %d",
                          typecode, col);

                psTrace("psFits.psFitsReadTableRow", 2,
                        "Column %d or row %d was of a non primitive type, %d",
                        col, row, typecode);
            }
        }

        if ( status != 0) {
            char fitsErr[MAX_STRING_LENGTH];
            (void)fits_get_errstatus(status, fitsErr);
            psError(PS_ERR_IO, true,
                    PS_ERRORTEXT_psFits_GET_TABLE_ELEMENT,
                    col,row,fitsErr);
            psFree(data);
            return NULL;
        }

    }

    return data;
}

psArray* psFitsReadTableColumn(const psFits* fits,
                               const char* colname)
{
    int colnum = 0;
    int status = 0;

    if (fits == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_NULL);
        return NULL;
    }

    // check to see if we even are positioned on a table HDU
    int hdutype;
    if ( fits_get_hdu_type(fits->fd,&hdutype, &status) != 0) {
        char fitsErr[MAX_STRING_LENGTH];
        (void)fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_GET_HDU_TYPE_FAILED,
                fitsErr);
        return NULL;
    }
    if (hdutype != ASCII_TBL && hdutype != BINARY_TBL) {
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_NOT_TABLE_TYPE);
        return NULL;
    }

    // find the column by name
    if ( fits_get_colnum(fits->fd, CASESEN, (char*)colname, &colnum, &status) != 0) {
        char fitsErr[MAX_STRING_LENGTH];
        (void)fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_FIND_COLUMN,
                colname, fitsErr);
        return NULL;
    }

    // get the number of rows
    long numRows = 0;
    fits_get_num_rows(fits->fd, &numRows, &status);

    // get the column length.
    int width;
    if ( fits_get_col_display_width(fits->fd, colnum, &width, &status) != 0) {
        char fitsErr[MAX_STRING_LENGTH];
        (void)fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_GET_COLTYPE,
                fitsErr);
        return NULL;
    }

    // allocate the buffers
    psArray* result = psArrayAlloc(numRows);
    for (int row = 0; row < numRows; row++) {
        result->data[row] = psAlloc((width+1)*sizeof(char));
    }
    result->n = numRows;

    fits_read_col_str(fits->fd,
                      colnum,
                      1, // firstrow
                      1, // firestelem
                      numRows,
                      "", // nulstr
                      (char**)result->data,
                      NULL,
                      &status);

    if ( status != 0) {
        char fitsErr[MAX_STRING_LENGTH];
        (void)fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_TABLE_READ_COL,
                fitsErr);
        return NULL;
    }

    return result;
}

psVector* psFitsReadTableColumnNum(const psFits* fits,
                                   const char* colname)
{
    int status = 0;
    int colnum = 0;

    if (fits == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_NULL);
        return NULL;
    }

    // check to see if we even are positioned on a table HDU
    int hdutype;
    if ( fits_get_hdu_type(fits->fd,&hdutype, &status) != 0) {
        char fitsErr[MAX_STRING_LENGTH];
        (void)fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_GET_HDU_TYPE_FAILED,
                fitsErr);
        return NULL;
    }
    if (hdutype != ASCII_TBL && hdutype != BINARY_TBL) {
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_NOT_TABLE_TYPE);
        return NULL;
    }

    // find the column by name
    if ( fits_get_colnum(fits->fd, CASESEN, (char*)colname, &colnum, &status) != 0) {
        char fitsErr[MAX_STRING_LENGTH];
        (void)fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_FIND_COLUMN,
                colname, fitsErr);
        return NULL;
    }

    // get the number of rows
    long numRows = 0;
    fits_get_num_rows(fits->fd,
                      &numRows,
                      &status);

    // get the column datatype.
    int typecode;
    long repeat;
    long width;
    if ( fits_get_eqcoltype(fits->fd, colnum, &typecode, &repeat, &width, &status) != 0) {
        char fitsErr[MAX_STRING_LENGTH];
        (void)fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_GET_COLTYPE,
                fitsErr);
        return NULL;
    }

    psVector* result = psVectorAlloc(numRows, convertFitsToPsType(typecode));

    fits_read_col(fits->fd,
                  typecode,
                  colnum,
                  1 /* firstrow */,
                  1 /* firstelem */,
                  numRows,
                  NULL,
                  (psPtr)(result->data.U8),
                  NULL,
                  &status);

    if ( status != 0) {
        char fitsErr[MAX_STRING_LENGTH];
        (void)fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_TABLE_READ_COL,
                fitsErr);
        return NULL;
    }

    return result;
}


psArray* psFitsReadTable(psFits* fits)
{
    int status = 0;

    if (fits == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_NULL);
        return NULL;
    }

    // check to see if we even are positioned on a table HDU
    int hdutype;
    if ( fits_get_hdu_type(fits->fd,&hdutype, &status) != 0) {
        char fitsErr[MAX_STRING_LENGTH];
        (void)fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_GET_HDU_TYPE_FAILED,
                fitsErr);
        return NULL;
    }
    if (hdutype != ASCII_TBL && hdutype != BINARY_TBL) {
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_NOT_TABLE_TYPE);
        return NULL;
    }

    // get the size of the FITS table
    long numRows = 0;
    fits_get_num_rows(fits->fd, &numRows, &status);
    if ( status != 0) {
        char fitsErr[MAX_STRING_LENGTH];
        (void)fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_GET_TABLE_SIZE_FAILED,
                fitsErr);
        return NULL;
    }


    psArray* table = psArrayAlloc(numRows);

    for (int row = 0; row < numRows; row++) {
        psTrace(".psFits.psFitsReadTable",5,"Reading row %i of %i\n",row, numRows);
        table->data[row] = psFitsReadTableRow(fits,row);
    }

    return table;
}

bool psFitsWriteTable(psFits* fits,
                      const psMetadata* header,
                      const psArray* table)
{
    int status = 0;

    if (fits == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_NULL);
        return false;
    }

    if (table == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_IMAGE_NULL);
        return false;
    }

    int rows = table->n;
    if (rows < 1) {
        // no table data, what can I do?
        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
                PS_ERRORTEXT_psFits_TABLE_EMPTY);
        return false;
    }

    // find all the columns needed
    psArray* columns = psArrayAlloc(((psMetadata*)table->data[0])->list->n);
    columns->n=0;

    // find the unique items in the array of metadata 'rows'
    psMetadataItem* item;
    for (int row=0; row < rows; row++) {
        psMetadata* rowMeta = table->data[row];
        if (rowMeta != NULL) {
            psListIterator* iter = psListIteratorAlloc(rowMeta->list,
                                   PS_LIST_HEAD,true);
            while ( (item=psListGetAndIncrement(iter)) != NULL) {
                if (PS_DATA_IS_PRIMITIVE(item->type) ||
                        item->type == PS_DATA_STRING ||
                        item->type == PS_DATA_VECTOR) {
                    bool found = false;
                    psMetadataItem* fItem = NULL;
                    int n;
                    for (n=0; n < columns->n && ! found; n++) {
                        fItem = (psMetadataItem*)(columns->data[n]);
                        if (strcmp(item->name, fItem->name) == 0) {
                            found = true;
                            break;
                        }
                    }
                    if (! found) {
                        psArrayAdd(columns, columns->nalloc, item);
                    } else if (item->type == PS_DATA_STRING &&
                               strlen(fItem->data.V) < strlen(item->data.V)) {
                        // got to keep the longest string value as to know what size to create the table column
                        psMemDecrRefCounter(fItem);
                        columns->data[n] = psMemIncrRefCounter(item);

                    } else if (item->type == PS_DATA_VECTOR &&
                               ((psVector*)(fItem->data.V))->n < ((psVector*)(item->data.V))->n) {
                        psMemDecrRefCounter(fItem);
                        columns->data[n] = psMemIncrRefCounter(item);
                    }
                } else {
                    // unsupported type -- treating as an error
                    psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                            "Unsupported data type (%d) for Metadata Item '%s' in row %d.",
                            item->type, item->name, row);
                    psFree(iter);
                    psFree(columns);
                    return false;
                }
            }
            psFree(iter);
        }
    }

    if (columns->n == 0) { // no table columns found
        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
                "Did not find any column data to write to a table.");
        psFree(columns);
        return false;
    }

    //create list of column names and types.
    psArray* columnNames = psArrayAlloc(columns->n);
    psArray* columnTypes = psArrayAlloc(columns->n);
    for (int n=0; n < columns->n; n++) {
        columnNames->data[n] = psMemIncrRefCounter(((psMetadataItem*)columns->data[n])->name);
        if ( ! getMetadataTForm((psMetadataItem*)columns->data[n],
                                (char**)&(columnTypes->data[n]),1) ) {
            psError(PS_ERR_UNKNOWN, true,
                    "Failed to determine the FITS data type of '%s' (type=%d).",
                    columnNames->data[n],
                    ((psMetadataItem*)columns->data[n])->type);
            psFree(columns);
            psFree(columnNames);
            psFree(columnTypes);
            return false;
        }
    }

    char* extname = NULL;
    if (header != NULL) {
        extname = psMetadataLookupPtr(NULL, header, "EXTNAME");
        if ( extname == NULL) {
            extname = psMetadataLookupPtr(NULL, header, "HDUNAME");
        }
    }

    fits_create_tbl(fits->fd,
                    BINARY_TBL,
                    table->n, // number of rows in table
                    columns->n, // number of columns in table
                    (char**)columnNames->data, // names of the columns
                    (char**)columnTypes->data, // format of the columns
                    NULL, // physical unit of columns
                    extname, // extension name
                    &status);

    psFree(columnNames);
    psFree(columnTypes);

    // fill in the table elements with data
    for (int n = 0; n < columns->n; n++) {
        int row;
        item = columns->data[n];
        if (PS_DATA_IS_PRIMITIVE(item->type) ) {
            psVector* col = psVectorAlloc(table->n, item->type);
            int dataSize = PSELEMTYPE_SIZEOF(item->type);
            int dataType;
            convertPsTypeToFits(item->type, NULL,NULL,&dataType);
            for (row = 0; row < table->n; row++) {
                psMetadataItem* dataItem = psMetadataLookup(
                                               table->data[row],
                                               item->name);
                memcpy(&col->data.U8[row*dataSize],&dataItem->data,dataSize);
            }
            fits_write_col(fits->fd,
                           dataType,
                           n+1, // column number
                           1, // firstrow
                           1, // firstelem
                           table->n, // nelements
                           col->data.U8,
                           &status);
            psFree(col);
        } else if (item->type == PS_DATA_STRING) {
            psArray* colArray = psArrayAlloc(table->n);
            for (row = 0; row < table->n; row++) {
                colArray->data[row] = psMemIncrRefCounter(
                                          psMetadataLookupStr(NULL,
                                                              table->data[row],
                                                              item->name));
            }
            fits_write_col_str(fits->fd,
                               n+1,
                               1,
                               1,
                               table->n,
                               (char**)colArray->data,
                               &status);
            psFree(colArray);
        } else if (item->type == PS_DATA_VECTOR) {
            psVector* vec = (psVector*)(item->data.V);
            psVector* col = psVectorAlloc(table->n*vec->n, vec->type.type);
            int dataSize = PSELEMTYPE_SIZEOF(vec->type.type)*vec->n;
            int dataType;
            convertPsTypeToFits(vec->type.type, NULL,NULL,&dataType);
            for (row = 0; row < table->n; row++) {
                psMetadataItem* dataItem = psMetadataLookup(
                                               table->data[row],
                                               item->name);
                memcpy(&col->data.U8[row*dataSize],
                       ((psVector*)(dataItem->data.V))->data.U8,
                       dataSize);
            }
            fits_write_col(fits->fd,
                           dataType,
                           n+1, // column number
                           1, // firstrow
                           1, // firstelem
                           table->n*vec->n, // nelements
                           col->data.U8,
                           &status);
            psFree(col);
        }
        if (status != 0) {
            char fitsErr[MAX_STRING_LENGTH];
            (void)fits_get_errstatus(status, fitsErr);
            psError(PS_ERR_IO, true,
                    PS_ERRORTEXT_psFits_WRITE_FAILED,
                    fits->filename, fitsErr);
            psFree(columns);
            return false;
        }

    }

    psFree(columns);

    return true;
}

bool psFitsUpdateTable(psFits* fits,
                       const psMetadata* data,
                       int row)
{
    int status = 0;

    if (fits == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_NULL);
        return false;
    }

    if (data == NULL) {
        psError(PS_ERR_BAD_PARAMETER_NULL, true,
                PS_ERRORTEXT_psFits_IMAGE_NULL);
        return false;
    }

    // check to see if we even are positioned on a table HDU
    int hdutype;
    if ( fits_get_hdu_type(fits->fd,&hdutype, &status) != 0) {
        char fitsErr[MAX_STRING_LENGTH];
        (void)fits_get_errstatus(status, fitsErr);
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_GET_HDU_TYPE_FAILED,
                fitsErr);
        return false;
    }
    if (hdutype != ASCII_TBL && hdutype != BINARY_TBL) {
        psError(PS_ERR_IO, true,
                PS_ERRORTEXT_psFits_NOT_TABLE_TYPE);
        return false;
    }

    psMetadataIterator* iter = psMetadataIteratorAlloc((psPtr)data,PS_LIST_HEAD,NULL);

    psMetadataItem* item;

    while ( (item=psMetadataGetAndIncrement(iter)) != NULL) {
        if (PS_DATA_IS_PRIMITIVE(item->type) ||
                item->type == PS_DATA_BOOL ||
                item->type == PS_DATA_STRING) {
            // operating on primitive data type or string, i.e., not a complex object
            int colnum = 0;

            if ( fits_get_colnum(fits->fd, CASESEN, item->name, &colnum, &status) == 0) {
                // cooresponding column found in table
                int dataType;
                convertPsTypeToFits(item->type, NULL, NULL, &dataType);

                if (fits_write_col(fits->fd, dataType, colnum, row+1, 1, 1, &item->data,&status) != 0) {
                    char fitsErr[MAX_STRING_LENGTH];
                    (void)fits_get_errstatus(status, fitsErr);
                    psError(PS_ERR_IO, true,
                            PS_ERRORTEXT_psFits_WRITE_FAILED,
                            fits->filename, fitsErr);
                    psFree(iter);
                    return false;
                }
            } else {
                // the column was not found.
                psWarning("No column with the name '%s' exists in the table.",
                          item->name);
            }
        }
    }

    psFree(iter);

    return true;
}
