Initial commit

1 files changed, 923 insertions, 0 deletions

diff --git a/vendor/x11iraf/ximtool/fitsio.c b/vendor/x11iraf/ximtool/fitsio.c
new file mode 100644
index 00000000..9bbd50c1
--- /dev/null
+++ b/vendor/x11iraf/ximtool/fitsio.c
@@ -0,0 +1,923 @@
+/*
+ * FITSIO.C -- Routines to load and save simple FITS files.
+ *
+ *       ival = isFITS  (fname)
+ *            loadFITS  (fname, pixels, w, h, r, g, b, ncolors, 
+ *		  	    zsc, zr, &z1, &z2, nsample)
+ *           writeFITS  (fp, pixels, pixtype, w,h, r,g,b, ncolors)
+ *    str = getFITSHdr  (fname)
+
+ * LoadFITS(fname, numcols)  -  loads a FITS file
+ * WriteFITS(fp, pic, w, h, rmap, gmap, bmap, numcols)
+ *
+ * isFITS    -- returns nonzero if the named file is a FITS file.
+ * loadFITS  -- reads a FITS file and returns the decoded pixel array and gray-
+ * 	        scale 8 bit colormap.  The caller is responsible for freeing
+ * 	        the pixels buffer.  
+ * writeFITS -- performs the converse operation, writing the given pixel array
+ * 		and colormap to the output Sun rasterfile.
+ * Based on contributed FITS I/O software for XV by David Robinson.
+ */
+
+
+#include <stdio.h>
+#include <math.h>
+#include <ctype.h>
+
+#define NCARDS 		36
+#define BLOCKSIZE 	2880
+#define EPSILON 	1.192e-7
+
+
+#define CONTRAST        0.25            /* zscaling parameters */
+#define NSAMPLE         1000
+
+
+/* MONO returns total intensity of r,g,b components */
+#define MONO(rd,gn,bl) (((rd)*11 + (gn)*16 + (bl)*5) >> 5)  /*.33R+ .5G+ .17B*/
+
+/* data types */
+enum datatype { T_INT, T_LOG, T_REAL, T_NOVAL, T_STRING };
+
+typedef unsigned char byte;
+#ifndef AIXV3
+#ifndef OSF1
+typedef unsigned char uchar;
+#endif
+#endif
+
+typedef struct {
+	FILE      *fp;              /* file pointer 			*/
+	int	  bitpix, size;     /* bits per pixel, sizeof(unit) 	*/
+	int	  naxis;            /* number of axes 			*/
+	long int  axes[2];     	    /* size of each axis 		*/
+	long int  ndata;       	    /* number of elements in data       */
+	long int  cpos;        	    /* current position in data file    */
+	char      title[80];        /* image title 			*/
+	int	  extend;           /* image has extensions?		*/
+	int	  nextend;          /* number of extensions		*/
+	float     bscale, bzero;    /* scaling parameters 		*/
+} FITS;
+
+/* Function prototypes */
+#ifdef __STDC__
+
+#include <stddef.h>
+#include <stdlib.h>
+
+static char *ftopen2d (FITS *fs, char *file, int *nx, int *ny, int *bitpix);
+static void ftclose (FITS *fs);
+static char *ftgbyte (FITS *fs, uchar *buffer, int nelem, int zsc, int zr, float *z1, float *z2, int nsample);
+static char *rdheader (FITS *fs);
+static char *wrheader (FILE *fp, int nx, int ny);
+static char *rdcard (char *card, char *name, enum datatype dtype, long int *kvalue, float *rvalue);
+static void wrcard (char *card, char *name, enum datatype dtype, int kvalue);
+static char *ftgdata (FITS *fs, void *buffer, int nelem);
+static char *ftfixdata (FITS *fs, void *buffer, int nelem);
+#else
+static char *ftopen2d ();
+static void ftclose ();
+static char *ftgbyte ();
+static char *rdheader ();
+static char *wrheader ();
+static char *rdcard ();
+static void wrcard ();
+static char *ftgdata ();
+static char *ftfixdata ();
+#endif
+
+
+/* ----------------
+ * Public routines.
+ * ----------------*/
+
+
+/* loadFits - Load a simple FITS file.
+ */
+char *
+loadFITS (fname, pix, nx, ny, r,g,b, ncolors, zsc, zr, z1, z2, nsample)
+char 	*fname;                		/* input filename 	*/
+uchar   **pix;                 		/* output pixels 	*/
+int 	*nx, *ny;                       /* dimensions 		*/
+uchar 	*r, *g, *b;               	/* colormap 		*/
+int 	*ncolors;                       /* number of colors 	*/
+int	zsc, zr;			/* z-scaling flags	*/
+float	*z1, *z2;			/* zscale values	*/
+int	nsample;			/* number of sample pts */
+
+{
+	FITS 	fs;
+	int	i, w = 0, h = 0, bitpix, np;
+	byte 	*image;
+	char	*error;
+
+	error = ftopen2d (&fs, fname, &w, &h, &bitpix);
+	if (error)
+	    return error;
+	if (fs.extend)
+	    return "Load support for MEF files\nis not currently implemented";
+
+	/* allocate memory for image and read it in */
+	np = w * h;
+	image = (byte *) malloc (np);
+	if (image == NULL) 
+	    return "Insufficient memory";
+
+	error = ftgbyte (&fs, image, np, zsc, zr, z1, z2, nsample);
+	ftclose (&fs);
+	if (error) {
+	    free (image);
+	    return error;
+	}
+
+	/* There seems to be a convention that fits files be displayed using
+         * a cartesian coordinate system. Thus the first pixel is in the lower
+	 * left corner. Fix this by reflecting in the line y=h/2.
+ 	*/
+	flip (image, w, h);
+
+	/* sucess ! */
+	*pix = (unsigned char *) image;
+	*nx = w;
+	*ny = h;
+	*ncolors = 256;
+	for (i = 0; i < 256; i++) 
+	    r[i] = g[i] = b[i] = i;
+
+	return NULL;
+}
+
+
+/* writeFITS --  Write the current frame buffer out as a FITS image.
+ */
+
+char *
+writeFITS (fp, image, w, h, rmap, gmap, bmap, numcols)
+FILE 	*fp;
+byte 	*image;
+int	w, h;
+byte 	*rmap, *gmap, *bmap;
+int	numcols;
+{
+	register int	i, j, np, nend;
+	register byte *ptr;
+	char	*error;
+	byte rgb[256];
+
+	error = wrheader(fp, w, h);
+	if (error != NULL)
+	    return error;
+
+	for (i = 0; i < numcols; i++)
+	    rgb[i] = MONO(rmap[i], gmap[i], bmap[i]);
+
+	/* flip line ordering when writing out */
+	for (i = h - 1; i >= 0; i--) {
+	    ptr = &image[i*w];
+	    for (j = 0; j < w; j++, ptr++)
+	    	putc(rgb[*ptr], fp);
+	}
+
+	np = w * h;
+	/* nend is the number of padding characters at the end of the last
+         * block.
+	 */
+	nend = ((np + BLOCKSIZE - 1) / BLOCKSIZE) * BLOCKSIZE - np;
+
+	if (nend) 
+	    for (i = 0; i < nend; i++) 
+	    	putc('\0', fp);
+
+	return NULL;
+}
+
+
+/* IsFITS -- Test a file to see if it is a FITS file.
+ */
+int 
+isFITS (fname)
+char	*fname;				/* input filename */
+{
+	register FILE *fp;
+	int value = 0;
+	char keyw[8], val;
+
+	if (fp = fopen (fname, "r")) {
+	    fscanf (fp, "%6s = %c", keyw, &val);
+	    if (strcmp ("SIMPLE", keyw) == 0 && val == 'T')
+		value = 1;
+	    fclose (fp);
+	}
+	return value;
+}
+
+
+/* getFITSHdr -- Get some set of header information for the GUI.
+ */
+
+char *
+getFITSHdr (fname)
+char	*fname;
+{
+	FITS 	fs;
+	char 	*error, *title, *line;
+	int  	w, h, bitpix;
+
+	line = (char *) malloc (80);
+
+	error = ftopen2d (&fs, fname, &w, &h, &bitpix);
+	if (error) {
+	    strcpy (line, error);
+	    return (line);
+	}
+
+	if (fs.extend) {
+	    sprintf (line, "%-16.16s  %3d  (%2d extns)   %s", 
+	    	fname, bitpix, fs.nextend, fs.title);
+	} else {
+	    sprintf (line, "%-16.16s  %3d  %5dx%-5d  %s", 
+	    	fname, bitpix, w,h, fs.title);
+	}
+
+	ftclose (&fs);
+
+	return (line);
+}
+
+
+
+/*
+ * Private Procedures
+ * ------------------
+ */
+
+/* Writes a minimalist FITS file header */
+static char
+*wrheader (fp, nx, ny)
+FILE 	*fp;
+int	nx, ny;
+{
+	char	*block;
+	int	i;
+
+	block = (char *) malloc (BLOCKSIZE);
+	if (block == NULL) 
+	    return "Insufficient memory for workspace";
+	memset(block, ' ', BLOCKSIZE);
+
+	i = 0;
+	wrcard(&block[80*i++], "SIMPLE", T_LOG, 1);  	/* SIMPLE keyword */
+	wrcard(&block[80*i++], "BITPIX", T_INT, 8);  	/* BITPIX keyword */
+	wrcard(&block[80*i++], "NAXIS", T_INT, 2);   	/* NAXIS keyword  */
+	wrcard(&block[80*i++], "NAXIS1", T_INT, nx); 	/* NAXIS1 keyword */
+	wrcard(&block[80*i++], "NAXIS2", T_INT, ny); 	/* NAXIS2 keyword */
+	wrcard(&block[80*i++], "END", T_NOVAL, 0);   	/* END keyword    */
+	i = fwrite(block, sizeof(char), BLOCKSIZE, fp);
+	if (i != BLOCKSIZE) 
+	    return "Error writing FITS file";
+	return NULL;
+}
+
+
+/* open a 2-dimensional fits file.
+ * Stores the dimensions of the file in nx and ny, and updates the FITS
+ * structure passed in fs.
+ * If successful, returns NULL otherwise returns an error message.
+ * Will return an error message if the primary data unit is not a 2-dimensional
+ * array.
+ */
+static char	
+*ftopen2d(fs, file, nx, ny, bitpix)
+FITS 	*fs;
+char	*file;
+int	*nx, *ny, *bitpix;
+{
+	FILE 	*fp;
+	int	naxis, i;
+	char	*error;
+
+	fp = fopen(file, "rb");
+	if (fp == NULL) 
+	    return "Unable to open FITS file";
+
+	fs->fp = fp;
+	fs->bitpix = 0;
+	fs->naxis = 0;
+	fs->cpos = 0;
+
+	/* read header */
+	error = rdheader(fs);
+
+	if (error != NULL) {
+	    ftclose(fs);
+	    return error;
+	}
+
+	/* get number of data */
+	fs->ndata = 1;
+	for (i = 0; i < fs->naxis; i++) 
+	    fs->ndata = fs->ndata * fs->axes[i];
+
+	naxis = fs->naxis;
+
+	*nx = fs->axes[0];
+	*ny = fs->axes[1];
+	*bitpix = fs->bitpix;
+
+	return error;
+}
+
+
+/* closes a fits file */
+static void
+ftclose (fs)
+FITS 	*fs;
+{
+	if (fs == NULL) 
+	    return;
+	if (fs->fp != NULL) 
+	    fclose(fs->fp);
+}
+
+
+/* reads the fits header, and updates the FITS structure fs.
+ * Returns NULL on success, or an error message otherwise.
+ */
+static char *
+rdheader (fs)
+FITS 	*fs;
+{
+	int	i, j, res;
+	char	name[9];
+	char	*block;
+	char	*error;
+	long int	val;         /* the value */
+	float	rval;		     /* floating point value */
+
+	block = (char *) malloc(BLOCKSIZE);
+	if (block == NULL) 
+	    return "Insufficient memory for workspace";
+
+	res = fread(block, sizeof(char), BLOCKSIZE, fs->fp);
+	if (res != BLOCKSIZE) 
+	    return "Error reading FITS file";
+	i = 0;
+
+	/* read SIMPLE key */
+	error = rdcard(block, "SIMPLE", T_LOG, &val, &rval);
+	if (error != NULL) 
+	    return error;
+	/*
+	if (val == 0) 
+	    return "Not a SIMPLE FITS file";
+	*/
+	i++;
+
+	/* read BITPIX key */
+	error = rdcard(&block[80], "BITPIX", T_INT, &val, &rval);
+	if (error != NULL) 
+	    return error;
+	if (val != 8 && val != 16 && val != 32 && val != 64 && val != -32 && 
+	    val != -64)
+	    return "Bad BITPIX value in FITS file";
+	fs->bitpix = val;
+	j = fs->bitpix;
+	if (j < 0) 
+	    j = -j;
+	fs->size = j / 8;
+	i++;
+
+	/* read NAXIS key */
+	error = rdcard(&block[2*80], "NAXIS", T_INT, &val, &rval);
+	if (error != NULL) 
+	    return error;
+	if (val < 0 || val > 999)
+	    return "Bad NAXIS value in FITS file";
+	if (val == 1 || val > 2)
+	    return "FITS file is not a two-dimensional image";
+	fs->naxis = val;
+	i++;
+
+	/* Check for an EXTEND/NEXTEND keyword pair. */
+	error = rdcard(&block[3*80], "EXTEND", T_LOG, &val, &rval);
+	if (error == NULL) {
+	    fs->extend = val; i++;
+	} else 
+	    fs->extend = 0;
+	error = rdcard(&block[4*80], "NEXTEND", T_INT, &val, &rval);
+	if (error == NULL) {
+	    fs->nextend = val; i++;
+	} else 
+	    fs->nextend = 0;
+
+
+	/* read NAXISnnn keys.
+ 	 * We allow NAXIS to be > 2 iff the dimensions of the extra axes are 1
+ 	 */
+	for (j = 0; j < fs->naxis; j++) {
+	    if (i == NCARDS) {
+	    	res = fread(block, sizeof(char), BLOCKSIZE, fs->fp);
+	    	if (res != BLOCKSIZE) 
+	    	    return "Error reading FITS file";
+	    	i = 0;
+	    }
+
+	    sprintf(name, "NAXIS%d", j + 1);
+	    error = rdcard(&block[i*80], name, T_INT, &val, &rval);
+	    if (error != NULL) 
+	    	return error;
+	    if (val < 0) 
+	    	return "Bad NAXISn value in FITS file";
+	    if (j < 2) 
+	    	fs->axes[j] = val;
+	    else if (val != 1) 
+	    	return "FITS file is not a two-dimensional image";
+	    i++;
+	}
+	fs->naxis = 2;
+
+	/* do remainder */
+	fs->bscale  = 1.0;
+	fs->bzero   = 0.0;
+	memset(fs->title, '\0', 80);
+	strcpy (fs->title, "No Title");
+	for (; ; ) {
+
+	    if (block[i*80] == 'B') {
+		/* Try reading a BSCALE or BZERO keyword from this card. */
+	        error = rdcard(&block[i*80], "BSCALE", T_REAL, &val, &rval);
+	        if (error == NULL)
+	            fs->bscale = rval;
+	        error = rdcard(&block[i*80], "BZERO", T_REAL, &val, &rval);
+	        if (error == NULL)
+	            fs->bzero = rval;
+	    }
+
+	    if (i == NCARDS) {
+	    	res = fread(block, sizeof(char), BLOCKSIZE, fs->fp);
+	    	if (res != BLOCKSIZE) 
+	    	    return "Unexpected eof in FITS file";
+	    	i = 0;
+	    }
+
+	    if (strncmp(&block[i*80], "OBJECT  ", 8) == 0) {
+		char *ip = &block[i*80];
+		int i = 0, j = 0;
+
+		/* Skip ahead to opening quote. */
+		while (*ip != '\'' && i < 80) ip++, i++; 
+		ip++, i++;
+		for (j=0; j < 80 && *ip != '\''; j++) {
+		    fs->title[j] = *ip;
+		    ip++;
+		}
+		fs->title[j] = '\0';
+	    }
+
+	    if (strncmp(&block[i*80], "END     ", 8) == 0) 
+	    	break;
+	    i++;
+	}
+	free(block);
+	return NULL;
+}
+
+
+/* write a header record into the 80 byte buffer card.
+ * The keyword name is passed in name. The value type is in dtype; this
+ * can have the following values:
+ *    dtype = T_NOVAL
+ *         no keyword value is written
+ *    dtype = T_LOG
+ *         a logical value, either 'T' or 'F' in column 30 is written
+ *    dtype = T_INT
+ *         an integer is written, right justified in columns 11-30
+ */
+static void
+wrcard (card, name, dtype, kvalue)
+char	*card, *name;
+enum 	datatype dtype;   /* type of value */
+int	kvalue;
+{
+	int	l;
+	memset(card, ' ', 80);
+	l = strlen(name);
+	if (l) 
+	    memcpy(card, name, l);   /* copy name */
+
+	if (dtype == T_NOVAL) 
+	    return;
+
+	card[8] = '=';
+
+	if (dtype == T_LOG)
+	    card[29] = kvalue ? 'T' : 'F';
+	else /* an integer */    {
+	    sprintf(&card[10], "%20d", kvalue);
+	    card[30] = ' ';
+	}
+
+}
+
+
+/* Read a header record, from the 80 byte buffer card.
+ * the keyword name must match 'name'; and parse its value according to
+ * dtype. This can have the following values:
+ *    dtype = T_LOG
+ *        value is logical, either 'T' or 'F' in column 30.
+ *    dtype = T_INT
+ *        value is an integer, right justified in columns 11-30.
+ *    dtype = T_REAL
+ *        value is a real
+ *
+ * The value is stored in kvalue.
+ * It returns NULL on success, or an error message otherwise.
+ */
+static char *
+rdcard (card, name, dtype, kvalue, rvalue)
+char	*card, *name;
+enum 	datatype dtype;   			/* type of value */
+long 	int *kvalue;
+float 	*rvalue;
+{
+	int	i, ptr;
+	char	namestr[9];
+	static char	error[45];
+
+	memcpy(namestr, card, 8);
+	i = 8;
+	do 
+	    i--; 
+	while (i >= 0 && namestr[i] == ' ');
+	namestr[i+1] = '\0';
+
+	if (strcmp(namestr, name) != 0) {
+	    sprintf(error, "Keyword %s not found in FITS file", name);
+	    return error;
+	}
+
+	/* get start of value */
+	ptr = 10;
+	while (ptr < 80 && card[ptr] == ' ') 
+	    ptr++;
+	if (ptr == 80 || card[ptr] == '/') 
+	    return "FITS file has missing keyword value"; /* no value */
+
+	if (dtype == T_LOG) {
+	    if ((card[ptr] != 'T' && card[ptr] != 'F'))
+	    	return "Keyword has bad logical value in FITS file";
+	    *kvalue = (card[ptr] == 'T');
+
+	} else /* an integer or real */    {
+	    int		j, end;
+	    long int	ival;
+	    float	fval;
+	    char	num[21];
+
+	    /*
+	    if (ptr > 29) 
+	    	return "Keyword has bad integer value in FITS file";
+	    */
+	    end = ptr;
+	    while (end < 80 && (card[end] != ' ' && card[end] != '/'))
+	        end++;
+	    memcpy(num, &card[ptr], end - ptr);
+	    num[end-ptr] = '\0';
+	    if (dtype == T_INT) {
+	       j = sscanf(num, "%ld", &ival);
+	       if (j != 1) 
+	    	   return "Keyword has bad integer value in FITS file";
+	       *kvalue = ival;
+	       *rvalue = 0.0;
+	    } else if (dtype == T_REAL) {
+	       j = sscanf(num, "%g", &fval);
+	       if (j != 1) 
+	    	   return "Keyword has bad real value in FITS file";
+	       *kvalue = 0;
+	       *rvalue = fval;
+	    }
+	}
+
+	return NULL;
+}
+
+
+/* reads nelem values into the buffer.
+ * returns NULL for success or an error message.
+ * Copes with the fact that the last 2880 byte record of the FITS file
+ * may be truncated, and should be padded out with zeros.
+ */
+static char *
+ftgdata (fs, buffer, nelem)
+FITS 	*fs;
+void 	*buffer;
+int	nelem;
+{
+	int	res;
+
+	if (nelem == 0) 
+	    return NULL;
+
+	res = fread(buffer, fs->size, nelem, fs->fp);
+	/* if failed to read all the data */
+	if (res != nelem) {
+	    /* nblock is the number of elements in a record. size is
+	     * always a factor of BLOCKSIZE */
+	    int	loffs, nblock = BLOCKSIZE / fs->size;
+
+	    if (!feof(fs->fp)) 
+	    	return "I/O error reading FITS file";
+
+	    /* The last record might be short; check this.
+ 	     * loffs is the offset of the start of the last record from
+	     * the current position.
+ 	     */
+	    loffs = ((fs->ndata + nblock - 1) / nblock - 1) * nblock - fs->cpos;
+
+	    /* if we didn't read to the end of the penultimate record */
+	    if (res < loffs) 
+	    	return "Unexpected EOF reading FITS file";
+
+	    /* pad with zeros */
+	    memset((char *)buffer + res * fs->size, '\0',
+		(nelem - res) * fs->size);
+	}
+	fs->cpos += res;
+	return ftfixdata(fs, buffer, nelem);
+}
+
+
+/* convert the raw data, as stored in the FITS file, to the format
+ * appropiate for the data representation of the host computer.
+ * Assumes that
+ *  short int = 2 byte integer
+ *  int       = 4 byte integer
+ */
+static char *
+ftfixdata (fs, buffer, nelem)
+FITS 	*fs;
+void 	*buffer;
+int	nelem;
+{
+	register int	i, n = nelem;
+	register uchar	*ptr = buffer;
+
+	/* conversions. Although the data may be signed, reverse using unsigned 
+         * variables.  Convert from big-endian two-byte signed integer to
+	 * native form 
+	 */
+	if (fs->bitpix == 16)
+	    for (i = 0; i < n; i++, ptr += 2)
+	    	*(unsigned short int *)ptr = (((int)*ptr) << 8) | (int)(ptr[1]);
+
+	/* convert from big-endian four-byte signed integer to native form */
+	else if (fs->bitpix == 32)
+	    for (i = 0; i < n; i++, ptr += 4)
+	    	*(unsigned int *)ptr = (((unsigned int)*ptr) << 24) | 
+	    	    ((unsigned int)ptr[1] << 16) | 
+	    	    ((unsigned int)ptr[2] << 8)  | 
+	    	    (unsigned int)ptr[3];
+
+	/* convert from IEE 754 single precision to native form */
+	else if (fs->bitpix == -32) {
+	    register int	j, k, expo;
+	    static float	*exps = NULL;
+
+	    if (exps == NULL) {
+	    	exps = (float *)calloc(256, sizeof(float));
+	    	if (exps == NULL) 
+	    	    return "Insufficient memory for workspace";
+	    	exps[150] = 1.;
+	    	for (i = 151; i < 256; i++) 
+	    	    exps[i] = 2. * exps[i-1];
+	    	for (i = 149; i >= 0; i--) 
+	    	    exps[i] = 0.5 * exps[i+1];
+	    }
+
+
+	    for (i = 0; i < n; i++, ptr += 4) {
+	    	k = (int)*ptr;
+	    	j = ((int)ptr[1] << 16) | ((int)ptr[2] << 8) | (int)ptr[3];
+	    	expo = ((k & 127) << 1) | (j >> 23);
+	    	if ((expo | j) == 0) 
+	    	    *(float *)ptr = 0.;
+	    	else
+	    	    *(float *)ptr = exps[expo] * (float)(j | 0x800000);
+	    	if (k & 128) 
+	    	    *(float *)ptr = -*(float *)ptr;
+	    }
+
+	/* convert from IEE 754 double precision to native form */
+	} else if (fs->bitpix == -64) {
+	    register int	expo, k, l;
+	    register unsigned int	j;
+	    static double	*exps = NULL;
+
+	    if (exps == NULL) {
+	    	exps = (double *)calloc(2048, sizeof(double));
+	    	if (exps == NULL) 
+	    	    return "Insufficient memory for workspace";
+	    	exps[1075] = 1.;
+	    	for (i = 1076; i < 2048; i++) 
+	    	    exps[i] = 2. * exps[i-1];
+	    	for (i = 1074; i >= 0; i--) 
+	    	    exps[i] = 0.5 * exps[i+1];
+	    }
+
+	    for (i = 0; i < n; i++, ptr += 8) {
+	    	k = (int)*ptr;
+	    	j = ((unsigned int)ptr[1] << 24) | ((unsigned int)ptr[2] << 16) | 
+	    	    ((unsigned int)ptr[3] << 8) | (unsigned int)ptr[4];
+	    	l = ((int)ptr[5] << 16) | ((int)ptr[6] << 8) | (int)ptr[7];
+	    	expo = ((k & 127) << 4) | (j >> 28);
+	    	if ((expo | j | l) == 0) 
+	    	    *(double *)ptr = 0.;
+	    	else
+	    	    *(double *)ptr = exps[expo] * (16777216. * 
+	    	        (double)((j & 0x0FFFFFFF) | 0x10000000) + (double)l);
+	    	if (k & 128) 
+	    	    *(double *)ptr = -*(double *)ptr;
+	    }
+	}
+	return NULL;
+}
+
+
+#undef max
+#define f_max(a,b) ((a) > (b) ? (a) : (b))
+#undef min
+#define f_min(a,b) ((a) < (b) ? (a) : (b)) 
+
+
+/* Reads a byte image from the FITS file fs. The image contains nelem pixels.
+ * If bitpix = 8, then the image is loaded as stored in the file if not scaled.
+ * Otherwise, it is rescaled so that the minimum value is stored as 0, and
+ * the maximum is stored as 255
+ */
+static char *
+ftgbyte(fs, cbuff, nelem, zsc, zr, z1, z2, nsample)
+FITS 	*fs;
+uchar	*cbuff;
+int	nelem;
+int	zsc, zr;
+float 	*z1, *z2;
+int	nsample;
+{
+	char * voidbuff;
+	register int	i, n = nelem;
+	char	*error;
+	int	pmin = 0, pmax = 255;
+	int	npts, stdline;
+	extern	void zscale();
+
+	/* if the data is uchar, then read it directly */
+	if (fs->bitpix == 8 && (fs->bscale == 1.0 && fs->bzero == 0.0)) {
+	    *z1 = 0.0;
+	    *z2 = 255.0;
+	    return ftgdata(fs, cbuff, nelem);
+	}
+
+	/* allocate a buffer to store the image */
+	if (fabs((double)fs->bscale-1.0) > EPSILON || 
+	    fabs((double)fs->bzero) > EPSILON)
+	        voidbuff = (char * )malloc(nelem * f_max(fs->size,4));
+	else
+	    voidbuff = (char * )malloc(nelem * fs->size);
+	if (voidbuff == NULL) 
+	    return "Insufficient memory for workspace";
+	error = ftgdata(fs, voidbuff, nelem);
+	if (error != NULL) {
+	    printf ("ftgbyte: %s\n", error);
+	    return error;
+	}
+
+	/* If we've got BSCALE/BZERO values compute the original pixel values
+	 * and convert the buffer to floating point before processing it.
+	 * The voidbuff was allocated above with this in mind so it should be
+	 * large enough that we can fix the pixels in place.
+	 */
+	            
+	npts = fs->axes[0] * fs->axes[1];
+	stdline = (int)((float)fs->axes[1] / sqrt((float)npts/(float)nsample));
+	if (fs->bscale != 1.0 || fs->bzero != 0.0) {
+	    register float *buf;
+
+	    buf = (float *)voidbuff;
+
+	    if (fs->bitpix ==  8) {
+	        for (i=(nelem-1); i >= 0; i--)
+		    buf[i] = (float) voidbuff[i] * fs->bscale + fs->bzero;
+	    } else if (fs->bitpix ==  16) {
+	        register short *old;
+	        for (i=(nelem-1); i >= 0; i--) {
+		    old = (short *) &voidbuff[i * 2];
+		    buf[i] = (float) *old * fs->bscale + fs->bzero;
+	        }
+	    } else if (fs->bitpix ==  32) {
+	        register int *old;
+	        for (i=(nelem-1); i >= 0; i--) {
+		    old = (int *) &voidbuff[i * 4];
+		    buf[i] = (float) *old * fs->bscale + fs->bzero;
+	        }
+	    } else if (fs->bitpix ==  -32) {
+	        register float *old;
+	        for (i=(nelem-1); i >= 0; i--) {
+		    old = (float *) &voidbuff[i * 4];
+		    buf[i] = (float) *old * fs->bscale + fs->bzero;
+	        }
+	    } else if (fs->bitpix ==  -64) {
+	        register double *old, *dbuf;
+		register float *fpix;
+
+	 	dbuf = (double *) malloc (nelem * sizeof(double));
+	        for (i=(nelem-1); i >= 0; i--) {
+		    old = (double *) &voidbuff[i * 8];
+		    dbuf[i] = (float) *old * fs->bscale + fs->bzero;
+	        }
+		fpix = (float *) voidbuff;
+		for (i=0; i<nelem; i++)
+		    buf[i] = fpix[i] = (float) dbuf[i];
+		free ((double *) dbuf);
+	    }
+	
+	    fs->size   = 4;
+	    fs->bitpix = -32;
+
+	    /* compute the optimal zscale values */
+	    if (zsc)
+	        zscale ((char *)buf, fs->axes[0], fs->axes[1], fs->bitpix,
+	            z1, z2, CONTRAST, nsample, stdline); 
+	    else if (zr)
+		min_max ((float *)buf, nelem, fs->bitpix, z1, z2);
+
+	} else {
+	    /* compute the optimal zscale values */
+	    if (zsc)
+	        zscale (voidbuff, fs->axes[0], fs->axes[1], fs->bitpix,
+	            z1, z2, CONTRAST, nsample, stdline); 
+	    else if (zr)
+		min_max ((float *)voidbuff, nelem, fs->bitpix, z1, z2);
+	}
+
+	/* convert short int to uchar */
+	if (fs->bitpix == 16) {
+	    register short int	*buffer = (short *)voidbuff;
+	    register int	max, min;
+	    register float	scale;
+
+	    min = (int) *z1;
+	    max = (int) *z2;
+	    scale = (max == min) ? 0. : 255. / (*z2 -*z1);
+
+	    /* rescale and convert */
+	    for (i = 0, buffer = (short *)voidbuff; i < n; i++)
+	    	cbuff[i] = f_max (pmin, f_min (pmax, 
+		    (int)(scale * (float)((int)buffer[i] - min)) ));
+
+	/* convert long int to uchar */
+	} else if (fs->bitpix == 32) {
+	    register int	*buffer = (int *)voidbuff;
+	    register int	max, min;
+	    register float	scale;
+
+	    min = (int) *z1;
+	    max = (int) *z2;
+	    scale = (max == min) ? 0. : 255. / (*z2 - *z1);
+
+	    /* rescale and convert */
+	    for (i = 0, buffer = (int *)voidbuff; i < n; i++)
+	    	cbuff[i] = f_max (pmin, f_min (pmax, 
+		    (int)(scale * (float)((int)buffer[i] - min)) ));
+
+	/* convert float to uchar */
+	} else if (fs->bitpix == -32) {
+	    register float	*buffer = (float *)voidbuff;
+	    register float	max, min, scale;
+
+	    min = *z1;
+	    max = *z2;
+	    scale = (max == min) ? 0. : 255. / (*z2 - *z1);
+
+	    /* rescale and convert */
+	    for (i = 0, buffer = (float *)voidbuff; i < n; i++)
+	    	cbuff[i] = f_max (pmin, f_min (pmax, 
+		    (int)(scale * ((float)buffer[i] - min)) ));
+
+	/* convert double to uchar */
+	} else if (fs->bitpix == -64) {
+	    register double	*buffer = (double *)voidbuff;
+	    register double	max, min, scale;
+
+	    min = (double) *z1;
+	    max = (double) *z2;
+	    scale = (max == min) ? 0. : 255. / (*z2 - *z1);
+
+	    /* rescale and convert */
+	    for (i = 0, buffer = (double *)voidbuff; i < n; i++)
+	    	cbuff[i] = f_max (pmin, f_min (pmax, 
+		    (int)(scale * ((double)buffer[i] - min)) ));
+	}
+
+	free( (char *)voidbuff);
+	return NULL;
+}