# include "dvoshell.h"
# define NCHUNK 10000

int pmeasure (int argc, char **argv) {
  
  FILE *f;
  int i, j, k, m, kapa, Narg, Npts, NPTS, status, VERBOSE;
  double Mz, Mr, mag;
  double Radius, Rmin, Rmax, R, D;
  unsigned IDclip, IDchoice, LimExclude;
  unsigned FlagChoice, FlagClip, flags;
  int PhotcodeClip;
  float *Xvec, *Yvec, *Zvec;

  SkyTable *sky;
  SkyList *skylist;
  Catalog catalog;
  Graphdata graphmode;

  if (!InitPhotcodes ()) return (FALSE);
  if (!style_args (&graphmode, &argc, argv, &kapa)) return FALSE;

  f = (FILE *) NULL;
  Mz = 17.0;
  Mr = -5.0;
  Rmin = graphmode.coords.crval1 - 182.0;
  Rmax = graphmode.coords.crval1 + 182.0;

  PhotcodeClip = -1;
  if ((Narg = get_argument (argc, argv, "-p"))) {
    remove_argument (Narg, &argc, argv);
    PhotcodeClip = GetPhotcodeCodebyName (argv[Narg]);
    remove_argument (Narg, &argc, argv);
  }
  IDchoice = 0;
  IDclip = FALSE;
  if ((Narg = get_argument (argc, argv, "-ID"))) {
    remove_argument (Narg, &argc, argv);
    IDchoice  = atoi(argv[Narg]);
    remove_argument (Narg, &argc, argv);
    IDclip = TRUE;
  }
  FlagChoice = 0;
  FlagClip = FALSE;
  if ((Narg = get_argument (argc, argv, "-flag"))) {
    remove_argument (Narg, &argc, argv);
    FlagChoice  = atoi(argv[Narg]);
    remove_argument (Narg, &argc, argv);
    FlagClip = TRUE;
  }

  // XXX Add time restrictions

  LimExclude = FALSE;
  if ((Narg = get_argument (argc, argv, "-x"))) {
    remove_argument (Narg, &argc, argv);
    LimExclude = TRUE;
  }

  VERBOSE = FALSE;
  if ((Narg = get_argument (argc, argv, "-v"))) {
    remove_argument (Narg, &argc, argv);
    VERBOSE = TRUE;
  }

  if ((Narg = get_argument (argc, argv, "-m"))) {
    remove_argument (Narg, &argc, argv);
    Mr  = atof(argv[Narg]);
    remove_argument (Narg, &argc, argv);
    Mz = atof(argv[Narg]);
    Mr = Mr - Mz;
    remove_argument (Narg, &argc, argv);
  }

  if (argc != 2) {
    gprint (GP_ERR, "USAGE: pmeasure (-all) [-m M M] [-p photcode] [-ID ID] [-flag value] [-x]\n");
    return (FALSE);
  }
  graphmode.style = 2; /* set style to points */
  graphmode.size = -1; /* point size determined by Zvec */
  graphmode.etype = 0; /* no errorbars */

  Radius = MAX (fabs(graphmode.xmax), fabs(graphmode.ymax));

  /* load sky from correct table */
  sky = GetSkyTable ();
  skylist = SkyListByRadius (sky, -1, graphmode.coords.crval1, graphmode.coords.crval2, Radius);
  
  /* storage for plotting the points */
  Npts = 0;
  NPTS = 1000;
  ALLOCATE (Xvec, float, NPTS);
  ALLOCATE (Yvec, float, NPTS);
  ALLOCATE (Zvec, float, NPTS);

  for (j = 0; j < skylist[0].Nregions; j++) {
    catalog.filename = skylist[0].filename[j];
    catalog.catflags = LOAD_AVES | LOAD_MEAS;
    catalog.Nsecfilt = 0;

    // an error exit status here is a significant error
    if (!dvo_catalog_open (&catalog, NULL, VERBOSE, "r")) {
      fprintf (stderr, "ERROR: failure to open catalog file %s\n", catalog.filename);
      exit (2);
    }
    dvo_catalog_unlock (&catalog);

    /* project stars to screen display coords */
    for (i = 0; i < catalog.Naverage; i++) {
      if (IDclip && (catalog.average[i].code != IDchoice)) continue;
      while (catalog.average[i].R < Rmin) catalog.average[i].R += 360.0;
      while (catalog.average[i].R > Rmax) catalog.average[i].R -= 360.0;
      m = catalog.average[i].offset;
      for (k = 0; k < catalog.average[i].Nm; k++) {
	if (FlagClip) {
	  flags = catalog.measure[m+k].flags;
	  if (!(flags & FlagChoice)) continue;
	}
	if ((PhotcodeClip != -1) && (catalog.measure[m+k].photcode != PhotcodeClip)) continue;
	mag = PhotCat (&catalog.measure[m+k]);
	Zvec[Npts] = MIN (1.0, MAX (0.01, (mag - Mz) / Mr));
	if (LimExclude && (Zvec[Npts] > 0.99)) continue;
	if (Zvec[Npts] < 0.011) continue;
	R = catalog.average[i].R - catalog.measure[m+k].dR/3600.0;
	D = catalog.average[i].D - catalog.measure[m+k].dD/3600.0;
	// XXX drop this check
	if ((R < Rmin) || (R > Rmax) || (D < -90.0) || (D > 90.0)) {
	  char *date;
	  date = sec_to_date (catalog.measure[m+k].t);
	  gprint (GP_LOG, "out: %f, %f : %s : (%f, %f) + (%f, %f)\n", R, D, date, catalog.average[i].R, catalog.average[i].D, catalog.measure[m+k].dR/3600.0, catalog.measure[m+k].dD/3600.0);
	  free (date);
	}
	status = fRD_to_XY (&Xvec[Npts], &Yvec[Npts], R, D, &graphmode.coords);
	if (!status) {
	  fprintf (stderr, ".");
	  continue;
	}
	Npts ++;

	if (Npts == NPTS - 1) {
	    NPTS += 1000;
	    REALLOCATE (Xvec, float, NPTS);
	    REALLOCATE (Yvec, float, NPTS);
	    REALLOCATE (Zvec, float, NPTS);
	}
	if (Npts > NCHUNK) {
	    PlotVectorTriplet (kapa, Npts, Xvec, Yvec, Zvec, &graphmode);
	    Npts = 0;
	}
      }
    }
    dvo_catalog_free (&catalog);
  }
  if (Npts > 0) PlotVectorTriplet (kapa, Npts, Xvec, Yvec, Zvec, &graphmode);

  free (Xvec);
  free (Yvec);
  free (Zvec);

  return (TRUE);

}
  
