/******************************************************************
 * DPUSER - The Next Generation
 *
 * File:     libfits/fits_file.cpp
 * Purpose:  Fits class methods to read and write fits files
 * Author:   Thomas Ott
 *
 * Revision history
 * ================
 * 05.08.1999: file created
 ******************************************************************/
 
#include "platform.h"
#include "fits.h"

#define SIMPLE_T "SIMPLE  =                    T                                                  "
#define FITS_EMPTY "                                                                                "
#define FITS_END "END                                                                             "
#define FITS_NAXIS0 "NAXIS   = 0                                                                     "
#define FITS_BITPIX8 "BITPIX  = 8                                                                     "
#define EXTENT_T "XTENSION= 'IMAGE   '                                                            "

/*
 * Utility function to correct byte order
 * if we are on a machine with little-endian byte order, we have to swap
 */

bool getByteOrder( void )
{
// Are we little or big endian?  From Harbison&Steele.
  union
  {
    long l;
    char c[sizeof (long)];
  } u;
  u.l = 1;
  return(u.c[sizeof (long) - 1] == 1);
}

void swapBytes(void *dataptr, int b, dpint64 n_elements)
{
    dpint64 i;
    char tempr;
    char *_tmp;

	if (!getByteOrder()) {
        _tmp = (char *)dataptr;
        if (b == C16) {
            for (i = 0; i < n_elements * 8 * 2; i += 8) {
                SWAP(_tmp[i], _tmp[i+7]);
                SWAP(_tmp[i+1], _tmp[i+6]);
                SWAP(_tmp[i+2], _tmp[i+5]);
                SWAP(_tmp[i+3], _tmp[i+4]);
			}
        } else if ((b == R8) || (b == I8)) {
            for (i = 0; i < n_elements * 8; i += 8) {
                SWAP(_tmp[i], _tmp[i+7]);
                SWAP(_tmp[i+1], _tmp[i+6]);
                SWAP(_tmp[i+2], _tmp[i+5]);
                SWAP(_tmp[i+3], _tmp[i+4]);
			}		
        } else if ((b == R4) || (b == I4)) {
            for (i = 0; i < n_elements * 4; i += 4) {
                SWAP(_tmp[i], _tmp[i+3]);
                SWAP(_tmp[i+1], _tmp[i+2]);
			}
        } else if (b == I2) {
            for (i = 0; i < n_elements * 2; i += 2) {
                SWAP(_tmp[i], _tmp[i+1]);
			}
		}
	}
}

/*!
our own fseek function
*/

int Fits::dpseek(long long offset, int whence) {
    if (gz) {
        return gzseek(gfd, offset, whence);
    } else {
#ifdef WIN
        return fseeko64(fd, offset, whence);
#else
        return fseek(fd, offset, whence);
#endif
    }
}

/*!
our own ftell function
*/

long long Fits::dptell() {
    if (gz) {
        return gztell(gfd);
    } else {
#ifdef WIN
        return ftello64(fd);
#else
        return ftell(fd);
#endif
    }
}

/*!
our own fread function
*/

dpint64 Fits::dpread(void *ptr, dpint64 nmemb) {
    if (gz) {
        return gzread(gfd, ptr, nmemb);
    } else {
        return fread(ptr, 1, nmemb, fd);
    }
}

/*!
Open a FITS file for READ ONLY
*/

bool Fits::OpenFITS(const char *fname) {
    char *ff = strdup(fname);
    unsigned char gz_magic_number[2];

    gz = FALSE;

    if ((fd = fopen(ff, "rb")) == NULL) {
        ff = (char*)realloc(ff, (strlen(ff)+4) * sizeof(char));
        strcat(ff, ".gz");
        if ((fd = fopen(ff, "rb")) == NULL) {
            free(ff);
            return fits_error("Could not open file for reading");
        }
    }
    if (fread(gz_magic_number, 1, 2, fd) != 2) {
        free(ff);
        fclose(fd);
        return fits_error("Could not open file for reading");
    }

    if ((gz_magic_number[0] == 0x1f) && (gz_magic_number[1] == 0x8b)) {
        fclose(fd);
        if ((gfd = gzopen(ff, "rb")) == NULL) {
            free(ff);
            return fits_error("Could not open file for reading");
        }
        gz = TRUE;
    } else {
        fseek(fd, 0L, SEEK_SET);
    }

    ro = TRUE;

	free(ff);
	return TRUE;
}

/*!
Open a FITS file for WRITING
*/

bool Fits::CreateFITS(const char *fname) {
    if ((fd = fopen(fname, "wb")) == NULL) {
		return fits_error("Could not open file for writing");
	}
	ro = FALSE;
	return TRUE;
}

/*!
Close a previously opened FITS file. If the file was opened for read/write, append enough 0's so
its total length is an integer multiple of 2880 bytes (FITS standard).
*/

void Fits::CloseFITS() {
    if (ro) {
        if (gz) {
            gzclose(gfd);
        } else {
            fclose(fd);
        }
    } else {
		long pos;
//		char w = (char)0;
		char w[2880];
		
		fseek(fd, 0L, SEEK_END);
		pos = ftell(fd);
//		fwrite(&w, 1, 2880 - (pos % 2880), fd);
        if ((pos % 2880) != 0) {
            fwrite(w, 2880 - (pos % 2880), 1, fd);
        }
		fclose(fd);
		fd = NULL;
	}
}

