# include "dvo.h"

int ddmags (int argc, char **argv) {
  
  Catalog catalog;
  Graphdata graphsky;
  RegionFile *regions;
  PhotCode *code[4];
  Vector *xvec, *yvec;

  char filename[256], catdir[256];
  double Radius;
  double *M1, *M2;
  int i, m, k, Npts, NPTS;
  int N1, N2, i1, i2, mode[4];
  int Ngraph, Nsec, Nregions, UseAverages;

  /* defaults */
  regions = NULL;
  catalog.average = NULL; 
  catalog.secfilt = NULL;
  catalog.measure = NULL;

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

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

  /* interpret command-line options */
  if (!SetPhotSelections (&argc, argv)) goto usage;
  if (argc != 8) goto usage;
  if (strcmp (argv[2], "-")) goto usage;
  if (strcmp (argv[4], ":")) goto usage;
  if (strcmp (argv[6], "-")) goto usage;

  /* interpret command-line options */
  if (!GetPhotcodeInfo (argv[1], &code[0], &mode[0])) return (FALSE);
  if (!GetPhotcodeInfo (argv[3], &code[1], &mode[1])) return (FALSE);
  if (!GetPhotcodeInfo (argv[5], &code[2], &mode[2])) return (FALSE);
  if (!GetPhotcodeInfo (argv[7], &code[3], &mode[3])) return (FALSE);
  if ((mode[2] == MAG_AVE) || (mode[2] == MAG_REF)) UseAverages = TRUE;

  /* find catalog files which overlap this region */
  Radius = MAX (fabs(graphsky.xmax), fabs(graphsky.ymax));
  regions = find_regions (graphsky.coords.crval1, graphsky.coords.crval2, Radius, &Nregions);

  /* init vectors to save data */
  Npts = 0;
  NPTS = 1;
  if ((xvec = SelectVector ("xv", ANYVECTOR, TRUE)) == NULL) goto escape;
  if ((yvec = SelectVector ("yv", ANYVECTOR, TRUE)) == NULL) goto escape;

  /* loop over regions, extract data for each region */
  for (k = 0; k < Nregions; k++) {
    /* lock, load, unlock catalog */
    sprintf (filename, "%s/%s", catdir, regions[k].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);

    /* get correct mags, convert to X,Y */
    for (i = 0; i < catalog.Naverage; i++) {
      M1 = M2 = NULL;
      m = catalog.average[i].offset;

      M1 = ExtractDMag (&code[0], &mode[0], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], &N1);
      if (N1 == 0) goto skip;

      M2 = ExtractDMag (&code[2], &mode[2], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], &N2);
      if (N2 == 0) goto skip;

      for (i1 = 0; i1 < N1; i1++) {
	for (i2 = 0; i2 < N2; i2++) {
	  xvec[0].elements[Npts] = M1[i1];
	  yvec[0].elements[Npts] = M2[i2];
	  Npts++;
	  if (Npts == NPTS) {
	    NPTS += 2000;
	    REALLOCATE (xvec[0].elements, float, NPTS);
	    REALLOCATE (yvec[0].elements, float, NPTS);
	  }
	}
      }
    skip:
      if (M1 != NULL) free (M1);
      if (M2 != NULL) free (M2);
    }
    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;
  }
  if (regions != NULL) free (regions);
  xvec[0].Nelements = yvec[0].Nelements = Npts;
  return (TRUE);

usage:
  fprintf (stderr, "USAGE: ddmags F - F : measure.param\n");
  return (FALSE);

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

