path: root/vendor/x11iraf/cdl/cdlfits.c

#include <stdio.h>
#include <math.h>
#include <ctype.h>
#define  CDL_LIBRARY_SOURCE
#include "cdl.h"


/*
 * CDLFITS.C -- Routines to load simple FITS files.
 *
 *          cdl_displayFITS  (cdl, fname, frame, fbconfig, zscale)
 *               cdl_isFITS  (fname)
 *             cdl_readFITS  (fname, &pix, &nx, &ny, &bitpix, title);
 */


#define NCARDS 		36
#define BLOCKSIZE 	2880
#define SZ_CARD 	80

/* data types */
enum datatype { T_INT, T_LOG, T_REAL, T_STR, T_NOVAL };

typedef struct {
	FILE      *fp;              /* file pointer                           */
	int	  bitpix, size;     /* number of 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          */
	float     bscale, bzero;    /* scaling parameters                     */
	char	  object[SZ_CARD];  /* object keyword name                    */
} FITS;

/* Function prototypes */
#ifdef __STDC__
#include <stddef.h>
#include <stdlib.h>
#endif

#undef max
#define f_max(a,b) ((a) > (b) ? (a) : (b))
#undef min
#define f_min(a,b) ((a) < (b) ? (a) : (b)) 

extern	int	cdl_debug;


#ifdef ANSI_FUNC

static int cdl_readFITSHdr(char *fname, float *bscale, float *bzero, char *obj);
static char *cdl_getFITSPixels(char *fname, uchar **pix, int *nx, int *ny, int *bitpix);
static char *cdl_openFITS(FITS *fs, char *file, int *nx, int *ny, int *bitpix);
static char *cdl_readHeader(FITS *fs);
static char *cdl_rdCard(char *card, char *name, enum datatype dtype, long *kvalue, float *rvalue, char *svalue);
static char *cdl_getData(FITS *fs, uchar *buffer, int nelem);
static char *cdl_fixData(FITS *fs, void *buffer, int nelem);

#else

static char *cdl_getFITSPixels();
static char *cdl_openFITS(), *cdl_readHeader ();
static char *cdl_rdCard(), *cdl_getData (), *cdl_fixData();
static int  cdl_readFITSHdr();

#endif


/*  CDL_DISPLAYFITS -- Display a simple FITS format image to the given frame, 
 *  optionally doing a Zscale transform of the image pixels.
 */

#ifdef ANSI_FUNC

int 
cdl_displayFITS (
    CDLPtr cdl,                         /* package ptr          */
    char *fname,			/* image name		*/
    int frame,				/* display frame	*/
    int fbconfig,			/* frame buffer config	*/
    int zscale				/* do zscale of image?	*/
)
#else

int
cdl_displayFITS (cdl, fname, frame, fbconfig, zscale)
CDLPtr  cdl;                            /* package ptr          */
char	*fname;				/* image name		*/
int	frame;				/* display frame	*/
int	fbconfig;			/* frame buffer config	*/
int	zscale;				/* do zscale of image?	*/
#endif
{
	char	obj[SZ_CARD];
	int	status, nx, ny, bitpix;
	uchar	*pix;
	float	z1 = 0.0, z2 = 0.0, bscale, bzero;

	/* See if this is a valid FITS file. */
	if (!cdl_isFITS (fname)) {
	    fprintf (stderr, "%s: not a simple FITS file or doesn't exist.\n",
		fname);
	    return (ERR);
	}

	/* Get the raw FITS image pixels. */
        if (cdl_readFITS (fname, &pix, &nx, &ny, &bitpix, obj))
	    return (ERR);

	if (cdl_debug) {
	    printf("[cdl_displayFITS] '%s' frame=%d zscale=%d\n",
		fname, frame, zscale);
	    printf("[cdl_displayFITS] %dx%d bitpix=%d pixels z1=%g z2=%g\n",
		nx, ny, bitpix, z1, z2);
	}

	/* Read the header. */
	if (cdl_readFITSHdr (fname, &bscale, &bzero, obj))
	    return (ERR);
	(void) cdl_setTitle (cdl, obj);

	(void) cdl_setName (cdl, fname);
	status =  cdl_displayPix (cdl, pix, nx, ny, bitpix, frame, fbconfig,
	    zscale);
	free ((char *) pix);
	return (status);
}


