# include "fits.h"

/** warning: this function requires 'dparse' which is included in libohana.a **/
int dparse (double *, int, char *region);

/*********************** fits read matrix ***********************************/
int fits_read_segment (char *filename, Matrix *matrix, char *region) {

  FILE *f;
  Header header;
  int status, i, nbytes, Nbytes, Nskip, NbytesData;
  int wantaxis[FT_MAX_NAXES][2];
  unsigned char *byte0, *byte1, *byte2, *byte3, *byte4, *byte5, *byte6, *byte7, ctmp;
  double tmp;
  int perpix;

  status = fits_read_header (filename, &header);
  if (!status) {
    fprintf (stderr, "error reading header of FITS file %s\n", filename);
    return (FALSE);
  }

  f = fopen (filename, "r");
  if (f == NULL) {
    return (FALSE);
  }

  for (i = 0; i < FT_MAX_NAXES; i++) {
    wantaxis[i][0] = wantaxis[i][1] = 0;
  }

  for (i = 0; i < header.Naxes; i++) {
    status = dparse (&tmp, 2*i + 1, region);
    if (!status || (tmp < 0))
      wantaxis[i][0] = 0;
    else 
      wantaxis[i][0] = (tmp > header.Naxis[i] - 1) ? header.Naxis[i] - 1 : tmp;
    status = dparse (&tmp, 2*i + 2, region);
    if (!status || (tmp < 0))
      wantaxis[i][1] = header.Naxis[i];
    else 
      wantaxis[i][1] = (tmp > header.Naxis[i]) ? header.Naxis[i] : tmp;
  }

  fseek (f, header.size, 0);

  matrix[0].bitpix = header.bitpix;
  matrix[0].unsign = header.unsign;
  matrix[0].bscale = header.bscale;
  matrix[0].bzero  = header.bzero;
  matrix[0].Naxes  = header.Naxes;
  for (i = 0; i < FT_MAX_NAXES; i++)
    matrix[0].Naxis[i] = wantaxis[i][1] - wantaxis[i][0];

  matrix[0].size = abs (matrix[0].bitpix) / 8;
  for (i = 0; i < matrix[0].Naxes; i++)
    matrix[0].size *= matrix[0].Naxis[i];

  Nskip = wantaxis[2][0]*matrix[0].Naxis[0]*matrix[0].Naxis[1]*abs(matrix[0].bitpix)/8;
  NbytesData = Nskip + matrix[0].size;

  if (NbytesData % FT_RECORD_SIZE) 
    Nbytes = FT_RECORD_SIZE * ((int) (NbytesData / FT_RECORD_SIZE) + 1) - Nskip;
  else 
    Nbytes = FT_RECORD_SIZE * ((int) (NbytesData / FT_RECORD_SIZE)) - Nskip;

  ALLOCATE (matrix[0].buffer, char, MAX (Nbytes, 1));

# ifndef SEEK_CUR 
# define SEEK_SET 0   
# define SEEK_CUR 1   
# define SEEK_END 2
# endif

  /* currently only good for reading in full planes in 3-D */
  fseek (f, Nskip, SEEK_CUR);

  nbytes = fread (matrix[0].buffer, sizeof(char), Nbytes, f);
  fclose (f);
  matrix[0].size = Nbytes;

# ifdef BYTE_SWAP  
  perpix = abs(header.bitpix) / 8;
  if (perpix > 1) {
    byte0 = (char *) matrix[0].buffer;
    byte1 = (char *) matrix[0].buffer + 1;
    byte2 = (char *) matrix[0].buffer + 2;
    byte3 = (char *) matrix[0].buffer + 3;
    byte4 = (char *) matrix[0].buffer + 4;
    byte5 = (char *) matrix[0].buffer + 5;
    byte6 = (char *) matrix[0].buffer + 6;
    byte7 = (char *) matrix[0].buffer + 7;
    if (perpix == 2) {
      for (i = 0; i < nbytes; i+=2, byte0 += 2, byte1 += 2) {
	ctmp = *byte0;
	*byte0 = *byte1;
	*byte1 = ctmp;
      }
    }
    if (perpix == 4) {
      for (i = 0; i < nbytes; i+=4, byte0 += 4, byte1 += 4, byte2 += 4, byte3 += 4) {
	ctmp = *byte0;
	*byte0 = *byte3;
	*byte3 = ctmp;
	ctmp = *byte1;
	*byte1 = *byte2;
	*byte2 = ctmp;
      }
    }
    if (perpix == 8) {
      for (i = 0; i < nbytes; i+=8, byte0 += 8, byte1 += 8, byte2 += 8, byte3 += 8, byte4 += 8, byte5 += 8, byte6 += 8, byte7 += 8) {
	ctmp = *byte0;
	*byte0 = *byte7;
	*byte7 = ctmp;
	ctmp = *byte1;
	*byte1 = *byte6;
	*byte6 = ctmp;
	ctmp = *byte1;
	*byte2 = *byte5;
	*byte5 = ctmp;
	ctmp = *byte1;
	*byte3 = *byte4;
	*byte4 = ctmp;
      }
    }
  }
# endif

  fits_free_header (&header);

  if (nbytes < Nbytes - 2880) {  /* this is a bad FITS error: image is not OK */
    fprintf (stderr, "error reading in matrix data from FITS file %s (%d < %d - 2880)\n", filename, nbytes, Nbytes);
    return (FALSE);
  }
  if (nbytes != Nbytes) {  /* this is a FITS error, but often the image is OK */
    fprintf (stderr, "incomplete block in %s: (%d, %d)\n", filename, nbytes, Nbytes);
    return (TRUE); 
  }

  return (TRUE);

}