/*!
Read in data in the FITS file
Partial loading is supported
*/

bool Fits::ReadFITSData(dpint64 x1, dpint64 x2, dpint64 y1, dpint64 y2, dpint64 z1, dpint64 z2) {
    dpint64 b;
    dpint64 n, count;
    dpint64 partial;

	b = abs(filebits) / 8;

// check if values are correct
	if ((x1 > x2) || (y1 > y2) || (z1 > z2)) {
		dp_output("Reading partial data: Wrong arguments.\n");
		return FALSE;
	}
    if ((x2 > (dpint64)naxis[1]) || (y2 > (dpint64)naxis[2]) || (z2 > (dpint64)naxis[3])) {
		dp_output("Reading partial data: Arguments out of bounds.\n");
		return FALSE;
	}

// check which case
	if ((x1 == 0) && (x2 == 0) && (y1 == 0) && (y2 == 0) && (z1 == 0) && (z2 == 0)) {

// read all data
		partial = b * n_elements / 100;
		if (partial < 10000) {
            if (dpread(dataptr, b * n_elements) != n_elements * b) {
				return fits_error("Premature end of data");
			}
		} else {
// read in in smaller chunks so we can display a progress bar
			char *_tmp = (char *)dataptr;
			dpProgress(-100, "");
			for (count = 0; count < 99; count++) {
                if (dpread(_tmp, partial) != partial) {
					dpProgress(0, "");
					return fits_error("Premature end of data");
				}
				_tmp += partial;
				dpProgress(count+1, "reading file...");
			}
			partial = b * n_elements - 99 * partial;
            if (dpread(_tmp, partial) != partial) {
				dpProgress(0, "");
				return fits_error("Premature end of data");
			}
			dpProgress(100, "");
		}
	} else if ((x1 == 0) && (x2 == 0) && (y1 == 0) && (y2 == 0)) {

// read complete images
        dpseek((naxis[1]) * naxis[2] * (z1 - 1) * b, SEEK_CUR);
        if (dpread(dataptr, b * n_elements) != n_elements * b) {
			return fits_error("Premature end of data");
		}
	} else if ((x1 == 0) && (x2 == 0)) {
		char *_tmp = (char *)dataptr;
		long pos = 0;
		
		if (z1 == 0) z1 = 1;
		if (z2 == 0) z2 = naxis[3];
        dpseek((naxis[1]) * naxis[2] * (z1 - 1) * b, SEEK_CUR);
        for (n = (int)z1; n <= (int)z2; n++) {
            dpseek((naxis[1]) * (y1 - 1) * b, SEEK_CUR);
            if (dpread(&_tmp[pos], b * naxis[1] * (y2 - y1 + 1)) != ((dpint64)(naxis[1]) * (y2 - y1 + 1) * b)) {
				return fits_error("Premature end of data");
			}	
			pos += naxis[1] * (y2 - y1 + 1) * b;
            dpseek((naxis[1]) * (naxis[2] - y2) * b, SEEK_CUR);
		}	
	} else {
		char *_tmp = (char *)dataptr;
		char *__tmp;
		long x, y, pos = 0;
		
		if (z1 == 0) z1 = 1;
		if (z2 == 0) z2 = naxis[3];
		if (y1 == 0) y1 = 1;
		if (y2 == 0) y2 = naxis[2];
		if ((__tmp = (char *)malloc(naxis[1] * (y2 - y1 + 1) * b)) == NULL) {
			return fits_error("Read partial: Could not allocate memory");
		}
        dpseek((naxis[1]) * naxis[2] * (z1 - 1) * b, SEEK_CUR);
		for (n = (int)z1; n <= (int)z2; n++) {
            dpseek((naxis[1]) * (y1 - 1) * b, SEEK_CUR);
            if (dpread(__tmp, b * naxis[1] * (y2 - y1 + 1)) != ((dpint64)(naxis[1]) * (y2 - y1 + 1) * b)) {
				return fits_error("Premature end of data");
			}
			for (y = 0; y <= (int)(y2 - y1); y++) {
				for (x = (int)(x1 - 1) * b; x < (int)x2 * b; x++) {
					_tmp[pos] = __tmp[y * naxis[1] * b + x];
					pos++;
				}
			}
            dpseek((naxis[1]) * (naxis[2] - y2) * b, SEEK_CUR);
		}
		free(__tmp);
	}

// swap bytes
	swapBytes(dataptr, filebits, n_elements);

// Deal with SHARP data
	if (isSHARPIdata()) convertSHARPIdata();

// everything ok
	return TRUE;
}

/*!
Read in one image of a FITS cube
*/