/*  CDL_ISFITS -- Test a file to see if it is a simple FITS file.
 */

#ifdef ANSI_FUNC

int 
cdl_isFITS (
    char *fname				/* input filename */
)
#else

int 
cdl_isFITS (fname)
char	*fname;				/* input filename */
#endif
{
	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;
}


/*  CDL_READFITS -- Read the pixels from a simple FITS format image, ret-
 *  urning the pixel array, dimensions, and image type.
 */

#ifdef ANSI_FUNC

int 
cdl_readFITS (
    char *fname,			/* image name		*/
    uchar **pix,			/* pixel array (output) */
    int *nx,
    int *ny,				/* dimensions (output)	*/
    int *bitpix,			/* pixel size (output)	*/
    char *title				/* image title (output) */
)
#else

int
cdl_readFITS (fname, pix, nx, ny, bitpix, title)
char	*fname;				/* image name		*/
uchar	**pix;				/* pixel array (output) */
int	*nx, *ny;			/* dimensions (output)	*/
int	*bitpix;			/* pixel size (output)	*/
char	*title;				/* image title (output) */
#endif
{
	char	*errstr, *cdl_getFITSPixels();
	float	bscale=1.0, bzero=0.0;

	/* See if this is a valid FITS file. */
	if (!cdl_isFITS (fname)) {
	    fprintf (stderr, "%s: not a simple FITS file or doesn't exist.\n",
		fname);
	    return (ERR);
	}

	if ((errstr = cdl_getFITSPixels (fname, pix, nx, ny, bitpix))) {
	    fprintf (stderr, "%s\n", errstr);
	    return (ERR);
	}
	if (cdl_readFITSHdr (fname, &bscale, &bzero, title))
	    return (ERR);

	if (cdl_debug)
	    printf ("[cdl_readFITS] '%s' nx=%d ny=%d bitpix=%d title='%s'\n",
		fname, *nx, *ny, *bitpix, title);

	return (OK);
}


/* ------------------
 * Private Procedures
 * ------------------*/


/*  CDL_READFITSHDR -- Read several header parameters from the image.
 */

#ifdef ANSI_FUNC

static int 
cdl_readFITSHdr (
    char *fname,			/* image name		*/
    float *bscale,
    float *bzero,			/* scaling (output)	*/
    char *obj				/* image title (output)	*/
)
#else

static int
cdl_readFITSHdr (fname, bscale, bzero, obj)
char	*fname;				/* image name		*/
float	*bscale, *bzero;		/* scaling (output)	*/
char	*obj;				/* image title (output)	*/
#endif
{
	char	*errstr;
	register FITS	*fs;
	register FILE	*fp;

	/* See if this is a valid FITS file. */
	if (!cdl_isFITS (fname)) {
	    fprintf (stderr, "%s: not a simple FITS file or doesn't exist.\n",
		fname);
	    return (ERR);
	}

	if ((fp = fopen(fname, "rb")) == NULL) 
	    return (ERR);

	fs = (FITS *) malloc (sizeof(FITS));
	fs->fp = fp;
	fs->bitpix = fs->naxis = fs->cpos = 0;

	/* Read the header. */
	if ((errstr = cdl_readHeader (fs))) {
	    fclose(fs->fp);
	    free ((char *) fs);
	    return (ERR);
	}

	*bscale = fs->bscale;
	*bzero = fs->bzero;
	if (obj == (char *)NULL)
	    obj = (char *) malloc (SZ_CARD+1);
	sprintf (obj, "%s", fs->object);

	if (cdl_debug)
	    printf ("[cdl_readFITSHdr] '%s' bscale=%g bzero=%g obj='%s'\n",
		fname, *bscale, *bzero, obj);

	fclose(fs->fp);
	free ((char *) fs);
	return (OK);
}


