# include "addstar.h"

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

  off_t i, j, n, N, J;
  double RADIUS, RADIUS2;
  double *X1, *Y1, *X2, *Y2;
  double dX, dY, dR;
  off_t *N1, *N2,  *next_meas;
  off_t Nave, NAVE, Nmeas, NMEAS, Nmatch;
  unsigned int objID, catID;
  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, off_t,  Nstars);

  /** allocate local arrays (catalog) **/
  Nave = catalog[0].Naverage;
  NAVE = Nave + 1000;
  ALLOCATE (X2, double, NAVE);
  ALLOCATE (Y2, double, NAVE);
  ALLOCATE (N2, off_t,  NAVE);
  ALLOCATE (catalog[0].found, off_t, 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 = Nmeas = catalog[0].Nmeasure;
  
  // 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].average.R, stars[i][0].average.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 */
  if (catalog[0].sorted && (catalog[0].Nmeasure == catalog[0].Nmeas_disk)) {
    // this version is only valid if we have done a full catalog load, and if the catalog
    // is sorted while processed
    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);
  }    

  /* 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); ) {
    if (!finite(X1[i]) || !finite(Y1[i])) { 
      i++; 
      continue;
    }
    if (!finite(X2[j]) || !finite(Y2[j])) { 
      j++; 
      continue;
    }
    
    dX = X1[i] - X2[j];
    if (dX <= -1.02*RADIUS) {
      i++;
      continue;
    }
    if (dX >= 1.02*RADIUS) {
      j++;
      continue;
    }

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

      /*** a match is found, add to average, measure ***/
      Nmatch ++;
      n = N2[J];
      N = N1[i];

      /** 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, stars[N])) continue;

      if (Nmeas == NMEAS) {
	NMEAS = Nmeas + 1000;
	REALLOCATE (next_meas, off_t, NMEAS);
	REALLOCATE (catalog[0].measure, Measure, NMEAS);
      }

      /* add to end of measurement list */
      add_meas_link (&catalog[0].average[n], next_meas, Nmeas, NMEAS);
	
      /** add measurements for this star **/

      // set the new measurements
      catalog[0].measure[Nmeas]          = stars[N][0].measure;

      /** *** dR,dD now in arcsec *** **/
      catalog[0].measure[Nmeas].dR       = 3600.0*(catalog[0].average[n].R - stars[N][0].average.R);
      catalog[0].measure[Nmeas].dD       = 3600.0*(catalog[0].average[n].D - stars[N][0].average.D);
      catalog[0].measure[Nmeas].dbFlags  = 0;
      catalog[0].measure[Nmeas].averef   = n;
      catalog[0].measure[Nmeas].objID    = catalog[0].average[n].objID;
      catalog[0].measure[Nmeas].catID    = catalog[0].catID;

      // rationalize dR:
      if (catalog[0].measure[Nmeas].dR > +180.0*3600.0) {
	  // average on high end of boundary, move star up
	  stars[N][0].average.R += 360.0;
	  catalog[0].measure[Nmeas].dR = 3600.0*(catalog[0].average[n].R - stars[N][0].average.R);
      }
      if (catalog[0].measure[Nmeas].dR < -180.0*3600.0) {
	  // average on low end of boundary, move star down
	  stars[N][0].average.R -= 360.0;
	  catalog[0].measure[Nmeas].dR = 3600.0*(catalog[0].average[n].R - stars[N][0].average.R);
      }

      catalog[0].measure[Nmeas].t        = (TIMEREF == 0) ? stars[N][0].measure.t      : TIMEREF; /** careful : time_t vs e_time **/
      catalog[0].measure[Nmeas].t_msec   = (TIMEREF == 0) ? stars[N][0].measure.t_msec : TIMEREF; /** careful : time_t vs e_time **/

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

      // update_coords (&catalog[0].average[n], &catalog[0].measure[0], next_meas);
    }
    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) {
    if (Nmeas >= NMEAS - NREFSTAR_GROUP) {
      NMEAS = Nmeas + 1000;
      REALLOCATE (next_meas, off_t, 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][0].found > -1) continue;
    if (!IN_REGION (stars[N][0].average.R, stars[N][0].average.D)) continue;

    catalog[0].average[Nave].R         	   = stars[N][0].average.R;
    catalog[0].average[Nave].D         	   = stars[N][0].average.D;

    catalog[0].average[Nave].Nmeasure      = NREFSTAR_GROUP;
    catalog[0].average[Nave].Nmissing      = 0;
    catalog[0].average[Nave].Nextend       = 0;

    catalog[0].average[Nave].measureOffset = Nmeas;
    catalog[0].average[Nave].missingOffset = -1;
    catalog[0].average[Nave].extendOffset  = -1;

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

    catalog[0].average[Nave].Xp            = 0;
    catalog[0].average[Nave].ChiSqAve  	   = 0.0;
    catalog[0].average[Nave].ChiSqPM   	   = 0.0;
    catalog[0].average[Nave].ChiSqPar  	   = 0.0;
    catalog[0].average[Nave].Tmean   	   = 0;
    catalog[0].average[Nave].Trange   	   = 0;
    catalog[0].average[Nave].Npos    	   = 0;

    catalog[0].average[Nave].objID     	   = objID;
    catalog[0].average[Nave].catID     	   = catID;
    catalog[0].average[Nave].flags     	   = 0;
    if (PSPS_ID) {
        catalog[0].average[Nave].extID = CreatePSPSObjectID(catalog[0].average[Nave].R,
                                                            catalog[0].average[Nave].D);
    } else {
        catalog[0].average[Nave].extID         = 0;
    }


    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;
      catalog[0].secfilt[Nave*Nsecfilt+j].M_20 	= NAN_S_SHORT;
      catalog[0].secfilt[Nave*Nsecfilt+j].M_80 	= NAN_S_SHORT;
      catalog[0].secfilt[Nave*Nsecfilt+j].Ncode = 0;
      catalog[0].secfilt[Nave*Nsecfilt+j].Nused = 0;
    }

    for (j = 0; j < NREFSTAR_GROUP; j++) {
      N = N1[i + j];

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

      catalog[0].measure[Nmeas].dR       = 0.0;
      catalog[0].measure[Nmeas].dD       = 0.0;

      catalog[0].measure[Nmeas].t        = (stars[N][0].measure.t == 0) ? TIMEREF : stars[N][0].measure.t;      /** careful : time_t vs e_time **/
      catalog[0].measure[Nmeas].t_msec   = (stars[N][0].measure.t == 0) ?       0 : stars[N][0].measure.t_msec; /** careful : time_t vs e_time **/
      catalog[0].measure[Nmeas].dbFlags  = 0;
      catalog[0].measure[Nmeas].averef   = Nave;
      catalog[0].measure[Nmeas].objID    = catalog[0].average[Nave].objID;
      catalog[0].measure[Nmeas].catID    = catalog[0].catID;

      stars[N][0].found = Nmeas;
      next_meas[Nmeas] = -1;  // initial value here, update below
      Nmeas ++; 
    }

    for (j = 0; j < NREFSTAR_GROUP - 1; j++) {
	next_meas[Nmeas - NREFSTAR_GROUP + j] = Nmeas - NREFSTAR_GROUP + j + 1;
    }
    Nave ++;
  }
      
  REALLOCATE (catalog[0].average, Average, Nave);

  // XXX allow for unsorted output?
  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 < Nstars; i++) {
    if (stars[i][0].found > -1) {
      stars[i][0].found = -2;
    } else {
      stars[i][0].found = -3;
    }
  }

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

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

  return (Nmatch);
}
