# include "dvoshell.h"

enum {ZERO, RA, DEC, Xm, AIRMASS, MCAL, dMCAL, PHOTCODE, TIME, FWHM, EXPTIME, NSTAR, SKY, NCAL, FLAG};

int imextract (int argc, char **argv) {
  
  int i, j, Nimage, mode, N, Ngraph, PhotcodeSelect;
  int TimeSelect, RegionSelect, *subset, Nsubset, TimeFormat, FlagSelect, FlagValue;
  double x, y, ra, dec, t, trange;
  time_t tzero, TimeReference;

  PhotCode *code;
  Graphdata graphsky;
  Image *image;
  Vector *vec;

  Ngraph = 0;
  if (!GetGraphData (&graphsky, NULL, &Ngraph)) return (FALSE);
  if (!InitPhotcodes ()) return (FALSE);

  /* check for time-based selection */
  TimeSelect = FALSE;
  if ((N = get_argument (argc, argv, "-time"))) {
    remove_argument (N, &argc, argv);
    if (!str_to_time (argv[N], &tzero)) { 
      gprint (GP_ERR, "syntax error\n");
      return (FALSE);
    }
    remove_argument (N, &argc, argv);
    if (!str_to_dtime (argv[N], &trange)) { 
      gprint (GP_ERR, "syntax error\n");
      return (FALSE);
    }
    remove_argument (N, &argc, argv);
    TimeSelect = TRUE;
  }

  /* check for region-based selection */
  RegionSelect = FALSE;
  if ((N = get_argument (argc, argv, "-region"))) {
    remove_argument (N, &argc, argv);
    RegionSelect = TRUE;
  }

  /* check for region-based selection */
  FlagValue = 0;
  FlagSelect = FALSE;
  if ((N = get_argument (argc, argv, "-flag"))) {
    remove_argument (N, &argc, argv);
    FlagValue = atoi (argv[N]);
    remove_argument (N, &argc, argv);
    FlagSelect = TRUE;
  }

  /* check for photcode-based selection */
  code = NULL;
  PhotcodeSelect = FALSE;
  if ((N = get_argument (argc, argv, "-photcode"))) {
    PhotcodeSelect = TRUE;
    remove_argument (N, &argc, argv);
    if ((code = GetPhotcodebyName (argv[N])) == NULL) {
      gprint (GP_ERR, "ERROR: photcode %s not found in photcode table\n", argv[N]);
      return (FALSE);
    }
    if ((code[0].type != PHOT_PRI) && (code[0].type != PHOT_SEC) && (code[0].type != PHOT_DEP)) {
      gprint (GP_ERR, "photcode must be primary, secondary, or dependent code\n");
      return (FALSE);
    }
    remove_argument (N, &argc, argv);
  }

  if (argc != 2) {
    gprint (GP_ERR, "USAGE: imextract (value) [-region] [-time start range] [-photcode photcode]\n");
    return (FALSE);
  }
  
  /* identify selection */
  mode = ZERO;
  if (!strcasecmp (argv[1], "ra")) mode = RA;
  if (!strcasecmp (argv[1], "dec")) mode = DEC;
  if (!strcasecmp (argv[1], "Xm")) mode = Xm;
  if (!strcasecmp (argv[1], "airmass")) mode = AIRMASS;
  if (!strcasecmp (argv[1], "Mcal")) mode = MCAL;
  if (!strcasecmp (argv[1], "dMcal")) mode = dMCAL;
  if (!strcasecmp (argv[1], "photcode")) mode = PHOTCODE;
  if (!strcasecmp (argv[1], "time")) mode = TIME;
  if (!strcasecmp (argv[1], "FWHM")) mode = FWHM;
  if (!strcasecmp (argv[1], "exptime")) mode = EXPTIME;
  if (!strcasecmp (argv[1], "nstar")) mode = NSTAR;
  if (!strcasecmp (argv[1], "ncal")) mode = NCAL;
  if (!strcasecmp (argv[1], "sky")) mode = SKY;
  if (!strcasecmp (argv[1], "flag")) mode = FLAG;
  if (mode == ZERO) {
    gprint (GP_ERR, "value may be one of the following:\n");
    gprint (GP_ERR, " ra dec airmass Mcal dMcal Xm photcode time fwhm exptime nstar ncal sky flag\n");
    return (FALSE);
  }
  if ((vec = SelectVector (argv[1], ANYVECTOR, TRUE)) == NULL) return (FALSE);

  if ((image = LoadImages (&Nimage)) == NULL) return (FALSE);
  image_subset (image, Nimage, &subset, &Nsubset, &graphsky, RegionSelect, tzero, trange, TimeSelect);
  if ((mode == RA) || (mode == DEC)) BuildChipMatch (image, Nimage);

  GetTimeFormat (&TimeReference, &TimeFormat);

  /* create storage vector */
  REALLOCATE (vec[0].elements, float, Nimage);
  vec[0].Nelements = Nimage;
  
  N = 0;
  /* assign vector values */
  for (i = 0; i < Nsubset; i++) {
    j = subset[i];
    if (PhotcodeSelect) {
      if (code[0].type == PHOT_DEP) {
	if (code[0].code != image[j].source) continue;
      } else {
	if (code[0].code != GetPhotcodeEquivCodebyCode (image[j].source)) continue;
      }
    }
    if (FlagSelect && (FlagValue != image[j].code)) continue;
    switch (mode) {
    case (RA):
      if (!FindMosaicForImage (image, Nimage, j)) continue;
      x = 0.5*image[j].NX;
      y = 0.5*image[j].NY;
      XY_to_RD (&ra, &dec, x, y, &image[j].coords);
      vec[0].elements[N] = ra;
      break;
    case (DEC):
      if (!FindMosaicForImage (image, Nimage, j)) continue;
      x = 0.5*image[j].NX;
      y = 0.5*image[j].NY;
      XY_to_RD (&ra, &dec, x, y, &image[j].coords);
      vec[0].elements[N] = dec;
      break;
    case (Xm):
      vec[0].elements[N] = pow(10.0, 0.01*image[j].Xm);
      break;
    case (AIRMASS):
      vec[0].elements[N] = image[j].secz_PS;
      break;
    case (MCAL):
      vec[0].elements[N] = image[j].Mcal_PS;
      break;
    case (dMCAL):
      vec[0].elements[N] = image[j].dMcal_PS;
      break;
    case (PHOTCODE):
      vec[0].elements[N] = image[j].source;
      break;
    case (TIME):
      t = image[j].tzero + 0.5*image[j].NY * image[j].trate / 10000;
      vec[0].elements[N] = TimeValue (t, TimeReference, TimeFormat);
      break;
    case (FWHM):
      vec[0].elements[N] = image[j].fwhm_x / 25.0;
      break;
    case (EXPTIME):
      vec[0].elements[N] = image[j].exptime;
      break;
    case (NSTAR):
      vec[0].elements[N] = image[j].nstar;
      break;
    case (NCAL):
      vec[0].elements[N] = image[j].Mxxxx;
      break;
    case (SKY):
      vec[0].elements[N] = image[j].Myyyy + 0x8000;
      break;
    case (FLAG):
      vec[0].elements[N] = image[j].code;
      break;
    }
    N++;
  }
  
  vec[0].Nelements = N;
  REALLOCATE (vec[0].elements, float, N);

  free (subset);
  free (image);
  return (TRUE);
  
}
  
