# include "dvo.h"

/* match code to measure  */
# define TESTCODE(C,M) \
  switch (C[0].type) { \
  case PHOT_DEP: \
  case PHOT_REF: \
    if (C[0].code != M.source) continue; \
    break; \
  case PHOT_PRI: \
  case PHOT_SEC: \
    if (C[0].code != GetPhotcodeEquivCodebyCode (M.source)) continue; \
    break; \
  default: \
    break; \
} 

/* exclusions based on measure.params  */
# define TESTMEASURE(M) \
  if (TimeSelect && (M.t < tzero)) continue; \
  if (TimeSelect && (M.t > tend)) continue; \
  if (ErrSelect  && (M.dM > ErrValue)) continue; \
  if (TypeSelect && (TypeValue != GetMeasureTypeCode (&M))) continue; \
  if (iMagSelect && (PhotInst (&M) < iMagMin)) continue; \
  if (FlagSelect && (M.flags != FlagValue)) continue;

# define SETMAG(MOUT,MEAS,MODE) \
  MOUT = NO_MAG; \
  if (MODE == MAG_INST) MOUT = PhotInst (&MEAS);  \
  if (MODE == MAG_CAT)  MOUT = PhotCat  (&MEAS); \
  if (MODE == MAG_SYS)  MOUT = PhotSys  (&MEAS, average, secfilt); \
  if (MODE == MAG_REL)  MOUT = PhotRel  (&MEAS, average, secfilt); \
  if (MODE == MAG_CAL)  MOUT = PhotCal  (&MEAS, average, secfilt, measure, GetPhotcodeEquivbyCode (MEAS.source)); \
  if (MagSelect && (MOUT > MagMax)) continue; \
  if (MagSelect && (MOUT < MagMin)) continue;

/* selection criteria */
/* selections based on Measure quantities */
static int TimeSelect;
static time_t tzero, tend;
static int MagSelect;
static double MagMax, MagMin;
static int TypeSelect, TypeValue;
static int ErrSelect, ErrValue;
static int iMagSelect;
static double iMagMin;
static int FlagSelect, FlagValue;

/* applied to Average quantities */
static int PhotcodeSelect;
static PhotCode *PhotcodeValue;

/* selections based on Average quantities */
static int ChiSelect;
static float ChiLimit;

/* selections based on ensemble quantities */
static int NphotSelect, NphotSign, NphotValue;
static int NcodeSelect, NcodeSign, NcodeValue;
static int FWHMSelect, FWHMsign;
static double FWHMvalue, FWHMfrac;

/* time concepts */
static time_t TimeReference;
static int TimeFormat;

/* db image table */
static Image *image = NULL;
static int *subset = NULL;
static int Nimage = 0;
static int Nsubset = 0;
static Coords mosaic;

int GetPhotcodeInfo (char *string, PhotCode **Code, int *Mode) {

  PhotCode *code;
  int mode, status;
  char *p, *tmpstring;

  /* save local copy */
  tmpstring = strcreate (string);

  /* check for code:mode in photcode name */
  mode = MAG_NONE;
  p = strchr (tmpstring, ':');
  if (p != NULL) {
    mode = GetMagMode (p + 1);
    if (mode == MAG_NONE) {
      fprintf (stderr, "syntax error in magnitude mode\n");
      free (tmpstring);
      return (FALSE);
    }
    *p = 0;
  }

  /* how do we handle this elsewhere? */
  if (!strcasecmp (tmpstring, "mag")) {
    /* need to validate mode */
    *Mode = mode;
    *Code = NULL;
    free (tmpstring);
    return (TRUE);
  }

  code = GetPhotcodebyName (tmpstring);
  if (code == NULL) {
    fprintf (stderr, "photcode not found in photcode table\n");
    free (tmpstring);
    return (FALSE);
  }

  /* test allowable cases and/or set default values */
  status = FALSE;
  if (code[0].type == PHOT_DEP) {
    if (mode == MAG_NONE) mode = MAG_REL;
    if (mode == MAG_INST) status = TRUE;
    if (mode == MAG_CAT)  status = TRUE;
    if (mode == MAG_SYS)  status = TRUE;
    if (mode == MAG_REL)  status = TRUE;
    if (mode == MAG_CAL)  status = TRUE;
  }  
  if ((code[0].type == PHOT_PRI) || (code[0].type == PHOT_SEC)) {
    if (mode == MAG_NONE) mode  = MAG_AVE;
    if (mode == MAG_INST) status = TRUE;
    if (mode == MAG_CAT)  status = TRUE;
    if (mode == MAG_SYS)  status = TRUE;
    if (mode == MAG_REL)  status = TRUE;
    if (mode == MAG_CAL)  status = TRUE;
    if (mode == MAG_AVE)  status = TRUE;
    if (mode == MAG_REF)  status = TRUE;
  }  
  if (code[0].type == PHOT_ALT) {
    if (mode == MAG_NONE) mode  = MAG_AVE;
    if (mode == MAG_AVE)  status = TRUE;
    if (mode == MAG_REF)  status = TRUE;
  }

  if (code[0].type == PHOT_REF) {
    if (mode == MAG_NONE) mode  = MAG_CAT;
    if (mode == MAG_CAT)  status = TRUE;
  }

  if (!status) {
    fprintf (stderr, "mismatch in photcode and magmode\n");
    free (tmpstring);
    return (FALSE);
  }
  *Code = code;
  *Mode = mode;
  free (tmpstring);
  return (TRUE);
}
 
