Initial commit

1 files changed, 2325 insertions, 0 deletions

diff --git a/tksao/wcssubs/fitsfile.c b/tksao/wcssubs/fitsfile.c
new file mode 100644
index 0000000..c832687
--- /dev/null
+++ b/tksao/wcssubs/fitsfile.c
@@ -0,0 +1,2325 @@
+/*** File libwcs/fitsfile.c
+ *** July 25, 2014
+ *** By Jessica Mink, jmink@cfa.harvard.edu
+ *** Harvard-Smithsonian Center for Astrophysics
+ *** Copyright (C) 1996-2014
+ *** Smithsonian Astrophysical Observatory, Cambridge, MA, USA
+
+    This library is free software; you can redistribute it and/or
+    modify it under the terms of the GNU Lesser General Public
+    License as published by the Free Software Foundation; either
+    version 2 of the License, or (at your option) any later version.
+
+    This library is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+    Lesser General Public License for more details.
+    
+    You should have received a copy of the GNU Lesser General Public
+    License along with this library; if not, write to the Free Software
+    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+    Correspondence concerning WCSTools should be addressed as follows:
+           Internet email: jmink@cfa.harvard.edu
+           Postal address: Jessica Mink
+                           Smithsonian Astrophysical Observatory
+                           60 Garden St.
+                           Cambridge, MA 02138 USA
+
+ * Module:      fitsfile.c (FITS file reading and writing)
+ * Purpose:     Read and write FITS image and table files
+ * fitsropen (inpath)
+ *		Open a FITS file for reading, returning a FILE pointer
+ * fitsrhead (filename, lhead, nbhead)
+ *		Read FITS header and return it
+ * fitsrtail (filename, lhead, nbhead)
+ *		Read appended FITS header and return it
+ * fitsrsect (filename, nbhead, header, fd, x0, y0, nx, ny)
+ *		Read section of a FITS image, having already read the header
+ * fitsrimage (filename, nbhead, header)
+ *		Read FITS image, having already ready the header
+ * fitsrfull (filename, nbhead, header)
+ *		Read a FITS image of any dimension
+ * fitsrtopen (inpath, nk, kw, nrows, nchar, nbhead)
+ *		Open a FITS table file for reading; return header information
+ * fitsrthead (header, nk, kw, nrows, nchar, nbhead)
+ *		Extract FITS table information from a FITS header
+ * fitsrtline (fd, nbhead, lbuff, tbuff, irow, nbline, line)
+ *		Read next line of FITS table file
+ * ftgetr8 (entry, kw)
+ *		Extract column from FITS table line as double
+ * ftgetr4 (entry, kw)
+ *		Extract column from FITS table line as float
+ * ftgeti4 (entry, kw)
+ *		Extract column from FITS table line as int
+ * ftgeti2 (entry, kw)
+ *		Extract column from FITS table line as short
+ * ftgetc (entry, kw, string, maxchar)
+ *		Extract column from FITS table line as a character string
+ * fitswimage (filename, header, image)
+ *		Write FITS header and image
+ * fitswext (filename, header, image)
+ *		Write FITS header and image as extension to existing FITS file
+ * fitswhdu (fd, filename, header, image)
+ *		Write FITS header and image as extension to file descriptor
+ * fitscimage (filename, header, filename0)
+ *		Write FITS header and copy FITS image
+ * fitswhead (filename, header)
+ *		Write FITS header and keep file open for further writing 
+ * fitswexhead (filename, header)
+ *		Write FITS header only to FITS extension without writing data
+ * isfits (filename)
+ *		Return 1 if file is a FITS file, else 0
+ * fitsheadsize (header)
+ *  		Return size of FITS header in bytes
+ */
+
+#include <stdlib.h>
+#ifndef VMS
+#include <unistd.h>
+#endif
+#include <stdio.h>
+#include <fcntl.h>
+#include <sys/file.h>
+#include <errno.h>
+#include <string.h>
+#include "fitsfile.h"
+
+static int verbose=0;		/* Print diagnostics */
+static char fitserrmsg[80];
+static int fitsinherit = 1;	/* Append primary header to extension header */
+void
+setfitsinherit (inh)
+int inh;
+{fitsinherit = inh; return;}
+
+static off_t ibhead = 0;	/* Number of bytes read before header starts */
+
+off_t
+getfitsskip()
+{return (ibhead);}
+
+/* FITSRHEAD -- Read a FITS header */
+
+char *
+fitsrhead (filename, lhead, nbhead)
+
+char	*filename;	/* Name of FITS image file */
+int	*lhead;		/* Allocated length of FITS header in bytes (returned) */
+int	*nbhead;	/* Number of bytes before start of data (returned) */
+			/* This includes all skipped image extensions */
+
+{
+    int fd;
+    char *header;	/* FITS image header (filled) */
+    int extend;
+    int nbytes,naxis, i;
+    int ntry,nbr,irec,nrec, nbh, ipos, npos, nbprim, lprim, lext;
+    int nax1, nax2, nax3, nax4, nbpix, ibpix, nblock, nbskip;
+    char fitsbuf[2884];
+    char *headend;	/* Pointer to last line of header */
+    char *headnext;	/* Pointer to next line of header to be added */
+    int hdu;		/* header/data unit counter */
+    int extnum;		/* desired header data number
+			   (0=primary -1=first with data -2=use EXTNAME) */
+    char extname[32];	/* FITS extension name */
+    char extnam[32];	/* Desired FITS extension name */
+    char *ext;		/* FITS extension name or number in header, if any */
+    char *pheader;	/* Primary header (naxis is 0) */
+    char cext = 0;
+    char *rbrac;	/* Pointer to right bracket if present in file name */
+    char *mwcs;		/* Pointer to WCS name separated by % */
+    char *newhead;	/* New larger header */
+    int nbh0;		/* Length of old too small header */
+    char *pheadend;
+    int inherit = 1;	/* Value of INHERIT keyword in FITS extension header */
+    int extfound = 0;	/* Set to one if desired FITS extension is found */
+    int npcount;
+
+    pheader = NULL;
+    lprim = 0;
+    header = NULL;
+
+    /* Check for FITS WCS specification and ignore for file opening */
+    mwcs = strchr (filename, '%');
+    if (mwcs != NULL)
+	*mwcs = (char) 0;
+
+    /* Check for FITS extension and ignore for file opening */
+    rbrac = NULL;
+    ext = strchr (filename, ',');
+    if (ext == NULL) {
+	ext = strchr (filename, '[');
+	if (ext != NULL) {
+	    rbrac = strchr (filename, ']');
+	    if (rbrac != NULL)
+		*rbrac = (char) 0;
+	    }
+	}
+    if (ext != NULL) {
+	cext = *ext;
+	*ext = (char) 0;
+	}
+
+    /* Open the image file and read the header */
+    if (strncasecmp (filename,"stdin",5)) {
+	fd = -1;
+	fd = fitsropen (filename);
+	}
+#ifndef VMS
+    else {
+	fd = STDIN_FILENO;
+	extnum = -1;
+	}
+#endif
+
+    if (ext != NULL) {
+	if (isnum (ext+1))
+	    extnum = atoi (ext+1);
+	else {
+	    extnum = -2;
+	    strcpy (extnam, ext+1);
+	    }
+	}
+    else
+	extnum = -1;
+
+    /* Repair the damage done to the file-name string during parsing */
+    if (ext != NULL)
+	*ext = cext;
+    if (rbrac != NULL)
+	*rbrac = ']';
+    if (mwcs != NULL)
+	*mwcs = '%';
+
+    if (fd < 0) {
+	fprintf (stderr,"FITSRHEAD:  cannot read file %s\n", filename);
+	return (NULL);
+	}
+
+    nbytes = FITSBLOCK;
+    *nbhead = 0;
+    headend = NULL;
+    nbh = FITSBLOCK * 20 + 4;
+    header = (char *) calloc ((unsigned int) nbh, 1);
+    (void) hlength (header, nbh);
+    headnext = header;
+    nrec = 1;
+    hdu = 0;
+    ibhead = 0;
+
+    /* Read FITS header from input file one FITS block at a time */
+    irec = 0;
+    ibhead = 0;
+    while (irec < 500) {
+	nbytes = FITSBLOCK;
+	for (ntry = 0; ntry < 10; ntry++) {
+	    for (i = 0; i < 2884; i++) fitsbuf[i] = 0;
+	    nbr = read (fd, fitsbuf, nbytes);
+	    if (verbose)
+		fprintf (stderr,"FITSRHEAD: %d header bytes read\n",nbr);
+
+	    /* Short records allowed only if they have the last header line */
+	    if (nbr < nbytes) {
+		headend = ksearch (fitsbuf,"END");
+		if (headend == NULL) {
+		    if (ntry < 9) {
+			if (verbose)
+			    fprintf (stderr,"FITSRHEAD: %d / %d bytes read %d\n",
+				     nbr,nbytes,ntry);
+			}
+		    else {
+			snprintf(fitserrmsg,79,"FITSRHEAD: '%d / %d bytes of header read from %s\n"
+				,nbr,nbytes,filename);
+#ifndef VMS
+			if (fd != STDIN_FILENO)
+#endif
+			    (void)close (fd);
+			free (header);
+			/* if (pheader != NULL)
+			    return (pheader); */
+    			if (extnum != -1 && !extfound) {
+			    *ext = (char) 0;
+			    if (extnum < 0) {
+	    			snprintf (fitserrmsg,79,
+				    "FITSRHEAD: Extension %s not found in file %s",
+				    extnam, filename);
+				}
+			    else {
+	    			snprintf (fitserrmsg,79,
+				    "FITSRHEAD: Extension %d not found in file %s",
+				    extnum, filename);
+				}
+			    *ext = cext;
+			    }
+			else if (hdu > 0) {
+	    		    snprintf (fitserrmsg,79,
+				"FITSRHEAD: No extensions found in file %s", filename);
+			    hdu = 0;
+			    if (pheader != NULL) {
+				*lhead = nbprim;
+				*nbhead = nbprim;
+				return (pheader);
+				}
+			    break;
+			    }
+			else {
+	    		    snprintf (fitserrmsg,79,
+				"FITSRHEAD: No header found in file %s", filename);
+			    }
+			return (NULL);
+			}
+		    }
+		else
+		    break;
+		}
+	    else
+		break;
+	    }
+
+	/* Replace control characters and nulls with spaces */
+	for (i = 0; i < 2880; i++)
+	    if (fitsbuf[i] < 32 || i > nbr) fitsbuf[i] = 32;
+	if (nbr < 2880)
+	    nbr = 2880;
+
+	/* Move current FITS record into header string */
+	strncpy (headnext, fitsbuf, nbr);
+	*nbhead = *nbhead + nbr;
+	nrec = nrec + 1;
+	*(headnext+nbr+1) = 0;
+	ibhead = ibhead + 2880;
+	if (verbose)
+	    fprintf (stderr,"FITSRHEAD: %d bytes in header\n",ibhead);
+
+	/* Check to see if this is the final record in this header */
+	headend = ksearch (fitsbuf,"END");
+	if (headend == NULL) {
+
+	    /* Double size of header buffer if too small */
+	    if (nrec * FITSBLOCK > nbh) {
+		nbh0 = nbh - 4;
+		nbh = (nrec * 2 * FITSBLOCK) + 4;
+		newhead = (char *) calloc (1,(unsigned int) nbh);
+		if (newhead) {
+		    for (i = 0; i < nbh0; i++)
+			newhead[i] = header[i];
+		    free (header);
+		    newhead[nbh-3] = (char) 0;
+		    header = newhead;
+		    (void) hlength (header, nbh);
+		    headnext = header + ((nrec-1) * FITSBLOCK);
+		    }
+		else {
+		    fprintf (stderr,"FITSRHEAD: %d bytes cannot be allocated for header\n",nbh);
+		    exit (1);
+		    }
+		}
+	    else
+		headnext = headnext + FITSBLOCK;
+	    }
+
+	else {
+	    naxis = 0;
+	    hgeti4 (header,"NAXIS",&naxis);
+
+	    /* If header has no data, save it for appending to desired header */
+	    if (naxis < 1) {
+		nbprim = nrec * FITSBLOCK;
+		headend = ksearch (header,"END");
+		lprim = headend + 80 - header;
+		pheader = (char *) calloc ((unsigned int) nbprim, 1);
+		for (i = 0; i < lprim; i++)
+		    pheader[i] = header[i];
+		for (i = lprim; i < nbprim; i++)
+		    pheader[i] = ' ';
+		}
+
+	    /* If header has no data, start with the next record */
+	    if (naxis < 1 && extnum == -1) {
+		extend = 0;
+		hgetl (header,"EXTEND",&extend);
+		if (naxis == 0 && extend) {
+		    headnext = header;
+		    *headend = ' ';
+		    headend = NULL;
+		    nrec = 1;
+		    hdu = hdu + 1;
+		    }
+		else {
+		    break;
+		    }
+		}
+
+	    /* If this is the desired header data unit, keep it */
+	    else if (extnum != -1) {
+		if (extnum > -1 && hdu == extnum) {
+		    extfound = 1;
+		    break;
+		    }
+		else if (extnum < 0) {
+		    extname[0] = 0;
+		    hgets (header, "EXTNAME", 32, extname);
+		    if (!strcmp (extnam,extname)) {
+			extfound = 1;
+			break;
+			}
+		    }
+
+		/* If this is not desired header data unit, skip over data */
+		hdu = hdu + 1;
+		nblock = 0;
+		ibhead = 0;
+		if (naxis > 0) {
+		    ibpix = 0;
+		    hgeti4 (header,"BITPIX",&ibpix);
+		    if (ibpix < 0) {
+			nbpix = -ibpix / 8;
+			}
+		    else {
+			nbpix = ibpix / 8;
+			}
+		    nax1 = 1;
+		    hgeti4 (header,"NAXIS1",&nax1);
+		    nax2 = 1;
+		    if (naxis > 1) {
+			hgeti4 (header,"NAXIS2",&nax2);
+			}
+		    nax3 = 1;
+		    if (naxis > 2) {
+			hgeti4 (header,"NAXIS3",&nax3);
+			}
+		    nax4 = 1;
+		    if (naxis > 3) {
+			hgeti4 (header,"NAXIS4",&nax4);
+			}
+		    nbskip = nax1 * nax2 * nax3 * nax4 * nbpix;
+		    nblock = nbskip / 2880;
+		    if (nblock*2880 < nbskip) {
+			nblock = nblock + 1;
+			}
+		    npcount = 0;
+		    hgeti4 (header,"PCOUNT", &npcount);
+		    if (npcount > 0) {
+			nbskip = nbskip + npcount;
+			nblock = nbskip / 2880;
+			if (nblock*2880 < nbskip)
+			    nblock = nblock + 1;
+			}
+		    }
+		else {
+		    nblock = 0;
+		    }
+		*nbhead = *nbhead + (nblock * 2880);
+
+		/* Set file pointer to beginning of next header/data unit */
+		if (nblock > 0) {
+#ifndef VMS
+		    if (fd != STDIN_FILENO) {
+			ipos = lseek (fd, *nbhead, SEEK_SET);
+			npos = *nbhead;
+			}
+		    else {
+#else
+			{
+#endif
+			ipos = 0;
+			for (i = 0; i < nblock; i++) {
+			    nbytes = FITSBLOCK;
+			    nbr = read (fd, fitsbuf, nbytes);
+			    if (nbr < nbytes) {
+				ipos = ipos + nbr;
+				break;
+				}
+			    else {
+				ipos = ipos + nbytes;
+				}
+			    }
+			npos = nblock * 2880;
+			}
+		    if (ipos < npos) {
+			snprintf (fitserrmsg,79,"FITSRHEAD: %d / %d bytes skipped\n",
+				 ipos,npos);
+			extfound = 0;
+			break;
+			}
+		    }
+		headnext = header;
+		headend = NULL;
+		nrec = 1;
+		}
+	    else {
+		break;
+		}
+	    }
+	}
+
+#ifndef VMS
+    if (fd != STDIN_FILENO)
+	(void)close (fd);
+#endif
+
+/* Print error message and return null if extension not found */
+    if (extnum != -1 && !extfound) {
+	if (extnum < 0)
+	    fprintf (stderr, "FITSRHEAD: Extension %s not found in file %s\n",extnam, filename);
+	else
+	    fprintf (stderr, "FITSRHEAD: Extension %d not found in file %s\n",extnum, filename);
+	if (pheader != NULL) {
+	    free (pheader);
+	    pheader = NULL;
+	    }
+	return (NULL);
+	}
+
+    /* Allocate an extra block for good measure */
+    *lhead = (nrec + 1) * FITSBLOCK;
+    if (*lhead > nbh) {
+	newhead = (char *) calloc (1,(unsigned int) *lhead);
+	for (i = 0; i < nbh; i++)
+	    newhead[i] = header[i];
+	free (header);
+	header = newhead;
+	(void) hlength (header, *lhead);
+	}
+    else
+	*lhead = nbh;
+
+    /* If INHERIT keyword is FALSE, never append primary header */
+    if (hgetl (header, "INHERIT", &inherit)) {
+	if (!inherit && fitsinherit)
+	    fitsinherit = 0;
+	}
+
+    /* Append primary data header to extension header */
+    if (pheader != NULL && extnum != 0 && fitsinherit && hdu > 0) {
+	extname[0] = 0;
+	hgets (header, "XTENSION", 32, extname);
+	if (!strcmp (extname,"IMAGE")) {
+	    strncpy (header, "SIMPLE  ", 8);
+	    hputl (header, "SIMPLE", 1);
+	    }
+	headend = blsearch (header,"END");
+	if (headend == NULL)
+	    headend = ksearch (header, "END");
+	lext = headend - header;
+
+	/* Update primary header for inclusion at end of extension header */
+	hchange (pheader, "SIMPLE", "ROOTHEAD");
+	hchange (pheader, "NEXTEND", "NUMEXT");
+	hdel (pheader, "BITPIX");
+	hdel (pheader, "NAXIS");
+	hdel (pheader, "EXTEND");
+	hputl (pheader, "ROOTEND",1);
+	pheadend = ksearch (pheader,"END");
+	lprim = pheadend + 320 - pheader;
+	if (lext + lprim > nbh) {
+	    nrec = (lext + lprim) / FITSBLOCK;
+	    if (FITSBLOCK*nrec < lext+lprim)
+		nrec = nrec + 1;
+	    *lhead = (nrec+1) * FITSBLOCK;
+	    newhead = (char *) calloc (1,(unsigned int) *lhead);
+	    for (i = 0; i < nbh; i++)
+		newhead[i] = header[i];
+	    free (header);
+	    header = newhead;
+	    headend = header + lext;
+	    (void) hlength (header, *lhead);
+	    }
+	hputs (header,"COMMENT","-------------------------------------------");
+	hputs (header,"COMMENT","Information from Primary Header");
+	hputs (header,"COMMENT","-------------------------------------------");
+	headend = blsearch (header,"END");
+	if (headend == NULL)
+	    headend = ksearch (header, "END");
+	pheader[lprim] = 0;
+	strncpy (headend, pheader, lprim);
+	if (pheader != NULL) {
+	    free (pheader);
+	    pheader = NULL;
+	    }
+	}
+
+    ibhead = *nbhead - ibhead;
+
+    return (header);
+}
+
+
+/* FITSRTAIL -- Read FITS header appended to graphics file */
+
+char *
+fitsrtail (filename, lhead, nbhead)
+
+char	*filename;	/* Name of image file */
+int	*lhead;		/* Allocated length of FITS header in bytes (returned) */
+int	*nbhead;	/* Number of bytes before start of data (returned) */
+			/* This includes all skipped image extensions */
+
+{
+    int fd;
+    char *header;	/* FITS image header (filled) */
+    int nbytes, i, ndiff;
+    int nbr, irec;
+    off_t offset;
+    char *mwcs;		/* Pointer to WCS name separated by % */
+    char *headstart;
+    char *newhead;
+
+    header = NULL;
+
+    /* Check for FITS WCS specification and ignore for file opening */
+    mwcs = strchr (filename, '%');
+    if (mwcs != NULL)
+	*mwcs = (char) 0;
+
+    /* Open the image file and read the header */
+    if (strncasecmp (filename,"stdin",5)) {
+	fd = -1;
+	fd = fitsropen (filename);
+	}
+#ifndef VMS
+    else {
+	fd = STDIN_FILENO;
+	}
+#endif
+
+    /* Repair the damage done to the file-name string during parsing */
+    if (mwcs != NULL)
+	*mwcs = '%';
+
+    if (fd < 0) {
+	fprintf (stderr,"FITSRTAIL:  cannot read file %s\n", filename);
+	return (NULL);
+	}
+
+    nbytes = FITSBLOCK;
+    *nbhead = 0;
+    *lhead = 0;
+
+    /* Read FITS header from end of input file one FITS block at a time */
+    irec = 0;
+    while (irec < 100) {
+	nbytes = FITSBLOCK * (irec + 2);
+	header = (char *) calloc ((unsigned int) nbytes, 1);
+	offset = lseek (fd, -nbytes, SEEK_END);
+	if (offset < 0) {
+	    free (header);
+	    header = NULL;
+	    nbytes = 0;
+	    break;
+	    }
+	for (i = 0; i < nbytes; i++) header[i] = 0;
+	nbr = read (fd, header, nbytes);
+
+	/* Check for SIMPLE at start of header */
+	for (i = 0; i < nbr; i++)
+	    if (header[i] < 32) header[i] = 32;
+	if ((headstart = ksearch (header,"SIMPLE"))) {
+	    if (headstart != header) {
+		ndiff = headstart - header;
+		newhead = (char *) calloc ((unsigned int) nbytes, 1);
+		for (i = 0; i < nbytes-ndiff; i++)
+		    newhead[i] = headstart[i];
+		free (header);
+		header = newhead;
+		}
+	    *lhead = nbytes;
+	    *nbhead = nbytes;
+	    break;
+	    }
+	free (header);
+	}
+    (void) hlength (header, nbytes);
+
+#ifndef VMS
+    if (fd != STDIN_FILENO)
+	(void)close (fd);
+#endif
+
+    return (header);
+}
+
+
+/* FITSRSECT -- Read a piece of a FITS image, having already read the header */
+
+char *
+fitsrsect (filename, header, nbhead, x0, y0, nx, ny, nlog)
+
+char	*filename;	/* Name of FITS image file */
+char	*header;	/* FITS header for image (previously read) */
+int	nbhead;		/* Actual length of image header(s) in bytes */
+int	x0, y0;		/* FITS image coordinate of first pixel */
+int	nx;		/* Number of columns to read (less than NAXIS1) */
+int	ny;		/* Number of rows to read (less than NAXIS2) */
+int	nlog;		/* Note progress mod this rows */
+{
+    int fd;		/* File descriptor */
+    int nbimage, naxis1, naxis2, bytepix, nbread;
+    int bitpix, naxis, nblocks, nbytes, nbr;
+    int x1, y1, nbline, nyleft;
+    off_t impos, nblin;
+    char *image, *imline, *imlast;
+    int ilog = 0;
+    int row;
+
+    /* Open the image file and read the header */
+    if (strncasecmp (filename,"stdin", 5)) {
+	fd = -1;
+
+	fd = fitsropen (filename);
+	if (fd < 0) {
+	    snprintf (fitserrmsg,79, "FITSRSECT:  cannot read file %s\n", filename);
+	    return (NULL);
+	    }
+
+	/* Skip over FITS header and whatever else needs to be skipped */
+	if (lseek (fd, nbhead, SEEK_SET) < 0) {
+	    (void)close (fd);
+	    snprintf (fitserrmsg,79, "FITSRSECT:  cannot skip header of file %s\n",
+		     filename);
+	    return (NULL);
+	    }
+	}
+#ifndef VMS
+    else
+	fd = STDIN_FILENO;
+#endif
+
+    /* Compute size of image in bytes using relevant header parameters */
+    naxis = 1;
+    hgeti4 (header,"NAXIS",&naxis);
+    naxis1 = 1;
+    hgeti4 (header,"NAXIS1",&naxis1);
+    naxis2 = 1;
+    hgeti4 (header,"NAXIS2",&naxis2);
+    bitpix = 0;
+    hgeti4 (header,"BITPIX",&bitpix);
+    if (bitpix == 0) {
+	/* snprintf (fitserrmsg,79, "FITSRSECT:  BITPIX is 0; image not read\n"); */
+	(void)close (fd);
+	return (NULL);
+	}
+    bytepix = bitpix / 8;
+    if (bytepix < 0) bytepix = -bytepix;
+
+    /* Keep X coordinates within image limits */
+    if (x0 < 1)
+	x0 = 1;
+    else if (x0 > naxis1)
+	x0 = naxis1;
+    x1 = x0 + nx - 1;
+    if (x1 < 1)
+	x1 = 1;
+    else if (x1 > naxis1)
+	x1 = naxis1;
+    nx = x1 - x0 + 1;
+
+    /* Keep Y coordinates within image limits */
+    if (y0 < 1)
+	y0 = 1;
+    else if (y0 > naxis2)
+	y0 = naxis2;
+    y1 = y0 + ny - 1;
+    if (y1 < 1)
+	y1 = 1;
+    else if (y1 > naxis2)
+	y1 = naxis2;
+    ny = y1 - y0 + 1;
+
+    /* Number of bytes in output image */
+    nbline = nx * bytepix;
+    nbimage = nbline * ny;
+
+    /* Set number of bytes to integral number of 2880-byte blocks */
+    nblocks = nbimage / FITSBLOCK;
+    if (nblocks * FITSBLOCK < nbimage)
+	nblocks = nblocks + 1;
+    nbytes = nblocks * FITSBLOCK;
+
+    /* Allocate image section to be read */
+    image = (char *) malloc (nbytes);
+    nyleft = ny;
+    imline = image;
+    nbr = 0;
+
+    /* Computer pointer to first byte of input image to read */
+    nblin = naxis1 * bytepix;
+    impos = ((y0 - 1) * nblin) + ((x0 - 1) * bytepix);
+    row = y0 - 1;
+
+    /* Read image section one line at a time */
+    while (nyleft-- > 0) {
+	if (lseek (fd, impos, SEEK_CUR) >= 0) {
+	    nbread = read (fd, imline, nbline);
+	    nbr = nbr + nbread;
+	    impos = nblin - nbread;
+	    imline = imline + nbline;
+	    row++;
+	    if (++ilog == nlog) {
+		ilog = 0;
+		fprintf (stderr, "Row %5d extracted   ", row);
+                (void) putc (13,stderr);
+		}
+	    }
+	}
+    if (nlog)
+	fprintf (stderr, "\n");
+
+    /* Fill rest of image with zeroes */
+    imline = image + nbimage;
+    imlast = image + nbytes;
+    while (imline++ < imlast)
+	*imline = (char) 0;
+
+    /* Byte-reverse image, if necessary */
+    if (imswapped ())
+	imswap (bitpix, image, nbytes);
+
+    return (image);
+}
+
+
+/* FITSRIMAGE -- Read a FITS image */
+
+char *
+fitsrimage (filename, nbhead, header)
+
+char	*filename;	/* Name of FITS image file */
+int	nbhead;		/* Actual length of image header(s) in bytes */
+char	*header;	/* FITS header for image (previously read) */
+{
+    int fd;
+    int nbimage, naxis1, naxis2, bytepix, nbread;
+    int bitpix, naxis, nblocks, nbytes, nbleft, nbr;
+    int simple;
+    char *image, *imleft;
+
+    /* Open the image file and read the header */
+    if (strncasecmp (filename,"stdin", 5)) {
+	fd = -1;
+
+	fd = fitsropen (filename);
+	if (fd < 0) {
+	    snprintf (fitserrmsg,79, "FITSRIMAGE:  cannot read file %s\n", filename);
+	    return (NULL);
+	    }
+
+	/* Skip over FITS header and whatever else needs to be skipped */
+	if (lseek (fd, nbhead, SEEK_SET) < 0) {
+	    (void)close (fd);
+	    snprintf (fitserrmsg,79, "FITSRIMAGE:  cannot skip header of file %s\n",
+		     filename);
+	    return (NULL);
+	    }
+	}
+#ifndef VMS
+    else
+	fd = STDIN_FILENO;
+#endif
+
+    /* If SIMPLE=F in header, simply put post-header part of file in buffer */
+    hgetl (header, "SIMPLE", &simple);
+    if (!simple) {
+	nbytes = getfilesize (filename) - nbhead;
+	if ((image = (char *) malloc (nbytes + 1)) == NULL) {
+	    /* snprintf (fitserrmsg,79, "FITSRIMAGE:  %d-byte image buffer cannot be allocated\n"); */
+	    (void)close (fd);
+	    return (NULL);
+	    }
+	hputi4 (header, "NBDATA", nbytes);
+	nbread = read (fd, image, nbytes);
+	return (image);
+	}
+
+    /* Compute size of image in bytes using relevant header parameters */
+    naxis = 1;
+    hgeti4 (header,"NAXIS",&naxis);
+    naxis1 = 1;
+    hgeti4 (header,"NAXIS1",&naxis1);
+    naxis2 = 1;
+    hgeti4 (header,"NAXIS2",&naxis2);
+    bitpix = 0;
+    hgeti4 (header,"BITPIX",&bitpix);
+    if (bitpix == 0) {
+	/* snprintf (fitserrmsg,79, "FITSRIMAGE:  BITPIX is 0; image not read\n"); */
+	(void)close (fd);
+	return (NULL);
+	}
+    bytepix = bitpix / 8;
+    if (bytepix < 0) bytepix = -bytepix;
+
+    /* If either dimension is one and image is 3-D, read all three dimensions */
+    if (naxis == 3 && (naxis1 ==1 || naxis2 == 1)) {
+	int naxis3;
+	hgeti4 (header,"NAXIS3",&naxis3);
+	nbimage = naxis1 * naxis2 * naxis3 * bytepix;
+	}
+    else
+	nbimage = naxis1 * naxis2 * bytepix;
+
+    /* Set number of bytes to integral number of 2880-byte blocks */
+    nblocks = nbimage / FITSBLOCK;
+    if (nblocks * FITSBLOCK < nbimage)
+	nblocks = nblocks + 1;
+    nbytes = nblocks * FITSBLOCK;
+
+    /* Allocate and read image */
+    image = (char *) malloc (nbytes);
+    nbleft = nbytes;
+    imleft = image;
+    nbr = 0;
+    while (nbleft > 0) {
+	nbread = read (fd, imleft, nbleft);
+	nbr = nbr + nbread;
+#ifndef VMS
+	if (fd == STDIN_FILENO && nbread < nbleft && nbread > 0) {
+	    nbleft = nbleft - nbread;
+	    imleft = imleft + nbread;
+	    }
+	else
+#endif
+	    nbleft = 0;
+	}
+#ifndef VMS
+    if (fd != STDIN_FILENO)
+	(void)close (fd);
+#endif
+    if (nbr < nbimage) {
+	snprintf (fitserrmsg,79, "FITSRIMAGE:  %d of %d bytes read from file %s\n",
+		 nbr, nbimage, filename);
+	return (NULL);
+	}
+
+    /* Byte-reverse image, if necessary */
+    if (imswapped ())
+	imswap (bitpix, image, nbytes);
+
+    return (image);
+}
+
+
+/* FITSRFULL -- Read a FITS image of any dimension */
+
+char *
+fitsrfull (filename, nbhead, header)
+
+char	*filename;	/* Name of FITS image file */
+int	nbhead;		/* Actual length of image header(s) in bytes */
+char	*header;	/* FITS header for image (previously read) */
+{
+    int fd;
+    int nbimage, naxisi, iaxis, bytepix, nbread;
+    int bitpix, naxis, nblocks, nbytes, nbleft, nbr, simple;
+    char keyword[16];
+    char *image, *imleft;
+
+    /* Open the image file and read the header */
+    if (strncasecmp (filename,"stdin", 5)) {
+	fd = -1;
+
+	fd = fitsropen (filename);
+	if (fd < 0) {
+	    snprintf (fitserrmsg,79, "FITSRFULL:  cannot read file %s\n", filename);
+	    return (NULL);
+	    }
+
+	/* Skip over FITS header and whatever else needs to be skipped */
+	if (lseek (fd, nbhead, SEEK_SET) < 0) {
+	    (void)close (fd);
+	    snprintf (fitserrmsg,79, "FITSRFULL:  cannot skip header of file %s\n",
+		     filename);
+	    return (NULL);
+	    }
+	}
+#ifndef VMS
+    else
+	fd = STDIN_FILENO;
+#endif
+
+    /* If SIMPLE=F in header, simply put post-header part of file in buffer */
+    hgetl (header, "SIMPLE", &simple);
+    if (!simple) {
+	nbytes = getfilesize (filename) - nbhead;
+	if ((image = (char *) malloc (nbytes + 1)) == NULL) {
+	    snprintf (fitserrmsg,79, "FITSRFULL:  %d-byte image buffer cannot be allocated\n",nbytes+1);
+	    (void)close (fd);
+	    return (NULL);
+	    }
+	hputi4 (header, "NBDATA", nbytes);
+	nbread = read (fd, image, nbytes);
+	return (image);
+	}
+
+    /* Find number of bytes per pixel */
+    bitpix = 0;
+    hgeti4 (header,"BITPIX",&bitpix);
+    if (bitpix == 0) {
+	snprintf (fitserrmsg,79, "FITSRFULL:  BITPIX is 0; image not read\n");
+	(void)close (fd);
+	return (NULL);
+	}
+    bytepix = bitpix / 8;
+    if (bytepix < 0) bytepix = -bytepix;
+    nbimage = bytepix;
+
+    /* Compute size of image in bytes using relevant header parameters */
+    naxis = 1;
+    hgeti4 (header,"NAXIS",&naxis);
+    for (iaxis = 1; iaxis <= naxis; iaxis++) {
+	sprintf (keyword, "NAXIS%d", iaxis);
+	naxisi = 1;
+	hgeti4 (header,keyword,&naxisi);
+	nbimage = nbimage * naxisi;
+	}
+
+    /* Set number of bytes to integral number of 2880-byte blocks */
+    nblocks = nbimage / FITSBLOCK;
+    if (nblocks * FITSBLOCK < nbimage)
+	nblocks = nblocks + 1;
+    nbytes = nblocks * FITSBLOCK;
+
+    /* Allocate and read image */
+    image = (char *) malloc (nbytes);
+    nbleft = nbytes;
+    imleft = image;
+    nbr = 0;
+    while (nbleft > 0) {
+	nbread = read (fd, imleft, nbleft);
+	nbr = nbr + nbread;
+#ifndef VMS
+	if (fd == STDIN_FILENO && nbread < nbleft && nbread > 0) {
+	    nbleft = nbleft - nbread;
+	    imleft = imleft + nbread;
+	    }
+	else
+#endif
+	    nbleft = 0;
+	}
+#ifndef VMS
+    if (fd != STDIN_FILENO)
+	(void)close (fd);
+#endif
+    if (nbr < nbimage) {
+	snprintf (fitserrmsg,79, "FITSRFULL:  %d of %d image bytes read from file %s\n",
+		 nbr, nbimage, filename);
+	return (NULL);
+	}
+
+    /* Byte-reverse image, if necessary */
+    if (imswapped ())
+	imswap (bitpix, image, nbytes);
+
+    return (image);
+}
+
+
+/* FITSROPEN -- Open a FITS file, returning the file descriptor */
+
+int
+fitsropen (inpath)
+
+char	*inpath;	/* Pathname for FITS tables file to read */
+
+{
+    int ntry;
+    int fd;		/* file descriptor for FITS tables file (returned) */
+    char *ext;		/* extension name or number */
+    char cext = 0;
+    char *rbrac;
+    char *mwcs;		/* Pointer to WCS name separated by % */
+
+/* Check for FITS WCS specification and ignore for file opening */
+    mwcs = strchr (inpath, '%');
+
+/* Check for FITS extension and ignore for file opening */
+    ext = strchr (inpath, ',');
+    rbrac = NULL;
+    if (ext == NULL) {
+	ext = strchr (inpath, '[');
+	if (ext != NULL) {
+	    rbrac = strchr (inpath, ']');
+	    }
+	}
+
+/* Open input file */
+    for (ntry = 0; ntry < 3; ntry++) {
+	if (ext != NULL) {
+	    cext = *ext;
+	    *ext = 0;
+	    }
+	if (rbrac != NULL)
+	    *rbrac = (char) 0;
+	if (mwcs != NULL)
+	    *mwcs = (char) 0;
+	fd = open (inpath, O_RDONLY);
+	if (ext != NULL)
+	    *ext = cext;
+	if (rbrac != NULL)
+	    *rbrac = ']';
+	if (mwcs != NULL)
+	    *mwcs = '%';
+	if (fd >= 0)
+	    break;
+	else if (ntry == 2) {
+	    snprintf (fitserrmsg,79, "FITSROPEN:  cannot read file %s\n", inpath);
+	    return (-1);
+	    }
+	}
+
+    if (verbose)
+	fprintf (stderr,"FITSROPEN:  input file %s opened\n",inpath);
+
+    return (fd);
+}
+
+
+static int offset1=0;
+static int offset2=0;
+
+/* FITSRTOPEN -- Open FITS table file and fill structure with
+ *		 pointers to selected keywords
+ *		 Return file descriptor (-1 if unsuccessful)
+ */
+
+int
+fitsrtopen (inpath, nk, kw, nrows, nchar, nbhead)
+
+char	*inpath;	/* Pathname for FITS tables file to read */
+int	*nk;		/* Number of keywords to use */
+struct Keyword	**kw;	/* Structure for desired entries */
+int	*nrows;		/* Number of rows in table (returned) */
+int	*nchar;		/* Number of characters in one table row (returned) */
+int	*nbhead;	/* Number of characters before table starts */
+
+{
+    char temp[16];
+    int fd;
+    int	lhead;		/* Maximum length in bytes of FITS header */
+    char *header;	/* Header for FITS tables file to read */
+
+/* Read FITS header from input file */
+    header = fitsrhead (inpath, &lhead, nbhead);
+    if (!header) {
+	snprintf (fitserrmsg,79,"FITSRTOPEN:  %s is not a FITS file\n",inpath);
+	return (0);
+	}
+
+/* Make sure this file is really a FITS table file */
+    temp[0] = 0;
+    (void) hgets (header,"XTENSION",16,temp);
+    if (strlen (temp) == 0) {
+	snprintf (fitserrmsg,79,
+		  "FITSRTOPEN:  %s is not a FITS table file\n",inpath);
+	free ((void *) header);
+	return (0);
+	}
+
+/* If it is a FITS file, get table information from the header */
+    else if (!strcmp (temp, "TABLE") || !strcmp (temp, "BINTABLE")) {
+	if (fitsrthead (header, nk, kw, nrows, nchar)) {
+	    snprintf (fitserrmsg,79,
+		      "FITSRTOPEN: Cannot read FITS table from %s\n",inpath);
+	    free ((void *) header);
+	    return (-1);
+	    }
+	else {
+	    fd = fitsropen (inpath);
+	    offset1 = 0;
+	    offset2 = 0;
+	    free ((void *) header);
+	    return (fd);
+	    }
+	}
+
+/* If it is another FITS extension note it and return */
+    else {
+	snprintf (fitserrmsg,79,
+		  "FITSRTOPEN:  %s is a %s extension, not table\n",
+		  inpath, temp);
+	free ((void *) header);
+	return (0);
+	}
+}
+
+static struct Keyword *pw;	/* Structure for all entries */
+static int *lpnam;		/* length of name for each field */
+static int bfields = 0;
+
+/* FITSRTHEAD -- From FITS table header, read pointers to selected keywords */
+
+int
+fitsrthead (header, nk, kw, nrows, nchar)
+
+char	*header;	/* Header for FITS tables file to read */
+int	*nk;		/* Number of keywords to use */
+struct Keyword	**kw;	/* Structure for desired entries */
+int	*nrows;		/* Number of rows in table (returned) */
+int	*nchar;		/* Number of characters in one table row (returned) */
+
+{
+    struct Keyword *rw;	/* Structure for desired entries */
+    int nfields;
+    int ifield, ik, i, ikf, ltform, kl;
+    char *h0, *h1, *tf1, *tf2;
+    char tname[12];
+    char temp[16];
+    char tform[16];
+    int tverb;
+    int bintable = 0;
+
+    h0 = header;
+
+/* Make sure this is really a FITS table file header */
+    temp[0] = 0;
+    hgets (header,"XTENSION",16,temp);
+    if (strlen (temp) == 0) {
+	snprintf (fitserrmsg,79, "FITSRTHEAD:  Not a FITS table header\n");
+	return (-1);
+	}
+    else if (!strcmp (temp, "BINTABLE")) {
+	bintable = 1;
+	}
+    else if (strcmp (temp, "TABLE")) {
+	snprintf (fitserrmsg,79, "FITSRTHEAD:  %s extension, not TABLE\n",temp);
+	return (-1);
+	}
+
+/* Get table size from FITS header */
+    *nchar = 0;
+    hgeti4 (header,"NAXIS1",nchar);
+    *nrows = 0;
+    hgeti4 (header,"NAXIS2", nrows);
+    if (*nrows <= 0 || *nchar <= 0) {
+	snprintf (fitserrmsg,79, "FITSRTHEAD: cannot read %d x %d table\n",
+		 *nrows,*nchar);
+	return (-1);
+	}
+
+/* Set up table for access to individual fields */
+    nfields = 0;
+    hgeti4 (header,"TFIELDS",&nfields);
+    if (verbose)
+	fprintf (stderr, "FITSRTHEAD: %d fields per table entry\n", nfields);
+    if (nfields > bfields) {
+	if (bfields > 0)
+	    free ((void *)pw);
+	pw = (struct Keyword *) calloc (nfields, sizeof(struct Keyword));
+	if (pw == NULL) {
+	    snprintf (fitserrmsg,79,"FITSRTHEAD: cannot allocate table structure\n");
+	    return (-1);
+	    }
+	if (bfields > 0)
+	    free ((void *)lpnam);
+	lpnam = (int *) calloc (nfields, sizeof(int));
+	if (lpnam == NULL) {
+	    snprintf (fitserrmsg,79,"FITSRTHEAD: cannot allocate length structure\n");
+	    return (-1);
+	    }
+	bfields = nfields;
+	}
+
+    tverb = verbose;
+    verbose = 0;
+    ikf = 0;
+
+    for (ifield = 0; ifield < nfields; ifield++) {
+
+    /* Name of field */
+	for (i = 0; i < 12; i++) tname[i] = 0;
+	sprintf (tname, "TTYPE%d", ifield+1);;
+	temp[0] = 0;
+	h1 = ksearch (h0,tname);
+	h0 = h1;
+	hgets (h0,tname,16,temp);
+	strcpy (pw[ifield].kname,temp);
+	pw[ifield].lname = strlen (pw[ifield].kname);
+
+    /* Sequence of field on line */
+	pw[ifield].kn = ifield + 1;
+
+    /* First column of field */
+	if (bintable)
+	    pw[ifield].kf = ikf;
+	else {
+	    for (i = 0; i < 12; i++) tname[i] = 0;
+	    sprintf (tname, "TBCOL%d", ifield+1);
+	    pw[ifield].kf = 0;
+	    hgeti4 (h0,tname, &pw[ifield].kf);
+	    }
+
+    /* Length of field */
+	for (i = 0; i < 12; i++) tname[i] = 0;
+	sprintf (tname, "TFORM%d", ifield+1);;
+	tform[0] = 0;
+	hgets (h0,tname,16,tform);
+	strcpy (pw[ifield].kform, tform);
+	ltform = strlen (tform);
+	if (tform[ltform-1] == 'A') {
+	    pw[ifield].kform[0] = 'A';
+	    for (i = 0; i < ltform-1; i++)
+		pw[ifield].kform[i+1] = tform[i];
+	    pw[ifield].kform[ltform] = (char) 0;
+	    tf1 = pw[ifield].kform + 1;
+	    kl = atof (tf1);
+	    }
+	else if (!strcmp (tform,"I"))
+	    kl = 2;
+	else if (!strcmp (tform, "J"))
+	    kl = 4;
+	else if (!strcmp (tform, "E"))
+	    kl = 4;
+	else if (!strcmp (tform, "D"))
+	    kl = 8;
+	else {
+	    tf1 = tform + 1;
+	    tf2 = strchr (tform,'.');
+	    if (tf2 != NULL)
+		*tf2 = ' ';
+	    kl = atoi (tf1);
+	    }
+	pw[ifield].kl = kl;
+	ikf = ikf + kl;
+	}
+
+/* Set up table for access to desired fields */
+    verbose = tverb;
+    if (verbose)
+	fprintf (stderr, "FITSRTHEAD: %d keywords read\n", *nk);
+
+/* If nk = 0, allocate and return structures for all table fields */
+    if (*nk <= 0) {
+	*kw = pw;
+	*nk = nfields;
+	return (0);
+	}
+    else
+	rw = *kw;
+
+/* Find each desired keyword in the header */
+    for (ik = 0; ik < *nk; ik++) {
+	if (rw[ik].kn <= 0) {
+	    for (ifield = 0; ifield < nfields; ifield++) {
+		if (rw[ik].lname != pw[ifield].lname)
+		    continue;
+		if (strcmp (pw[ifield].kname, rw[ik].kname) == 0) {
+		    break;
+		    }
+		}
+	    }
+	else
+	    ifield = rw[ik].kn - 1;
+
+/* Set pointer, lentth, and name in returned array of structures */
+	rw[ik].kn = ifield + 1;
+	rw[ik].kf = pw[ifield].kf - 1;
+	rw[ik].kl = pw[ifield].kl;
+	strcpy (rw[ik].kform, pw[ifield].kform);
+	strcpy (rw[ik].kname, pw[ifield].kname);
+	}
+
+    return (0);
+}
+
+
+int
+fitsrtline (fd, nbhead, lbuff, tbuff, irow, nbline, line)
+
+int	fd;		/* File descriptor for FITS file */
+int	nbhead;		/* Number of bytes in FITS header */
+int	lbuff;		/* Number of bytes in table buffer */
+char	*tbuff;		/* FITS table buffer */
+int	irow;		/* Number of table row to read */
+int	nbline;		/* Number of bytes to read for this line */
+char	*line;		/* One line of FITS table (returned) */
+
+{
+    int nbuff, nlbuff;
+    int nbr = 0;
+    int offset, offend, ntry, ioff;
+    char *tbuff1;
+
+    offset = nbhead + (nbline * irow);
+    offend = offset + nbline - 1;
+
+/* Read a new buffer of the FITS table into memory if needed */
+    if (offset < offset1 || offend > offset2) {
+	nlbuff = lbuff / nbline;
+	nbuff = nlbuff * nbline;
+	for (ntry = 0; ntry < 3; ntry++) {
+	    ioff = lseek (fd, offset, SEEK_SET);
+	    if (ioff < offset) {
+		if (ntry == 2)
+		    return (0);
+		else
+		    continue;
+		}
+	    nbr = read (fd, tbuff, nbuff);
+	    if (nbr < nbline) {
+		if (verbose)
+		    fprintf (stderr, "FITSRTLINE: %d / %d bytes read %d\n",
+				nbr,nbuff,ntry);
+		if (ntry == 2)
+		    return (nbr);
+		}
+	    else
+		break;
+	    }
+	offset1 = offset;
+	offset2 = offset + nbr - 1;
+	strncpy (line, tbuff, nbline);
+	return (nbline);
+	}
+    else {
+	tbuff1 = tbuff + (offset - offset1);
+	strncpy (line, tbuff1, nbline);
+	return (nbline);
+	}
+}
+
+
+void
+fitsrtlset ()
+{
+    offset1 = 0;
+    offset2 = 0;
+    return;
+}
+
+
+/* FTGETI2 -- Extract n'th column from FITS table line as short */
+
+short
+ftgeti2 (entry, kw)
+
+char	*entry;		/* Row or entry from table */
+struct Keyword *kw;	/* Table column information from FITS header */
+{
+    char temp[30];
+    short i;
+    int j;
+    float r;
+    double d;
+
+    if (ftgetc (entry, kw, temp, 30)) {
+	if (!strcmp (kw->kform, "I"))
+	    moveb (temp, (char *) &i, 2, 0, 0);
+	else if (!strcmp (kw->kform, "J")) {
+	    moveb (temp, (char *) &j, 4, 0, 0);
+	    i = (short) j;
+	    }
+	else if (!strcmp (kw->kform, "E")) {
+	    moveb (temp, (char *) &r, 4, 0, 0);
+	    i = (short) r;
+	    }
+	else if (!strcmp (kw->kform, "D")) {
+	    moveb (temp, (char *) &d, 8, 0, 0);
+	    i = (short) d;
+	    }
+	else
+	    i = (short) atof (temp);
+	return (i);
+	}
+    else
+	return ((short) 0);
+}
+
+
+/* FTGETI4 -- Extract n'th column from FITS table line as int */
+
+int
+ftgeti4 (entry, kw)
+
+char	*entry;		/* Row or entry from table */
+struct Keyword *kw;	/* Table column information from FITS header */
+{
+    char temp[30];
+    short i;
+    int j;
+    float r;
+    double d;
+
+    if (ftgetc (entry, kw, temp, 30)) {
+	if (!strcmp (kw->kform, "I")) {
+	    moveb (temp, (char *) &i, 2, 0, 0);
+	    j = (int) i;
+	    }
+	else if (!strcmp (kw->kform, "J"))
+	    moveb (temp, (char *) &j, 4, 0, 0);
+	else if (!strcmp (kw->kform, "E")) {
+	    moveb (temp, (char *) &r, 4, 0, 0);
+	    j = (int) r;
+	    }
+	else if (!strcmp (kw->kform, "D")) {
+	    moveb (temp, (char *) &d, 8, 0, 0);
+	    j = (int) d;
+	    }
+	else
+	    j = (int) atof (temp);
+	return (j);
+	}
+    else
+	return (0);
+}
+
+
+/* FTGETR4 -- Extract n'th column from FITS table line as float */
+
+float
+ftgetr4 (entry, kw)
+
+char	*entry;		/* Row or entry from table */
+struct Keyword *kw;	/* Table column information from FITS header */
+{
+    char temp[30];
+    short i;
+    int j;
+    float r;
+    double d;
+
+    if (ftgetc (entry, kw, temp, 30)) {
+	if (!strcmp (kw->kform, "I")) {
+	    moveb (temp, (char *) &i, 2, 0, 0);
+	    r = (float) i;
+	    }
+	else if (!strcmp (kw->kform, "J")) {
+	    moveb (temp, (char *) &j, 4, 0, 0);
+	    r = (float) j;
+	    }
+	else if (!strcmp (kw->kform, "E"))
+	    moveb (temp, (char *) &r, 4, 0, 0);
+	else if (!strcmp (kw->kform, "D")) {
+	    moveb (temp, (char *) &d, 8, 0, 0);
+	    r = (float) d;
+	    }
+	else
+	    r = (float) atof (temp);
+	return (r);
+	}
+    else
+	return ((float) 0.0);
+}
+
+
+/* FTGETR8 -- Extract n'th column from FITS table line as double */
+
+double
+ftgetr8 (entry, kw)
+
+char	*entry;		/* Row or entry from table */
+struct Keyword *kw;	/* Table column information from FITS header */
+{
+    char temp[30];
+    short i;
+    int j;
+    float r;
+    double d;
+
+    if (ftgetc (entry, kw, temp, 30)) {
+	if (!strcmp (kw->kform, "I")) {
+	    moveb (temp, (char *) &i, 2, 0, 0);
+	    d = (double) i;
+	    }
+	else if (!strcmp (kw->kform, "J")) {
+	    moveb (temp, (char *) &j, 4, 0, 0);
+	    d = (double) j;
+	    }
+	else if (!strcmp (kw->kform, "E")) {
+	    moveb (temp, (char *) &r, 4, 0, 0);
+	    d = (double) r;
+	    }
+	else if (!strcmp (kw->kform, "D"))
+	    moveb (temp, (char *) &d, 8, 0, 0);
+	else
+	    d = atof (temp);
+	return (d);
+	}
+    else
+	return ((double) 0.0);
+}
+
+
+/* FTGETC -- Extract n'th column from FITS table line as character string */
+
+int
+ftgetc (entry, kw, string, maxchar)
+
+char	*entry;		/* Row or entry from table */
+struct Keyword *kw;	/* Table column information from FITS header */
+char	*string;	/* Returned string */
+int	maxchar;	/* Maximum number of characters in returned string */
+{
+    int length = maxchar;
+
+    if (kw->kl < length)
+	length = kw->kl;
+    if (length > 0) {
+	strncpy (string, entry+kw->kf, length);
+	string[length] = 0;
+	return ( 1 );
+	}
+    else
+	return ( 0 );
+}
+
+extern int errno;
+
+
+/*FITSWIMAGE -- Write FITS header and image */
+
+int
+fitswimage (filename, header, image)
+
+char	*filename;	/* Name of FITS image file */
+char	*header;	/* FITS image header */
+char	*image;		/* FITS image pixels */
+
+{
+    int fd;
+
+    /* Open the output file */
+    if (strcasecmp (filename,"stdout") ) {
+
+	if (!access (filename, 0)) {
+	    fd = open (filename, O_WRONLY);
+	    if (fd < 3) {
+		snprintf (fitserrmsg,79, "FITSWIMAGE:  file %s not writeable\n", filename);
+		return (0);
+		}
+	    }
+	else {
+	    fd = open (filename, O_RDWR+O_CREAT, 0666);
+	    if (fd < 3) {
+		snprintf (fitserrmsg,79, "FITSWIMAGE:  cannot create file %s\n", filename);
+		return (0);
+		}
+	    }
+	}
+#ifndef VMS
+    else
+	fd = STDOUT_FILENO;
+#endif
+
+    return (fitswhdu (fd, filename, header, image));
+}
+
+
+/*FITSWEXT -- Write FITS header and image as extension to a file */
+
+int
+fitswext (filename, header, image)
+
+char	*filename;	/* Name of IFTS image file */
+char	*header;	/* FITS image header */
+char	*image;		/* FITS image pixels */
+
+{
+    int fd;
+
+    /* Open the output file */
+    if (strcasecmp (filename,"stdout") ) {
+
+	if (!access (filename, 0)) {
+	    fd = open (filename, O_WRONLY);
+	    if (fd < 3) {
+		snprintf (fitserrmsg,79, "FITSWEXT:  file %s not writeable\n",
+			 filename);
+		return (0);
+		}
+	    }
+	else {
+	    fd = open (filename, O_APPEND, 0666);
+	    if (fd < 3) {
+		snprintf (fitserrmsg,79, "FITSWEXT:  cannot append to file %s\n",
+			 filename);
+		return (0);
+		}
+	    }
+	}
+#ifndef VMS
+    else
+	fd = STDOUT_FILENO;
+#endif
+
+    return (fitswhdu (fd, filename, header, image));
+}
+
+
+/* FITSWHDU -- Write FITS head and image as extension */
+
+int
+fitswhdu (fd, filename, header, image)
+
+int	fd;		/* File descriptor */
+char	*filename;	/* Name of IFTS image file */
+char	*header;	/* FITS image header */
+char	*image;		/* FITS image pixels */
+{
+    int nbhead, nbimage, nblocks, bytepix, i, nbhw;
+    int bitpix, naxis, iaxis, naxisi, nbytes, nbw, nbpad, nbwp, simple;
+    char *endhead, *padding;
+    double bzero, bscale;
+    char keyword[32];
+
+    /* Change BITPIX=-16 files to BITPIX=16 with BZERO and BSCALE */
+    bitpix = 0;
+    hgeti4 (header,"BITPIX",&bitpix);
+    if (bitpix == -16) {
+	if (!hgetr8 (header, "BZERO", &bzero) &&
+	    !hgetr8 (header, "BSCALE", &bscale)) {
+	    bitpix = 16;
+	    hputi4 (header, "BITPIX", bitpix);
+	    hputr8 (header, "BZERO", 32768.0);
+	    hputr8 (header, "BSCALE", 1.0);
+	    }
+	}
+
+    /* Write header to file */
+    endhead = ksearch (header,"END") + 80;
+    nbhead = endhead - header;
+    nbhw = write (fd, header, nbhead);
+    if (nbhw < nbhead) {
+	snprintf (fitserrmsg,79, "FITSWHDU:  wrote %d / %d bytes of header to file %s\n",
+		 nbhw, nbhead, filename);
+	(void)close (fd);
+	return (0);
+	}
+
+    /* Write extra spaces to make an integral number of 2880-byte blocks */
+    nblocks = nbhead / FITSBLOCK;
+    if (nblocks * FITSBLOCK < nbhead)
+	nblocks = nblocks + 1;
+    nbytes = nblocks * FITSBLOCK;
+    nbpad = nbytes - nbhead;
+    padding = (char *)calloc (1, nbpad);
+    for (i = 0; i < nbpad; i++)
+	padding[i] = ' ';
+    nbwp = write (fd, padding, nbpad);
+    if (nbwp < nbpad) {
+	snprintf (fitserrmsg,79, "FITSWHDU:  wrote %d / %d bytes of header padding to file %s\n",
+		 nbwp, nbpad, filename);
+	(void)close (fd);
+	return (0);
+	}
+    nbhw = nbhw + nbwp;
+    free (padding);
+
+    /* Return if file has no data */
+    if (bitpix == 0 || image == NULL) {
+	/* snprintf (fitserrmsg,79, "FITSWHDU:  BITPIX is 0; image not written\n"); */
+	(void)close (fd);
+	return (0);
+	}
+
+    /* If SIMPLE=F in header, just write whatever is in the buffer */
+    hgetl (header, "SIMPLE", &simple);
+    if (!simple) {
+	hgeti4 (header, "NBDATA", &nbytes);
+	nbimage = nbytes;
+	}
+
+    else {
+
+	/* Compute size of pixel in bytes */
+	bytepix = bitpix / 8;
+	if (bytepix < 0) bytepix = -bytepix;
+	nbimage = bytepix;
+
+	/* Compute size of image in bytes using relevant header parameters */
+	naxis = 1;
+	hgeti4 (header,"NAXIS",&naxis);
+	for (iaxis = 1; iaxis <= naxis; iaxis++) {
+	    sprintf (keyword, "NAXIS%d", iaxis);
+	    naxisi = 1;
+	    hgeti4 (header,keyword,&naxisi);
+	    nbimage = nbimage * naxisi;
+	    }
+
+	/* Number of bytes to write is an integral number of FITS blocks */
+	nblocks = nbimage / FITSBLOCK;
+	if (nblocks * FITSBLOCK < nbimage)
+	    nblocks = nblocks + 1;
+	nbytes = nblocks * FITSBLOCK;
+
+	/* Byte-reverse image before writing, if necessary */
+	if (imswapped ())
+	    imswap (bitpix, image, nbimage);
+	}
+
+    /* Write image to file */
+    nbw = write (fd, image, nbimage);
+    if (nbw < nbimage) {
+	snprintf (fitserrmsg,79, "FITSWHDU:  wrote %d / %d bytes of image to file %s\n",
+		 nbw, nbimage, filename);
+	return (0);
+	}
+
+    /* Write extra zeroes to make an integral number of 2880-byte blocks */
+    nbpad = nbytes - nbimage;
+    if (nbpad > 0) {
+	padding = (char *)calloc (1, nbpad);
+	nbwp = write (fd, padding, nbpad);
+	if (nbwp < nbpad) {
+	    snprintf (fitserrmsg,79, "FITSWHDU:  wrote %d / %d bytes of image padding to file %s\n",
+		 nbwp, nbpad, filename);
+	    (void)close (fd);
+	    return (0);
+	    }
+	free (padding);
+	}
+    else
+	nbwp = 0;
+
+    (void)close (fd);
+
+    /* Byte-reverse image after writing, if necessary */
+    if (imswapped ())
+	imswap (bitpix, image, nbimage);
+
+    nbw = nbw + nbwp + nbhw;
+    return (nbw);
+}
+
+
+/*FITSCIMAGE -- Write FITS header and copy FITS image
+		Return number of bytes in output image, 0 if failure */
+
+int
+fitscimage (filename, header, filename0)
+
+char	*filename;	/* Name of output FITS image file */
+char	*header;	/* FITS image header */
+char	*filename0;	/* Name of input FITS image file */
+
+{
+    int fdout, fdin;
+    int nbhead, nbimage, nblocks, bytepix;
+    int bitpix, naxis, naxis1, naxis2, nbytes, nbw, nbpad, nbwp;
+    char *endhead, *lasthead, *padding;
+    char *image;	/* FITS image pixels */
+    char *oldhead;	/* Input file image header */
+    int nbhead0;	/* Length of input file image header */
+    int lhead0;
+    int nbbuff, nbuff, ibuff, nbr, nbdata;
+
+    /* Compute size of image in bytes using relevant header parameters */
+    naxis = 1;
+    hgeti4 (header, "NAXIS", &naxis);
+    naxis1 = 1;
+    hgeti4 (header, "NAXIS1", &naxis1);
+    naxis2 = 1;
+    hgeti4 (header, "NAXIS2", &naxis2);
+    hgeti4 (header, "BITPIX", &bitpix);
+    bytepix = bitpix / 8;
+    if (bytepix < 0) bytepix = -bytepix;
+
+    /* If either dimension is one and image is 3-D, read all three dimensions */
+    if (naxis == 3 && (naxis1 ==1 || naxis2 == 1)) {
+	int naxis3;
+	hgeti4 (header,"NAXIS3",&naxis3);
+	nbimage = naxis1 * naxis2 * naxis3 * bytepix;
+	}
+    else
+	nbimage = naxis1 * naxis2 * bytepix;
+
+    nblocks = nbimage / FITSBLOCK;
+    if (nblocks * FITSBLOCK < nbimage)
+	nblocks = nblocks + 1;
+    nbytes = nblocks * FITSBLOCK;
+
+    /* Allocate image buffer */
+    nbbuff = FITSBLOCK * 100;
+    if (nbytes < nbbuff)
+	nbbuff = nbytes;
+    image = (char *) calloc (1, nbbuff);
+    nbuff = nbytes / nbbuff;
+    if (nbytes > nbuff * nbbuff)
+	nbuff = nbuff + 1;
+
+    /* Read input file header */
+    if ((oldhead = fitsrhead (filename0, &lhead0, &nbhead0)) == NULL) {
+	snprintf (fitserrmsg, 79,"FITSCIMAGE: header of input file %s cannot be read\n",
+		 filename0);
+	return (0);
+	}
+
+    /* Find size of output header */
+    nbhead = fitsheadsize (header);
+
+    /* If overwriting, be more careful if new header is longer than old */
+    if (!strcmp (filename, filename0) && nbhead > nbhead0) {
+	if ((image = fitsrimage (filename0, nbhead0, oldhead)) == NULL) {
+	    snprintf (fitserrmsg,79, "FITSCIMAGE:  cannot read image from file %s\n",
+		     filename0);
+	    free (oldhead);
+	    return (0);
+	    }
+	return (fitswimage (filename, header, image));
+	}
+    free (oldhead);
+
+    /* Open the input file and skip over the header */
+    if (strcasecmp (filename0,"stdin")) {
+	fdin = -1;
+	fdin = fitsropen (filename0);
+	if (fdin < 0) {
+	    snprintf (fitserrmsg, 79,"FITSCIMAGE:  cannot read file %s\n", filename0);
+	    return (0);
+	    }
+
+	/* Skip over FITS header */
+	if (lseek (fdin, nbhead0, SEEK_SET) < 0) {
+	    (void)close (fdin);
+	    snprintf (fitserrmsg,79, "FITSCIMAGE:  cannot skip header of file %s\n",
+		     filename0);
+	    return (0);
+	    }
+	}
+#ifndef VMS
+    else
+	fdin = STDIN_FILENO;
+#endif
+
+    /* Open the output file */
+    if (!access (filename, 0)) {
+	fdout = open (filename, O_WRONLY);
+	if (fdout < 3) {
+	    snprintf (fitserrmsg,79, "FITSCIMAGE:  file %s not writeable\n", filename);
+	    return (0);
+	    }
+	}
+    else {
+	fdout = open (filename, O_RDWR+O_CREAT, 0666);
+	if (fdout < 3) {
+	    snprintf (fitserrmsg,79, "FITSCHEAD:  cannot create file %s\n", filename);
+	    return (0);
+	    }
+	}
+
+    /* Pad header with spaces */
+    endhead = ksearch (header,"END") + 80;
+    lasthead = header + nbhead;
+    while (endhead < lasthead)
+	*(endhead++) = ' ';
+
+    /* Write header to file */
+    nbw = write (fdout, header, nbhead);
+    if (nbw < nbhead) {
+	snprintf (fitserrmsg, 79,"FITSCIMAGE:  wrote %d / %d bytes of header to file %s\n",
+		 nbw, nbytes, filename);
+	(void)close (fdout);
+	(void)close (fdin);
+	return (0);
+	}
+
+    /* Return if no data */
+    if (bitpix == 0) {
+	(void)close (fdout);
+	(void)close (fdin);
+	return (nbhead);
+	}
+
+    nbdata = 0;
+    for (ibuff = 0; ibuff < nbuff; ibuff++) {
+	nbr = read (fdin, image, nbbuff);
+	if (nbr > 0) {
+	    nbw = write (fdout, image, nbr);
+	    nbdata = nbdata + nbw;
+	    }
+	}
+
+    /* Write extra to make integral number of 2880-byte blocks */
+    nblocks = nbdata / FITSBLOCK;
+    if (nblocks * FITSBLOCK < nbdata)
+	nblocks = nblocks + 1;
+    nbytes = nblocks * FITSBLOCK;
+    nbpad = nbytes - nbdata;
+    padding = (char *)calloc (1,nbpad);
+    nbwp = write (fdout, padding, nbpad);
+    nbw = nbdata + nbwp;
+    free (padding);
+
+    (void)close (fdout);
+    (void)close (fdin);
+
+    if (nbw < nbimage) {
+	snprintf (fitserrmsg, 79, "FITSWIMAGE:  wrote %d / %d bytes of image to file %s\n",
+		 nbw, nbimage, filename);
+	return (0);
+	}
+    else
+	return (nbw);
+}
+
+
+/* FITSWHEAD -- Write FITS header and keep file open for further writing */
+
+int
+fitswhead (filename, header)
+
+char	*filename;	/* Name of IFTS image file */
+char	*header;	/* FITS image header */
+
+{
+    int fd;
+    int nbhead, nblocks;
+    int nbytes, nbw;
+    char *endhead, *lasthead;
+
+    /* Open the output file */
+    if (!access (filename, 0)) {
+	fd = open (filename, O_WRONLY);
+	if (fd < 3) {
+	    snprintf (fitserrmsg, 79, "FITSWHEAD:  file %s not writeable\n", filename);
+	    return (0);
+	    }
+	}
+    else {
+	fd = open (filename, O_RDWR+O_CREAT, 0666);
+	if (fd < 3) {
+	    snprintf (fitserrmsg, 79, "FITSWHEAD:  cannot create file %s\n", filename);
+	    return (0);
+	    }
+	}
+
+    /* Write header to file */
+    endhead = ksearch (header,"END") + 80;
+    nbhead = endhead - header;
+    nblocks = nbhead / FITSBLOCK;
+    if (nblocks * FITSBLOCK < nbhead)
+	nblocks = nblocks + 1;
+    nbytes = nblocks * FITSBLOCK;
+
+    /* Pad header with spaces */
+    lasthead = header + nbytes;
+    while (endhead < lasthead)
+	*(endhead++) = ' ';
+    
+    nbw = write (fd, header, nbytes);
+    if (nbw < nbytes) {
+	fprintf (stderr, "FITSWHEAD:  wrote %d / %d bytes of header to file %s\n",
+		 nbw, nbytes, filename);
+	(void)close (fd);
+	return (0);
+	}
+    return (fd);
+}
+
+
+/* FITSWEXHEAD -- Write FITS header in place */
+
+int
+fitswexhead (filename, header)
+
+char	*filename;	/* Name of FITS image file with ,extension */
+char	*header;	/* FITS image header */
+
+{
+    int fd;
+    int nbhead, lhead;
+    int nbw, nbnew, nbold;
+    char *endhead, *lasthead, *oldheader;
+    char *ext, cext;
+
+    /* Compare size of existing header to size of new header */
+    fitsinherit = 0;
+    oldheader = fitsrhead (filename, &lhead, &nbhead);
+    if (oldheader == NULL) {
+	snprintf (fitserrmsg, 79, "FITSWEXHEAD:  file %s cannot be read\n", filename);
+	return (-1);
+	}
+    nbold = fitsheadsize (oldheader);
+    nbnew = fitsheadsize (header);
+
+    /* Return if the new header is bigger than the old header */
+    if (nbnew > nbold) {
+	snprintf (fitserrmsg, 79, "FITSWEXHEAD:  old header %d bytes, new header %d bytes\n", nbold,nbnew);
+	free (oldheader);
+	oldheader = NULL;
+	return (-1);
+	}
+
+    /* Add blank lines if new header is smaller than the old header */
+    else if (nbnew < nbold) {
+	strcpy (oldheader, header);
+	endhead = ksearch (oldheader,"END");
+	lasthead = oldheader + nbold;
+	while (endhead < lasthead)
+	    *(endhead++) = ' ';
+	strncpy (lasthead-80, "END", 3);
+	}
+
+    /* Pad header with spaces */
+    else {
+	endhead = ksearch (header,"END") + 80;
+	lasthead = header + nbnew;
+	while (endhead < lasthead)
+	    *(endhead++) = ' ';
+	strncpy (oldheader, header, nbnew);
+	}
+
+    /* Check for FITS extension and ignore for file opening */
+    ext = strchr (filename, ',');
+    if (ext == NULL)
+	ext = strchr (filename, '[');
+    if (ext != NULL) {
+	cext = *ext;
+	*ext = (char) 0;
+	}
+
+    /* Open the output file */
+    fd = open (filename, O_WRONLY);
+    if (ext != NULL)
+	*ext = cext;
+    if (fd < 3) {
+	snprintf (fitserrmsg, 79, "FITSWEXHEAD:  file %s not writeable\n", filename);
+	return (-1);
+	}
+
+    /* Skip to appropriate place in file */
+    (void) lseek (fd, ibhead, SEEK_SET);
+
+    /* Write header to file */
+    nbw = write (fd, oldheader, nbold);
+    (void)close (fd);
+    free (oldheader);
+    oldheader = NULL;
+    if (nbw < nbold) {
+	fprintf (stderr, "FITSWHEAD:  wrote %d / %d bytes of header to file %s\n",
+		 nbw, nbold, filename);
+	return (-1);
+	}
+    return (0);
+}
+
+
+/* ISFITS -- Return 1 if FITS file, else 0 */
+int
+isfits (filename)
+
+char    *filename;      /* Name of file for which to find size */
+{
+    int diskfile;
+    char keyword[16];
+    char *comma;
+    int nbr;
+
+    /* First check to see if this is an assignment */
+    if (strchr (filename, '='))
+	return (0);
+
+    /* Check for stdin (input from pipe) */
+    else if (!strcasecmp (filename,"stdin"))
+	return (1);
+
+    /* Then check file extension
+    else if (strsrch (filename, ".fit") ||
+	strsrch (filename, ".fits") ||
+	strsrch (filename, ".fts"))
+	return (1); */
+
+    /* If no FITS file extension, try opening the file */
+    else {
+	if ((comma = strchr (filename,',')))
+	    *comma = (char) 0;
+	if ((diskfile = open (filename, O_RDONLY)) < 0) {
+	    if (comma)
+		*comma = ',';
+	    return (0);
+	    }
+	else {
+	    nbr = read (diskfile, keyword, 8);
+	    if (comma)
+		*comma = ',';
+	    close (diskfile);
+	    if (nbr < 8)
+		return (0);
+	    else if (!strncmp (keyword, "SIMPLE", 6))
+		return (1);
+	    else
+		return (0);
+	    }
+	}
+}
+
+
+/* FITSHEADSIZE -- Find size of FITS header */
+
+int
+fitsheadsize (header)
+
+char	*header;	/* FITS header */
+{
+    char *endhead;
+    int nbhead, nblocks;
+
+    endhead = ksearch (header,"END") + 80;
+    nbhead = endhead - header;
+    nblocks = nbhead / FITSBLOCK;
+    if (nblocks * FITSBLOCK < nbhead)
+        nblocks = nblocks + 1;
+    return (nblocks * FITSBLOCK);
+}
+
+
+/* Print error message */
+void
+fitserr ()
+{   fprintf (stderr, "%s\n",fitserrmsg);
+    return; }
+
+
+/* MOVEB -- Copy nbytes bytes from source+offs to dest+offd (any data type) */
+
+void
+moveb (source, dest, nbytes, offs, offd)
+
+char *source;	/* Pointer to source */
+char *dest;	/* Pointer to destination */
+int nbytes;	/* Number of bytes to move */
+int offs;	/* Offset in bytes in source from which to start copying */
+int offd;	/* Offset in bytes in destination to which to start copying */
+{
+char *from, *last, *to;
+        from = source + offs;
+        to = dest + offd;
+        last = from + nbytes;
+        while (from < last) *(to++) = *(from++);
+        return;
+}
+
+/*
+ * Feb  8 1996	New subroutines
+ * Apr 10 1996	Add subroutine list at start of file
+ * Apr 17 1996	Print error message to stderr
+ * May  2 1996	Write using stream IO
+ * May 14 1996	If FITSRTOPEN NK is zero, return all keywords in header
+ * May 17 1996	Make header internal to FITSRTOPEN
+ * Jun  3 1996	Use stream I/O for input as well as output
+ * Jun 10 1996	Remove unused variables after running lint
+ * Jun 12 1996	Deal with byte-swapped images
+ * Jul 11 1996	Rewrite code to separate header and data reading
+ * Aug  6 1996  Fixed small defects after lint
+ * Aug  6 1996  Drop unused NBHEAD argument from FITSRTHEAD
+ * Aug 13 1996	If filename is stdin, read from standard input instead of file
+ * Aug 30 1996	Use write for output, not fwrite
+ * Sep  4 1996	Fix mode when file is created
+ * Oct 15 1996	Drop column argument from FGET* subroutines
+ * Oct 15 1996	Drop unused variable 
+ * Dec 17 1996	Add option to skip bytes in file before reading the header
+ * Dec 27 1996	Turn nonprinting header characters into spaces
+ *
+ * Oct  9 1997	Add FITS extension support as filename,extension
+ * Dec 15 1997	Fix minor bugs after lint
+ *
+ * Feb 23 1998	Do not append primary header if getting header for ext. 0
+ * Feb 23 1998	Accept either bracketed or comma extension
+ * Feb 24 1998	Add SIMPLE keyword to start of extracted extension
+ * Apr 30 1998	Fix error return if not table file after Allan Brighton
+ * May  4 1998	Fix error in argument sequence in HGETS call
+ * May 27 1998	Include fitsio.h and imio.h
+ * Jun  1 1998	Add VMS fixes from Harry Payne at STScI
+ * Jun  3 1998	Fix bug reading EXTNAME
+ * Jun 11 1998	Initialize all header parameters before reading them
+ * Jul 13 1998	Clarify argument definitions
+ * Aug  6 1998	Rename fitsio.c to fitsfile.c to avoid conflict with CFITSIO
+ * Aug 13 1998	Add FITSWHEAD to write only header
+ * Sep 25 1998	Allow STDIN or stdin for standard input reading
+ * Oct  5 1998	Add isfits() to decide whether a file is FITS
+ * Oct  9 1998	Assume stdin and STDIN to be FITS files in isfits()
+ * Nov 30 1998	Fix bug found by Andreas Wicenec when reading large headers
+ * Dec  8 1998	Fix bug introduced by previous bug fix
+ *
+ * Jan  4 1999	Do not print error message if BITPIX is 0
+ * Jan 27 1999	Read and write all of 3D images if one dimension is 1
+ * Jan 27 1999	Pad out data to integral number of 2880-byte blocks
+ * Apr 29 1999	Write BITPIX=-16 files as BITPIX=16 with BSCALE and BZERO
+ * Apr 30 1999	Add % as alternative to , to denote sub-images
+ * May 25 1999	Set buffer offsets to 0 when FITS table file is opened
+ * Jul 14 1999	Do not try to write image data if BITPIX is 0
+ * Sep 27 1999	Add STDOUT as output filename option in fitswimage()
+ * Oct  6 1999	Set header length global variable hget.lhead0 in fitsrhead()
+ * Oct 14 1999	Update header length as it is changed in fitsrhead()
+ * Oct 20 1999	Change | in if statements to ||
+ * Oct 25 1999	Change most malloc() calls to calloc()
+ * Nov 24 1999	Add fitscimage()
+ *
+ * Feb 23 2000	Fix problem with some error returns in fitscimage()
+ * Mar 17 2000	Drop unused variables after lint
+ * Jul 20 2000	Drop BITPIX and NAXIS from primary header if extension printerd
+ * Jul 20 2000	Start primary part of header with ROOTHEAD keyword
+ * Jul 28 2000	Add loop to deal with buffered stdin
+ *
+ * Jan 11 2001	Print all messages to stderr
+ * Jan 12 2001	Add extension back onto filename after fitsropen() (Guy Rixon)
+ * Jan 18 2001	Drop EXTEND keyword when extracting an extension
+ * Jan 18 2001	Add fitswext() to append HDU and fitswhdu() to do actual writing
+ * Jan 22 2001	Ignore WCS name or letter following a : in file name in fitsrhead()
+ * Jan 30 2001	Fix FITSCIMAGE so it doesn't overwrite data when overwriting a file
+ * Feb 20 2001	Ignore WCS name or letter following a : in file name in fitsropen()
+ * Feb 23 2001	Initialize rbrac in fitsropen()
+ * Mar  8 2001	Use % instead of : for WCS specification in file name
+ * Mar  9 2001	Fix bug so primary header is always appended to secondary header
+ * Mar  9 2001	Change NEXTEND to NUMEXT in appended primary header
+ * Mar 20 2001	Declare fitsheadsize() in fitschead()
+ * Apr 24 2001	When matching column names, use longest length
+ * Jun 27 2001	In fitsrthead(), allocate pw and lpnam only if more space needed
+ * Aug 24 2001	In isfits(), return 0 if argument contains an equal sign
+ *
+ * Jan 28 2002	In fitsrhead(), allow stdin to include extension and/or WCS selection
+ * Jun 18 2002	Save error messages as fitserrmsg and use fitserr() to print them
+ * Oct 21 2002	Add fitsrsect() to read a section of an image
+ *
+ * Feb  4 2003	Open catalog file rb instead of r (Martin Ploner, Bern)
+ * Apr  2 2003	Drop unused variable in fitsrsect()
+ * Jul 11 2003	Use strcasecmp() to check for stdout and stdin
+ * Aug  1 2003	If no other header, return root header from fitsrhead()
+ * Aug 20 2003	Add fitsrfull() to read n-dimensional FITS images
+ * Aug 21 2003	Modify fitswimage() to always write n-dimensional FITS images
+ * Nov 18 2003	Fix minor bug in fitswhdu()
+ * Dec  3 2003	Remove unused variable lasthead in fitswhdu()
+ *
+ * May  3 2004	Do not always append primary header to extension header
+ * May  3 2004	Add ibhead as position of header read in file
+ * May 19 2004	Do not reset ext if NULL in fitswexhead()
+ * Jul  1 2004	Initialize INHERIT to 1
+ * Aug 30 2004	Move fitsheadsize() declaration to fitsfile.h
+ * Aug 31 2004	If SIMPLE=F, put whatever is in file after header in image
+ *
+ * Mar 17 2005	Use unbuffered I/O in isfits() for robustness
+ * Jun 27 2005	Drop unused variable nblocks in fitswexhead()
+ * Aug  8 2005	Fix space-padding bug in fitswexhead() found by Armin Rest
+ * Sep 30 2005	Fix fitsrsect() to position relatively, not absolutely
+ * Oct 28 2005	Add error message if desired FITS extension is not found
+ * Oct 28 2005	Fix initialization problem found by Sergey Koposov
+ *
+ * Feb 23 2006	Add fitsrtail() to read appended FITS headers
+ * Feb 27 2006	Add file name to header-reading error messages
+ * May  3 2006	Remove declarations of unused variables
+ * Jun 20 2006	Initialize uninitialized variables
+ * Nov  2 2006	Change all realloc() calls to calloc()
+ *
+ * Jan  5 2007	In fitsrtail(), change control characters in header to spaces
+ * Apr 30 2007	Improve error reporting in FITSRFULL
+ * Nov 28 2007	Add support to BINTABLE in ftget*() and fitsrthead()
+ * Dec 20 2007	Add data heap numerated by PCOUNT when skipping HDU in fitsrhead()
+ * Dec 20 2007	Return NULL pointer if fitsrhead() cannot find requested HDU
+ *
+ * Apr  7 2008	Drop comma from name when reading file in isfits()
+ * Jun 27 2008	Do not append primary data header if it is the only header
+ * Nov 21 2008	In fitswhead(), print message if too few bytes written
+ *
+ * Sep 18 2009	In fitswexhead() write to error string instead of stderr
+ * Sep 22 2009	In fitsrthead(), fix lengths for ASCII numeric table entries
+ * Sep 25 2009	Add subroutine moveb() and fix calls to it
+ * Sep 25 2009	Fix several small errors found by Jessicalas Burke
+ *
+ * Mar 29 2010	In fitswhead(), always pad blocks to 2880 bytes with spaces
+ * Mar 31 2010	In fitsrhead(), fix bug reading long primary headers
+ *
+ * Sep 15 2011	In fitsrsect() declare impos and nblin off_t
+ * Sep 15 2011	In fitsrtail() declare offset off_t
+ * Sep 15 2011	Declare global variable ibhead off_t
+ *
+ * Jul 25 2014	Fix bug when reallocating buffer for long headers
+ */