# include "dvoshell.h"

// need to upgrade this to support multiple fields and WHERE clauses

enum {ZERO, 
      RA, 
      DEC, 
      Xm, 
      AIRMASS, 
      MCAL, 
      dMCAL, 
      PHOTCODE, 
      TIME, 
      FWHM, 
      EXPTIME, 
      NSTAR, 
      NCAL, 
      SKY, 
      FLAG, 
      NX_PIX, 
      NY_PIX, 
      THETA, 
      SKEW, 
      SCALE, 
      DSCALE, 
      IMAGE_ID,
      X_LL_CHIP,
      X_LR_CHIP,
      X_UL_CHIP,
      X_UR_CHIP,
      Y_LL_CHIP,
      Y_LR_CHIP,
      Y_UL_CHIP,
      Y_UR_CHIP,
      X_LL_FP,
      X_LR_FP,
      X_UL_FP,
      X_UR_FP,
      Y_LL_FP,
      Y_LR_FP,
      Y_UL_FP,
      Y_UR_FP,
};

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

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

  if (!InitPhotcodes ()) return (FALSE);

  // parse skyregion options
  if ((selection = SetRegionSelection (&argc, argv)) == NULL) {
    gprint (GP_ERR, "invalid sky region selection\n");
    return FALSE;
  }

  /* check for time-based selection */
  TimeSelect = FALSE;
  if ((N = get_argument (argc, argv, "-time"))) {
    remove_argument (N, &argc, argv);
    if (!ohana_str_to_time (argv[N], &tzero)) { 
      gprint (GP_ERR, "syntax error\n");
      return (FALSE);
    }
    remove_argument (N, &argc, argv);
    if (!ohana_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 */
  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 (!strcasecmp (argv[1], "imageID"  )) mode = IMAGE_ID;
  if (!strcasecmp (argv[1], "X_LL_CHIP")) mode = X_LL_CHIP;
  if (!strcasecmp (argv[1], "X_LR_CHIP")) mode = X_LR_CHIP;
  if (!strcasecmp (argv[1], "X_UL_CHIP")) mode = X_UL_CHIP;
  if (!strcasecmp (argv[1], "X_UR_CHIP")) mode = X_UR_CHIP;
  if (!strcasecmp (argv[1], "Y_LL_CHIP")) mode = Y_LL_CHIP;
  if (!strcasecmp (argv[1], "Y_LR_CHIP")) mode = Y_LR_CHIP;
  if (!strcasecmp (argv[1], "Y_UL_CHIP")) mode = Y_UL_CHIP;
  if (!strcasecmp (argv[1], "Y_UR_CHIP")) mode = Y_UR_CHIP;
  if (!strcasecmp (argv[1], "X_LL_FP"  )) mode = X_LL_FP;
  if (!strcasecmp (argv[1], "X_LR_FP"  )) mode = X_LR_FP;
  if (!strcasecmp (argv[1], "X_UL_FP"  )) mode = X_UL_FP;
  if (!strcasecmp (argv[1], "X_UR_FP"  )) mode = X_UR_FP;
  if (!strcasecmp (argv[1], "Y_LL_FP"  )) mode = Y_LL_FP;
  if (!strcasecmp (argv[1], "Y_LR_FP"  )) mode = Y_LR_FP;
  if (!strcasecmp (argv[1], "Y_UL_FP"  )) mode = Y_UL_FP;
  if (!strcasecmp (argv[1], "Y_UR_FP"  )) mode = Y_UR_FP;
  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, selection, 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 IMAGE_ID:
	vec[0].elements[N] = image[j].imageID;
	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; }

    // reference pixel extractions
    case X_LL_CHIP:
    case Y_LL_CHIP:
    case Y_LR_CHIP:
    case X_UL_CHIP:
      vec[0].elements[N] = 0.0;
      break;
    case X_LR_CHIP:
    case X_UR_CHIP:
      vec[0].elements[N] = image[j].NX;
      break;
    case Y_UL_CHIP:
    case Y_UR_CHIP:
      vec[0].elements[N] = image[j].NX;
      break;

    case X_LL_FP:
    case Y_LL_FP:
      XY_to_LM (&x, &y, 0.0, 0.0, &image[j].coords);
      vec[0].elements[N] = (mode == X_LL_FP) ? x : y;
      break;
    case X_LR_FP:
    case Y_LR_FP:
      XY_to_LM (&x, &y, image[j].NX, 0.0, &image[j].coords);
      vec[0].elements[N] = (mode == X_LR_FP) ? x : y;
      break;
    case X_UL_FP:
    case Y_UL_FP:
      XY_to_LM (&x, &y, 0.0, image[j].NY, &image[j].coords);
      vec[0].elements[N] = (mode == X_UL_FP) ? x : y;
      break;
    case X_UR_FP:
    case Y_UR_FP:
      XY_to_LM (&x, &y, image[j].NX, image[j].NY, &image[j].coords);
      vec[0].elements[N] = (mode == X_UR_FP) ? x : y;
      break;
    }
    N++;
  }
  
  vec[0].Nelements = N;
  REALLOCATE (vec[0].elements, float, N);

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