bool Fits::ReadFITSCubeImage(int which) {
    dpint64 n, b;

	b = abs(filebits) / 8;
	if (!getHeaderInformation()) return FALSE;
	if (!recreate(naxis[1], naxis[2])) return FALSE;
// seek file position
    if (dpseek(hdr + b * (which - 1) * n_elements, SEEK_SET) == -1) {
		return fits_error("Cannot acces cube image");
	}

// read data
    if ((n = dpread(dataptr, b * n_elements)) != n_elements * b) {
		dp_output("Could only read %li elements.", n);
		return fits_error("Premature end of data");
	}

// swap bytes
	swapBytes(dataptr, filebits, n_elements);
	membits = filebits;

// deal with special SHARP I case
	if (isSHARPIdata()) convertSHARPIdata();

// everything ok
	return TRUE;
}

/*!
Write the data to file.
*/

bool Fits::WriteFITSData(void) {
	if (fd == NULL) return FALSE;
	swapBytes(dataptr, membits, n_elements);
    if (fwrite(dataptr, abs((int)membits / 8), n_elements, fd) != n_elements) {
		swapBytes(dataptr, membits, n_elements);
		CloseFITS();
		return FALSE;
	}
	swapBytes(dataptr, membits, n_elements);
	return TRUE;
}

/*!
Read in a complete FITS file
by default, the image(s) is stored in a float array
*/

bool Fits::ReadFITS(const char *fname, dpint64 x1, dpint64 x2, dpint64 y1, dpint64 y2, dpint64 z1, dpint64 z2) {
    dpint64 mem;
	bool hasrefpix;

// open file
	if (!OpenFITS(fname)) return FALSE;
//	if (!OpenFITS(fname)) throw dpFitsException(3);

// read header
    if ((hdr = ReadFitsHeader()) < 0) return FALSE;

// get Header info and allocate memory
	if (!getHeaderInformation()) return FALSE;
	if ((x1) || (x2)) mem = x2 - x1 + 1;
	else mem = naxis[1];
	if ((y1) || (y2)) mem *= y2 - y1 + 1;
	else mem *= naxis[2];
	if ((z1) || (z2)) mem *= z2 - z1 + 1;
	else mem *= naxis[3];
	mem *= abs(membits / 8);
	if (mem <= 0) return fits_error("ReadFits: invalid range");
	if (!allocateMemory(mem)) return FALSE;

// read data
	if (!ReadFITSData(x1, x2, y1, y2, z1, z2)) return FALSE;

	hasrefpix = hasRefPix();
	if ((x1) || (x2)) {
		naxis[1] = x2 - x1 + 1;
		if ((x1 > 1) && hasrefpix) SetFloatKey("CRPIX1", getCRPIX(1) - (double)(x1 - 1));
	}
	if ((y1) || (y2)) {
		naxis[2] = y2 - y1 + 1;
		if ((y1 > 1) && hasrefpix) SetFloatKey("CRPIX2", getCRPIX(2) - (double)(y1 - 1));
	}
	if ((z1) || (z2)) {
		naxis[3] = z2 - z1 + 1;
		if ((z1 > 1) && hasrefpix) SetFloatKey("CRPIX3", getCRPIX(3) - (double)(z1 - 1));
	}
	deflate();
//	if (naxis[3] < 2) naxis[0] = 2;

// close file
	CloseFITS();

// everything ok
	strncpy(FileName, fname, 255);
    extensionNumber = 0;
	
	return TRUE;
}	


/*!
Read in a complete FITS file
by default, the image(s) is stored in a float array
*/

bool Fits::ReadFITSExtension(const char *fname, int ext, dpint64 x1, dpint64 x2, dpint64 y1, dpint64 y2, dpint64 z1, dpint64 z2)
{
    dpint64 mem;
    bool hasrefpix;

    // open file
    if (!OpenFITS(fname))
        return FALSE;
    //	if (!OpenFITS(fname))
    //      throw dpFitsException(3);

    // read header
    if ((hdr = ReadFitsExtensionHeader(ext)) < 0)
        return FALSE;

    // get Header info and allocate memory
    if (!getHeaderInformation())
        return FALSE;
    if ((x1) || (x2))
        mem = x2 - x1 + 1;
    else
        mem = naxis[1];

    if ((y1) || (y2))
        mem *= y2 - y1 + 1;
    else
        mem *= naxis[2];

    if ((z1) || (z2))
        mem *= z2 - z1 + 1;
    else
        mem *= naxis[3];

    mem *= abs(membits / 8);

    if (mem <= 0)
        return fits_error("ReadFits: invalid range");

    if (!allocateMemory(mem))
        return FALSE;

    // read data only if NAXIS > 0 (no empty data section)
    int val = 0;
    GetIntKey("NAXIS", &val);
    if(val > 0) {
        if (!ReadFITSData(x1, x2, y1, y2, z1, z2)) return FALSE;
    }

    hasrefpix = hasRefPix();
    if ((x1) || (x2)) {
        naxis[1] = x2 - x1 + 1;
        if ((x1 > 1) && hasrefpix)
            SetFloatKey("CRPIX1", getCRPIX(1) - (double)(x1 - 1));
    }
    if ((y1) || (y2)) {
        naxis[2] = y2 - y1 + 1;
        if ((y1 > 1) && hasrefpix)
            SetFloatKey("CRPIX2", getCRPIX(2) - (double)(y1 - 1));
    }
    if ((z1) || (z2)) {
        naxis[3] = z2 - z1 + 1;
        if ((z1 > 1) && hasrefpix)
            SetFloatKey("CRPIX3", getCRPIX(3) - (double)(z1 - 1));
    }
    deflate();
//    if (naxis[3] < 2) naxis[0] = 2;

    // close file
    CloseFITS();

    // everything ok
    strncpy(FileName, fname, 255);
    char tmp[20];
    sprintf(tmp, ", ext #%d", ext);
    strncat(FileName, tmp, 20);
    return TRUE;
    extensionNumber = ext;
}

