# include "dvo.h"

int bracket (double *list, double value, int Nlist, int mode);

int subpix (int argc, char **argv) {
  
  int i, j, I, Nlo, Nhi, Nentry, Nstars, Nregions, Nimage, Nmeasure;
  int *index, *entry;
  int Nmin, Nsub, NSUB, status;
  int TimeFormat;
  time_t Timage, TimeReference;
  double X, Y, Mabs, t;
  double Ra, Dec, Radius, Radius2, r, Rmin;
  double *RA, *DEC;
  char catdir[256], filename[256];
  
  Measure *measure;
  Image *image;
  RegionFile *regions;
  Catalog catalog;

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

  VarConfig ("CATDIR", "%s", catdir);
  GetTimeFormat (&TimeReference, &TimeFormat);

  if (argc != 4) {
    fprintf (stderr, "USAGE: subpix ra dec radius\n");
    return (FALSE);
  }
  if (!str_to_radec (&Ra, &Dec, argv[1], argv[2])) return (FALSE);
  if (Ra < 0) Ra += 360.0;
  if (Ra > 360.0) Ra -= 360.0;
  Radius = atof (argv[3]);

  /* load star nearest position */
  regions = find_regions (Ra, Dec, Radius, &Nregions);
  if (Nregions > 1) {
    fprintf (stderr, "warning, radius overlaps region boundary, not yet implemented\n");
  }

  /* lock, load, unlock catalog */
  sprintf (filename, "%s/%s", catdir, regions[0].name);
  catalog.filename = filename;
  switch (lock_catalog (&catalog, LCK_SOFT)) {
  case 2:
    unlock_catalog (&catalog);
  case 0:
    return (FALSE);
  }
  if (!load_catalog (&catalog, LOAD_AVES | LOAD_MEAS, TRUE)) {
    unlock_catalog (&catalog);
    return (FALSE);
  }
  unlock_catalog (&catalog);

  /* quick search of star list for Ra, Dec */
  Nstars = catalog.Naverage;
  ALLOCATE (RA, double, Nstars);
  ALLOCATE (DEC, double, Nstars);
  ALLOCATE (index, int, Nstars);
  for (i = 0; i < Nstars; i++) {
    RA[i] = catalog.average[i].R;
    DEC[i] = catalog.average[i].D;
    index[i] = i;
  }
  if (Nstars > 1) dsort_lists (DEC, RA, index, Nstars);

  /* bracket the DEC range of interest */
  Nlo = bracket (DEC, Dec - Radius, Nstars, FALSE);
  Nhi = bracket (DEC, Dec + Radius, Nstars, TRUE);
  ALLOCATE (entry, int, MAX (Nhi - Nlo, 1));
  Nentry = 0;

  /* find the list of stars */
  Radius2 = Radius*Radius;
  for (i = Nlo; i < Nhi; i++) {
    r = SQ(Dec - DEC[i]) + SQ(Ra - RA[i]);
    if (r < Radius2) {
      entry[Nentry] = i;
      Nentry ++;
    }
  }
  if (!Nentry) {
    fprintf (stderr, "no stars found\n");
    free (RA);
    free (DEC);
    free (entry);
    free (index);
    if (catalog.average != 0) free (catalog.average);
    if (catalog.measure != 0) free (catalog.measure);
    free (regions);
    return (TRUE);
  }

  /* find the closest star */
  Nmin = 0;
  Rmin = SQ(Dec - DEC[entry[0]]) + SQ(Ra - RA[entry[0]]);
  for (i = 1; i < Nentry; i++) {
    r = SQ(Dec - DEC[entry[i]]) + SQ(Ra - RA[entry[i]]);
    if (r < Rmin) {
      Rmin = r;
      Nmin = i;
    }
  }
  Nentry = index[entry[Nmin]];
  Ra = RA[entry[Nmin]];
  Dec = DEC[entry[Nmin]];
  fprintf (stderr, "finding subpix values for star @ %f %f\n", Ra, Dec);

  free (RA);
  free (DEC);
  free (entry);
  free (index);

  /* storage for the image references */
  Nsub = 0;
  NSUB = 100;
  ALLOCATE (index, int, NSUB);

  /* load all images, extract those touching Ra, Dec */
  image = LoadImages (&Nimage);
  for (i = 0; i < Nimage; i++) {
    status = RD_to_XY (&X, &Y, Ra, Dec, &image[i].coords);
    if (!status || (X < 0) || (X > image[i].NX) || (Y < 0) || (Y > image[i].NY)) continue;
    index[Nsub] = i;
    Nsub ++;
    if (Nsub == NSUB - 1) {
      NSUB += 100;
      REALLOCATE (index, int, NSUB);
    }
  }

  /* only print the entries for existing measurements of this star */ 
  measure = &catalog.measure[catalog.average[Nentry].offset];
  Nmeasure = catalog.average[Nentry].Nm;
  for (i = 0; i < Nsub; i++) {
    I = index[i];
    Timage = image[I].tzero;
    for (j = 0; j < Nmeasure; j++) {
      if (measure[j].t == Timage) { 
	Mabs = PhotCat (&measure[j]);
	RD_to_XY (&X, &Y, Ra, Dec, &image[I].coords);
	t = TimeValue (measure[j].t, TimeReference, TimeFormat);
	fprintf (GetOutfile(), "%f %6.3f %7.2f %7.2f %5.3f\n", t, Mabs, X, Y, 0.001*image[I].secz);
      } 
    }
  }

  if (catalog.average != 0) free (catalog.average);
  if (catalog.measure != 0) free (catalog.measure);
  free (image);
  free (index);
  free (regions);

  return (TRUE);
}

/* fast operation to find an entry just below (0) or above (1) value */
int bracket (double *list, double value, int Nlist, int mode) {

  int Nlo, Nhi, N;

  if (mode == 0) {
    Nlo = 0; Nhi = Nlist;
    while (Nhi - Nlo > 10) {
      N = 0.5*(Nlo + Nhi);
      if (list[N] < value) {
	Nlo = N;
      } else {
	Nhi = N + 1;
      }
    }
    return (Nlo);
  }
  if (mode == 1) {
    Nlo = 0; Nhi = Nlist;
    while (Nhi - Nlo > 10) {
      N = 0.5*(Nlo + Nhi);
      if (list[N] > value) {
	Nhi = N;
      } else {
	Nlo = N - 1;
      }
    }
    return (Nhi);
  }
  return (0);
}
