# include "dvo.h"

int mextract (int argc, char **argv) {
  
  int i, j, k, m, N, N1, NPTS;
  int param, mode, Nregions, Ngraph, Nsec;
  char filename[128], catdir[256], *p;
  double *M1;

  PhotCode *code;
  Catalog catalog;
  RegionFile *regions;
  Graphdata graphsky;
  Vector *vec;

  regions = NULL;
  catalog.average = (Average *) NULL; 
  catalog.secfilt = (SecFilt *) NULL;
  catalog.measure = (Measure *) NULL;

  /* load photcode information */
  if (!InitPhotcodes ()) goto escape;
  Nsec = GetPhotcodeNsecfilt ();

  /* load data about plot windows */
  Ngraph = 0;
  if (!GetGraphData (&graphsky, NULL, &Ngraph)) goto escape;
 
  /* find CATDIR in config system */
  VarConfig ("CATDIR", "%s", catdir);

  /* specify catalog files
  RegionList = FALSE;
  if (N = get_argument (argc, argv, "-list")) {
    RegionList = TRUE;
    remove_argument (N, &argc, argv);
    strcpy (RegionListFile, argv[N]);
    remove_argument (N, &argc, argv);
  } */

  /* interpret command-line options */
  if (!SetPhotSelections (&argc, argv)) goto usage;
  if (argc != 3) goto usage;

  /* identify selection */
  code = NULL;
  mode = MAG_REL;
  param = GetMeasureParam (argv[2]);
  if (param == MEAS_ZERO) {
    if (!GetPhotcodeInfo (argv[2], &code, &mode)) {
      GetMeasureParam ("help");
      fprintf (stderr, "value may also be a valid photcode\n");
      fprintf (stderr, "photcodes or 'mag' may have optional magnitude type: mag,[Minst, Mcat, Msys, Mrel, Mcal]\n");
      goto escape;
    }
    param = MEAS_MAG;
    for (p = argv[2]; *p != 0; p++) { 
      if (*p == '.') *p = ':';
    }
  } 
  /* test PhotSelections: is photcode specified if needed? */
  if (!TestPhotSelections (&code)) {
    fprintf (stderr, "photcode selection rules violated\n");
    fprintf (stderr, "average and ensemble restrictions require a corresponding photcode\n");
    goto escape;
  }

  if ((vec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) goto escape;

  /* determine region-file names */
  if (!strcmp (argv[1], "all")) {
    float Radius;
    Radius = MAX (fabs(graphsky.xmax), fabs(graphsky.ymax));
    regions = find_regions (graphsky.coords.crval1, graphsky.coords.crval2, Radius, &Nregions);
    /* I don't have a tool to select images by the overlapping region */
    if (!SetImageSelection (param, &graphsky, TRUE)) goto escape;
  } else {
    Nregions = 1;
    ALLOCATE (regions, RegionFile, 1);
    strcpy (regions[0].name, argv[1]);
    if (!SetImageSelection (param, NULL, FALSE)) goto escape;
  }

  /* create storage vector */
  NPTS = 1000;
  REALLOCATE (vec[0].elements, float, NPTS);
  vec[0].Nelements = N = 0;

  for (i = 0; i < Nregions; i++) {
    /* lock, load, unlock catalog */
    sprintf (filename, "%s/%s", catdir, regions[i].name);
    catalog.filename = filename;
    switch (lock_catalog (&catalog, LCK_SOFT)) {
    case 2:
      unlock_catalog (&catalog);
    case 0:
      continue;
    }
    if (!load_catalog (&catalog, LOAD_AVES | LOAD_MEAS | LOAD_SECF, FALSE)) {
      unlock_catalog (&catalog);
      continue;
    }
    unlock_catalog (&catalog);

    for (j = 0; j < catalog.Naverage; j++) {
      M1 = NULL;
      m = catalog.average[j].offset;
      M1 = ExtractMeasures (code, mode, &catalog.average[j], &catalog.secfilt[j*Nsec], &catalog.measure[m], &N1, param);
      for (k = 0; k < N1; k++) {
	vec[0].elements[N] = M1[k];
	N++;
	if (N == NPTS - 1) {
	  NPTS += 2000;
	  REALLOCATE (vec[0].elements, float, NPTS);
	}
      }
      if (M1 != NULL) free (M1);
    }
    if (catalog.average != NULL) free (catalog.average);
    if (catalog.secfilt != NULL) free (catalog.secfilt);
    if (catalog.measure != NULL) free (catalog.measure);
    catalog.average = (Average *) NULL; 
    catalog.secfilt = (SecFilt *) NULL;
    catalog.measure = (Measure *) NULL;
  }
  vec[0].Nelements = N;
  REALLOCATE (vec[0].elements, float, MAX(1,N));

  if (regions != NULL) free (regions);
  FreeImageSelection ();
  return (TRUE);

usage:
  fprintf (stderr, "USAGE: mextract (from) (value) [options]\n");
  fprintf (stderr, "  from: cpt name or 'all'\n");
  fprintf (stderr, "  value: measure.parameter or photcode\n");
  return (FALSE);

escape:
  if (regions != NULL) free (regions);
  FreeImageSelection ();
  if (catalog.average != NULL) free (catalog.average);
  if (catalog.secfilt != NULL) free (catalog.secfilt);
  if (catalog.measure != NULL) free (catalog.measure);
  return (FALSE);
}