/*  CDL_GETFITSPIXELS -- Given an FITS filename return a pointer to the pixel
 *  array and the image dimensions.
 */

#ifdef ANSI_FUNC

static char *
cdl_getFITSPixels (
    char *fname,               		/* input filename 	*/
    uchar **pix,               		/* output pixels 	*/
    int *nx,
    int *ny,                       	/* dimensions 		*/
    int *bitpix				/* pixel size		*/
)
#else

static char *
cdl_getFITSPixels (fname, pix, nx, ny, bitpix)
char 	*fname;                		/* input filename 	*/
uchar   **pix;                 		/* output pixels 	*/
int 	*nx, *ny;                       /* dimensions 		*/
int	*bitpix;			/* pixel size		*/
#endif
{
	FITS 	fs;
	register int	i, nelem;
	char	*error = NULL;


	if ((error = cdl_openFITS (&fs, fname, nx, ny, bitpix)))
	    return error;

        /* Allocate a buffer to store the image. */
	nelem = (*nx) * (*ny);
        if (fs.bscale != 1.0 || fs.bzero != 0.0)
            *pix = (uchar *) malloc (nelem * f_max(fs.size,4));
        else
            *pix = (uchar *) malloc (nelem * fs.size);
        if (*pix == NULL) 
            return "Insufficient memory for workspace";

        /* If the data is uchar, then read it directly */
        if (fs.bitpix == 8 && (fs.bscale == 1.0 || fs.bzero == 0.0))
            return cdl_getData (&fs, *pix, nelem);
	else if ((error = cdl_getData (&fs, *pix, nelem)))
            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 pix buffer was allocated above with this in mind so it should
         * be large enough that we can fix the pixels in place.
         */
        if (fs.bscale != 1.0 || fs.bzero != 0.0) {
            register float *buf;

            buf = (float *)*pix;

            if (fs.bitpix ==  8) {
                for (i=(nelem-1); i >= 0; i--)
                    buf[i] = (float) (*pix)[i] * fs.bscale + fs.bzero;
            } else if (fs.bitpix ==  16) {
                register short *old = (short *) *pix;
                for (i=(nelem-1); i >= 0; i--)
                    buf[i] = (float) old[i] * fs.bscale + fs.bzero;
            } else if (fs.bitpix ==  32) {
                register int *old = (int *) *pix;
                for (i=(nelem-1); i >= 0; i--)
                    buf[i] = (float) old[i] * fs.bscale + fs.bzero;
            } else if (fs.bitpix ==  -32) {
                register float *old = (float *) *pix;
                for (i=(nelem-1); i >= 0; i--)
                    buf[i] = (float) old[i] * fs.bscale + fs.bzero;
            } else if (fs.bitpix ==  -64) {
                register double *old = (double *) *pix, *dbuf;
                register float *fpix;

                dbuf = (double *) malloc (nelem * sizeof(double));
                for (i=(nelem-1); i >= 0; i--)
                    dbuf[i] = (float) old[i] * fs.bscale + fs.bzero;
                fpix = (float *) *pix;
                for (i=0; i<nelem; i++)
                    fpix[i] = (float) dbuf[i];
                free ((double *) dbuf);
            }
        
            fs.size   = 4;
            fs.bitpix = -32;
            *bitpix = -32;
        }

	/* sucess ! */
	fclose(fs.fp);
	return NULL;
}


/*  CDL_OPENFITS -- 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.
 */

#ifdef ANSI_FUNC

static char *
cdl_openFITS (FITS *fs, char *file, int *nx, int *ny, int *bitpix)
#else