int GetMagMode (char *string) {

  int mode;

  mode = MAG_NONE;
  if (!strcasecmp (string, "inst"))   mode = MAG_INST;
  if (!strcasecmp (string, "cat"))    mode = MAG_CAT;
  if (!strcasecmp (string, "sys"))    mode = MAG_SYS;
  if (!strcasecmp (string, "rel"))    mode = MAG_REL;
  if (!strcasecmp (string, "cal"))    mode = MAG_CAL;
  if (!strcasecmp (string, "ave"))    mode = MAG_AVE;
  if (!strcasecmp (string, "ref"))    mode = MAG_REF;
  return (mode);
}

int GetMeasureParam (char *parname) {

  int param;

  param = MEAS_ZERO;
  if (!strcasecmp (parname, "ra"))   	 param = MEAS_RA;
  if (!strcasecmp (parname, "dec"))  	 param = MEAS_DEC;
  if (!strcasecmp (parname, "mag")) 	 param = MEAS_MAG;
  if (!strcasecmp (parname, "dmag")) 	 param = MEAS_dMAG;
  if (!strcasecmp (parname, "airmass"))  param = MEAS_AIRMASS;
  if (!strcasecmp (parname, "exptime"))  param = MEAS_EXPTIME;
  if (!strcasecmp (parname, "photcode")) param = MEAS_PHOTCODE;
  if (!strcasecmp (parname, "time"))     param = MEAS_TIME;
  if (!strcasecmp (parname, "dR")) 	 param = MEAS_dR;
  if (!strcasecmp (parname, "dD")) 	 param = MEAS_dD;
  if (!strcasecmp (parname, "fwhm"))   	 param = MEAS_FWHM;
  if (!strcasecmp (parname, "dophot")) 	 param = MEAS_DOPHOT;
  if (!strcasecmp (parname, "FLAGS"))    param = MEAS_FLAGS;
  if (!strcasecmp (parname, "XCCD"))   	 param = MEAS_XCCD;
  if (!strcasecmp (parname, "YCCD"))   	 param = MEAS_YCCD;
  if (!strcasecmp (parname, "XMOSAIC"))  param = MEAS_XMOSAIC;
  if (!strcasecmp (parname, "YMOSAIC"))  param = MEAS_YMOSAIC;
  if (!strcasecmp (parname, "help")) {
    fprintf (stderr, "value may be one of the following:\n");
    fprintf (stderr, " ra dR dec dD mag dmag Mrel Mcal photcode time fwhm dophot xccd yccd xmosaic ymosaic flags\n");
  }
  return (param);
}
  
int TestPhotSelections (PhotCode **code) {

  int NeedPhotcode;

  /* if i've supplied a photcode (code != NULL), i'm not allowed to restrict it */
  if ((code[0] != NULL) && PhotcodeSelect) return (FALSE);

  /* if I have an average or ensemble restriction, I must have a photcode */
  if (ChiSelect || NphotSelect || NcodeSelect || FWHMSelect) {
    if (code[0] != NULL) return (TRUE);
    if (PhotcodeSelect) {
      code[0] = PhotcodeValue;
      return (TRUE);
    }
    return (FALSE);
  }

  /* if I have a photcode restriction, apply it to the code */
  if (PhotcodeSelect) {
    code[0] = PhotcodeValue;
    return (TRUE);
  }
  return (TRUE);
}

int GetAverageParam (char *parname) {

  int param;

  param = AVE_ZERO;
  if (!strcasecmp (parname, "ra"))    param = AVE_RA;
  if (!strcasecmp (parname, "dec"))   param = AVE_DEC;
  if (!strcasecmp (parname, "dmag"))  param = AVE_dMAG;
  if (!strcasecmp (parname, "Nmeas")) param = AVE_NMEAS;
  if (!strcasecmp (parname, "Nmiss")) param = AVE_NMISS;
  if (!strcasecmp (parname, "Xp"))    param = AVE_Xp;
  if (!strcasecmp (parname, "Xm"))    param = AVE_Xm;
  if (!strcasecmp (parname, "flag"))  param = AVE_FLAG;
  if (!strcasecmp (parname, "type"))  param = AVE_TYPE;
  if (!strcasecmp (parname, "Nphot")) param = AVE_NPHOT;
  if (!strcasecmp (parname, "Ncode")) param = AVE_NCODE;
  if (!strcasecmp (parname, "help")) {
    fprintf (stderr, "value may be one of the following:\n");
    fprintf (stderr, " ra dec dmag Nmeas Nmiss Xm Xp Nphot Ncode flag type\n");
  }
  return (param);
}

