#ifndef lint
static char SccsId[] = "%W%  %G%";
#endif

/* Module:	readfits.c (Read FITS)
 * Purpose:	Read a FITS file header and prepare file for basic array read
 * Subroutine:	int init_fits()
 * Note:	Only SIMPLE = T is permitted, but will handle BITPIX = 8, 16,
 *		32, -16(u_short), -32(float), -64(double).  If img->nimage
 *		is set to a value greater than 1, file will be lseeked to the
 *		nimage'th image (assuming there are several images stored as
 *		a stack (on 3rd dimension).
 * Copyright:	1989 Smithsonian Astrophysical Observatory
 *		You may do anything you like with this file except remove
 *		this copyright.  The Smithsonian Astrophysical Observatory
 *		makes no representations about the suitability of this
 *		software for any purpose.  It is provided "as is" without
 *		express or implied warranty.
 * Modified:	{0} Michael VanHilst	initial version	     23 January 1989
 *              {1} Doug Mink           implement nimage     11 October 1990
 *              {2} Doug Mink           fix dimensions        3 May     1995
 *              {3} Doug Mink           move memset          18 October 1995
 *              {4} Doug Mink           set 3rd dimension   14 February 1997
 *		{n} <who> -- <does what> -- <when>
 */

#include <stdio.h>		/* define stderr, NULL, etc. */
#include <X11/Xlib.h>		/* needed for control.h */
#include "hfiles/constant.h"
#include "hfiles/control.h"	/* define IOP codes */
#include "hfiles/image.h"

#define FITSBLOCK 2880

/*
 * Subroutine:	init_fits
 * Purpose:	Open fits file, get parameters and get ready to read data.
 * Returns:	-1 if failure, else fd of file, positioned at start of data
 * Note:	Data read will be handled as standard array read.
 */
int init_fits ( img )
     struct imageRec *img;
{
  float scale, bias;
  int fd;
  int headlen, image_start;
  int bitpix, naxes, naxis[3];
  char *header;
  static char *load_fitsheader();
  int read_fitsheader(), open_disk(), lseek_disk();
  void close_disk();

  /* open the image file */
  if( (fd = open_disk(img->filename, IOP_Read, 0)) <= 0 )
    return( -1 );
  headlen = img->headersize;
  /* read and interpret the header */
  if( (header = load_fitsheader(fd, &headlen, img->filename)) == NULL ) {
    close_disk(fd, img->filename);
    return( -1 );
  }
  if( read_fitsheader(header, headlen, &bitpix, &naxes, naxis, &scale, &bias)
      != 0 ) {
    free(header);
    switch( bitpix ) {
    case 8:
      img->storage_type = ARR_U1;
      img->bytepix = 1;
      break;
    case 16:
      img->storage_type = ARR_I2;
      img->bytepix = 2;
      break;
    case 32:
      img->storage_type = ARR_I4;
      img->bytepix = 4;
      break;
    case -16:
      img->storage_type = ARR_U2;
      img->bytepix = 2;
      break;
    case -32:
      img->storage_type = ARR_R4;
      img->bytepix = 4;
      break;
    case -64:
      img->storage_type = ARR_R8;
      img->bytepix = 8;
      break;
    default:
      (void)fprintf(stderr,"Illegal FITS BITPIX: %d\n",bitpix);
      close_disk(fd, img->filename);
      return( -1 );
    }
    if( img->nimage > 1 ) {
      if( (naxes <= 2) || (naxis[2] < img->nimage) ) {
	(void)fprintf(stderr,
		      "Only %d images in file %s\n", naxis[2], img->filename);
	if (naxis[2] > 0)
	    img->nimage = naxis[2];
	else
	    img->nimage = 1;

	/* close_disk(fd, img->filename);
	return( -1 ); */
      }
      image_start = headlen + 
			((img->nimage-1) * naxis[0] * naxis[1] * img->bytepix);
      if( lseek_disk(fd, image_start, img->filename) < 0 ) {
	img->headersize = 0;
	close_disk(fd, img->filename);
	return( -1 );
      }
    } else {
      img->headersize = headlen;
    }
    img->filecols = naxis[0];
    img->filerows = naxis[1];
    img->filenimg = naxis[2];
    if( (scale != 1.0) || (bias != 0.0) ) {
      img->fscale = scale;
      img->fbias = bias;
      img->fscaled = 1;
    } else
      img->fscaled = 0;
    return( fd );
  } else {
    free(header);
    close_disk(fd, img->filename);
    return( -1 );
  }
}