bool Fits::tfieldWidth(int column, int *bytes, int *repeat, char *type) {
	char key[9];
	char tform[81];
	int c;

	sprintf(key, "TFORM%i", column);
	if (!GetStringKey(key, tform)) return FALSE;
	*repeat = 1;
	c = 0;
	if (isdigit(tform[c])) {
		*repeat = tform[c] - '0';
		c++;
        while (isdigit(tform[c])) {
			*repeat = 10 * (*repeat) + (tform[c] - '0');
			c++;
		}
	}
	switch (tform[c]) {
		case 'I':
			*bytes = 2;
			break;
		case 'J':
		case 'E':
			*bytes = 4;
			break;
		case 'D':
		case 'C':
		case 'P':
        case 'K':
			*bytes = 8;
			break;
		case 'M':
			*bytes = 16;
			break;
		default:
			*bytes = 1;
			break;
	}
	*type = tform[c];
	
	return TRUE;
}
	

/*
read the nth FITS Bintable column and return Fits
*/
bool Fits::GetBintableColumn(int column, Fits &result) {
    int tfields;
    if ((column < 1) ||
        (!GetIntKey("TFIELDS", &tfields)) ||
        (tfields < column))
    {
        return false;
    }
	
    // go through the TFORMi keywords and count the bytes
    int bytes  = 0,
        width  = 0,
        repeat = 0;
    char fieldtype;
    for (int i = 1; i < column; i++) {
        if (!tfieldWidth(i, &width, &repeat, &fieldtype)) {
            return false;
        }
		bytes += repeat * width;
	}
    if (!tfieldWidth(column, &width, &repeat, &fieldtype)) {
        return false;
    }
	width *= repeat;
	
    FitsBitpix bits = NA;
	switch (fieldtype) {
		case 'A':
		case 'B':
		case 'L': bits = I1; break;
		case 'I': bits = I2; break;
		case 'J': bits = I4; break;
        case 'K': bits = I8; break;
		case 'E': bits = R4; break;
		case 'D': bits = R8; break;
        case 'C':
        case 'M': bits = C16; break;
		default: bits = NA; break;
	}
    if (bits == NA) {
        return false;
    }

    // take care that vectors are setup correctly here
    int xDim = 0, yDim = 0;
    if (repeat == 1) {
        xDim = Naxis(2);
        yDim = repeat;
    } else {
        xDim = repeat;
        yDim = Naxis(2);
    }

    result.create(xDim, yDim, bits);
	
    int c = 0;
    for (int y = 0; y < Naxis(2); y++) {
        for (int x = bytes; x < bytes + width; x++) {
			result.i1data[c] = i1data[C_I(x, y)];
			c++;
		}
	}

    if (fieldtype != 'C') {
        swapBytes(result.dataptr, result.membits, result.Nelements());
    } else {
        swapBytes(result.dataptr, R4, result.Nelements() * 2);
        int i, j;
        for (i = result.Nelements() * 2 - 1, j = result.Nelements() - 1; i >= 0; i -= 2, j--) {
            result.cdata[j].i = (double)result.r4data[i];
            result.cdata[j].r = (double)result.r4data[i-1];
        }
    }
	
    double value;
    char key[9];
	sprintf(key, "TSCAL%i", column);
    if (GetFloatKey(key, &value)) {
        result.bscale = value;
    }
	sprintf(key, "TZERO%i", column);
    if (GetFloatKey(key, &value)) {
        result.bzero = value;
    }

// check for TDIMi
    char tdim[81];
    sprintf(key, "TDIM%i", column);
    if (GetStringKey(key, tdim)) {
        if (tdim[0] == '(') {
            int xdim = 1, ydim = 1, zdim = 1;
            int count = 0;
            count = sscanf(tdim, "(%i,%i", &xdim, &ydim);
            if (count == 2) {
                if (result.Naxis(2) > 1) zdim = result.Naxis(2);
                if (result.Nelements() == xdim * ydim * zdim) {
                    result.setNaxis(0, (zdim > 1 ? 3 : 2));
                    result.setNaxis(1, xdim);
                    result.setNaxis(2, ydim);
                    result.setNaxis(3, zdim);

                    result.extensionType = BINTABLEIMAGE;
                } else {
                    result.extensionType = BINTABLE;
                }
            }
        }
    } else {
        result.extensionType = BINTABLE;
    }

	result.updateHeader();

    // get column name
    sprintf(key, "TTYPE%i", column);
    GetStringKey(key, result.columnName);
	
    return true;
}