/* load images based on parameters and region, etc */
int SetImageSelection (int mode, Graphdata *graphsky, int RegionSelect) {

  image = NULL;
  subset = NULL;
  
  switch (mode) {
  case MEAS_XCCD:
  case MEAS_YCCD:
    break;
  case MEAS_XMOSAIC: 
  case MEAS_YMOSAIC:
    /* mosaic defines a frame with 0,0 at the mosaic center, and 1 arcsec / pixel */
    mosaic.crpix1 = mosaic.crpix2 = 0.0;
    mosaic.cdelt1 = mosaic.cdelt2 = 1.0 / 3600;
    mosaic.pc1_1  = mosaic.pc2_2  = 1.0;
    mosaic.pc1_2  = mosaic.pc2_1  = 0.0;
    mosaic.Npolyterms = 0;
    strcpy (mosaic.ctype, "RA---SIN");
    break;
  default:
    return (TRUE);
  }

  image = LoadImages (&Nimage);
  if (Nimage == 0) {
    return (FALSE);
  }
  image_subset (image, Nimage, &subset, &Nsubset, graphsky, RegionSelect, tzero, (double) tend - tzero, TimeSelect);

  sort_image_subset (image, subset, Nsubset);
  return (TRUE);
}

/* free loaded images */
void FreeImageSelection () {
  if (image != NULL) free (image);
  if (subset != NULL) free (subset);
  return;
}

/* (re)load photcodes from photcode table */
int InitPhotcodes () {

  double ZERO_POINT;
  char PhotCodeFile[256];

  if (VarConfig ("ZERO_PT", "%lf", &ZERO_POINT) == (char *) NULL) return (FALSE);
  SetZeroPoint (ZERO_POINT);

  if (VarConfig ("PHOTCODE_FILE", "%s", PhotCodeFile) == (char *) NULL) return (FALSE);
  if (!LoadPhotcodes (PhotCodeFile)) {
    fprintf (stderr, "error loading photcodes\n");
    return (FALSE);
  }
  return (TRUE);
}

