# include "dvo.h"

enum {Nd, Nm, NC, NR, ND, Np, Nc, Nt, Nx, Nd1, Nd2, NVEC};

int calextract (int argc, char **argv) {
  
  Catalog catalog;
  Graphdata graphsky;
  RegionFile *regions;
  PhotCode *code[2];
  Vector **vec;

  int i, j, m, N, Nr, Nregions, mode[2];
  int RegionList, Ngraph, Nsec, NSTAR;
  char filename[256], catdir[256], RegionListFile[256];
  double Radius, M1, M2, dM2, color;

  /* these need to be freed in the end */
  regions = NULL;
  catalog.average = NULL; 
  catalog.secfilt = NULL;
  catalog.measure = NULL;

  /* load photcode information */
  if (!InitPhotcodes ()) return (FALSE);
  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);
  }

  /* command line arguments */
  if (!SetPhotSelections (&argc, argv)) goto usage;
  if (argc != 4) goto usage;
  if (strcmp (argv[2], "-")) goto usage;
  if (!GetPhotcodeInfo (argv[1], &code[0], &mode[0])) goto usage;
  if (!GetPhotcodeInfo (argv[3], &code[1], &mode[1])) goto usage;
  /* code.type must be PHOT_REF */

  /* one unique value per star */
  N = 0;
  NSTAR = 1;
  ALLOCATE (vec, Vector *, NVEC);
  if ((vec[Nd] 	= SelectVector ("cal:dmag",     ANYVECTOR, TRUE)) == NULL) goto escape;
  if ((vec[Nm] 	= SelectVector ("cal:mag",      ANYVECTOR, TRUE)) == NULL) goto escape;
  if ((vec[NC] 	= SelectVector ("cal:color",    ANYVECTOR, TRUE)) == NULL) goto escape;
  if ((vec[NR] 	= SelectVector ("cal:ra",       ANYVECTOR, TRUE)) == NULL) goto escape;
  if ((vec[ND] 	= SelectVector ("cal:dec",      ANYVECTOR, TRUE)) == NULL) goto escape;
  if ((vec[Np] 	= SelectVector ("cal:nphot",    ANYVECTOR, TRUE)) == NULL) goto escape;
  if ((vec[Nc] 	= SelectVector ("cal:ncode",    ANYVECTOR, TRUE)) == NULL) goto escape;
  if ((vec[Nt] 	= SelectVector ("cal:ncrit",    ANYVECTOR, TRUE)) == NULL) goto escape;
  if ((vec[Nx] 	= SelectVector ("cal:chisq",    ANYVECTOR, TRUE)) == NULL) goto escape;
  if ((vec[Nd1] = SelectVector ("cal:dm1",      ANYVECTOR, TRUE)) == NULL) goto escape;
  if ((vec[Nd2] = SelectVector ("cal:dm2",      ANYVECTOR, TRUE)) == NULL) goto escape;

  /* select region files */
  if (RegionList) {
    regions = region_list (RegionListFile, &Nregions);
  } else {
    Radius = MAX (fabs(graphsky.xmax), fabs(graphsky.ymax));
    regions = find_regions (graphsky.coords.crval1, graphsky.coords.crval2, Radius, &Nregions);
  }

  for (Nr = 0; Nr < Nregions; Nr++) {
    if (Nr && !(Nr % 500)) { fprintf (stderr, "."); }

    /* lock, load, unlock catalog */
    sprintf (filename, "%s/%s", catdir, regions[Nr].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);

    /* extract values, assign to vectors */
    for (i = 0; i < catalog.Naverage; i++) {
      m = catalog.average[i].offset;

      if (code[0][0].c1 && code[0][0].c2 && !PhotColor (&catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], code[0][0].c1, code[0][0].c2, &color)) continue;

      /* find data for filter 2 (PHOT_REF) */
      M2 = NO_MAG;
      for (j = 0; j < catalog.average[i].Nm; j++) {
	if (catalog.measure[m+j].source != code[1][0].code) continue;
	M2 = PhotCat  (&catalog.measure[m+j]);
	dM2 = 0.001*catalog.measure[m+j].dM;
      }	
      if (M2 == NO_MAG) continue;

      /* find data for filter 1 */
      M1 = ExtractAverages (code[0], mode[0], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], AVE_MAG);
      if (M1 == NO_MAG) continue;

      vec[Nd ][0].elements[N] = M1 - M2;
      vec[Nm ][0].elements[N] = M2;
      vec[NC ][0].elements[N] = color;
      vec[NR ][0].elements[N] = catalog.average[i].R;
      vec[ND ][0].elements[N] = catalog.average[i].D;
      vec[Nd1][0].elements[N] = ExtractAverages (code[0], mode[0], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], AVE_dMAG);
      vec[Nd2][0].elements[N] = dM2;
      vec[Nx ][0].elements[N] = ExtractAverages (code[0], mode[0], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], AVE_Xm);
      vec[Nc ][0].elements[N] = ExtractAverages (code[0], mode[0], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], AVE_NCODE);
      vec[Np ][0].elements[N] = ExtractAverages (code[0], mode[0], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], AVE_NPHOT);
      vec[Nt ][0].elements[N] = ExtractAverages (code[0], mode[0], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], AVE_NCRIT);
      N ++;
      if (N == NSTAR) {
	NSTAR += 100;
	for (i = 0; i < NVEC; i++) {
	  REALLOCATE (vec[i][0].elements, float, NSTAR);
	}
      }
    }
    if (catalog.average != (Average *) NULL) free (catalog.average);
    if (catalog.measure != (Measure *) NULL) free (catalog.measure);
    if (catalog.secfilt != (SecFilt *) NULL) free (catalog.secfilt);
    catalog.average = (Average *) NULL; 
    catalog.secfilt = (SecFilt *) NULL;
    catalog.measure = (Measure *) NULL;
  }
  if (regions != NULL) free (regions);  
  for (i = 0; i < NVEC; i++) {
    vec[i][0].Nelements = N;
  }
  return (TRUE);
  
usage:
  fprintf (stderr, "USAGE: calextract F - F\n");
  return (FALSE);

 escape:
  
  if (regions != NULL) free (regions);  
  for (i = 0; i < NVEC; i++) {
    DeleteVector (vec[i]);
  }
  free (vec);

  return (FALSE);
}
