# include <ohana.h>
# include <gfitsio.h>

int gfits_fread_matrix (FILE *f, Matrix *matrix, Header *header) {
  
  int status;
  
  status = gfits_load_matrix (f, matrix, header);
  return (status);
}

/*********************** fits read matrix ***********************************/
int gfits_load_matrix (FILE *f, Matrix *matrix, Header *header) {

  int i, nbytes, Nbytes;
  unsigned char *byte0, *byte1, *byte2, *byte3, *byte4, *byte5, *byte6, *byte7, tmp;
  int perpix;

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

  Nbytes = gfits_matrix_size (header);
  ALLOCATE (matrix[0].buffer, char, Nbytes);
  matrix[0].size = Nbytes;

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

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

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

  return (TRUE);

}