/* remove standard photometry filtering options, set selections */
/* not all functions respect all selections... */
int SetPhotSelections (int *argc, char **argv) {

  int N;
  double trange;

  if (N = get_argument (*argc, argv, "-phothelp")) {
    fprintf (stderr, "optional photometry selection criteria:\n");
    fprintf (stderr, " -magrange min max\n");
    fprintf (stderr, " -imaglim min\n");
    fprintf (stderr, " -flag value\n");
    fprintf (stderr, " -chisq value\n");
    fprintf (stderr, " -photcode code\n");
    fprintf (stderr, " -time start range\n");
    fprintf (stderr, " -errorlim value\n");
    fprintf (stderr, " -type type\n");
    fprintf (stderr, " -nmeas [+/-]N\n");
    fprintf (stderr, " -fwhm [+/-]fraction\n");
    return (FALSE);
  }

  /* select based on measured mag (MEASURE ONLY) */
  MagSelect = FALSE;
  if (N = get_argument (*argc, argv, "-magrange")) {
    MagSelect = TRUE;
    remove_argument (N, argc, argv);
    MagMin = atof (argv[N]);
    remove_argument (N, argc, argv);
    MagMax = atof (argv[N]);
    remove_argument (N, argc, argv);
  }

  /* select based on instrument mag (MEASURE ONLY) */
  iMagSelect = FALSE;
  if (N = get_argument (*argc, argv, "-imaglim")) {
    iMagSelect = TRUE;
    remove_argument (N, argc, argv);
    iMagMin = atof (argv[N]);
    remove_argument (N, argc, argv);
  }

  /* select on value of flag (MEASURE ONLY) */
  FlagSelect = FALSE;
  if (N = get_argument (*argc, argv, "-flag")) {
    FlagSelect = TRUE;
    remove_argument (N, argc, argv);
    FlagValue = atof (argv[N]);
    remove_argument (N, argc, argv);
  }

  /* select on value of flag (MEASURE ONLY) */
  ChiSelect = FALSE;
  if (N = get_argument (*argc, argv, "-chisq")) {
    ChiSelect = TRUE;
    remove_argument (N, argc, argv);
    ChiLimit = atof (argv[N]);
    remove_argument (N, argc, argv);
  }

  /* select on value of photcode */
  PhotcodeSelect = FALSE;
  if (N = get_argument (*argc, argv, "-photcode")) {
    PhotcodeSelect = TRUE;
    remove_argument (N, argc, argv);
    PhotcodeValue = GetPhotcodebyName (argv[N]);
    if (PhotcodeValue == NULL) {
      fprintf (stderr, "photcode not found in photcode table\n");
      PhotcodeSelect = FALSE;
    }
    remove_argument (N, argc, argv);
  }

  /* selection on basis of time range (MEASURE only) */
  TimeSelect = FALSE;
  if (N = get_argument (*argc, argv, "-time")) {
    remove_argument (N, argc, argv);
    if (!str_to_time (argv[N], &tzero)) { 
      fprintf (stderr, "syntax error\n");
      return (FALSE);
    }
    remove_argument (N, argc, argv);
    if (!str_to_dtime (argv[N], &trange)) { 
      fprintf (stderr, "syntax error\n");
      return (FALSE);
    }
    remove_argument (N, argc, argv);
    TimeSelect = TRUE;
    if (trange < 0) {
      trange = fabs (trange);
      tend = tzero;
      tzero -= trange;
    } else {
      tend = tzero + trange;
    }
  }

  /* select by error (on measure or average) */
  ErrSelect = FALSE;
  if (N = get_argument (*argc, argv, "-errorlim")) {
    remove_argument (N, argc, argv);
    ErrValue = 1000*atof (argv[N]);
    remove_argument (N, argc, argv);
    ErrSelect = TRUE;
  }

  /* select on measurement type: 1,2,3 (AVERAGE ONLY) */
  TypeSelect = FALSE;
  if (N = get_argument (*argc, argv, "-type")) {
    remove_argument (N, argc, argv);
    TypeValue = atoi (argv[N]);
    remove_argument (N, argc, argv);
    TypeSelect = TRUE;
  }

  /* select by number of measurements (AVERAGE ONLY) */
  NphotSelect = FALSE;
  if (N = get_argument (*argc, argv, "-nphot")) {
    remove_argument (N, argc, argv);
    NphotValue = abs (atoi (argv[N]));
    NphotSign = 0;
    if (argv[N][0] == '-') NphotSign = -1;
    if (argv[N][0] == '+') NphotSign = +1;
    remove_argument (N, argc, argv);
    NphotSelect = TRUE;
  }

  /* select by number of measurements (AVERAGE ONLY) */
  NcodeSelect = FALSE;
  if (N = get_argument (*argc, argv, "-ncode")) {
    remove_argument (N, argc, argv);
    NcodeValue = abs (atoi (argv[N]));
    NcodeSign = 0;
    if (argv[N][0] == '-') NcodeSign = -1;
    if (argv[N][0] == '+') NcodeSign = +1;
    remove_argument (N, argc, argv);
    NcodeSelect = TRUE;
  }

  /* -fwhm value frac (AVERAGE ONLY) */
  FWHMSelect = FALSE;
  if (N = get_argument (*argc, argv, "-fwhm")) {
    remove_argument (N, argc, argv);
    FWHMvalue = abs (atof (argv[N]));
    FWHMsign = 0;
    if (argv[N][0] == '-') FWHMsign = -1;
    if (argv[N][0] == '+') FWHMsign = +1;
    remove_argument (N, argc, argv);
    FWHMSelect = TRUE;
    remove_argument (N, argc, argv);
    FWHMfrac = atof (argv[N]);
  }

  GetTimeFormat (&TimeReference, &TimeFormat);
  return (TRUE);
}

