#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:	1998 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.
 */

#include <stdio.h>		/* define stderr, NULL, etc. */

#ifndef VMS
#ifdef SYSV
#include <string.h>
#else
#include <strings.h>          /* strlen, etc. for unenlightened BSD's */
#endif
#else
#include <string.h>
#endif
#include <X11/Xlib.h>		/* needed for control.h */
#include "hfiles/constant.h"
#include "hfiles/image.h"
#include "hfiles/control.h"	/* define IOP codes */
#include "hfiles/extern.h"
#include "wcs.h"
#include "fitsfile.h"

/*
 * 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, nbhead;
  int bitpix, naxes, naxis[3];
  char *header, *ext, *mwcs, cext, *rbrac;
  int lseek_disk();
  struct WorldCoor *wcsinit();

  /* Read the FITS header from the file */
  header = fitsrhead (img->filename, &headlen, &nbhead);
  if (header == NULL)
    return (-1);

  headlen = img->headersize;

  /* read and interpret the header */
  hgeti4 (header, "BITPIX", &bitpix);
  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);
      return (-1);
    }
  hgeti4 (header, "NAXIS", &naxes);
  naxis[0] = 1;
  if (naxes > 0)
    hgeti4 (header, "NAXIS1", &naxis[0]);
  naxis[1] = 1;
  if (naxes > 1)
    hgeti4 (header, "NAXIS2", &naxis[1]);
  naxis[2] = 1;
  if (naxes > 2)
    hgeti4 (header, "NAXIS3", &naxis[2]);
  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;
      }
    }
  else {
    img->headersize = headlen;
    }
  img->filecols = naxis[0];
  img->filerows = naxis[1];
  img->filenimg = naxis[2];

  scale = 1.0;
  hgetr4 (header, "BSCALE", &scale);
  bias = 0.0;
  hgetr4 (header, "BZERO", &bias);

  if (scale != 1.0 || bias != 0.0) {
    img->fscale = scale;
    img->fbias = bias;
    img->fscaled = 1;
    }
  else
    img->fscaled = 0;

  /* Check for FITS WCS specification and ignore for file opening */
    mwcs = strchr (img->filename, '%');
    if (mwcs != NULL)
      *mwcs = (char) 0;

  /* Set file pointer to start of data */
  rbrac = NULL;
  ext = strchr (img->filename, ',');
  if (ext != NULL) {
    cext = *ext;
    *ext = 0;
    }
  else {
    ext = strchr (img->filename, '[');
    if (ext != NULL) {
      cext = *ext;
      rbrac = strchr (img->filename, ']');
      if (rbrac != NULL)
        *rbrac = (char) 0;
      }
    }

  fd = open_disk(img->filename, IOP_Read, 0);
  if (ext != NULL)
    *ext = cext;
  if (rbrac != NULL)
    *rbrac = ']';
  if (mwcs != NULL)
    *mwcs = '%';
  if (fd <= 0 )
    return( -1 );

  /* Initialize WCS structure */
  if (mwcs != NULL)
    wcs = wcsinitn (header, mwcs+1);
  else
    wcs = wcsinit (header);
  (void)setwcslin (wcs, 1);
  if (control.verbose) {
    if (wcs)
      wcscent (wcs);
    else
      wcserr ();
    }
  wcsoutinit (wcs, getwcscoor());

  if (img->nimage > 1 )
    nbhead = nbhead + ((img->nimage - 1) * naxis[0] * naxis[1] * img->bytepix);
  if (lseek_disk (fd, nbhead, img->filename)) {
    free (header);
    return (-1);
    }
  else {
    img->header = header;
    return (fd);
    }
}
/* Jan 23 1989	Michael VanHilst Initial version
 *
 * Oct 11 1990  Doug Mink Implement nimage
 *
 * May  3 1995  Doug Mink Fix dimensions
 * Oct 18 1995  Doug Mink Move memset
 *
 * Feb 14 1997  Doug Mink Set 3rd dimension
 * Oct  9 1997  Doug Mink Use WCS FITSIO subs
 *
/* Jul  7 1998	Change setlinmode() to setwcslin()
 * Jul 13 1998  Doug Mink Set 3rd dimension
 * Jul 13 1998	Move file pointer to appropriate frame in 3-D image
 * Jul 17 1998  Doug Mink Change wcs setting
 * Jul 17 1998	Use getwcscom() and getwcscoor() instead of external variables
 * Aug 12 1998  Doug Mink Include fitsfile.h instead of fitsio.h
 * Aug 14 1998  Doug Mink Change WCS command initialization
 * Aug 14 1998	Use setwcscom() instead of wcscominit()
 * Sep 29 1998  Doug Mink Drop setwcscom() as it is in wcsinit()
 * Oct 14 1998  Doug Mink Save FITS header in image struct
 *
 * Feb 21 2001  Doug Mink Specify multiple WCS's by :name or :char after filename
 * Mar  8 2001	Doug Mink Change WCS specification from : to %c
 * Nov 14 2001	Initialize uninitialized parameter rbrac (from David Berry)
 */