# include "dvoshell.h"
# define D_NSTARS 1000
# define BYTES_STAR 31
# define BLOCK 1000

int simage (int argc, char **argv) {

  char *buffer;
  Vector Xvec, Yvec, Zvec;
  double R, D, X, Y, M, zero, range;
  FILE *f;
  Header header;
  Coords coords;
  int i, j, Nstars, nstars, Nbytes, nbytes, Npts, N;
  Graphdata graphmode;
  int Ngraph;

  Ngraph = 0;
  if (!GetGraph (&graphmode, NULL, &Ngraph)) return (FALSE);

  zero = 17.0;
  range = -5.0;
  if ((N = get_argument (argc, argv, "-m"))) {
    remove_argument (N, &argc, argv);
    range = atof(argv[N]);
    remove_argument (N, &argc, argv);
    zero  = atof(argv[N]);
    range = range - zero;
    remove_argument (N, &argc, argv);
  }

  if (argc != 2) {
    gprint (GP_ERR, "USAGE: image (filename)\n");
    return (FALSE);
  }

  gprint (GP_ERR, "not working at the moment (cmp format)\n");
  return (FALSE);
  
  /* read header */
  if (!gfits_read_header (argv[1], &header)) {
    gprint (GP_ERR, "ERROR: can't find image file %s\n", argv[1]);
    return (FALSE);
  }
  /* get astrometry information */
  strcpy (coords.ctype, "NONE");
  gfits_scan (&header, "CTYPE1",   "%s",  1, coords.ctype);
  if (strcmp (coords.ctype, "RA---PLY")) {
    gprint (GP_ERR, "ERROR: wrong astrometric info in header\n");
    return (FALSE);
  }
  gfits_scan (&header, "CDELT1",   "%f", 1, &coords.cdelt1); 
  gfits_scan (&header, "CDELT2",   "%f", 1, &coords.cdelt2);
  gfits_scan (&header, "CRVAL1",   "%lf", 1, &coords.crval1);
  gfits_scan (&header, "CRVAL2",   "%lf", 1, &coords.crval2);  
  gfits_scan (&header, "CRPIX1",   "%f", 1, &coords.crpix1);
  gfits_scan (&header, "CRPIX2",   "%f", 1, &coords.crpix2);
  gfits_scan (&header, "PC001001", "%f", 1, &coords.pc1_1);
  gfits_scan (&header, "PC001002", "%f", 1, &coords.pc1_2);
  gfits_scan (&header, "PC002001", "%f", 1, &coords.pc2_1);
  gfits_scan (&header, "PC002002", "%f", 1, &coords.pc2_2);
  /* RA Terms */
  gfits_scan (&header, "PCA1X2Y0", "%f", 1, &coords.polyterms[0][0]);
  gfits_scan (&header, "PCA1X1Y1", "%f", 1, &coords.polyterms[1][0]);
  gfits_scan (&header, "PCA1X0Y2", "%f", 1, &coords.polyterms[2][0]);
  gfits_scan (&header, "PCA1X3Y0", "%f", 1, &coords.polyterms[3][0]);
  gfits_scan (&header, "PCA1X2Y1", "%f", 1, &coords.polyterms[4][0]);
  gfits_scan (&header, "PCA1X1Y2", "%f", 1, &coords.polyterms[5][0]);
  gfits_scan (&header, "PCA1X0Y3", "%f", 1, &coords.polyterms[6][0]);
  /* Dec Terms */			    
  gfits_scan (&header, "PCA2X2Y0", "%f", 1, &coords.polyterms[0][1]);
  gfits_scan (&header, "PCA2X1Y1", "%f", 1, &coords.polyterms[1][1]);
  gfits_scan (&header, "PCA2X0Y2", "%f", 1, &coords.polyterms[2][1]);
  gfits_scan (&header, "PCA2X3Y0", "%f", 1, &coords.polyterms[3][1]);
  gfits_scan (&header, "PCA2X2Y1", "%f", 1, &coords.polyterms[4][1]);
  gfits_scan (&header, "PCA2X1Y2", "%f", 1, &coords.polyterms[5][1]);
  gfits_scan (&header, "PCA2X0Y3", "%f", 1, &coords.polyterms[6][1]);
  coords.Npolyterms = 2; /* how many do we use? */

  /* find number of stars */
  gfits_scan (&header, "NSTARS", "%d", 1, &Nstars);
  if (Nstars == 0) {
    gprint (GP_ERR, "no stars in file\n");
    return (FALSE);
  }

  /* open file data */
  f = fopen (argv[1], "r");
  if (f == NULL) {
    gprint (GP_ERR, "can't find data in file %s\n", argv[1]);
    return (FALSE);
  }
  fseek (f, header.size, SEEK_SET); 

  /* set up storage buffers */
  Xvec.Nelements = Nstars;
  Yvec.Nelements = Nstars;
  Zvec.Nelements = Nstars;
  ALLOCATE (Xvec.elements, float, Xvec.Nelements);
  ALLOCATE (Yvec.elements, float, Yvec.Nelements);
  ALLOCATE (Zvec.elements, float, Zvec.Nelements);
  ALLOCATE (buffer, char, (BLOCK*BYTES_STAR));

  /* load in stars by blocks of 1000 */
  nstars = 0;
  Nbytes = Nstars*BYTES_STAR;
  for (i = 0; i < (int)(Nbytes / (BLOCK*BYTES_STAR)); i++) {
    nbytes = fread (buffer, 1, (BLOCK*BYTES_STAR), f);
    if (nbytes != BLOCK*BYTES_STAR) {
      gprint (GP_ERR, "failed to read in stars (1)\n");
      free (Xvec.elements);
      free (Yvec.elements);
      free (Zvec.elements);
      free (buffer);
      return (FALSE);
    }
    for (j = 0; j < BLOCK; j++, nstars++) {
      dparse (&X,  1, &buffer[j*BYTES_STAR]);
      dparse (&Y,  2, &buffer[j*BYTES_STAR]);
      dparse (&M,  3, &buffer[j*BYTES_STAR]);
      XY_to_RD (&R, &D, X, Y, &coords);
      fRD_to_XY (&Xvec.elements[nstars], &Yvec.elements[nstars], R, D, &graphmode.coords);
      Zvec.elements[nstars] = MIN (1.0, MAX (0.01, (M - zero) / range));
    }
  }
  /* left over fraction of a block */
  nbytes = fread (buffer, 1, (Nbytes % (BLOCK*BYTES_STAR)), f);
  if (nbytes != (Nbytes % (BLOCK*BYTES_STAR))) {
    gprint (GP_ERR, "ERROR: failed to read in stars (2)\n");
    free (Xvec.elements);
    free (Yvec.elements);
    free (Zvec.elements);
    free (buffer);
    return (FALSE);
  }
  for (j = 0; j < nbytes / BYTES_STAR; j++, nstars++) {
    dparse (&X,  1, &buffer[j*BYTES_STAR]);
    dparse (&Y,  2, &buffer[j*BYTES_STAR]);
    dparse (&M,  3, &buffer[j*BYTES_STAR]);
    XY_to_RD (&R, &D, X, Y, &coords);
    fRD_to_XY (&Xvec.elements[nstars], &Yvec.elements[nstars], R, D, &graphmode.coords);
    Zvec.elements[nstars] = MIN (1.0, MAX (0.01, (M - zero) / range));
  }
  
  if (nstars != Nstars) {
    gprint (GP_ERR, "ERROR: failed to read in all stars (%d of %d)\n", nstars, Nstars);
    free (Xvec.elements);
    free (Yvec.elements);
    free (Zvec.elements);
    free (buffer);
    return (FALSE);
  }

  graphmode.style = 2;
  graphmode.size = -1;
  graphmode.etype = 0;
  Npts = Xvec.Nelements;
  PrepPlotting (Npts, &graphmode);

  PlotVector (Npts, Xvec.elements);
  PlotVector (Npts, Yvec.elements);
  PlotVector (Npts, Zvec.elements);

  free (Xvec.elements);
  free (Yvec.elements);
  free (Zvec.elements);
  free (buffer);

  return (TRUE);

}