/* extract a list of measure parameters from the specified average entry based on the pre-set selections */
double *ExtractMeasures (PhotCode *code, int mode, Average *average, SecFilt *secfilt, Measure *measure, int *nlist, int param) {

  int i, k, m, Nlist, NLIST, np, found, Nm;
  time_t t;
  double ra, dec, x, y, value, M, fwhm;
  double *list;
  
  *nlist = 0; 
  Nlist = 0;
  NLIST = MAX (1, average[0].Nm);
  ALLOCATE (list, double, NLIST);

  /* chisq, fwhm, Nphot, Ncode */
  if (!TestAverage (code, average, secfilt, measure)) return (list);

  /* look for measures */
  for (i = 0; i < average[0].Nm; i++) {
    if (code != NULL) {
      TESTCODE (code, measure[i]);
    }
    /* exclusions based on measure.params  */
    TESTMEASURE (measure[i]);

    /* need to have the mag value to test MagSelect */ 
    SETMAG (M, measure[i], mode);

    /* assign value */
    switch (param) {
    case MEAS_MAG: /* magnitudes are already determined above */
      value = M;
      break;
    case MEAS_RA: /* OK */
      value = average[0].R - measure[i].dR / 360000.0;
      break;
    case MEAS_DEC: /* OK */
      value = average[0].D - measure[i].dD / 360000.0;
      break;
    case MEAS_dMAG: /* OK */
      value = 0.001*measure[i].dM;
      break;
    case MEAS_AIRMASS: /* OK */
      value = 0.001*measure[i].airmass;
      break;
    case MEAS_EXPTIME: /* OK */
      value = pow (10.0, measure[i].dt * 0.0004);
      break;
    case MEAS_PHOTCODE: /* OK */
      value = measure[i].source;
      break;
    case MEAS_TIME: /* OK */
      value = TimeValue (measure[i].t, TimeReference, TimeFormat);
      break;
    case MEAS_dR: /* OK */
      value = 0.01*measure[i].dR;
      break;
    case MEAS_dD: /* OK */
      value = 0.01*measure[i].dD;
      break;
    case MEAS_FWHM: /* OK */
      value = 0.01*measure[i].FWx;
      break;
    case MEAS_DOPHOT: /* OK */
      value = measure[i].dophot;
      break;
    case MEAS_FLAGS: /* ? */
      value = measure[i].flags;
      break;
    case MEAS_XCCD: /* OK */
      ra  = average[0].R - measure[i].dR / 360000.0;
      dec = average[0].D - measure[i].dD / 360000.0;
      m = match_image_subset (image, subset, Nsubset, measure[i].t, measure[i].source);
      if (m == -1) {
	value = 0.0;
	break;
      }
      RD_to_XY (&x, &y, ra, dec, &image[m].coords);
      value = x;
      break;
    case MEAS_YCCD: /* OK */
      ra  = average[0].R - measure[i].dR / 360000.0;
      dec = average[0].D - measure[i].dD / 360000.0;
      m = match_image_subset (image, subset, Nsubset, measure[i].t, measure[i].source);
      if (m == -1) {
	value = 0.0;
	break;
      }
      RD_to_XY (&x, &y, ra, dec, &image[m].coords);
      value = y;
      break;
    case MEAS_XMOSAIC: /* OK */
      ra  = average[0].R - measure[i].dR / 360000.0;
      dec = average[0].D - measure[i].dD / 360000.0;
      m = match_image_subset (image, subset, Nsubset, measure[i].t, measure[i].source);
      if (m == -1) {
	value = 0.0;
	break;
      }
      mosaic.crval1 = image[m].coords.crval1;
      mosaic.crval2 = image[m].coords.crval2;
      RD_to_XY (&x, &y, ra, dec, &mosaic);
      value = x;
      break;
    case MEAS_YMOSAIC: /* OK */
      ra  = average[0].R - measure[i].dR / 360000.0;
      dec = average[0].D - measure[i].dD / 360000.0;
      m = match_image_subset (image, subset, Nsubset, measure[i].t, measure[i].source);
      if (m == -1) {
	value = 0.0;
	break;
      }
      mosaic.crval1 = image[m].coords.crval1;
      mosaic.crval2 = image[m].coords.crval2;
      RD_to_XY (&x, &y, ra, dec, &mosaic);
      value = y;
      break;
    }
    list[Nlist] = value;
    Nlist ++;
  }
  *nlist = Nlist;
  return (list);
}

/* return average.param based on the selection */
double ExtractAverages (PhotCode *code, int mode, Average *average, SecFilt *secfilt, Measure *measure, int param) {

  int i, k, Nm, Type;
  double fwhm, value, dM;

  value = NO_MAG;

  /* need valid code if param is: */
  switch (param) {
  case AVE_MAG:
  case AVE_dMAG:
  case AVE_Xm:
  case AVE_TYPE:
  case AVE_NCODE:
  case AVE_NPHOT: 
    if (code == NULL) return (NO_MAG);
    break;
  default:
    break;
  }

  /* chisq, fwhm, Nphot, Ncode */
  if (!TestAverage (code, average, secfilt, measure)) return (NO_MAG);

  /* for ErrSelect, check average errors */
  if (ErrSelect) {
    dM = iPhotdM (code, average, secfilt);
    if (dM > ErrValue) return (NO_MAG);
  }
  
  /* for TypeSelect, check on TypeCode for this object */
  if (TypeSelect) {
    Type = DetermineTypeCode (average, measure, code[0].code);
    if (Type != TypeValue) return (NO_MAG);
  }

  /* assign vector values */
  switch (param) {
  case AVE_RA:
    value = average[0].R;
    break;
  case AVE_DEC:
    value = average[0].D;
    break;
  case AVE_NMEAS:
    value = average[0].Nm;
    break;
  case AVE_NMISS:
    value = average[0].Nn;
    break;
  case AVE_Xp:
    value = 0.01*average[0].Xp;
    break;
  case AVE_FLAG:
    value = average[0].code;
    break;
  case AVE_MAG:
    switch (mode) {
    case MAG_AVE:
      value = PhotAve  (code, average, secfilt);
      break;
    case MAG_REF:
      value = PhotRef  (code, average, secfilt, measure);
      break;
    }
    break;
  case AVE_dMAG:
    value = PhotdM (code, average, secfilt);
    break;
  case AVE_Xm:
    value = PhotXm (code, average, secfilt);
    break;
  case AVE_TYPE:
    value = DetermineTypeCode (average, measure, code[0].code);
    break;
  case AVE_NCODE:
    value = 0;
    for (i = 0; i < average[0].Nm; i++) {
      if (code[0].code != GetPhotcodeEquivCodebyCode (measure[i].source)) continue;
      value ++;
    }
    break;
  case AVE_NPHOT:
    value = 0;
    for (i = 0; i < average[0].Nm; i++) {
      if (code[0].code != GetPhotcodeEquivCodebyCode (measure[i].source)) continue;
      if (measure[i].flags & (ID_MEAS_POOR | ID_MEAS_SKIP)) continue;
      value ++;
    }
    break;
  case AVE_NCRIT:
    value = 0;
    for (i = 0; i < average[0].Nm; i++) {
      if (code[0].code != GetPhotcodeEquivCodebyCode (measure[i].source)) continue;
      if (ErrSelect && (measure[i].dM > ErrValue)) continue;
      if (FlagSelect && (measure[i].flags != FlagValue)) continue;
      if (TypeSelect && (TypeValue != GetMeasureTypeCode (&measure[i]))) continue;
      if (iMagSelect && (PhotInst (&measure[i]) < iMagMin)) continue;
      value ++;
    }
    break;
  }
  return (value);
}  