/*
read the nth FITS Bintable column and return dpStringList
*/
bool Fits::GetBintableColumn(int column, dpStringList &result) {
    bool retValue = false;
    int tfields;
    if ((column < 1) ||
        (!GetIntKey("TFIELDS", &tfields)) ||
        (tfields < column))
    {
        return false;
    }

    // go through the TFORMi keywords and count the bytes
    int bytes  = 0,
        width  = 0,
        repeat = 0;
    char fieldtype;
    for (int i = 1; i < column; i++) {
        if (!tfieldWidth(i, &width, &repeat, &fieldtype)) {
            return false;
        }
        bytes += repeat * width;
    }
    if (!tfieldWidth(column, &width, &repeat, &fieldtype)) {
        return false;
    }
    width *= repeat;

    if (fieldtype == 'A') {
        for (int y = 0; y < Naxis(2); y++) {
            dpString t;
            for (int x = bytes; x < bytes + width; x++) {
                t.append(1, (char)i1data[C_I(x, y)]);
            }
            result.push_back(t);
        }
        retValue = true;
    }

    // get column name
    char key[9],
         value[72];
    sprintf(key, "TTYPE%i", column);

    GetStringKey(key, value);
    result.setColumnName(value);

    return retValue;
}

nodeEnum Fits::GetBintableType(int column) {
    int tfields;
    if ((column < 1) ||
        (!GetIntKey("TFIELDS", &tfields)) ||
        (tfields < column))
    {
        return typeUnknown;
    }

    // go through the TFORMi keywords and count the bytes
    int bytes  = 0,
        width  = 0,
        repeat = 0;
    char fieldtype;
    for (int i = 1; i < column; i++) {
        if (!tfieldWidth(i, &width, &repeat, &fieldtype)) {
            return typeUnknown;
        }
        bytes += repeat * width;
    }
    if (!tfieldWidth(column, &width, &repeat, &fieldtype)) {
        return typeUnknown;
    }
    width *= repeat;

    nodeEnum type = typeUnknown;
    switch (fieldtype) {
        case 'A':
            type = typeStrarr;
            break;
        case 'B':
        case 'L':
        case 'I':
        case 'J':
        case 'K':
        case 'E':
        case 'D':
        case 'C':
        case 'M':
            type = typeFits;
            break;
        default:
            type = typeUnknown;
            break;
    }

    return type;
}

/*
When using this function, create a new Fits and delete it afterwards!
Otherwise header information for the first (primary) header could be corrupted.
*/
int  Fits::CountExtensions(const char *fname) {
    int ext = 1;

    // open file
    if (!OpenFITS(fname))
        return -1;

    // read header
    while (ReadFitsExtensionHeader(ext, TRUE) > 0) {
        ext++;
    }


    // close file
    CloseFITS();

    return ext-1;
}

/*
When using this function, create a new Fits and delete it afterwards!
Otherwise header information for the first (primary) header could be corrupted.
*/
int Fits::FindExtensionByName(const char *fname, const char *extname) {
	int ext = 1;
	char extensionName[81];

	if (!OpenFITS(fname))
		return -1;

// read header
    while (ReadFitsExtensionHeader(ext, TRUE) > 0) {
		if (GetStringKey("EXTNAME", extensionName)) {
			dpString s(extensionName);
			s.strtrim(0);
			if (strcasecmp(extname, s.c_str()) == 0) {
				CloseFITS();
				return ext;
			}
		}
		ext++;
	}
	CloseFITS();

	return -1;
}

/*
When using this function, create a new Fits and delete it afterwards!
Otherwise header information for the first (primary) header could be corrupted.
*/
dpStringList Fits::ListFitsExtensions(const char *fname) {
	int ext = 1;
	char extensionName[81];
	dpStringList result;

	if (!OpenFITS(fname))
		return result;

// read header
    while (ReadFitsExtensionHeader(ext, TRUE) > 0) {
		if (GetStringKey("EXTNAME", extensionName)) {
			dpString s(extensionName);
			s.strtrim(0);
			result.append(s);
        } else {
            result.append(dpString::number(ext));
        }
		ext++;
	}
	CloseFITS();

	return result;
}

int Fits::FindColumnByName(const char *colname) {
	char columnName[81], columnKey[9];
	int rv = -1;
	int i = 1;
	bool found = FALSE;

	i = 1;
	while (!found) {
		sprintf(columnKey, "TTYPE%i", i);
		if (!GetStringKey(columnKey, columnName)) {
				return -1;
		}
		dpString s(columnName);
		s.strtrim(0);
		found = (strcasecmp(colname, s.c_str()) == 0);
		if (found) {
			return i;
		}
		i++;
		if (i > 999) {
			return -1;
		}
	}

	return -1;
}

dpStringList Fits::ListTableColumns(const char *fname, const int extension) {
	int column = 1;
	char columnKey[9];
	char columnName[81];
	dpStringList result;
	bool found = TRUE;

	if (!OpenFITS(fname))
		return result;

// read header
    if (ReadFitsExtensionHeader(extension, TRUE) > 0) {
		while (found) {
			sprintf(columnKey, "TTYPE%i", column);
			if (GetStringKey(columnKey, columnName)) {
				dpString s(columnName);
				s.strtrim(0);
				result.append(s);
			} else {
				found = FALSE;
			}
			column++;
			if (column > 999) {
				found = FALSE;
			}
		}
	}
	CloseFITS();

	return result;
}

