# 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;
  float *X1, *Y1, *X2, *Y2;
  float dX, dY, dR;
  int *N1, *N2,  *next, *next_miss, last, last_miss;
  int Nave, NAVE, Nmeas, NMEAS, Nmiss, NMISS, Nmatch;
  Measure *tmpmeasure;
  Missing *tmpmissing;
  Coords tcoords;
  int Nsecfilt;

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

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

  /** allocate local arrays (catalog) **/
  Nave = catalog[0].Naverage;
  NAVE = Nave + 1000;
  ALLOCATE (X2, float, NAVE);
  ALLOCATE (Y2, float, 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);
  
  /* 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++) {
    fRD_to_XY (&X1[i], &Y1[i], stars[i][0].R, stars[i][0].D, &tcoords);
    N1[i] = i;
  }
  if (Nstars > 1) sort_lists (X1, Y1, N1, Nstars);
  
  for (i = 0; i < Nave; i++) {
    fRD_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_lists (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 **/
  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].offset;  

      /** 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].Nm - 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_PS       = 3600.0*(catalog[0].average[n].R - stars[N][0].R);
      catalog[0].measure[Nmeas].dD_PS       = 3600.0*(catalog[0].average[n].D - stars[N][0].D);
      catalog[0].measure[Nmeas].M_PS        = MIN (stars[N][0].M,  NO_MAG);
      catalog[0].measure[Nmeas].dM_PS       = MIN (stars[N][0].dM, NO_ERR);
      catalog[0].measure[Nmeas].Mcal_PS     = 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].source      = stars[N][0].code;
      catalog[0].measure[Nmeas].dophot      = 0;
      catalog[0].measure[Nmeas].flags       = 0;
      catalog[0].measure[Nmeas].dt_PS       = 0xffff;

      catalog[0].measure[Nmeas].Mgal_PS     = NO_MAG;
      catalog[0].measure[Nmeas].airmass_PS  = 0;
      catalog[0].measure[Nmeas].FWx         = NO_MAG;
      catalog[0].measure[Nmeas].FWy         = NO_MAG;
      catalog[0].measure[Nmeas].theta       = NO_ERR;
	
      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].flags |= BLEND_IMAGE;
	catalog[0].measure[Nmeas].flags |= BLEND_IMAGE;
      } 
      if (stars[N][0].found == -2) { /* this image star matches a catalog star on a neighboring catalog */
	catalog[0].measure[Nmeas].flags |= BLEND_IMAGE_NEIGHBOR;
      } 
      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]].flags |= BLEND_CATALOG;
	catalog[0].measure[Nmeas].flags |= BLEND_CATALOG;
      } else {
	catalog[0].found[n] = Nmeas;
      }

      catalog[0].average[n].Nm ++;
      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++) {
    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;

    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].M         = NO_MAG;
    catalog[0].average[Nave].dM        = NO_MAG;
    catalog[0].average[Nave].Nm        = 1;
    catalog[0].average[Nave].Nn        = 0;
    catalog[0].average[Nave].Xp        = NO_MAG;
    catalog[0].average[Nave].Xm        = NO_MAG;
    catalog[0].average[Nave].Xg        = NO_MAG;
    catalog[0].average[Nave].offset    = Nmeas;
    catalog[0].average[Nave].missing   = -1;
    catalog[0].average[Nave].code      = 0;

    for (j = 0; j < Nsecfilt; j++) {
      catalog[0].secfilt[Nave*Nsecfilt+j].M_PS  = NO_MAG;
      catalog[0].secfilt[Nave*Nsecfilt+j].dM_PS = NO_MAG;
      catalog[0].secfilt[Nave*Nsecfilt+j].Xm    = NO_MAG;
    }

    catalog[0].measure[Nmeas].dR_PS       = 0.0;
    catalog[0].measure[Nmeas].dD_PS       = 0.0;
    catalog[0].measure[Nmeas].M_PS        = MIN (stars[N][0].M,  NO_MAG);
    catalog[0].measure[Nmeas].dM_PS       = MIN (stars[N][0].dM, NO_ERR);
    catalog[0].measure[Nmeas].Mcal_PS     = 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].source   	  = stars[N][0].code;
    catalog[0].measure[Nmeas].dophot   	  = 0;
    catalog[0].measure[Nmeas].flags    	  = 0;
    catalog[0].measure[Nmeas].dt_PS       = 0xffff;

    catalog[0].measure[Nmeas].Mgal_PS     = NO_MAG;
    catalog[0].measure[Nmeas].airmass_PS  = 0;
    catalog[0].measure[Nmeas].FWx      	  = NO_MAG;
    catalog[0].measure[Nmeas].FWy      	  = NO_ERR;
    catalog[0].measure[Nmeas].theta    	  = NO_ERR;

    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].offset;
    catalog[0].average[i].offset = N;
    for (k = 0; k < catalog[0].average[i].Nm; 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].Nn > 0) {
      n = catalog[0].average[i].missing;
      catalog[0].average[i].missing = N;
      for (k = 0; k < catalog[0].average[i].Nn; 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].Naverage = Nave;
  catalog[0].Nmeasure = Nmeas;
  catalog[0].Nmissing = Nmiss;
  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas, Nmiss: %d %d %d %d, (%d matches)\n", Nstars, Nave, Nmeas, Nmiss, Nmatch);
  return (Nmatch);
}