/* determine the representative dophot type for this photcode */
int DetermineTypeCode (Average *average, Measure *measure, int code) {

  int k, N, Nt[3];
  
  Nt[0] = Nt[1] = Nt[2] = 0;
  for (k = 0; k < average[0].Nm; k++) {
    if (code != GetPhotcodeEquivCodebyCode (measure[k].source)) continue;
    N = GetMeasureTypeCode (&measure[k]);
    Nt[N] ++;
  }
  if (Nt[0]) return (0);
  if (Nt[1]) return (1);
  if (Nt[2]) return (2);
  return (3);
}

int GetMeasureTypeCode (Measure *measure) {
  switch (measure[0].dophot) {
  case 0:
  case 1:
  case 2:
    return (0);
    break;
  case 3:
  case 4:
  case 5:
  case 7:
  case 9:
    return (1);
    break;
  case 10:
  default:
    return (2);
  }
  return (2);
}  

int Quality (Measure *measure, int IsDophot) {

  return (TRUE);
  
  if (IsDophot) {
    
    if (measure[0].dophot == 4) return (FALSE);
    return (TRUE);
  
  } else {
    
    if (0.01 * measure[0].FWx < 3.0) return (FALSE);
    if (0.01 * measure[0].FWx > 10.0) return (FALSE);

    return (TRUE);

  }
}

int TestAverage (PhotCode *code, Average *average, SecFilt *secfilt, Measure *measure) {

  int i, Nm;
  double value, fwhm;

  /* exclusions based on average.params  */
  if (ChiSelect) {
    /* some average.params need a valid photcode for reference */
    if (code == NULL) return (FALSE);
    value = PhotXm (code, average, secfilt);
    if (value == -1) return (FALSE);
    if (value > ChiLimit) return (FALSE);
  }
  
  /* pass objects with more than FWHMfrac points with FWHM above / below FWHMvalue */ 
  if (FWHMSelect) {
    Nm = 0;
    for (i = 0; i < average[0].Nm; i++) {
      fwhm = measure[i].FWx / 100.0;
      switch (FWHMsign) {
      case 0:
	if (fwhm == FWHMvalue) break;
	continue;
      case +1:
	if (fwhm >= FWHMvalue) break;
	continue;
      case -1:
	if (fwhm <= FWHMvalue) break;
	continue;
      }
      Nm++;
    }
    if (average[0].Nm * FWHMfrac > Nm) return (FALSE);
  }

  /* for NcodeSelect, count Nmeas for appropriate photcode */
  if (NcodeSelect) {
    Nm = 0;
    for (i = 0; i < average[0].Nm; i++) {
      TESTCODE (code, measure[i]);
      Nm++;
    }
    switch (NcodeSign) {
    case 0:
      if (Nm == NcodeValue) break;
      return (FALSE);
    case 1:
      if (Nm >= NcodeValue) break;
      return (FALSE);
    case -1:
      if (Nm <= NcodeValue) break;
      return (FALSE);
    default:
      return (FALSE);
    }
  }

  /* for NphotSelect, count Nmeas for appropriate photcode */
  if (NphotSelect) {
    Nm = 0;
    for (i = 0; i < average[0].Nm; i++) {
      TESTCODE (code, measure[i]);
      if (measure[i].flags && ID_MEAS_SKIP) continue;
      Nm++;
    }
    switch (NphotSign) {
    case 0:
      if (Nm == NphotValue) break;
      return (FALSE);
    case 1:
      if (Nm >= NphotValue) break;
      return (FALSE);
    case -1:
      if (Nm <= NphotValue) break;
      return (FALSE);
    default:
      return (FALSE);
    }
  }
  return (TRUE);
}

/* for this function, we don't need to call PhotRel, etc, but we
   do need to multiply by 0.001:
   average[].M is stored as 1000*mag where mag is PhotAbs
   measure[].M for PHOT_REL is the same 
*/ 