int Fits::CountTableColumns(const char *fname, const int extension) {
    int  column = 1;
    char columnKey[9],
         columnName[81];
    bool found = TRUE;

    if (!OpenFITS(fname))
        return -1;

    // read header
    if (ReadFitsExtensionHeader(extension, TRUE) > 0) {
        while (found) {
            sprintf(columnKey, "TTYPE%i", column);
            if (!GetStringKey(columnKey, columnName)) {
                found = FALSE;
            }
            column++;
            if (column > 999) {
                found = FALSE;
            }
        }
    }
    CloseFITS();

    return column-1;
}

/*!
Write a fits file
*/

bool Fits::WriteFITS(const char *fname) {
// open file
	if (!CreateFITS(fname)) return FALSE;

// write header
	if (!updateHeader()) return FALSE;
    if (!WriteFitsHeader(fd)) return FALSE;

// write data
    if (!WriteFITSData()) return FALSE;

// close file
	CloseFITS();

// everything ok
	strncpy(FileName, fname, 255);
	return TRUE;
}

/*!
Write a FITS extension
*/
bool Fits::WriteFITSExtension(const char *fname) {
    char s[7];
    dpint64 pos;
    char w[2880];
    int i;

// Check of file exists already, and is a uncompressed FITS file
    if ((fd = fopen(fname, "rb")) != NULL) {
        // Yes. Is this a FITS file?
        fread(s, 1, 6, fd);
        s[6] = 0;
        fclose(fd);
        if (strncmp(s, "SIMPLE", 6) != 0) {
            return FALSE;
        }
    }
    if ((fd = fopen(fname, "ab")) == NULL) {
        return FALSE;
    }
    ro = FALSE;
    fseek(fd, 0L, SEEK_END);

// we are at the end of the file already. See if we have to pad 0's...
    pos = ftell(fd);
    if (pos == 0) {
        fwrite(SIMPLE_T, 80, 1, fd);
        fwrite(FITS_NAXIS0, 80, 1, fd);
        fwrite(FITS_BITPIX8, 80, 1, fd);
        fwrite(FITS_END, 80, 1, fd);
        for (i = 0; i < 32; i++) fwrite(FITS_EMPTY, 80, 1, fd);
// write an empty primary FITS header
    } else if ((pos % 2880) != 0) {
        fwrite(w, 2880 - (pos % 2880), 1, fd);
    }
// Set PCOUNT and GCOUNT header keys (required by the FITS standard)
    SetIntKey("PCOUNT", 0, "size of special data area");
    SetIntKey("GCOUNT", 1, "one data group (required keyword)");

// Write the new header
    fwrite(EXTENT_T, 1, 80, fd);
    fwrite(header+80, 1, HeaderLength-80, fd);

//...and the data
    if (!WriteFITSData()) {
        return FALSE;
    }

    CloseFITS();

    return TRUE;
}

/*!
Update Header values
*/

bool Fits::updateHeader() {
	SetIntKey("BITPIX", (int)membits);
	SetIntKey("NAXIS", naxis[0]);
	if (naxis[0] > 0) {
		SetIntKey("NAXIS1", naxis[1]);
	} else DeleteKey("NAXIS1");
	if (naxis[0] > 1) {
		SetIntKey("NAXIS2", naxis[2]);
	} else DeleteKey("NAXIS2");
	if (naxis[0] > 2)
		SetIntKey("NAXIS3", naxis[3]);
	else DeleteKey("NAXIS3");
	if (membits > 0) {
		SetFloatKey("BSCALE", bscale);
		SetFloatKey("BZERO", bzero);
	} else {
		DeleteKey("BSCALE");
		DeleteKey("BZERO");
	}
/*
	if ((naxis[0] > 1) && (hasRefPix())) {
		char value[10];
		DeleteKey("CD1_1");
		DeleteKey("CD1_2");
		DeleteKey("CD2_1");
		DeleteKey("CD2_2");
		SetFloatKey("CROTA1", crot);
		SetFloatKey("CROTA2", crot);
		if (strlen(crtype) > 0) {
			sprintf(value, "RA--%s", crtype);
			SetStringKey("CTYPE1", value);
			sprintf(value, "DEC-%s", crtype);
			SetStringKey("CTYPE2", value);
		}
	} else {
		char key[10];
		for (int i = naxis[0]+1; i <= MAXNAXIS; i++) {
			sprintf(key, "CRPIX%i", i);
			DeleteKey(key);
			sprintf(key, "CRVAL%i", i);
			DeleteKey(key);
			sprintf(key, "CDELT%i", i);
			DeleteKey(key);
			sprintf(key, "CTYPE%i", i);
			DeleteKey(key);
			sprintf(key, "CUNIT%i", i);
			DeleteKey(key);
		}
		DeleteKey("CD1_1");
		DeleteKey("CD1_2");
		DeleteKey("CD2_1");
		DeleteKey("CD2_2");
	}
*/
	return TRUE;
}

