# include "dvoshell.h"

// define the sky region for which extractions are limited
int skyregion (int argc, char **argv) {
  
  int N;

  // dvo_client should have 2 standard arguments: -hostID and -hostdir
  int SaveRegion = FALSE;
  if ((N = get_argument (argc, argv, "-save"))) {
    remove_argument (N, &argc, argv);
    SaveRegion = TRUE;
  }

  if (argc == 1) {
    double Rmin, Rmax, Dmin, Dmax;
    get_skyregion(&Rmin, &Rmax, &Dmin, &Dmax);

    if (SaveRegion) {
      set_variable ("Rmin", Rmin);
      set_variable ("Rmax", Rmax);
      set_variable ("Dmin", Dmin);
      set_variable ("Dmax", Dmax);
      return TRUE;
    } else {
      gprint (GP_ERR, "current skyregion: %f - %f : %f - %f\n", Rmin, Rmax, Dmin, Dmax);
      gprint (GP_ERR, "USAGE:  skyregion (min RA) (max RA) (min DEC) (max DEC)\n");
      return FALSE;
    }
  }

  if (argc != 5) {
    gprint (GP_ERR, "USAGE: skyregion (min RA) (max RA) (min DEC) (max DEC)\n");
    return (FALSE);
  }

  set_skyregion (atof(argv[1]), atof(argv[2]), atof(argv[3]), atof(argv[4]));

  return (TRUE);
}

/* find region which overlaps c at given depth (-1 : populated ) */
int SkyRegionByPoint_r (SkyTable *table, SkyList *list, int depth, double ra, double dec) {
  
  int i, Ns, Ne, No;
  SkyRegion *skyregion;

  list[0].Nregions = 0;
  list[0].ownElements = FALSE; // this list is only holding a view to the elements

  skyregion = table[0].regions;

  Ns = 0;
  Ne = 1;

  while (1) {
    No = -1;
    for (i = Ns; (i < Ne) && (i < table[0].Nregions); i++) {
      if (ra  < skyregion[i].Rmin) continue;
      if (ra  > skyregion[i].Rmax) continue;
      if (dec < skyregion[i].Dmin) continue;
      if (dec > skyregion[i].Dmax) continue;
      No = i;
      break;
    }
    if (No == -1) return (FALSE);
    if ((depth == -1) && (skyregion[No].table)) break;
    if (depth == skyregion[No].depth) break;
    if ((depth > skyregion[No].depth) && !skyregion[No].child) return (FALSE);

    /* need to check Ns, Ne, or guarantee valid range */
    Ns = skyregion[No].childS;
    Ne = skyregion[No].childE;
  }

  list[0].regions[0] = &skyregion[No];
  list[0].filename[0] = table[0].filename[No];
  list[0].Nregions = 1;
  return (TRUE);
}

SkyList *SelectRegionsByCoordVectors (Vector *RA, Vector *DEC) {
  
  int i, j, Npts, Nout, NOUT, *found;
  double ra, dec;
  SkyList *new, *list;
  SkyTable *sky;

  sky = GetSkyTable ();

  Npts = RA->Nelements;
  ALLOCATE (found, int, Npts);
  memset (found, 0, Npts*sizeof(int));

  ALLOCATE (new, SkyList, 1);
  ALLOCATE (new[0].regions,  SkyRegion *, 1);
  ALLOCATE (new[0].filename,  char *, 1);
  new[0].Nregions = 0;
  new[0].ownElements = FALSE; // this list is only holding a view to the elements
  strcpy (new[0].hosts, sky[0].hosts);

  // output list
  Nout = 0;
  NOUT = 100;
  ALLOCATE (list, SkyList, 1);
  ALLOCATE (list[0].regions,  SkyRegion *, NOUT);
  ALLOCATE (list[0].filename,  char *, NOUT);
  list[0].Nregions = 0;
  list[0].ownElements = FALSE; // this list is only holding a view to the elements
  strcpy (list[0].hosts, sky[0].hosts);

  for (i = 0; i < Npts; i++) {
    if (found[i]) continue;

    SkyRegionByPoint_r (sky, new, -1, RA->elements.Flt[i], DEC->elements.Flt[i]);
    if (new->Nregions == 0) continue;
    assert (new->Nregions == 1);
    // append new region to output list
    list[0].regions[Nout] = new[0].regions[0];
    list[0].filename[Nout] = new[0].filename[0];
    Nout ++;
    if (Nout >= NOUT) {
	NOUT += 100;
	REALLOCATE (list[0].regions, SkyRegion *, NOUT);
	REALLOCATE (list[0].filename, char *, NOUT);
    }
    found[i] = TRUE;

    // scan over the remaining points to find any that lie within this region
    // if we sorted the ra,dec inputs we could break out of this more quickly
    for (j = i + 1; j < Npts; j++) {
      ra = RA->elements.Flt[j];
      dec = DEC->elements.Flt[j];
      if (ra  < new[0].regions[0][0].Rmin) continue;
      if (ra  > new[0].regions[0][0].Rmax) continue;
      if (dec < new[0].regions[0][0].Dmin) continue;
      if (dec > new[0].regions[0][0].Dmax) continue;
      found[j] = TRUE;
    }
  }
  list[0].Nregions = Nout;

  free (found);
  SkyListFree (new);

  return (list);
}