/* send in:
   Nphot - photcode number
   Tphot - photcode type
   Ns    - secfilt entry (-1 for PRI)
   &catalog.average[i], 
   &catalog.measure[catalog.average[i].offset], 
   &catalog.secfilt[i*Nsec] 
*/


double *ExtractMagnitudes (PhotCode *code, int mode, Average *average, SecFilt *secfilt, Measure *measure, int *n) {
  
  double *M, mag;
  int N;

  if ((mode == MAG_AVE) || (mode == MAG_REF)) {
    ALLOCATE (M, double, 1);
    mag = ExtractAverages (code, mode, average, secfilt, measure, AVE_MAG);
    if (mag == NO_MAG) {
      N = 0;
    } else {
      N = 1;
      M[0] = mag;
    }
  } else {
    M = ExtractMeasures (code, mode, average, secfilt, measure, &N, MEAS_MAG);
  }
  
  *n = N;
  return (M);
}

double *ExtractDMag (PhotCode **code, int *mode, Average *average, SecFilt *secfilt, Measure *measure, int *nlist) {

  int i, j, A1, A2, N1, N2, Nlist, NLIST;
  double *M1, *M2, *list;

  /* check for special case of measure-measure */
  A1 = ((mode[0] == MAG_AVE) || (mode[0] == MAG_REF));
  A2 = ((mode[1] == MAG_AVE) || (mode[1] == MAG_REF));
  if (!A1 && !A2) {
    list = ExtractMeasuresDMag (code, mode, average, secfilt, measure, nlist);
    return (list);
  }

  *nlist = 0; 
  Nlist = 0;
  NLIST = MAX (1, average[0].Nm*average[0].Nm);
  ALLOCATE (list, double, NLIST);
  M1 = M2 = NULL;

  /* one of the two is an average, must do independently */
  M1 = ExtractMagnitudes (code[0], mode[0], average, secfilt, measure, &N1);
  if (N1 == 0) goto skip;
  
  M2 = ExtractMagnitudes (code[1], mode[1], average, secfilt, measure, &N2);
  if (N2 == 0) goto skip;

  /* average mags can return NO_MAG : skip them? */
  for (i = 0; i < N1; i++) {
    for (j = 0; j < N2; j++) {
      list[Nlist] = M1[i] - M2[j];
      Nlist ++;
    }
  }
  *nlist = Nlist;
  return (list);

 skip: 
  if (M1 != NULL) free (M1);
  if (M2 != NULL) free (M2);
  M1 = M2 = NULL;
  return (list);
}
  
/* extract a list of delta-measure-mags from the specified average entry based on the pre-set selections */
double *ExtractMeasuresDMag (PhotCode **code, int *mode, Average *average, SecFilt *secfilt, Measure *measure, int *nlist) {

  int i, j, Nlist, NLIST;
  double *list, M1, M2, value;
  
  *nlist = 0; 
  Nlist = 0;
  NLIST = MAX (1, average[0].Nm*average[0].Nm);
  ALLOCATE (list, double, NLIST);

  /* must have two code values */
  if (code == NULL) return (list);
  if (code[0] == NULL) return (list);
  if (code[1] == NULL) return (list);

  /* exclude based on chisq for both codes  */
  if (ChiSelect) {
    value = PhotXm (code[0], average, secfilt);
    if (value == -1) return (list);
    if (value > ChiLimit) return (list);
    value = PhotXm (code[1], average, secfilt);
    if (value == -1) return (list);
    if (value > ChiLimit) return (list);
  }

  /* loop twice over all measures */
  for (i = 0; i < average[0].Nm; i++) {
    TESTCODE (code[0], measure[i]);
    TESTMEASURE (measure[i]);
    SETMAG(M1, measure[i], mode[0]);
    for (j = 0; j < average[0].Nm; j++) {
      if (i == j) continue;
      TESTCODE (code[1], measure[j]);
      TESTMEASURE (measure[j]);
      SETMAG(M2, measure[j], mode[1]);
      list[Nlist] = M1 - M2;
      Nlist ++;
    }
  }
  *nlist = Nlist;
  return (list);
}

double *ExtractByDMag (PhotCode **code, int *mode, Average *average, SecFilt *secfilt, Measure *measure, int *nlist, int param) {

  int A1, A2, N1;
  double *list;

  /* check for special case of measure-measure */
  A1 = ((mode[0] == MAG_AVE) || (mode[0] == MAG_REF));
  A2 = ((mode[1] == MAG_AVE) || (mode[1] == MAG_REF));
  if (!A1 && !A2) {
    list = ExtractMeasuresByDMag (code, mode, 1, average, secfilt, measure, nlist, param);
    return (list);
  }

  if (A1) {
    list = ExtractMeasures (code[1], mode[1], average, secfilt, measure, &N1, param);
  } else {
    list = ExtractMeasures (code[0], mode[0], average, secfilt, measure, &N1, param);
  }
  *nlist = N1;
  return (list);
}
  