/*!
Update local variables to that of the FITS header
*/

bool Fits::getHeaderInformation() {
	int iv, i, intvalue;
	double v;
	
	if (!GetIntKey("BITPIX", &iv)) return FALSE;
	filebits = membits = (FitsBitpix)iv;
	if (!GetIntKey("NAXIS", &intvalue)) return FALSE;
    naxis[0] = (dpint64)intvalue;
	for (i = 0; i < MAXNAXIS; i++) {
		naxis[i + 1] = 1;
	}
	intvalue = 0;
	GetIntKey("NAXIS1", &intvalue);
    if (intvalue > 1) naxis[1] = (dpint64)intvalue;
	intvalue = 0;
	GetIntKey("NAXIS2", &intvalue);
    if (intvalue > 1) naxis[2] = (dpint64)intvalue;
	intvalue = 0;
	GetIntKey("NAXIS3", &intvalue);
    if (intvalue > 1) naxis[3] = (dpint64)intvalue;
	if (!GetFloatKey("BSCALE", &v)) bscale = 1.0;
	else bscale = v;
	if (!GetFloatKey("BZERO", &v)) bzero = 0.0;
	else bzero = v;

// Paranoia
	if (naxis[3] < 2) {
		naxis[3] = 1;
		naxis[0] = 2;
	}
// even more paranoia
	if (naxis[0] > MAXNAXIS) naxis[0] = MAXNAXIS;

	if (naxis[0] == 2) naxis[3] = 1;
	if (bscale == 0.0) {
		bscale = 1.0;
		bzero = 0.0;
	}

// more precise, at least for the first two axis
	wcsinfo info;
	info = readWCSinfo(*this);
//	crot = info.rot;
	
//	if (info.xinc != getCDELT(1)) SetFloatKey("CDELT1", info.xinc);
//	if (info.yinc != getCDELT(2)) SetFloatKey("CDELT2", info.yinc);
	
	strncpy(crtype, info.type, 9);

	return TRUE;
}

/*!
Find out if the header represents SHARP I data: This is the case, if:
1> The keys INSTRUME and BITPIX do exist and have "SHARPI  " and 16, respectively
2> The keys BSCALE and BZERO don't exist.
*/

bool Fits::isSHARPIdata() {
	char instrument[81], imagetyp[81];
	double bscale, bzero;
	int bitpix;

	if ((GetStringKey("INSTRUME", instrument)) && (GetIntKey("BITPIX", &bitpix)) && (GetStringKey("IMAGETYP", imagetyp))) {
		if ((!GetFloatKey("BSCALE", &bscale)) && (!GetFloatKey("BZERO", &bzero))) {
			instrument[5] = 0;
			if ((strcmp(instrument, "SHARP") == 0) && (strcmp(imagetyp, "INTEGER ") == 0) && (bitpix == 16)) return TRUE;
		}
	}
	return FALSE;
}

/*!
Convert SHARP I data to valid FITS format:
SHARP I data is unsigned short, which does not exist in FITS standard. So
we subtract 32767 from each pixel and typecast to signed short. The data values
represent the correct values if we also set bzero to 32767 and bscale to 1.0.
*/

void Fits::convertSHARPIdata() {
	unsigned short *usdata = (unsigned short *)i2data;
    dpint64 i;

	for (i = 0; i < Nelements(); i++) i2data[i] = (signed short)((long)usdata[i] - 32767);
	bzero = 32767.0;
	bscale = 1.0;
}

struct bmpfileheader {
    char bfType[2];
    unsigned long   bfSize; 
    unsigned short    bfReserved1; 
    unsigned short    bfReserved2; 
    unsigned long   bfOffBits;
};

typedef struct {
	unsigned long   biSize; 
    long    biWidth; 
    long    biHeight; 
    unsigned short    biPlanes; 
    unsigned short    biBitCount; 
    unsigned long  biCompression; 
    unsigned long  biSizeImage; 
    long   biXPelsPerMeter; 
    long   biYPelsPerMeter; 
    unsigned long  biClrUsed; 
    unsigned long  biClrImportant;
} bmpinfo;

/*!
Utility function to swap the bytes in a BMP header.
*/
void swapBMPBytes(bmpfileheader *h, bmpinfo *b) {
    char *_tmp;
    char tempr;

    _tmp = (char *)h;
    SWAP(_tmp[4], _tmp[7]);
    SWAP(_tmp[5], _tmp[6]);
    SWAP(_tmp[8], _tmp[9]);
    SWAP(_tmp[10], _tmp[11]);
    SWAP(_tmp[12], _tmp[15]);
    SWAP(_tmp[13], _tmp[14]);

	_tmp = (char *)b;
    SWAP(_tmp[0], _tmp[3]);
    SWAP(_tmp[1], _tmp[2]);
    SWAP(_tmp[4], _tmp[7]);
    SWAP(_tmp[5], _tmp[6]);
    SWAP(_tmp[8], _tmp[11]);
    SWAP(_tmp[9], _tmp[10]);
    SWAP(_tmp[12], _tmp[13]);
    SWAP(_tmp[14], _tmp[15]);
    SWAP(_tmp[16], _tmp[19]);
    SWAP(_tmp[17], _tmp[18]);
    SWAP(_tmp[20], _tmp[23]);
    SWAP(_tmp[21], _tmp[22]);
    SWAP(_tmp[24], _tmp[27]);
    SWAP(_tmp[25], _tmp[26]);
    SWAP(_tmp[28], _tmp[31]);
    SWAP(_tmp[29], _tmp[30]);
    SWAP(_tmp[32], _tmp[35]);
    SWAP(_tmp[33], _tmp[34]);
    SWAP(_tmp[36], _tmp[39]);
    SWAP(_tmp[37], _tmp[38]);
}

