# include "dvoshell.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, 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;
  
  SkyTable *sky;
  SkyList *skylist;
  Measure *measure;
  Image *image;
  Catalog catalog;

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

  GetTimeFormat (&TimeReference, &TimeFormat);

  if (argc != 4) {
    gprint (GP_ERR, "USAGE: subpix ra dec radius\n");
    return (FALSE);
  }
  if (!ohana_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 */
  sky = GetSkyTable ();
  skylist = SkyListByRadius (sky, -1, Ra, Dec, Radius);
  if (skylist[0].Nregions > 1) {
    gprint (GP_ERR, "warning, radius overlaps region boundary, not yet implemented\n");
  }

  /* lock, load, unlock catalog */
  catalog.filename = skylist[0].filename[0];
  catalog.catflags = LOAD_AVES | LOAD_MEAS;
  catalog.Nsecfilt = 0;

  // an error exit status here is a significant error
  if (!dvo_catalog_open (&catalog, NULL, FALSE, "r")) {
      fprintf (stderr, "ERROR: failure to open catalog file %s\n", catalog.filename);
      exit (2);
  }
  dvo_catalog_unlock (&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) {
    gprint (GP_ERR, "no stars found\n");
    free (RA);
    free (DEC);
    free (entry);
    free (index);
    dvo_catalog_free (&catalog);
    SkyListFree (skylist);
    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]];
  gprint (GP_ERR, "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 */
  if ((image = LoadImages (&Nimage)) == NULL) return (FALSE);
  BuildChipMatch (image, Nimage);

  for (i = 0; i < Nimage; i++) {
    if (!FindMosaicForImage (image, Nimage, i)) continue;
    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];
    if (!FindMosaicForImage (image, Nimage, I)) continue;
    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);
	gprint (GP_LOG, "%f %6.3f %7.2f %7.2f %5.3f\n", t, Mabs, X, Y, image[I].secz);
      } 
    }
  }

  dvo_catalog_free (&catalog);
  free (image);
  free (index);
  SkyListFree (skylist);

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