# include "addstar.h"

int find_matches_refstars (SkyRegion *region, Stars **stars, int Nstars, Catalog *catalog, AddstarClientOptions options) {

  int i, j, k, n, m, N, J;
  double RADIUS, RADIUS2;
  double *X1, *Y1, *X2, *Y2;
  double dX, dY, dR;
  int *N1, *N2,  *next, *next_miss, last, last_miss;
  int Nave, NAVE, Nmeas, NMEAS, Nmiss, NMISS, Nmatch;
  unsigned int objID, catID;
  Measure *tmpmeasure;
  Missing *tmpmissing;
  Coords tcoords;
  int Nsecfilt;

  if ((NREFSTAR_GROUP != 1) && (NREFSTAR_GROUP != 3)) {
      fprintf (stderr, "ERROR: NREFSTAR_GROUP NOT SET!\n");
      exit (1);
  }

  /* photcode data -- should not have to modify secfilt / average */
  Nsecfilt = GetPhotcodeNsecfilt ();

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

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

  Nmatch = 0;
  Nmeas = catalog[0].Nmeasure;
  NMEAS = Nmeas + 1000;
  ALLOCATE (next, int, NMEAS);
  REALLOCATE (catalog[0].measure, Measure, NMEAS);
  
  Nmiss = catalog[0].Nmissing;
  NMISS = Nmiss + 1000;
  if ((NMISS < 1) || (NMISS > 1e10)) {
      fprintf (stderr, "weird value for NMISS: %d\n", NMISS);
  }
  ALLOCATE (next_miss, int, NMISS);
  REALLOCATE (catalog[0].missing, Missing, NMISS);
  
  // current max obj ID for this catalog
  objID = catalog[0].objID;
  catID = catalog[0].catID;

  /* project onto rectilinear grid with 1 arcsec pixels, sort by X */
  /* reference for coords is catalog center */
  tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax);
  tcoords.crval2 = 0.5*(region[0].Dmin + region[0].Dmax);
  tcoords.crpix1 = tcoords.crpix2 = 0.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;
  strcpy (tcoords.ctype, "RA---ZEA");
  tcoords.Npolyterms = 0;
  
  for (i = 0; i < Nstars; i++) {
    RD_to_XY (&X1[i], &Y1[i], stars[i][0].R, stars[i][0].D, &tcoords);
    N1[i] = i;
  }
  if (Nstars > 1) sort_coords_index (X1, Y1, N1, Nstars);
  
  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 measurements */
  for (i = 0; i < Nmeas - 1; i++) {
    next[i] = i+1;
  }
  next[i] = -1;
  last = i;

  for (i = 0; i < Nmiss - 1; i++) {
    next_miss[i] = i+1;
  }
  next_miss[i] = -1;
  last_miss = i;

  /* choose a radius for matches */
  if (options.radius == 0) {
    RADIUS = 2.0; /* hardwired default for refstars */
  } else {
    RADIUS = options.radius; /* provided by config */
  }
  RADIUS2 = RADIUS*RADIUS;

  /** find matched stars **/
  // XXX could use NREFSTAR_GROUP to do this match more quicky
  for (i = j = 0; (i < Nstars) && (j < Nave); ) {
    
    dX = X1[i] - X2[j];
    if (dX <= -2*RADIUS) {
      i++;
      continue;
    }
    if (dX >= 2*RADIUS) {
      j++;
      continue;
    }

    /* negative dX: j is too large, positive dX, i is too large */
    for (J = j; (dX > -2*RADIUS) && (J < Nave); J++) {
      dX = X1[i] - X2[J];
      dY = Y1[i] - Y2[J];
      dR = dX*dX + dY*dY;
      if (dR > RADIUS2) continue;

      Nmatch ++;
      n = N2[J];
      N = N1[i];
      m = catalog[0].average[n].measureOffset;  

      /** in replace mode, search for entry and replace values M, dM, R, D */
      if (options.replace && replace_match (&catalog[0].average[n], &catalog[0].measure[m], stars[N])) continue;

      /** insert star in measurement list */
      /* find last measurement of this star */
      for (k = 0; k < catalog[0].average[n].Nmeasure - 1; k++) m = next[m];
      /* set up references */
      next[Nmeas] = next[m];
      next[m] = Nmeas;
      /* last just was moved */
      if (next[Nmeas] == -1) last = Nmeas;
	
      /** add measurements for this star **/
      /** *** dR,dD now in arcsec *** **/
      catalog[0].measure[Nmeas].dR       = 3600.0*(catalog[0].average[n].R - stars[N][0].R);
      catalog[0].measure[Nmeas].dD       = 3600.0*(catalog[0].average[n].D - stars[N][0].D);
      catalog[0].measure[Nmeas].M        = stars[N][0].M;
      catalog[0].measure[Nmeas].dM       = stars[N][0].dM;
      catalog[0].measure[Nmeas].Mcal     = 0;
      catalog[0].measure[Nmeas].t        = (TIMEREF == 0) ? stars[N][0].t : TIMEREF; /** careful : time_t vs e_time **/
      catalog[0].measure[Nmeas].averef   = n;
      catalog[0].measure[Nmeas].photcode = stars[N][0].code;
      catalog[0].measure[Nmeas].dophot   = 0;
      catalog[0].measure[Nmeas].dbFlags  = 0;
      catalog[0].measure[Nmeas].dt       = 0xffff;

      catalog[0].measure[Nmeas].Map      = NAN;
      catalog[0].measure[Nmeas].airmass  = 0;
      catalog[0].measure[Nmeas].FWx      = stars[N][0].fx; // XXX make sure these are zero'ed as needed
      catalog[0].measure[Nmeas].FWy      = stars[N][0].fy; // XXX make sure these are zero'ed as needed
      catalog[0].measure[Nmeas].theta    = stars[N][0].df; // XXX make sure these are zero'ed as needed
	
      catalog[0].measure[Nmeas].photFlags = stars[N][0].flags; // XXX make sure these are zero'ed as needed
      catalog[0].measure[Nmeas].qPSF      = 0;
      catalog[0].measure[Nmeas].psfChisq  = 0;
      catalog[0].measure[Nmeas].crNsigma  = 0;
      catalog[0].measure[Nmeas].extNsigma = 0;
      catalog[0].measure[Nmeas].Sky       = 0;
      catalog[0].measure[Nmeas].dSky      = 0;

      catalog[0].measure[Nmeas].stargal   = 0; // XXX not yet set

      catalog[0].measure[Nmeas].detID     = 0;
      catalog[0].measure[Nmeas].imageID   = 0;

      catalog[0].measure[Nmeas].dXccd     = 0;
      catalog[0].measure[Nmeas].dYccd     = 0;

      catalog[0].measure[Nmeas].Xccd     = stars[N][0].X; // XXX make sure these are zero'ed as needed
      catalog[0].measure[Nmeas].Yccd     = stars[N][0].Y; // XXX make sure these are zero'ed as needed

      if (ACCEPT_MOTION) {
	catalog[0].average[n].uR         = stars[N][0].uR;
	catalog[0].average[n].uD         = stars[N][0].uD;
	catalog[0].average[n].duR        = stars[N][0].duR;
	catalog[0].average[n].duD        = stars[N][0].duD;
	catalog[0].average[n].P          = stars[N][0].P;
	catalog[0].average[n].dP         = stars[N][0].dP;
      }

      /** don't update average / secfilt values for REF photcodes **/

      /*** handle multiple stars */
      /* this image star matches more than one catalog star */
      if (stars[N][0].found > -1) {
	catalog[0].measure[stars[N][0].found].dbFlags |= ID_MEAS_BLEND_MEAS;
	catalog[0].measure[Nmeas].dbFlags |= ID_MEAS_BLEND_MEAS;
      } 
      if (stars[N][0].found == -2) { /* this image star matches a catalog star on a neighboring catalog */
	catalog[0].measure[Nmeas].dbFlags |= ID_MEAS_BLEND_MEAS_X;
      } 
      if (stars[N][0].found == -1) { /* this image star matches only this catalog star */
	stars[N][0].found = Nmeas;  /* save first match, in case coincidences are found */
      }
      /* this catalog star matches more than one image star */
      if (catalog[0].found[n] > -1) {
	catalog[0].measure[catalog[0].found[n]].dbFlags |= ID_MEAS_BLEND_OBJ;
	catalog[0].measure[Nmeas].dbFlags |= ID_MEAS_BLEND_OBJ;
      } else {
	catalog[0].found[n] = Nmeas;
      }

      catalog[0].average[n].Nmeasure ++;
      Nmeas ++;
      if (Nmeas == NMEAS) {
	NMEAS = Nmeas + 1000;
	REALLOCATE (next, int, NMEAS);
	REALLOCATE (catalog[0].measure, Measure, NMEAS);
      }

      update_coords (&catalog[0].average[n], &catalog[0].measure[0], next);
    }
    i++;
  }

  /* we don't add missed entries for refcat
     (already in database, not refcat) */

  /* incorporate unmatched refcat stars */
  /* skip if we want to require matches
     combined with -replace, this lets us keep 
     the reference up-to-date with known stars only */

  for (i = 0; (i < Nstars) && !options.only_match; i+=NREFSTAR_GROUP) {
    N = N1[i];
    if (stars[N][0].found >= 0) continue;

    catalog[0].average[Nave].R         	   = stars[N][0].R;
    catalog[0].average[Nave].D         	   = stars[N][0].D;
    catalog[0].average[Nave].Xp        	   = 0;
    catalog[0].average[Nave].Nmeasure      = NREFSTAR_GROUP;
    catalog[0].average[Nave].Nmissing      = 0;
    catalog[0].average[Nave].measureOffset = Nmeas;
    catalog[0].average[Nave].missingOffset = -1;
    catalog[0].average[Nave].code          = 0;

    if (ACCEPT_MOTION) {
      catalog[0].average[Nave].dR    = stars[N][0].dR;
      catalog[0].average[Nave].dD    = stars[N][0].dD;
      catalog[0].average[Nave].uR    = stars[N][0].uR;
      catalog[0].average[Nave].uD    = stars[N][0].uD;
      catalog[0].average[Nave].duR   = stars[N][0].duR;
      catalog[0].average[Nave].duD   = stars[N][0].duD;
      catalog[0].average[Nave].P     = stars[N][0].P;
      catalog[0].average[Nave].dP    = stars[N][0].dP;
    } else {
      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;
    }

    catalog[0].average[Nave].objID     = objID;
    catalog[0].average[Nave].catID     = catID;
    objID ++;

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

    for (j = 0; j < NREFSTAR_GROUP; j++) {
      N = N1[i + j];
      catalog[0].measure[Nmeas].dR       = 0.0;
      catalog[0].measure[Nmeas].dD       = 0.0;
      catalog[0].measure[Nmeas].M        = stars[N][0].M;
      catalog[0].measure[Nmeas].dM       = stars[N][0].dM;
      catalog[0].measure[Nmeas].Mcal     = 0;
      catalog[0].measure[Nmeas].t        = (stars[N][0].t == 0) ? TIMEREF : stars[N][0].t; /** careful : time_t vs e_time **/
      catalog[0].measure[Nmeas].averef   = Nave;
      catalog[0].measure[Nmeas].photcode = stars[N][0].code;
      catalog[0].measure[Nmeas].dophot   = 0;
      catalog[0].measure[Nmeas].dbFlags  = 0;
      catalog[0].measure[Nmeas].dt       = 0xffff;

      catalog[0].measure[Nmeas].photFlags = stars[N][0].flags; // XXX make sure these are zero'ed as needed
      catalog[0].measure[Nmeas].qPSF      = 0;
      catalog[0].measure[Nmeas].psfChisq  = 0;
      catalog[0].measure[Nmeas].crNsigma  = 0;
      catalog[0].measure[Nmeas].extNsigma = 0;
      catalog[0].measure[Nmeas].Sky       = 0;
      catalog[0].measure[Nmeas].dSky      = 0;

      catalog[0].measure[Nmeas].stargal   = 0; // XXX not yet set

      catalog[0].measure[Nmeas].detID     = 0;
      catalog[0].measure[Nmeas].imageID   = 0;

      catalog[0].measure[Nmeas].dXccd     = 0;
      catalog[0].measure[Nmeas].dYccd     = 0;

      catalog[0].measure[Nmeas].Xccd     = stars[N][0].X; // XXX make sure these are zero'ed as needed
      catalog[0].measure[Nmeas].Yccd     = stars[N][0].Y; // XXX make sure these are zero'ed as needed

      catalog[0].measure[Nmeas].airmass  = 0;
      catalog[0].measure[Nmeas].Map      = NAN;
      catalog[0].measure[Nmeas].FWx      = stars[N][0].fx; // XXX make sure these are zero'ed as needed
      catalog[0].measure[Nmeas].FWy      = stars[N][0].fy; // XXX make sure these are zero'ed as needed
      catalog[0].measure[Nmeas].theta    = stars[N][0].df; // XXX make sure these are zero'ed as needed

      catalog[0].measure[Nmeas].Xccd     = 0.0;
      catalog[0].measure[Nmeas].Yccd     = 0.0;

      stars[N][0].found = Nmeas;
      next[last] = Nmeas;
      next[Nmeas] = -1;
      last = Nmeas;
      Nmeas ++;
      if (Nmeas == NMEAS) {
	NMEAS = Nmeas + 1000;
	REALLOCATE (next, int, NMEAS);
	REALLOCATE (catalog[0].measure, Measure, NMEAS);
      }
    }

    Nave ++;
    if (Nave == NAVE) {
      NAVE = Nave + 1000;
      REALLOCATE (catalog[0].average, Average, NAVE);
      REALLOCATE (catalog[0].secfilt, SecFilt, NAVE*catalog[0].Nsecfilt);
    }
  }
      
  free (catalog[0].found);
  REALLOCATE (catalog[0].average, Average, Nave);

  /* fix order of Measure (memory intensive, but fast) */
  N = 0; 
  ALLOCATE (tmpmeasure, Measure, Nmeas);
  for (i = 0; i < Nave; i++) {
    n = catalog[0].average[i].measureOffset;
    catalog[0].average[i].measureOffset = N;
    for (k = 0; k < catalog[0].average[i].Nmeasure; k++, N++) {
      tmpmeasure[N] = catalog[0].measure[n]; 
      tmpmeasure[N].averef = i;
      n = next[n];
    }
  }
  free (catalog[0].measure);
  catalog[0].measure = tmpmeasure;
    
  /* fix order of Missing (memory intensive, but fast) */
  N = 0; 
  ALLOCATE (tmpmissing, Missing, Nmiss);
  for (i = 0; i < Nave; i++) {
    if (catalog[0].average[i].Nmissing > 0) {
      n = catalog[0].average[i].missingOffset;
      catalog[0].average[i].missingOffset = N;
      for (k = 0; k < catalog[0].average[i].Nmissing; k++, N++) {
	tmpmissing[N] = catalog[0].missing[n]; 
	n = next_miss[n];
      }
    }
  }
  free (catalog[0].missing);
  catalog[0].missing = tmpmissing;

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

  free (X1);
  free (Y1);
  free (N1);
  free (X2);
  free (Y2);
  free (N2);
  free (next);
  free (next_miss);

  catalog[0].objID    = objID; // new max value, save on catalog close
  catalog[0].Naverage = Nave;
  catalog[0].Nmeasure = Nmeas;
  catalog[0].Nmissing = Nmiss;
  catalog[0].Nsecf_mem = Nave*Nsecfilt;
  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas, Nmiss: %d %d %d %d, (%d matches)\n", Nstars, Nave, Nmeas, Nmiss, Nmatch);
  return (Nmatch);
}