static char *
cdl_openFITS(fs, file, nx, ny, bitpix)
FITS 	*fs;
char	*file;
int	*nx, *ny, *bitpix;
#endif
{
	register int	i;
	FILE 	*fp;
	char	*error = NULL;

	if ((fp = fopen(file, "rb")) == NULL) 
	    return "Unable to open FITS file";

	fs->fp = fp;
	fs->bitpix = 0;
	fs->naxis = 0;
	fs->cpos = 0;

	/* Read the header. */
	if ((error = cdl_readHeader (fs))) {
	    fclose(fs->fp);
	    return error;
	}

	if (fs->naxis > 2) {
	    fclose(fs->fp);
	    return "Not a 2-D FITS image.";
	}

	/* Get number of pixel. */
	fs->ndata = 1;
	for (i = 0; i < fs->naxis; i++) 
	    fs->ndata = fs->ndata * fs->axes[i];

	*nx = fs->axes[0];
	*ny = fs->axes[1];
	*bitpix = fs->bitpix;

	return NULL;
}


/*  CDL_READHEADER -- Reads the fits header, and updates the FITS structure fs.
 *  Returns NULL on success, or an error message otherwise.
 */

#ifdef ANSI_FUNC

static char *
cdl_readHeader (FITS *fs)
#else

static char *
cdl_readHeader (fs)
FITS *fs;
#endif
{
	register int	 i, j, res;
	char	 name[9];
	char	 *block;
	char	 *error = NULL;
	long int val;         	     	/* the value 		*/
	float	 rval;		     	/* floating point value */
	char	 sval[SZ_FNAME];	/* string value		*/


	if ((block = (char *) malloc(BLOCKSIZE)) == NULL) 
	    return "Insufficient memory for workspace";

	if ((res = fread (block, sizeof(char), BLOCKSIZE, fs->fp)) != BLOCKSIZE)
	    return "Error reading FITS file";
	i = 0;

	/* Read SIMPLE keyword. */
	if ((error = cdl_rdCard (block, "SIMPLE", T_LOG, &val, &rval, sval)))
	    return error;
	if (val == 0) 
	    return "Not a SIMPLE FITS file";
	i++;

	/* Read BITPIX keyword. */
	if ((error = cdl_rdCard(&block[SZ_CARD], "BITPIX", T_INT, &val, &rval,
	    sval)))
	        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 keyword. */
	if ((error=cdl_rdCard(&block[2*SZ_CARD], "NAXIS", T_INT, &val, &rval,
            sval)))
	        return error;
	if (val < 0 || val > 999)
	    return "Bad NAXIS value in FITS file";
	if (val < 2) 
	    return "FITS file is not a two-dimensional image";
	fs->naxis = val;
	i++;

	/* 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);
	    if ((error=cdl_rdCard (&block[i*SZ_CARD], name, T_INT, &val, &rval,
		sval)))
	    	    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 the remainder. */
	fs->bscale  = 1.0;
	fs->bzero   = 0.0;
	fs->object[0] = '\0';
	while (1) {

	    /* Try reading a BSCALE or BZERO keyword from this card. */
	    if (block[i*SZ_CARD] == 'B') {
	        error = cdl_rdCard(&block[i*SZ_CARD], "BSCALE", T_REAL, &val,
		    &rval, sval);
	        if (error == NULL)
	            fs->bscale = rval;
	        error = cdl_rdCard (&block[i*SZ_CARD], "BZERO", T_REAL, &val,
		    &rval, sval);
	        if (error == NULL)
	            fs->bzero = rval;
	    } else if (strncmp("OBJECT", &block[i*SZ_CARD], 6) == 0) {
	        /* Read OBJECT keyword.  */
	        if ((error = cdl_rdCard (&block[i*SZ_CARD], "OBJECT", T_STR,
		    &val, &rval, sval))) {
	                strcpy (fs->object, "no title");
			break;
		}
	        strncpy (fs->object, sval, SZ_CARD);
	        i++;
	    }

	    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*SZ_CARD], "END     ", 8) == 0) 
	    	break;
	    i++;
	}

	free (block);
	return NULL;
}


