# include "dvoshell.h"

static double RAs = 0.0;
static double RAe = 0.0;
static double DECs = 0.0;
static double DECe = 0.0;

// define the sky region for which extractions are limited
int skyregion (int argc, char **argv) {
  
  if (argc == 1) {
    gprint (GP_ERR, "current skyregion: %f - %f : %f - %f\n", RAs, RAe, DECs, DECe);
    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);
}

int get_skyregion (double *Rs, double *Re, double *Ds, double *De) {

  *Rs = RAs;
  *Re = RAe;
  *Ds = DECs;
  *De = DECe;

  return TRUE;
}

int set_skyregion (double Rs, double Re, double Ds, double De) {

  RAs  = Rs;
  RAe  = Re;
  DECs = Ds;
  DECe = De;

  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);

  // 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

  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);
}

