# include "addstar.h"

void find_matches (GSCRegion *region, Stars *stars, int Nstars, Catalog *catalog) {

  int i, j, k, n, m, N, first_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, Nsec;

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

  Nsecfilt = GetPhotcodeNsecfilt ();
  Nsec     = GetPhotcodeNsec (thiscode[0].code);

  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;
  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].RA[0] + region[0].RA[1]);
  tcoords.crval2 = 0.5*(region[0].DEC[0] + region[0].DEC[1]);
  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---TAN");
  tcoords.Npolyterms = 0;
  
  for (i = 0; i < Nstars; i++) {
    fRD_to_XY (&X1[i], &Y1[i], stars[i].R, stars[i].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 */
  /* 0.02 corrects cerror to arcsec from storage units */
  if (DEFAULT_RADIUS == 0) {
    RADIUS = 2.0;
  } else {
    RADIUS = DEFAULT_RADIUS;
  }
  RADIUS2 = RADIUS*RADIUS;

  /** find matched stars **/
  for (i = j = 0; (i < Nstars) && (j < Nave); ) {
    
    dX = X1[i] - X2[j];
    /* fprintf (stderr, "%f   %f %f  %f %f    %f %f  %f %f\n", dX, X1[i], Y1[i], X2[j], Y2[j], stars[i].R, stars[i].D, catalog[0].average[j].R, catalog[0].average[j].D); */

    if (dX <= -2*RADIUS) {
      i++;
      continue;
    }
    if (dX >= 2*RADIUS) {
      j++;
      continue;
    }

    /* negative dX: j is too large, positive dX, i is too large */
    first_j = j;
    for (; (dX > -2*RADIUS) && (j < Nave); j++) {
      dX = X1[i] - X2[j];
      dY = Y1[i] - Y2[j];
      dR = dX*dX + dY*dY;
      if (dR < RADIUS2) {  /* new measurement of this star */
	Nmatch ++;
	n = N2[j];
	N = N1[i];

	/** in replace mode, search for entry and replace values M, dM, R, D */
	if (REPLACE) {
	  int found;
	  found = FALSE;
	  m = catalog[0].average[n].offset;  
	  for (k = 0; !found && (k < catalog[0].average[n].Nm); k++) {
	    if (catalog[0].measure[m+k].source == thiscode[0].code) {
	      catalog[0].measure[m+k].dR       = 360000.0*(catalog[0].average[n].R - stars[N].R);
	      catalog[0].measure[m+k].dD       = 360000.0*(catalog[0].average[n].D - stars[N].D);
	      catalog[0].measure[m+k].M        = 1000.0*stars[N].M;
	      catalog[0].measure[m+k].dM       = 1000.0*MIN (255, stars[N].dM);  /* error in input files stored in thousandths of mag */
	      found = TRUE;
	      stars[N].found = m+k;  /* save first match, in case coincidences are found */
	    }
	  }
	  if (found) continue;
	}

	/** insert star in measurement list */
	/* find last measurement of this star */
	m = catalog[0].average[n].offset;  
	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 **/
	catalog[0].measure[Nmeas].dR       = 360000.0*(catalog[0].average[n].R - stars[N].R);
	catalog[0].measure[Nmeas].dD       = 360000.0*(catalog[0].average[n].D - stars[N].D);
	catalog[0].measure[Nmeas].M        = 1000.0*stars[N].M;
	catalog[0].measure[Nmeas].Mcal     = 0;
	catalog[0].measure[Nmeas].dM       = 1000.0*MIN (255, stars[N].dM);  /* error in input files stored in thousandths of mag */
	catalog[0].measure[Nmeas].t        = TIMEREF;
	catalog[0].measure[Nmeas].averef   = n;
	catalog[0].measure[Nmeas].source   = thiscode[0].code;
	catalog[0].measure[Nmeas].dophot   = stars[N].dophot;  /* error in input files stored in thousandths of mag */
	catalog[0].measure[Nmeas].flags    = 0;
	catalog[0].measure[Nmeas].dt       = 0xffff;

	catalog[0].measure[Nmeas].Mgal     = NO_MAG;
	catalog[0].measure[Nmeas].airmass  = 0;
	catalog[0].measure[Nmeas].FWx      = NO_MAG;
	catalog[0].measure[Nmeas].fwy      = 0xff;
	catalog[0].measure[Nmeas].theta    = 0xff;
	
	/* check for entries in the secfilt lists */
	if ((thiscode[0].type == PHOT_PRI) && (catalog[0].average[n].M == NO_MAG)) {
	  catalog[0].average[n].M = 1000.0*stars[N].M;
	}
	if ((thiscode[0].type == PHOT_SEC) && (catalog[0].secfilt[n*Nsecfilt+Nsec].M == NO_MAG)) {
	  catalog[0].secfilt[n*Nsecfilt+Nsec].M = 1000.0*stars[N].M;
	}

	/*** handle multiple stars */
	/* this image star matches more than one catalog star */
	if (stars[N].found > -1) {
	  catalog[0].measure[stars[N].found].flags |= BLEND_IMAGE;
	  catalog[0].measure[Nmeas].flags |= BLEND_IMAGE;
	} 
	if (stars[N].found == -2) { /* this image star matches a catalog star on a neighboring catalog */
	  catalog[0].measure[Nmeas].flags |= BLEND_IMAGE_NEIGHBOR;
	} 
	if (stars[N].found == -1) { /* this image star matches only this catalog star */
	  stars[N].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 values of Ro, Do */
	if (catalog[0].average[n].Nm > 2) {
	  double R, D, r, d, r2, d2, Npt;
	  Npt = r = d = r2 = d2 = 0;
	  m = catalog[0].average[n].offset;  /* first measurement of this star */
	  for (k = 0; k < catalog[0].average[n].Nm; k++) {
	    if (catalog[0].measure[m].t == 0) {
	      m = next[m];
	      continue;
	    }
	    R = catalog[0].measure[m].dR;
	    D = catalog[0].measure[m].dD;
	    r += R;
	    d += D;
	    r2 += R*R;
	    d2 += D*D;
	    m = next[m];
	    Npt += 1.0;
	  }
	  if (Npt > 2) {
	    r = r / Npt;  /* these are corrections in 1/100 arcsec to RA and DEC */
	    d = d / Npt;
	    R = r2 / Npt - r*r;
	    D = d2 / Npt - d*d;
	    /* Xp is scatter in position in hundredths of arcsec */
	    catalog[0].average[n].Xp = sqrt (D + R / SQ(cos(catalog[0].average[n].D*RAD_DEG)));
	    m = catalog[0].average[n].offset;  /* first measurement of this star */
	    for (k = 0; k < catalog[0].average[n].Nm; k++) {
	      catalog[0].measure[m].dR = catalog[0].measure[m].dR - r;
	      catalog[0].measure[m].dD = catalog[0].measure[m].dD - d;
	      m = next[m];
	    }
	    catalog[0].average[n].R = catalog[0].average[n].R - r / 360000.0;
	    catalog[0].average[n].D = catalog[0].average[n].D - d / 360000.0;
	  }
	}
      }
    }
    j = first_j;
    i++;
  }

  /* incorporate unmatched image stars, if this star is in field of this catalog */
  /* skip if we want to force matches - combined with -replace, this lets us keep 
     the reference up-to-date with known stars only */
  for (i = 0; (i < Nstars) && !ONLY_MATCH; i++) {
    N = N1[i];
    if (stars[N].found < 0) {
      catalog[0].average[Nave].R = stars[N].R;
      catalog[0].average[Nave].D = stars[N].D;
      catalog[0].average[Nave].M = NO_MAG;
      for (j = 0; j < Nsecfilt; j++) {
	catalog[0].secfilt[Nave*Nsecfilt+j].M  = NO_MAG;
	catalog[0].secfilt[Nave*Nsecfilt+j].Xm = NO_MAG;
      }

      if (thiscode[0].type == PHOT_PRI) {
	catalog[0].average[Nave].M = 1000.0*stars[N].M;
      }
      if (thiscode[0].type == PHOT_SEC) {
	catalog[0].secfilt[Nave*Nsecfilt+Nsec].M = 1000.0*stars[N].M;
      }
      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].dM        = NO_MAG;
      catalog[0].average[Nave].offset    = Nmeas;
      catalog[0].average[Nave].missing   = -1;
      catalog[0].average[Nave].code      = 0;

      catalog[0].measure[Nmeas].dR       = 0.0;
      catalog[0].measure[Nmeas].dD       = 0.0;
      catalog[0].measure[Nmeas].M        = 1000.0*stars[N].M;
      catalog[0].measure[Nmeas].Mcal     = 0;
      catalog[0].measure[Nmeas].dM       = MIN (255, stars[N].dM);
      catalog[0].measure[Nmeas].t        = 0;
      catalog[0].measure[Nmeas].averef   = Nave;
      catalog[0].measure[Nmeas].source   = thiscode[0].code;
      catalog[0].measure[Nmeas].dophot   = stars[N].dophot;  /* error in input files stored in thousandths of mag */
      catalog[0].measure[Nmeas].flags    = 0;
      catalog[0].measure[Nmeas].dt       = 0xffff;

      catalog[0].measure[Nmeas].Mgal     = NO_MAG;
      catalog[0].measure[Nmeas].airmass  = 0;
      catalog[0].measure[Nmeas].FWx      = NO_MAG;
      catalog[0].measure[Nmeas].fwy      = 0xff;
      catalog[0].measure[Nmeas].theta    = 0xff;

      stars[N].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);
  REALLOCATE (catalog[0].measure, Measure, Nmeas);

  /* 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].found > -1) {
      stars[i].found = -2;
    } else {
      stars[i].found = -3;
    }
  }

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


# if (0)
      /* this is a bit of an issue.  i don't want to include missing obs if the
	 ref star has not been previously detected (why?), but if i do this, it won't 
	 get all missed entries if ref star is detected in the future... */
      /** now add references from all previous non-detection observations of this spot on the sky */
      for (j = 0; j < Nimage; j++) {
	if (in_image (catalog[0].average[Nave].R, catalog[0].average[Nave].D, &pimage[j])) {
	  if (catalog[0].average[Nave].Nn < 1) {
	    catalog[0].average[Nave].missing = Nmiss;
	  }
	  next_miss[last_miss] = Nmiss;
	  next_miss[Nmiss] = -1;
	  last_miss = Nmiss;
	  /* this can now be done exactly */
	  RD_to_XY (&X, &Y, catalog[0].average[Nave].R, catalog[0].average[Nave].D, &pimage[j].coords);	  
	  catalog[0].missing[Nmiss].t  = pimage[j].tzero + 1e-4*Y*pimage[j].trate;  /* rough guess at time */
	  catalog[0].average[Nave].Nn ++;
	  Nmiss ++;
	  if (Nmiss == NMISS) {
	    NMISS = Nmiss + 1000;
	    REALLOCATE (next_miss, int, NMISS);
	    REALLOCATE (catalog[0].missing, Missing, NMISS);
	  }
  
	}
      }
# endif