/*  CDL_RDCARD -- 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 int, right justified in columns 11-30
 *      dtype = T_REAL     # value is a real
 *      dtype = T_STR      # value is a string
 *
 *  The value is stored in kvalue.  Returns NULL on success, or an error
 *  message otherwise.
 */

#ifdef ANSI_FUNC

static char *
cdl_rdCard (
    char *card,				/* FITS keyword card 	*/
    char *name,				/* keyword name		*/
    enum datatype dtype,		/* type of value 	*/
    long *kvalue,			/* integer value 	*/
    float *rvalue,			/* real value 		*/
    char *svalue			/* string value 	*/
)
#else

static char *
cdl_rdCard (card, name, dtype, kvalue, rvalue, svalue)
char	*card;				/* FITS keyword card 	*/
char	*name;				/* keyword name		*/
enum datatype dtype;   			/* type of value 	*/
long 	*kvalue;			/* integer value 	*/
float 	*rvalue;			/* real value 		*/
char  	*svalue;			/* string value 	*/
#endif
{
	register int	i, ptr;
	char	namestr[9], *ip;

	/* Get the keyword from the card.  */
	bcopy (card, namestr, 8);
        for (i=7; i >= 0 && namestr[i] == ' '; i--)
	    ;
	namestr[i+1] = '\0';
	if (strcmp(namestr, name) != 0)
	    return "Keyword not found in FITS file.";

	/* Get start of value. */
	ptr = 10;
	while (ptr < SZ_CARD && card[ptr] == ' ') 
	    ptr++;
	if (ptr == SZ_CARD) 
	    return "FITS file has missing keyword value"; /* no value */

	if (dtype == T_LOG) {
	    *kvalue = (card[ptr] == 'T');

	} else if (dtype == T_STR) {
	    for (ip=&card[9]; *ip != '\''; ip++)
		;
	    for (i=0, ip=&card[ptr+1]; *ip != '\'' && i < SZ_FNAME; i++)
		svalue[i] = *ip++;
	    svalue[i] = '\0';
	} else {
	    /* An integer or real value. */
	    long int	ival;
	    float	fval;
	    char	num[21];

	    if (ptr > 29) 
	    	return "Keyword has bad integer value in FITS file";
	    memcpy (num, &card[ptr], 30 - ptr);
	    num[30-ptr] = '\0';
	    if (dtype == T_INT) {
	       i = sscanf (num, "%ld", &ival);
	       if (i != 1) 
	    	   return "Keyword has bad integer value in FITS file";
	       *kvalue = ival;
	       *rvalue = 0.0;
	    } else if (dtype == T_REAL) {
	       i = sscanf (num, "%g", &fval);
	       if (i != 1) 
	    	   return "Keyword has bad real value in FITS file";
	       *kvalue = 0;
	       *rvalue = fval;
	    }
	}

	return NULL;
}


/*  CDL_GETDATA -- Reads nelem values into the buffer.  Copes with the fact
 *  that the last 2880 byte record of the FITS file may be truncated, and
 *  should be padded out with zeros.  Returns NULL for success or an error
 *  message.
 */

#ifdef ANSI_FUNC

static char *
cdl_getData (FITS *fs, uchar *buffer, int nelem)
#else

static char *
cdl_getData (fs, buffer, nelem)
FITS 	*fs;
uchar 	*buffer;
int	nelem;
#endif
{
	int	res;

	if (nelem == 0) 
	    return NULL;

	if ((res = fread ((void *)buffer, fs->size, nelem, fs->fp)) != 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 cdl_fixData(fs, buffer, nelem);
}


/*  CDL_FIXDATA -- 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
 */

#ifdef ANSI_FUNC

static char *
cdl_fixData (FITS *fs, void *buffer, int nelem)
#else

static char *
cdl_fixData (fs, buffer, nelem)
FITS 	*fs;
void 	*buffer;
int	nelem;
#endif
{
	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;
}