# include "data.h"

/* there is some confusion in this function with the several possible options */
int rd (int argc, char **argv) {
  
  int i, N, status, plane, Nplane, extend, Nextend, Nskip, JustHead, blank;
  int ccdsel, done, Nword, IsCompressed;
  char region[512], *ccdid, *filename;
  FILE *f;
  Buffer *buf;

  JustHead = FALSE;
  if ((N = get_argument (argc, argv, "-head"))) {
    remove_argument (N, &argc, argv);
    JustHead = TRUE;
  }

  plane = 1;
  if ((N = get_argument (argc, argv, "-plane"))) {
    remove_argument (N, &argc, argv);
    plane  = atof(argv[N]);
    remove_argument (N, &argc, argv);
  }

  extend = FALSE;
  Nextend = -1;
  if ((N = get_argument (argc, argv, "-x"))) {
    remove_argument (N, &argc, argv);
    Nextend  = atof(argv[N]);
    remove_argument (N, &argc, argv);
    extend = TRUE;
  }

  ccdsel = FALSE;
  ccdid = (char *) NULL;
  if ((N = get_argument (argc, argv, "-n"))) {
    remove_argument (N, &argc, argv);
    ccdid  = strcreate (argv[N]);
    remove_argument (N, &argc, argv);
    ccdsel = TRUE;
  }

  if (extend && ccdsel) {
    gprint (GP_ERR, "only specify one of -n and -x\n");
    return (FALSE);
  }

  if (argc != 3) {
    gprint (GP_ERR, "USAGE: rd <buffer> <filename> [-head] [-plane N] [-n ccdid] [-x extnum]\n");
    return (FALSE);
  }

  /* test if file exists */
  f = fopen (argv[2], "r");
  if (f == (FILE *) NULL) {
    gprint (GP_ERR, "file %s not found\n", argv[2]);
    return (FALSE);
  }

  /* find matrix, free old data */
  if ((buf = SelectBuffer (argv[1], ANYBUFFER, TRUE)) == NULL) {
      fclose (f);
      return (FALSE);
  }
  gfits_free_matrix (&buf[0].matrix);
  gfits_free_header (&buf[0].header);

  /* save file name */
  filename = filebasename (argv[2]);
  strcpy (buf[0].file, filename);
  free (filename);

  status = FALSE;
  IsCompressed = FALSE;

  /*** advance to the correct FITS extension ***/

  /* FITS extension by number */
  if (extend) {
    /* load in appropriate header */
    Nskip = gfits_fread_Xheader (f, &buf[0].header, Nextend);
    if (!Nskip) {
      gprint (GP_ERR, "entry in %s not found\n", argv[2]);
      DeleteBuffer (buf);
      fclose (f);
      return (FALSE);
    }
    if (gfits_extension_is_compressed (&buf[0].header)) {
	IsCompressed = TRUE;
    }
  } 

  /* FITS extension by name */
  if (ccdsel) {
    char *CCDKeyword, ID[64];

    CCDKeyword = get_variable ("CCDKEYWORD");
    if (CCDKeyword == (char *) NULL) {
      // gprint (GP_ERR, "CCDKEYWORD variable is not set; ");
      // gprint (GP_ERR, "using EXTNAME as default\n");
      CCDKeyword = strcreate ("EXTNAME");
    }
    done = FALSE;
    for (i = 0; !done; i++) {
      Nskip = gfits_fread_Xheader (f, &buf[0].header, i);
      if (!Nskip) {
	gprint (GP_ERR, "extension %s in %s not found\n", ccdid, argv[2]);
	DeleteBuffer (buf);
	free (CCDKeyword);
	fclose (f);
	return (FALSE);
      }

      // for compressed data tables, EXTNAME may be duplicated, with the first one containing the
      // word 'COMPRESSED_IMAGE'.  in this case, check the second EXTNAME, if CCDKeyword == EXTNAME 
      // this may have to be a more obscure test specifically for 'imcopy' data...
      // need to check each header, since file may contain a mix
      
      Nword = 1;
      IsCompressed = FALSE;
      if (gfits_extension_is_compressed (&buf[0].header)) {
	if (!strcmp (CCDKeyword, "EXTNAME")) Nword = 1;
	IsCompressed = TRUE;
      }
      if (!gfits_scan (&buf[0].header, CCDKeyword, "%s", Nword, ID)) {
	gprint (GP_ERR, "%s not in header\n", CCDKeyword);
	DeleteBuffer (buf);
	free (CCDKeyword);
	fclose (f);
	return (FALSE);
      }

      /* compare as numbers if both are pure numbers, else as strings */
      done = strnumcmp (ccdid, ID);
      if (!done) gfits_free_header (&buf[0].header);
    }
    free (CCDKeyword);
  }

  /* fix up header, if needed */
  if (extend || ccdsel) {
    if (!IsCompressed) {
      gfits_extended_to_primary (&buf[0].header, TRUE, "Standard FITS");
    }
  } else {
      gfits_fread_header (f, &buf[0].header);
  }

  /* for JustHead, we skip reading the data segment */ 
  // XXX for compressed data, we need to convert the header to the equivalent uncompressed version
  if (JustHead) {
    // XXX what are we doing here exactly?
    buf[0].header.Naxes = 0;
    ALLOCATE (buf[0].matrix.buffer, char, 1);
    buf[0].matrix.datasize = 0;
    buf[0].bitpix = 16;
    buf[0].bzero = 0;
    buf[0].bscale = 1;
    buf[0].matrix.bitpix = 16;
    buf[0].matrix.bzero = 0;
    buf[0].matrix.bscale = 1;
    buf[0].header.bitpix = 16;
    buf[0].header.bzero = 0;
    buf[0].header.bscale = 1;
    fclose (f);
    return (TRUE);
  }

  /* check for valid plane */
  Nplane = buf[0].header.Naxis[2];
  if (Nplane == 0) Nplane = 1;
  if (plane > Nplane) {
    gprint (GP_ERR, "-plane is too large: %d total planes\n", Nplane);
    DeleteBuffer (buf);
    fclose (f);
    return (FALSE);
  }

  /* load matrix data */
  if (IsCompressed) {
    FTable ftable;
    Header theader;
    ftable.header = &theader;
    ftable.header[0].buffer = NULL;
    gfits_copy_header (&buf[0].header, ftable.header);
    status = gfits_fread_ftable_data (f, &ftable);  // this just reads the bytes (not even a SWAP)
    status = gfits_uncompress_image (&buf[0].header, &buf[0].matrix, &ftable);
    // uncompressing the image leaves the format as an extension
    gfits_extended_to_primary (&buf[0].header, TRUE, "Standard FITS");
    gfits_free_table (&ftable);
    // XXX this currently does not work for a cube (we get a cube back, not a specific plane)
  } else {
    sprintf (region, "-1 -1 -1 -1 %d %d", (plane - 1), plane);
    status = gfits_fread_matrix_segment (f, &buf[0].matrix, &buf[0].header, region);
  }
  fclose (f);

  if (!status) {
    gprint (GP_ERR, "problem reading file, buffer not opened\n");
    DeleteBuffer (buf);
    return (FALSE);
  }

  /* adjust buffer to represent 2D data */
  if (Nplane > 1) {
    buf[0].header.Naxis[2] = 0;
    buf[0].header.Naxes = 2;
    gfits_modify (&buf[0].header, "NAXIS", "%d", 1, 2);
    gfits_delete (&buf[0].header, "NAXIS3", 1);
  }

  /* we need to return a 2D array, convert 1D images to 2D (Naxis[1] = 1) */
  if (buf[0].header.Naxes == 1) {
    buf[0].header.Naxes = 2;
    buf[0].header.Naxis[1] = 1;
    buf[0].matrix.Naxis[1] = 1;
  }    

  buf[0].bitpix = buf[0].header.bitpix;    /* store the original values */
  buf[0].bscale = buf[0].header.bscale;    /* store the original values */
  buf[0].bzero  = buf[0].header.bzero;     /* store the original values */
  buf[0].unsign = buf[0].header.unsign;

  gprint (GP_LOG, "read %lld bytes from %s into buffer %s\n", 
	  (long long) buf[0].header.datasize + buf[0].matrix.datasize, argv[2], argv[1]);

  blank = 0xffff;
  gfits_scan (&buf[0].header, "BLANK", "%d", 1, &blank);

  /** now - convert the matrix values to floats for internal use **/
  gfits_convert_format (&buf[0].header, &buf[0].matrix, -32, 1.0, 0.0, blank, gfits_get_unsign_mode());

  return (TRUE);
}