/*
 * Subroutine:	load_fitsheader
 * Purpose:	Load the FITS header from the disk file into memory.
 * Notes:	FITS headers are made in block units (2880 bytes each).
 *		If length is not initially zero, it is taken as a directive
 *		to override FITS standard (header of given size is loaded).
 *		FITS header must have "SIMPLE =" in first 10 bytes (and 'T'
 *		in next 20 for ximage use).
 */
static char *load_fitsheader ( fd, length, filename )
     int fd;
     int *length;
     char *filename;
{
  int headlen, nbytes;
  char *header;
  char TorF[4];
  static char *get_keyfield();
  char *realloc(), *calloc_errchk();
  int read_disk();
  void no_fitscomment();

  if( *length > 0 )
    headlen = *length;
  else
    headlen = FITSBLOCK;
  /* allocate initial sized buffer */
  header = calloc_errchk (headlen+1, 1, "FITS header");
  /* read in first block */
  if( (nbytes = read_disk(fd, header, headlen, 1, filename, "FITS header"))
     != headlen ) {
    free (header);
    return( NULL );
  }
  /* consistency check for FITS header */
  if( strncmp(header, "SIMPLE  =", 9) != 0) {
    (void)fprintf(stderr,"No SIMPLE keyword in FITS header");
    free(header);
    return( NULL );
  }
  /* we only support SIMPLE = T */
  no_fitscomment(header+10, 20);
  if( (sscanf(header+10, "%1s", TorF) <= 0) || (TorF[0] != 'T') ) {
    (void)fprintf(stderr,"SIMPLE = %c not supported\n", TorF[0]);
    free(header);
    return( NULL );
  }
  if( *length != 0 ) {
    /* if header length was given, it must check out with an END */
    if( get_keyfield(header, "END     ", headlen, 0) == NULL ) {
      (void)fprintf(stderr,"%d byte FITS header has no END card!\n", headlen);
      free(header);
      return( NULL );
    }
  } else {
    /* if END keyword not found, read in more header until END is found */
    while( get_keyfield(header, "END     ", headlen, 0) == NULL ) {
      if ((header = realloc(header,(unsigned)(headlen+FITSBLOCK+1))) == NULL) {
	fputs("Reallocation failure reading header\n", stderr);

	return( NULL );
      }
      /* initialize the newly allocated memory */
      memset (header + headlen, 0, FITSBLOCK + 1);
      if( (nbytes = read_disk(fd, header+headlen, FITSBLOCK, 1,
			      filename, "extended header")) != FITSBLOCK ) {
	(void)fprintf(stderr, "Read %d bytes of header with no END card!\n",
		      headlen + nbytes);
	(void)fflush(stderr);
	free( header );
	return( NULL );
      }
      headlen += FITSBLOCK;
    }
    *length = headlen;
  }
  return( header );
}


/*
 * Subroutine:	get_keyfield
 * Purpose:	Return the data field for a given FITS header keyword
 * Returns:	If key not found, return NULL (0).
 */
static char *get_keyfield ( header, keyword, length, report_error )
     char *header;	/* buffer start */
     char *keyword;	/* keyword to match */
     int length;	/* if zero, search up to "END" keyword */
     int report_error;	/* if > 0, fatal if key not found */
{
  int key_not_end, i;
  void no_fitscomment();

  key_not_end = (strncmp(keyword, "END     ", 8) != 0);
  for( i=0; i<length; i+=80 ) {
    /* check for END keyword marking end of header, unless END is the key */
    if( key_not_end && strncmp(header+i,"END     ",8) == 0 )
      break;
    /* check for desired keyword */
    if( strncmp(header+i,keyword,8) == 0 ) {
      no_fitscomment (header+i+10, 20);
      return( header+i+10 );
    }
  }
  if( report_error )
    (void)fprintf(stderr, "No `%s' keyword in FITS header\n", keyword);
  return( (char *)NULL );
}