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

#include <unistd.h>

#include "psFits.h"
#include "string.h"
#include "psError.h"
#include "psFileUtilsErrors.h"
#include "psImageExtraction.h"
#include "psMemory.h"
#include "psString.h"
#include "psLogMsg.h"
#include "psTrace.h"

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

static char* standardFitsKeys[] = {
                                      "SIMPLE",
                                      "BITPIX",
                                      "NAXIS",
                                      "NAXIS1",
                                      "NAXIS2",
                                      "NAXIS3",
                                      "EXTEND",
                                      "COMMENT",
                                      "EXTNAME",
                                      "HDUNAME",
                                      "XTENSION",
                                      "TFIELDS",
                                      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;
    default:
        return PS_TYPE_PTR;
    }
}


static bool convertPsTypeToFits(psElemType 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;

    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 bool convertPsTypeToBinaryTForm(psElemType type, char** fitsType)
{
    switch (type) {
    case PS_TYPE_BOOL:
        *fitsType = psStringCopy("1L");
        break;
    case PS_TYPE_S32:
        *fitsType = psStringCopy("1J");
        break;
    case PS_TYPE_F32:
        *fitsType = psStringCopy("1E");
        break;
    case PS_TYPE_F64:
        *fitsType = psStringCopy("1D");
        break;
        // XXX: Handle other types, e.g., Vectors, etc.
    default:
        return false;
    }

    return true;
}

static bool isHDUEmpty(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->p_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->p_fd, &status);
        psFree((void*)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->p_fd = fptr;
    strcpy((char*)fits->filename,name);
    p_psMemSetDeallocator(fits,(psFreeFcn)fitsFree);

    return fits;
}

bool psFitsMoveExtName(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->p_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(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->p_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->p_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(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->p_fd,&hdunum) - 1;
}

char* psFitsGetExtName(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->p_fd, "EXTNAME", name, NULL, &status) != 0) {
        status = 0;
        if (fits_read_key_str(fits->p_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->p_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(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->p_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(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->p_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->p_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->p_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_TYPE_S32,
                                        PS_META_PRIMITIVE,
                                        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_TYPE_F64,
                                        PS_META_PRIMITIVE,
                                        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_TYPE_PTR,
                                            PS_META_STR,
                                            keyComment,
                                            keyValue+1);
                } else {
                    success = psMetadataAdd(out,
                                            PS_LIST_TAIL,
                                            keyName,
                                            PS_TYPE_PTR,
                                            PS_META_STR,
                                            keyComment,
                                            keyValue);
                }
                break;
            case 'L':
                tempBool = (keyValue[0] == 'T') ? 1 : 0;
                success = psMetadataAdd(out,
                                        PS_LIST_TAIL,
                                        keyName,
                                        PS_TYPE_BOOL,
                                        PS_META_PRIMITIVE,
                                        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;
}

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

    if (out == NULL) {
        out = psHashAlloc(10);
        // XXX: what is the appropriate number of buckets? Got 10 from psMetadataAlloc.
        if (out == NULL) {
            psError(PS_ERR_UNKNOWN, false,
                    "Failed to allocate a new psHash 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) {
            psHashAdd(out, name, header);
        } else { // XXX: is this a warning or error?
            psLogMsg(__func__, PS_LOG_WARN,
                     "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.
                         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->p_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->p_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->p_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->p_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->p_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,
                      const psMetadata* header,
                      const psImage* input,
                      int numZPlanes,
                      char* extname)
{

    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->p_fd, bitPix, naxis, naxes, &status);

    if (extname != NULL) {
        fits_update_key_str(fits->p_fd, "EXTNAME", (char*)extname, NULL, &status);
    }

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

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

    if (input->parent == NULL) { // if no parent, assume that the image data is contiguous
        fits_write_img(fits->p_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->p_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->p_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->p_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,0,0,dx+1,dy+1);
    } else {
        subset = psMemIncrRefCounter((psImage*)input);
    }

    fits_write_subset(fits->p_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* header,
                       psFits* fits)
{

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

    if (header == 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(header->list,PS_LIST_HEAD,true);
    psMetadataItem* item;
    while ( (item=psListGetAndIncrement(iter)) != NULL ) {
        if (item->type == PS_META_PRIMITIVE) { // we only write primitive data types.
            switch (item->pType) {
            case PS_TYPE_BOOL: {
                    int value = item->data.B;
                    fits_update_key(fits->p_fd,
                                    TLOGICAL,
                                    item->name,
                                    &value,
                                    item->comment,
                                    &status);
                    break;
                }
            case PS_TYPE_S32:
                fits_update_key(fits->p_fd,
                                TINT,
                                item->name,
                                &item->data.S32,
                                item->comment,
                                &status);
                break;
            case PS_TYPE_F32:
                fits_update_key(fits->p_fd,
                                TFLOAT,
                                item->name,
                                &item->data.F32,
                                item->comment,
                                &status);
                break;
            case PS_TYPE_F64:
                fits_update_key(fits->p_fd,
                                TDOUBLE,
                                item->name,
                                &item->data.F64,
                                item->comment,
                                &status);
                break;
            default:
                psError(PS_ERR_BAD_PARAMETER_TYPE, true,
                        PS_ERRORTEXT_psFits_METADATA_PTYPE_UNSUPPORTED,
                        item->pType);
                break;
            }
        } else if (item->type == PS_META_STR) {
            fits_update_key(fits->p_fd,
                            TSTRING,
                            item->name,
                            item->data.V,
                            item->comment,
                            &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;
}

psMetadata* psFitsReadTableRow(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->p_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->p_fd, &numRows, &status);
    fits_get_num_cols(fits->p_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;
    }

    // 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->p_fd, col, name,
                               NULL, NULL, NULL, NULL, NULL, NULL, NULL, &status);
        } else {
            fits_get_acolparms(fits->p_fd, col, name,
                               NULL, NULL, NULL, NULL, NULL, NULL, NULL, &status);
        }

        // get the column type
        fits_get_coltype(fits->p_fd, col, &typecode, &repeat, &width, &status);

        if (status == 0) {

            #define READ_TABLE_ROW_CASE(FITSTYPE, NATIVETYPE, TYPE) \
        case FITSTYPE: { \
                NATIVETYPE value = 0; \
                int anynul = 0; \
                fits_read_col(fits->p_fd, FITSTYPE, col,row+1, \
                              1, 1, NULL, &value, &anynul, &status); \
                psMetadataAdd(data,PS_LIST_TAIL, name, \
                              PS_TYPE_##TYPE, PS_META_PRIMITIVE, \
                              "", (ps##TYPE)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;
                    int anynul = 0;
                    fits_read_col(fits->p_fd, TSTRING, col,row+1,
                                  1, 1, NULL, &value, &anynul, &status);
                    if (anynul == 0) {
                        psMetadataAdd(data,PS_LIST_TAIL, name,
                                      PS_TYPE_PTR, PS_META_STR,
                                      "", value);
                    }
                    break;
                }
            default:
                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(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->p_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->p_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->p_fd, &numRows, &status);

    psArray* result = psArrayAlloc(numRows);

    fits_read_col_str(fits->p_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(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->p_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->p_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->p_fd, &numRows, &status);

    // get the column datatype.
    int typecode;
    long repeat;
    long width;
    if ( fits_get_eqcoltype(fits->p_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->p_fd, typecode, colnum, 1 /* firstrow */,
                  1 /* firstelem */, numRows, NULL, result->data.V,
                  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->p_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->p_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++) {
        table->data[row] = psFitsReadTableRow(fits,row);
    }

    return table;
}

bool psFitsWriteTable(psFits* fits,
                      psMetadata* header,
                      psArray* table,
                      char* extname)
{
    int status = 0;
    psMetadataItem* item;

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

    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->size);
    columns->n=0;

    // find the unique items in the array of metadata 'rows'
    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 (item->type == PS_META_PRIMITIVE) {
                    bool found = false;
                    for (int n=0; n < columns->n && ! found; n++) {
                        if (strcmp(item->name,
                                   ((psMetadataItem*)(columns->data[n]))->name) == 0) {
                            found = true;
                        }
                    }
                    if (! found) {
                        psArrayAdd(columns, columns->nalloc, psMemIncrRefCounter(item));
                    }
                }
            }
            psFree(iter);
        }
    }

    if (columns->n == 0) { // no table columns found
        // XXX: Error?
        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++) {
        char* fitsType;
        columnNames->data[n] = psMemIncrRefCounter(((psMetadataItem*)columns->data[n])->name);
        if ( ! convertPsTypeToBinaryTForm(((psMetadataItem*)columns->data[n])->pType,
                                          &fitsType)) {
            // XXX: error message
            return false;
        }
        columnTypes->data[n] = fitsType;
    }

    fits_create_tbl(fits->p_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 (item->type == PS_META_PRIMITIVE) {
            psVector* col = psVectorAlloc(table->n, item->pType);
            switch (item->pType) {
            case PS_TYPE_S32:
                for (row = 0; row < table->n; row++) {
                    col->data.S32[row] = psMetadataLookupS32(
                                             table->data[row],item->name,NULL);
                }
                fits_write_col_int(fits->p_fd,
                                   n+1, // column number
                                   1, // firstrow
                                   1, // firstelem
                                   table->n, // nelements
                                   col->data.S32,
                                   &status);
                break;
            case PS_TYPE_F32:
                for (row = 0; row < table->n; row++) {
                    col->data.F32[row] = psMetadataLookupF32(
                                             table->data[row],item->name,NULL);
                }
                fits_write_col_flt(fits->p_fd,
                                   n+1, // column number
                                   1, // firstrow
                                   1, // firstelem
                                   table->n, // nelements
                                   col->data.F32,
                                   &status);
                break;
            case PS_TYPE_F64:
                for (row = 0; row < table->n; row++) {
                    col->data.F64[row] = psMetadataLookupF64(
                                             table->data[row],item->name,NULL);
                }
                fits_write_col_dbl(fits->p_fd,
                                   n+1, // column number
                                   1, // firstrow
                                   1, // firstelem
                                   table->n, // nelements
                                   col->data.F64,
                                   &status);
                break;
            case PS_TYPE_BOOL:
                for (row = 0; row < table->n; row++) {
                    col->data.S8[row] = psMetadataLookupBool(
                                            table->data[row],item->name,NULL);
                }
                fits_write_col_log(fits->p_fd,
                                   n+1, // column number
                                   1, // firstrow
                                   1, // firstelem
                                   table->n, // nelements
                                   col->data.S8,
                                   &status);
                break;
            default:
                // XXX: error message
                break;
            }
            psFree(col);
        } else if (item->type == PS_META_STR) {
            psArray* col = psArrayAlloc(table->n);
            for (row = 0; row < table->n; row++) {
                col->data[row] = item->data.V;
            }
            fits_write_col_str(fits->p_fd,
                               n, // column number
                               1, // firstrow
                               1, // firstelem
                               table->n, // nelements
                               (char**)col->data,
                               &status);
            psFree(col);
        }
    }

    psFree(columns);

    return true;
}