/* extract a list of delta-measure-mags from the specified average entry based on the pre-set selections */
double *ExtractMeasuresByDMag (PhotCode **code, int *mode, int use_first, Average *average, SecFilt *secfilt, Measure *measure, int *nlist, int param) {

  int i, j, k, m, n, Nlist, NLIST;
  time_t t;
  int np, found;
  double *list, M1, M2, value;
  double ra, dec, x, y;
  
  *nlist = 0; 
  Nlist = 0;
  NLIST = MAX (1, average[0].Nm*average[0].Nm);
  ALLOCATE (list, double, NLIST);

  /* must have two code values */
  if (code == NULL) return (list);
  if (code[0] == NULL) return (list);
  if (code[1] == NULL) return (list);

  /* chisq, fwhm, Nphot, Ncode */
  if (!TestAverage (code[0], average, secfilt, measure)) return (list);
  if (!TestAverage (code[1], average, secfilt, measure)) return (list);

  /* loop twice over all measures */
  for (i = 0; i < average[0].Nm; i++) {
    TESTCODE (code[0], measure[i]);
    TESTMEASURE (measure[i]);
    SETMAG(M1, measure[i], mode[0]);
    for (j = 0; j < average[0].Nm; j++) {
      if (i == j) continue;
      TESTCODE (code[1], measure[j]);
      TESTMEASURE (measure[j]);
      SETMAG(M2, measure[j], mode[1]);
      if (use_first) 
	n = i;
      else 
	n = j;
      /* assign value */
      switch (param) {
      case MEAS_MAG: /* magnitudes are already determined above */
	value = use_first ? M1 : M2;
	break;
      case MEAS_RA: /* OK */
	value = average[0].R - measure[n].dR / 360000.0;
	break;
      case MEAS_DEC: /* OK */
	value = average[0].D - measure[n].dD / 360000.0;
	break;
      case MEAS_dMAG: /* OK */
	value = 0.001*measure[n].dM;
	break;
      case MEAS_AIRMASS: /* OK */
	value = 0.001*measure[n].airmass;
	break;
      case MEAS_EXPTIME: /* OK */
	value = pow (10.0, measure[n].dt * 0.0004);
	break;
      case MEAS_PHOTCODE: /* OK */
	value = measure[n].source;
	break;
      case MEAS_TIME: /* OK */
	value = TimeValue (measure[n].t, TimeReference, TimeFormat);
	break;
      case MEAS_dR: /* OK */
	value = 0.01*measure[n].dR;
	break;
      case MEAS_dD: /* OK */
	value = 0.01*measure[n].dD;
	break;
      case MEAS_FWHM: /* OK */
	value = 0.01*measure[n].FWx;
	break;
      case MEAS_DOPHOT: /* OK */
	value = measure[n].dophot;
	break;
      case MEAS_FLAGS: /* ? */
	value = measure[n].flags;
	break;
    case MEAS_XCCD: /* OK */
      ra  = average[0].R - measure[n].dR / 360000.0;
      dec = average[0].D - measure[n].dD / 360000.0;
      m = match_image_subset (image, subset, Nsubset, measure[n].t, measure[n].source);
      if (m == -1) {
	value = 0.0;
	break;
      }
      RD_to_XY (&x, &y, ra, dec, &image[m].coords);
      value = x;
      break;
    case MEAS_YCCD: /* OK */
      ra  = average[0].R - measure[n].dR / 360000.0;
      dec = average[0].D - measure[n].dD / 360000.0;
      m = match_image_subset (image, subset, Nsubset, measure[n].t, measure[n].source);
      if (m == -1) {
	value = 0.0;
	break;
      }
      RD_to_XY (&x, &y, ra, dec, &image[m].coords);
      value = y;
      break;
    case MEAS_XMOSAIC: /* OK */
      ra  = average[0].R - measure[n].dR / 360000.0;
      dec = average[0].D - measure[n].dD / 360000.0;
      m = match_image_subset (image, subset, Nsubset, measure[n].t, measure[n].source);
      if (m == -1) {
	value = 0.0;
	break;
      }
      mosaic.crval1 = image[m].coords.crval1;
      mosaic.crval2 = image[m].coords.crval2;
      RD_to_XY (&x, &y, ra, dec, &mosaic);
      value = x;
      break;
    case MEAS_YMOSAIC: /* OK */
      ra  = average[0].R - measure[n].dR / 360000.0;
      dec = average[0].D - measure[n].dD / 360000.0;
      m = match_image_subset (image, subset, Nsubset, measure[n].t, measure[n].source);
      if (m == -1) {
	value = 0.0;
	break;
      }
      mosaic.crval1 = image[m].coords.crval1;
      mosaic.crval2 = image[m].coords.crval2;
      RD_to_XY (&x, &y, ra, dec, &mosaic);
      value = y;
      break;
      }
      list[Nlist] = value;
      Nlist ++;
    }
  }
  *nlist = Nlist;
  return (list);
}