/*!
Write a BMP file.
*/

bool Fits::WriteBMP(const char *fname) {
	FILE *fd;
	bmpfileheader h;
	bmpinfo b;
	long x, y;
	unsigned char color[3];
	dpint64 off;

	off = Naxis(1) % 4;

	h.bfType[0] = 'B';
	h.bfType[1] = 'M';
	h.bfSize = 14 + 40 + Nelements() * 3 + off * Naxis(2);
	h.bfReserved1 = 0;
	h.bfReserved2 = 0;
	h.bfOffBits = 14 + 40;

	b.biSize = 40;
	b.biWidth = Naxis(1);
	b.biHeight = Naxis(2);
	b.biPlanes = 1;
	b.biBitCount = 24;
	b.biCompression = 0;
	b.biSizeImage = Nelements() * 3 + off * Naxis(2);
	b.biXPelsPerMeter = 2834;
	b.biYPelsPerMeter = 2834;
	b.biClrUsed = 0;
	b.biClrImportant = 0;

/* Swap bytes if necessary */
	if (getByteOrder()) swapBMPBytes(&h, &b);

	if ((fd = fopen(fname, "w+b")) == NULL) {
		dp_output("Fits::WriteBMP: Could not open file %s for writing\n", fname);
		return FALSE;
	}

	fwrite(h.bfType, 2, 1, fd);
	fwrite(&h.bfSize, 4, 1, fd);
	fwrite(&h.bfReserved1, 2, 1, fd);
	fwrite(&h.bfReserved2, 2, 1, fd);
	fwrite(&h.bfOffBits, 4, 1, fd);
	fwrite(&b, 40, 1, fd);

	for (y = 0; y < Naxis(2); y++) {
		for (x = 0; x < Naxis(1); x++) {
			color[0] = color[1] = color[2] = (unsigned char)ValueAt(C_I(x, y));
			fwrite(color, sizeof(color), 1, fd);
		}
		if (off) fwrite(color, off, 1, fd);
	}
	fclose(fd);

	return TRUE;
}
/*!
Write a color BMP file.
*/

bool Fits::WriteBMP(const char *fname, Fits &red, Fits &green, Fits &blue) {
	FILE *fd;
	bmpfileheader h;
	bmpinfo b;
	long x, y, n;
	unsigned char color[3];
	dpint64 off;

/* Consistency check */
	if ((!red.matches(green)) || (!red.matches(blue))) {
		dp_output("Fits::WriteBMP: The arrays don't match in size\n");
		return FALSE;
	}
	off = red.Naxis(1) % 4;

	h.bfType[0] = 'B';
	h.bfType[1] = 'M';
	h.bfSize = 14 + 40 + red.Nelements() * 3 + off * red.Naxis(2);
	h.bfReserved1 = 0;
	h.bfReserved2 = 0;
	h.bfOffBits = 14 + 40;

	b.biSize = 40;
	b.biWidth = red.Naxis(1);
	b.biHeight = red.Naxis(2);
	b.biPlanes = 1;
	b.biBitCount = 24;
	b.biCompression = 0;
	b.biSizeImage = red.Nelements() * 3 + off * red.Naxis(2);
	b.biXPelsPerMeter = 2834;
	b.biYPelsPerMeter = 2834;
	b.biClrUsed = 0;
	b.biClrImportant = 0;

/* Swap bytes if necessary */
	if (getByteOrder()) swapBMPBytes(&h, &b);

	if ((fd = fopen(fname, "w+b")) == NULL) {
		dp_output("Fits::WriteBMP: Could not open file %s for writing\n", fname);
		return FALSE;
	}

	fwrite(h.bfType, 2, 1, fd);
	fwrite(&h.bfSize, 4, 1, fd);
	fwrite(&h.bfReserved1, 2, 1, fd);
	fwrite(&h.bfReserved2, 2, 1, fd);
	fwrite(&h.bfOffBits, 4, 1, fd);
	fwrite(&b, 40, 1, fd);

	for (y = 0; y < red.Naxis(2); y++) {
		for (x = 0; x < red.Naxis(1); x++) {
			n = red.C_I(x, y);
			color[2] = (unsigned char)red.ValueAt(n);
			color[1] = (unsigned char)green.ValueAt(n);
			color[0] = (unsigned char)blue.ValueAt(n);
			fwrite(color, sizeof(color), 1, fd);
		}
		if (off) fwrite(color, off, 1, fd);
	}
	fclose(fd);

	return TRUE;
}