# include "dvoshell.h"

enum {ZERO, RA, DEC, Xm, AIRMASS, MCAL, dMCAL, PHOTCODE, TIME, FWHM, EXPTIME, NSTAR, SKY, NCAL, FLAG, NX_PIX, NY_PIX, THETA, SKEW, SCALE, DSCALE};

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_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 (!strcasecmp (argv[1], "NX")) mode = NX_PIX;
  if (!strcasecmp (argv[1], "NY")) mode = NY_PIX;
  if (!strcasecmp (argv[1], "theta")) mode = THETA;
  if (!strcasecmp (argv[1], "skew")) mode = SKEW;
  if (!strcasecmp (argv[1], "scale")) mode = SCALE;
  if (!strcasecmp (argv[1], "dscale")) mode = DSCALE;
  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].photcode) continue;
      } else {
	if (code[0].code != GetPhotcodeEquivCodebyCode (image[j].photcode)) 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;
	break;
      case MCAL:
	vec[0].elements[N] = image[j].Mcal;
	break;
      case dMCAL:
	vec[0].elements[N] = image[j].dMcal;
	break;
      case PHOTCODE:
	vec[0].elements[N] = image[j].photcode;
	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;
      case NX_PIX:
	vec[0].elements[N] = image[j].NX;
	break;
      case NY_PIX:
	vec[0].elements[N] = image[j].NY;
	break;
      case THETA: {
	double theta1, theta2, s1, s2;
	s1 = SIGN(image[j].coords.pc1_1);
	s2 = SIGN(image[j].coords.pc2_2);
	theta1 = DEG_RAD*atan2 (+s1*image[j].coords.pc1_2, s1*image[j].coords.pc1_1);
	theta2 = DEG_RAD*atan2 (-s2*image[j].coords.pc2_1, s2*image[j].coords.pc2_2);
	vec[0].elements[N] = 0.5*(theta1+theta2);
	break; }
      case SKEW: {
	double theta1, theta2, s1, s2;
	s1 = SIGN(image[j].coords.pc1_1);
	s2 = SIGN(image[j].coords.pc2_2);
	theta1 = DEG_RAD*atan2 (+s1*image[j].coords.pc1_2, s1*image[j].coords.pc1_1);
	theta2 = DEG_RAD*atan2 (-s2*image[j].coords.pc2_1, s2*image[j].coords.pc2_2);
	vec[0].elements[N] = (theta1-theta2);
	break; }
      case SCALE: {
	double scale1, scale2;
	scale1 = fabs(image[j].coords.cdelt1);
	scale2 = fabs(image[j].coords.cdelt2);
	vec[0].elements[N] = 0.5*(scale1+scale2);
	break; }
      case DSCALE: {
	double scale1, scale2;
	scale1 = fabs(image[j].coords.cdelt1);
	scale2 = fabs(image[j].coords.cdelt2);
	vec[0].elements[N] = (scale1-scale2);
	break; }
	  }
    N++;
  }
  
  vec[0].Nelements = N;
  REALLOCATE (vec[0].elements, float, N);

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