# include "addstar.h"

int find_matches_closest (SkyRegion *region, Stars *stars, int NstarsIn, Catalog *catalog, AddstarClientOptions options) {

  int i, j, n, N, J, Jmin, status, Nstars;
  double RADIUS, RADIUS2, Rmin;
  double *X1, *Y1, *X2, *Y2;
  double dX, dY, dR;
  int *N1, *N2,  *next_meas, *next_miss;
  int Nave, NAVE, Nmeas, NMEAS, Nmiss, NMISS, Nmatch;
  Coords tcoords;
  int Nsecfilt, Nsec;

  /* photcode data - must by of type DEP; options.photcode is equiv photcode for all input
     images this function requires incoming stars to have the same photcode.equiv value.  if
     this value is not a valid photcode (ie, 0), then no modification is made to the average 
     magnitudes (Nsec will be -1) */

  Nsecfilt = GetPhotcodeNsecfilt ();
  Nsec     = GetPhotcodeNsec (options.photcode);

  /** allocate local arrays (stars) **/
  ALLOCATE (X1, double, NstarsIn);
  ALLOCATE (Y1, double, NstarsIn);
  ALLOCATE (N1, int,   NstarsIn);

  /** allocate local arrays (catalog) **/
  NAVE = Nave = catalog[0].Naverage;
  ALLOCATE (X2, double, NAVE);
  ALLOCATE (Y2, double, NAVE);
  ALLOCATE (N2, int,   NAVE);
  ALLOCATE (catalog[0].found, int, NAVE);
  /* for secfilt j and star i, secfilt[i*Nsecfilt+j] */

  /* internal counters */
  Nmatch = 0;
  NMEAS = Nmeas = catalog[0].Nmeasure;
  NMISS = Nmiss = catalog[0].Nmissing;

  /* project onto rectilinear grid with 1 arcsec pixels. the choice of ZEA projection has the
   * advantage that every point in R,D has a mapping to a unique X,Y.  However, note that not all
   * possible X,Y points map back to R,D and the local plate scale changes substantially far from
   * the projection pole.  a better mapping might be ARC, not yet implemented (see
   * coordops.update.c).  We use the center of the region (catalog) for crval1,2. 
   */
  tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax);
  if (region[0].Dmax < 90) {
    tcoords.crval2 = 0.5*(region[0].Dmin + region[0].Dmax);
  } else {
    tcoords.crval2 = 90.0;
  }
  tcoords.crpix1 = 0;
  tcoords.crpix2 = 0;
  tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
  tcoords.pc1_1 = tcoords.pc2_2 = 1.0;
  tcoords.pc1_2 = tcoords.pc2_1 = 0.0;
  tcoords.Npolyterms = 1;
  strcpy (tcoords.ctype, "RA---ARC");

  /* build spatial index (RA sort) referencing input array sequence */
  Nstars = 0;
  for (i = 0; i < NstarsIn; i++) {
    status = RD_to_XY (&X1[Nstars], &Y1[Nstars], stars[i].R, stars[i].D, &tcoords);
    if (!status) continue;
    N1[Nstars] = i;
    Nstars ++;
  }
  if (Nstars < 1) {
    if (VERBOSE) fprintf (stderr, "skipping %s, no overlapping stars\n", catalog[0].filename);
    free (catalog[0].found);
    free (X1);
    free (Y1);
    free (N1);
    free (X2);
    free (Y2);
    free (N2);
    return (0);
  }
  if (Nstars > 1) sort_coords_index (X1, Y1, N1, Nstars);

  /* build spatial index (RA sort) */
  for (i = 0; i < Nave; i++) {
    RD_to_XY (&X2[i], &Y2[i], catalog[0].average[i].R, catalog[0].average[i].D, &tcoords);
    N2[i] = i;
    catalog[0].found[N2[i]] = -1;
  }
  if (Nave > 1) sort_coords_index (X2, Y2, N2, Nave);

  /* set up pointers for linked list of measure, missing */
  if (catalog[0].sorted) {
    next_meas = init_measure_links (catalog[0].average, Nave, catalog[0].measure, Nmeas);
  } else {
    next_meas = build_measure_links (catalog[0].average, Nave, catalog[0].measure, Nmeas);
  }    
  next_miss = init_missing_links (catalog[0].average, Nave, catalog[0].missing, Nmiss);
  /* missing MUST be written 'sorted', or not at all */

  /* choose a radius for matches */
  RADIUS = options.radius; /* provided by config */
  RADIUS2 = RADIUS*RADIUS;

  /** find matched stars **/
  for (i = j = 0; (i < Nstars) && (j < Nave); ) {
    if (!finite(X1[i]) || !finite(Y1[i])) { 
      i++; 
      continue;
    }
    if (!finite(X2[j]) || !finite(Y2[j])) { 
      j++; 
      continue;
    }
    
    /* negative dX: j is too large */
    dX = X1[i] - X2[j];
    if (dX <= -2*RADIUS) {
      i++;
      continue;
    }
    /* positive dX, i is too large */
    if (dX >= 2*RADIUS) {
      j++;
      continue;
    }

    Jmin = -1;
    Rmin = RADIUS2;
    for (J = j; (dX > -2*RADIUS) && (J < Nave); J++) {
      /* find closest match for this detection */
      dX = X1[i] - X2[J];
      dY = Y1[i] - Y2[J];
      dR = dX*dX + dY*dY;
      if (dR > RADIUS2) continue;
      if (dR < Rmin) {
	Rmin = dR;
	Jmin  = J;
      }
    }

    /* no match, try next detection */ 
    if (Jmin == -1) {
      i++;
      continue;
    }

    /*** a match is found, add to average, measure ***/

    /* make sure there is space for next entry */
    if (Nmeas >= NMEAS) {
      NMEAS = Nmeas + 1000;
      REALLOCATE (next_meas, int, NMEAS);
      REALLOCATE (catalog[0].measure, Measure, NMEAS);
    }

    Nmatch ++;
    n = N2[Jmin];
    N = N1[i];

    /* add to end of measurement list */
    add_meas_link (&catalog[0].average[n], next_meas, Nmeas, NMEAS);
	
    /** add measurements for this star **/
    /** dR,dD now represent arcsec **/
    catalog[0].measure[Nmeas].dR       = 3600.0*(catalog[0].average[n].R - stars[N].R);
    if (catalog[0].measure[Nmeas].dR > +180.0*3600.0) {
      // average on high end of boundary, move star up
      stars[N].R += 360.0;
      catalog[0].measure[Nmeas].dR = 3600.0*(catalog[0].average[n].R - stars[N].R);
    }
    if (catalog[0].measure[Nmeas].dR < -180.0*3600.0) {
      // average on low end of boundary, move star down
      stars[N].R -= 360.0;
      catalog[0].measure[Nmeas].dR = 3600.0*(catalog[0].average[n].R - stars[N].R);
    }
    if (fabs(catalog[0].measure[Nmeas].dR) > 10*RADIUS) {
      fprintf (stderr, "error: %10.6f,%10.6f vs %10.6f,%10.6f (%f,%f vs %f,%f)\n", 
	       catalog[0].average[n].R, catalog[0].average[n].D, 
	       stars[N].R, stars[N].D,
	       X1[i], X2[Jmin], 
	       Y1[i], Y2[Jmin]);
    }
    catalog[0].measure[Nmeas].dD       = 3600.0*(catalog[0].average[n].D - stars[N].D);

    /* XXX need to add dX, dY : need to load into stars[N].dX,dY */
    /* XXX need to add stargal, Sky, dSky, qPSF, detID, imageID */
    catalog[0].measure[Nmeas].Xccd     = stars[N].X;
    catalog[0].measure[Nmeas].Yccd     = stars[N].Y;

    catalog[0].measure[Nmeas].M        = stars[N].M;
    catalog[0].measure[Nmeas].dM       = stars[N].dM;  /* error in input files stored in thousandths of mag */
    catalog[0].measure[Nmeas].Mcal     = stars[N].Mcal;
    catalog[0].measure[Nmeas].t        = stars[N].t;
    catalog[0].measure[Nmeas].averef   = n;
    catalog[0].measure[Nmeas].photcode = stars[N].code;  /* photcode */
    catalog[0].measure[Nmeas].dophot   = stars[N].dophot;  
    catalog[0].measure[Nmeas].dbFlags  = 0;
    catalog[0].measure[Nmeas].dt       = stars[N].dt;
    catalog[0].measure[Nmeas].airmass  = stars[N].airmass;

    catalog[0].measure[Nmeas].photFlags = stars[N].flags;
    catalog[0].measure[Nmeas].qPSF      = stars[N].psfQual;
    catalog[0].measure[Nmeas].psfProb   = stars[N].psfProb;
    catalog[0].measure[Nmeas].crNsigma  = stars[N].crNsigma;
    catalog[0].measure[Nmeas].extNsigma = stars[N].extNsigma;
    catalog[0].measure[Nmeas].Sky       = stars[N].sky;
    catalog[0].measure[Nmeas].dSky      = stars[N].dsky;

    catalog[0].measure[Nmeas].stargal   = 0;

    catalog[0].measure[Nmeas].detID     = stars[N].detID;
    catalog[0].measure[Nmeas].imageID   = options.imageID;

    catalog[0].measure[Nmeas].dXccd     = stars[N].dX;
    catalog[0].measure[Nmeas].dYccd     = stars[N].dY;

    catalog[0].measure[Nmeas].Mgal     = stars[N].Mgal;

    // XXX saturate range for FWx, FWy, theta
    // XXX convert to unsigned int for these...
    catalog[0].measure[Nmeas].FWx      = 100*stars[N].fx;
    catalog[0].measure[Nmeas].FWy      = 100*stars[N].fy;
    catalog[0].measure[Nmeas].theta    = (0xffff/360.0)*stars[N].df;
	
    /* set the average magnitude if not already set and the photcode.equiv is not 0 */
    /* in UPDATE mode, this value is not saved; use relphot to recalculate */
    if (Nsec > -1) { 
      if (isnan(catalog[0].secfilt[n*Nsecfilt+Nsec].M)) {
	catalog[0].secfilt[n*Nsecfilt+Nsec].M = PhotCat (&catalog[0].measure[Nmeas]);
      }
    }

    /* adds the measurement to the calibration if appropriate color terms are found */
    if (options.calibrate) {
      AddToCalibration (&catalog[0].average[n], &catalog[0].secfilt[n*Nsecfilt], catalog[0].measure, &catalog[0].measure[Nmeas], next_meas, N);
    }

    /* if we choose to flag close encounters, see find_matches.c */
    /* if we choose to calculate RA,DEC averages, see update_coords.c */

    /* Nm is updated, but not written out in -update mode (for existing entries)
       Nm is recalculated in build_meas_links if loaded table is not sorted */
    stars[N].found = Nmeas;
    catalog[0].found[n] = Nmeas;
    catalog[0].average[n].Nm ++;
    Nmeas ++;
    i++;
  }

  /** incorporate unmatched image stars, if this star is in field of this catalog **/
  /* these new entries are all written out in UPDATE mode */ 
  for (i = 0; (i < Nstars) && !options.only_match; i++) {
    /* make sure there is space for next entry */
    if (Nmeas >= NMEAS) {
      NMEAS = Nmeas + 1000;
      REALLOCATE (next_meas, int, NMEAS);
      REALLOCATE (catalog[0].measure, Measure, NMEAS);
    }
    if (Nave >= NAVE) {
      NAVE = Nave + 1000;
      REALLOCATE (catalog[0].average, Average, NAVE);
      REALLOCATE (catalog[0].secfilt, SecFilt, NAVE*catalog[0].Nsecfilt);
    }

    N = N1[i];
    if (stars[N].found >= 0) continue;
    if (!IN_REGION (stars[N].R, stars[N].D)) continue;

    catalog[0].average[Nave].R         = stars[N].R;
    catalog[0].average[Nave].D         = stars[N].D;
    catalog[0].average[Nave].Xp        = 0;
    catalog[0].average[Nave].Nm        = 1;
    catalog[0].average[Nave].Nn        = 0;
    catalog[0].average[Nave].offset    = Nmeas;
    catalog[0].average[Nave].missing   = -1;
    catalog[0].average[Nave].code      = 0;

    catalog[0].average[Nave].dR        = 0;
    catalog[0].average[Nave].dD        = 0;
    catalog[0].average[Nave].uR        = 0;
    catalog[0].average[Nave].uD        = 0;
    catalog[0].average[Nave].duR       = 0;
    catalog[0].average[Nave].duD       = 0;
    catalog[0].average[Nave].P         = 0;
    catalog[0].average[Nave].dP        = 0;

    for (j = 0; j < Nsecfilt; j++) {
      catalog[0].secfilt[Nave*Nsecfilt+j].M  = NAN;
      catalog[0].secfilt[Nave*Nsecfilt+j].dM = NAN;
      catalog[0].secfilt[Nave*Nsecfilt+j].Xm = NAN_S_SHORT;
    }

    catalog[0].measure[Nmeas].Xccd     = stars[N].X;
    catalog[0].measure[Nmeas].Yccd     = stars[N].Y;

    catalog[0].measure[Nmeas].dR       = 0.0;
    catalog[0].measure[Nmeas].dD       = 0.0;
    catalog[0].measure[Nmeas].M        = stars[N].M;
    catalog[0].measure[Nmeas].dM       = stars[N].dM;
    catalog[0].measure[Nmeas].Mcal     = stars[N].Mcal;
    catalog[0].measure[Nmeas].t        = stars[N].t;
    catalog[0].measure[Nmeas].averef   = Nave;
    catalog[0].measure[Nmeas].photcode = stars[N].code;  /* photcode */
    catalog[0].measure[Nmeas].dophot   = stars[N].dophot;  
    catalog[0].measure[Nmeas].dbFlags  = 0;
    catalog[0].measure[Nmeas].dt       = stars[N].dt;
    catalog[0].measure[Nmeas].airmass  = stars[N].airmass;

    catalog[0].measure[Nmeas].photFlags = stars[N].flags;
    catalog[0].measure[Nmeas].qPSF      = stars[N].psfQual;
    catalog[0].measure[Nmeas].psfProb   = stars[N].psfProb;
    catalog[0].measure[Nmeas].crNsigma  = stars[N].crNsigma;
    catalog[0].measure[Nmeas].extNsigma = stars[N].extNsigma;
    catalog[0].measure[Nmeas].Sky       = stars[N].sky;
    catalog[0].measure[Nmeas].dSky      = stars[N].dsky;

    catalog[0].measure[Nmeas].stargal   = 0;

    catalog[0].measure[Nmeas].detID     = stars[N].detID;
    catalog[0].measure[Nmeas].imageID   = options.imageID;

    catalog[0].measure[Nmeas].dXccd     = stars[N].dX;
    catalog[0].measure[Nmeas].dYccd     = stars[N].dY;

    catalog[0].measure[Nmeas].Mgal     = stars[N].Mgal;
    catalog[0].measure[Nmeas].FWx      = 100*stars[N].fx;
    catalog[0].measure[Nmeas].FWy      = 100*stars[N].fy;
    catalog[0].measure[Nmeas].theta    = (0xffff/360.0)*stars[N].df;

    /* set the average magnitude if not already set and the photcode.equiv is not 0 */
    /* in UPDATE mode, this value is not saved; use relphot to recalculate */
    if (Nsec > -1) { 
	catalog[0].secfilt[Nave*Nsecfilt+Nsec].M = PhotCat (&catalog[0].measure[Nmeas]);
    }

    /* next[Nmeas] should always be -1 in this context (it is always the only
       measurement for the star) */
    stars[N].found = Nmeas;
    next_meas[Nmeas] = -1;
    Nmeas ++;
    Nave ++;
  }
      
  REALLOCATE (catalog[0].average, Average, Nave);
  REALLOCATE (catalog[0].measure, Measure, Nmeas);
 
  if (options.nosort) {
    catalog[0].sorted = FALSE;
  } else {
    catalog[0].sorted = TRUE;
    catalog[0].measure = sort_measure (catalog[0].average, Nave, catalog[0].measure, Nmeas, next_meas);
  }

  /* note stars which have been found in this catalog */
  for (i = 0; i < NstarsIn; i++) {
    if (stars[i].found > -1) {
      stars[i].found = -2;
    } else {
      stars[i].found = -3;
    }
  }

  /* check if the catalog has changed?  if no change, no need to write */
  catalog[0].Naverage = Nave;
  catalog[0].Nmeasure = Nmeas;
  catalog[0].Nsecf_mem = Nave*Nsecfilt;
  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas: %d %d %d, (%d matches)\n", Nstars, Nave, Nmeas, Nmatch);

  free (catalog[0].found);
  free (X1);
  free (Y1);
  free (N1);
  free (X2);
  free (Y2);
  free (N2);
  return (Nmatch);
}

/* 
   notes:
   
   for finding if a catalog star is in an image or an image star is in the catalog:
   
   catalogs have boundaries defined by RA and DEC, but they may curve in projection
   images have boundaries which are lines in pixels coords, but curve in RA and DEC
   
   catalog[0].found[Ncat] but stars[Nstar].found
   
*/
