# include "dvo.h"

enum {Nd, Nm1, Nm2, Nc, Ns, Nt, Nz, NR, ND, Nxc, Nyc, Nxm, Nym, NT, NP, Nd1, Nd2, NVEC};
int ConcatMeasures (Vector *vec, PhotCode *code, int mode, Average *average, SecFilt *secfilt, Measure *measure, int param, int Nin);

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

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

  /* 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 */

  /* 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;
  }


  /* returned vectors are dmag, mag, color, time, airmass, ra, dec, x, y, exptime */
  N = 0;
  Nstar = 0;
  NSTAR = 100;
  ALLOCATE (vec, Vector *, NVEC);
  if ((vec[Nd ] = SelectVector ("cal:dmag",     ANYVECTOR, TRUE)) == NULL) goto escape;
  if ((vec[Nm1] = SelectVector ("cal:mag1",     ANYVECTOR, TRUE)) == NULL) goto escape;
  if ((vec[Nm2] = SelectVector ("cal:mag2",     ANYVECTOR, TRUE)) == NULL) goto escape;
  if ((vec[Nc ] = SelectVector ("cal:color",    ANYVECTOR, TRUE)) == NULL) goto escape;
  if ((vec[Ns ] = SelectVector ("cal:star",     ANYVECTOR, TRUE)) == NULL) goto escape;
  if ((vec[Nt ] = SelectVector ("cal:time",     ANYVECTOR, TRUE)) == NULL) goto escape;
  if ((vec[Nz ] = SelectVector ("cal:airmass",  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[Nxc] = SelectVector ("cal:xccd",     ANYVECTOR, TRUE)) == NULL) goto escape;
  if ((vec[Nyc] = SelectVector ("cal:yccd",     ANYVECTOR, TRUE)) == NULL) goto escape;
  if ((vec[Nxm] = SelectVector ("cal:xmosaic",  ANYVECTOR, TRUE)) == NULL) goto escape;
  if ((vec[Nym] = SelectVector ("cal:ymosaic",  ANYVECTOR, TRUE)) == NULL) goto escape;
  if ((vec[NT ] = SelectVector ("cal:exptime",  ANYVECTOR, TRUE)) == NULL) goto escape;
  if ((vec[NP ] = SelectVector ("cal:photcode", 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;
  for (k = 0; k < NVEC; k++) {
    REALLOCATE (vec[k][0].elements, float, NSTAR);
    vec[k][0].Nelements = 0;
  }

  /* select region files */
  if (RegionList) {
    regions = region_list (RegionListFile, &Nregions);
    if (!SetImageSelection (MEAS_XMOSAIC, NULL, FALSE)) goto escape;
  } else {
    Radius = MAX (fabs(graphsky.xmax), fabs(graphsky.ymax));
    regions = find_regions (graphsky.coords.crval1, graphsky.coords.crval2, Radius, &Nregions);
    if (!SetImageSelection (MEAS_XMOSAIC, &graphsky, TRUE)) goto escape;
  }

  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;

      /* PRI/SEC must have data for color term */
      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 (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 = ExtractMeasures (code[0], mode[0], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], &N1, MEAS_MAG);
      if (N1 == 0) goto skip;

      /* extend storage vectors to take new data, if needed */
      if (N + N1 >= NSTAR) {
	NSTAR += N1 + 100;
	for (k = 0; k < NVEC; k++) {
	  REALLOCATE (vec[k][0].elements, float, NSTAR);
	}
      }

      ConcatMeasures (vec[Nt ], code[0], mode[0], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], N1, MEAS_TIME); 
      ConcatMeasures (vec[Nz ], code[0], mode[0], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], N1, MEAS_AIRMASS); 
      ConcatMeasures (vec[NT ], code[0], mode[0], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], N1, MEAS_EXPTIME); 
      ConcatMeasures (vec[NP ], code[0], mode[0], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], N1, MEAS_PHOTCODE); 
      ConcatMeasures (vec[Nd1], code[0], mode[0], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], N1, MEAS_dMAG); 
      ConcatMeasures (vec[Nxc], code[0], mode[0], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], N1, MEAS_XCCD); 
      ConcatMeasures (vec[Nyc], code[0], mode[0], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], N1, MEAS_YCCD); 
      ConcatMeasures (vec[Nxm], code[0], mode[0], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], N1, MEAS_XMOSAIC); 
      ConcatMeasures (vec[Nym], code[0], mode[0], &catalog.average[i], &catalog.secfilt[i*Nsec], &catalog.measure[m], N1, MEAS_YMOSAIC); 

      for (j = 0; j < N1; j++, N++) {
	vec[Nd ][0].elements[N] = M1[j] - M2;
	vec[Nm1][0].elements[N] = M1[j];
	vec[Nm2][0].elements[N] = M2;
	vec[Nd2][0].elements[N] = dM2;
	vec[Nc ][0].elements[N] = color;
	vec[Ns ][0].elements[N] = NStar;
	vec[NR ][0].elements[N] = catalog.average[i].R;
	vec[ND ][0].elements[N] = catalog.average[i].D;
      }
      NStar ++; 
    skip:
      if (M1 != NULL) free (M1);
    }
    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: dmags F - F : measure.param\n");
  return (FALSE);

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

int ConcatMeasures (Vector *vec, PhotCode *code, int mode, Average *average, SecFilt *secfilt, Measure *measure, int Nin, int param) {

  int i, Ns, N;
  double *value;

  value = ExtractMeasures (code, mode, average, secfilt, measure, &N, param); 
  if (N != Nin) {
    fprintf (stderr, "error!\n");
    return (FALSE);
  }

  Ns = vec[0].Nelements;
  for (i = 0; i < N; i++) {
    vec[0].elements[Ns+i] = value[i];
  }
  vec[0].Nelements = Ns + N;

  free (value);
  return (TRUE);
}
