# 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;
  int ccdsel, done;
  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) {
    fprintf (stderr, "only specify one of -n and -x\n");
    return (FALSE);
  }

  if (argc != 3) {
    fprintf (stderr, "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) {
    fprintf (stderr, "file %s not found\n", argv[2]);
    return (FALSE);
  }
  fclose (f);

  /* find matrix, free old data */
  if ((buf = SelectBuffer (argv[1], ANYBUFFER, TRUE)) == NULL) return (FALSE);
  fits_free_matrix (&buf[0].matrix);
  fits_free_header (&buf[0].header);

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

  if (extend) {
    /* load in appropriate header */
    Nskip = fits_read_Xheader (argv[2], &buf[0].header, Nextend);
    if (!Nskip) {
      fprintf (stderr, "entry in %s not found\n", argv[2]);
      DeleteBuffer (buf);
      return (FALSE);
    }
    /* fix up header */
    {
      static char simple[] = "SIMPLE  =                    T / Standard FITS";
      int Ns, No;

      Ns = strlen (simple);
      No = 80 - Ns;
      strncpy (buf[0].header.buffer, simple, Ns);
      memset (&buf[0].header.buffer[Ns], ' ', No);
    }
    if (JustHead) {
      /* should be in CreateMatrix / CreateBuffer function */
      buf[0].header.Naxes = 0;
      ALLOCATE (buf[0].matrix.buffer, char, 1);
      buf[0].matrix.size = 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;
      return (TRUE);
    }
    /* load matrix data */
    f = fopen (argv[2], "r");
    if (f == (FILE *) NULL) {
      fprintf (stderr, "file %s not found\n", argv[2]);
      DeleteBuffer (buf);
      return (FALSE);
    }
    fseek (f, Nskip, SEEK_SET);
    status = fits_load_matrix (f, &buf[0].matrix, &buf[0].header);
    fclose (f);
  } 

  if (ccdsel) {
    char *CCDKeyword, ID[64];

    CCDKeyword = get_variable ("CCDKEYWORD");
    if (CCDKeyword == (char *) NULL) {
      fprintf (stderr, "CCDKEYWORD variable is not set; ");
      fprintf (stderr, "using IMAGEID as default\n");
      CCDKeyword = strcreate ("IMAGEID");
    }
    done = FALSE;
    for (i = 0; !done; i++) {
      Nskip = fits_read_Xheader (argv[2], &buf[0].header, i);
      if (!Nskip) {
	fprintf (stderr, "extension %s in %s not found\n", ccdid, argv[2]);
	DeleteBuffer (buf);
	free (CCDKeyword);
	return (FALSE);
      }
      if (!fits_scan (&buf[0].header, CCDKeyword, "%s", 1, ID)) {
	fprintf (stderr, "%s not in header\n", CCDKeyword);
	DeleteBuffer (buf);
	free (CCDKeyword);
	return (FALSE);
      }

      /* compare as numbers if both are pure numbers, else as strings */
      done = strnumcmp (ccdid, ID);

      if (!done) fits_free_header (&buf[0].header);
    }
    free (CCDKeyword);

    /* fix up header */
    {
      static char simple[] = "SIMPLE  =                    T / Standard FITS";
      int Ns, No;

      Ns = strlen (simple);
      No = 80 - Ns;
      strncpy (buf[0].header.buffer, simple, Ns);
      memset (&buf[0].header.buffer[Ns], ' ', No);
    }
    if (JustHead) {
      buf[0].header.Naxes = 0;
      ALLOCATE (buf[0].matrix.buffer, char, 1);
      buf[0].matrix.size = 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;
      return (TRUE);
    }
    /* load matrix data */
    f = fopen (argv[2], "r");
    if (f == (FILE *) NULL) {
      fprintf (stderr, "file %s not found\n", argv[2]);
      DeleteBuffer (buf);
      return (FALSE);
    }
    fseek (f, Nskip, SEEK_SET);
    status = fits_load_matrix (f, &buf[0].matrix, &buf[0].header);
    fclose (f);
  }

  if (!ccdsel && !extend) {
    status = fits_read_header (argv[2], &buf[0].header);
    if (JustHead) {
      buf[0].header.Naxes = 0;
      ALLOCATE (buf[0].matrix.buffer, char, 1);
      buf[0].matrix.size = 0;
      buf[0].bitpix = 16;
      buf[0].bzero = 0;
      buf[0].matrix.bscale = 1;
      buf[0].matrix.bitpix = 16;
      buf[0].matrix.bzero = 0;
      buf[0].matrix.bscale = 1;
      return (TRUE);
    }
    Nplane = buf[0].header.Naxis[2];
    if (Nplane > 0) {
      if (plane > Nplane) {
	fprintf (stderr, "-plane is too large: %d total planes\n", Nplane);
	DeleteBuffer (buf);
	return (FALSE);
      }
      buf[0].header.Naxis[2] = 0;
      buf[0].header.Naxes = 2;
      fits_modify (&buf[0].header, "NAXIS", "%d", 1, 2);
      fits_delete (&buf[0].header, "NAXIS3", 1);
    }
    sprintf (region, "-1 -1 -1 -1 %d %d", (plane - 1), plane);
    status = fits_read_segment (argv[2], &buf[0].matrix, region);
  }

  if (!status) {
    fprintf (stderr, "problem reading file, buffer not opened\n");
    DeleteBuffer (buf);
    return (FALSE);
  }
  if (buf[0].header.Naxes == 1) {
    /* we need to return an array, so make Naxis[1] = 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;
  fprintf (GetOutfile(), "read %d bytes from %s into buffer %s\n", 
	   buf[0].header.size + buf[0].matrix.size, argv[2], argv[1]);

  /** now - convert the matrix values to floats for internal use **/
  fits_convert_format (&buf[0].header, &buf[0].matrix, -32, 1.0, 0.0, FT_UNSIGN_MODE);
  
  return (TRUE);

}
