# include "addstar.h"

# define IN_CATALOG(R,D) ( \
((D) >= region[0].DEC[0]) && ((D) < region[0].DEC[1]) && \
((R) >= region[0].RA[0])  && ((R) < region[0].RA[1]))

void find_matches (GSCRegion *region, Stars *stars, int Nstars, Catalog *catalog, Image *image, Image *pimage, int Nimage) {

  int i, j, k, n, m, N, first_j, found, found0, found1;
  double X, Y, RADIUS, RADIUS2;
  float *X1, *Y1, *X2, *Y2, CalM, CalC, *Xs, *Ys;
  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, MTIME;
  short int Mrel, *Mval;

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

  /** allocate local arrays **/

  Nsecfilt = GetPhotcodeNsecfilt ();

  Nave = catalog[0].Naverage;
  NAVE = Nave + 1000;
  ALLOCATE (X2, float, NAVE);
  ALLOCATE (Y2, float, NAVE);
  ALLOCATE (N2, int, NAVE);
  ALLOCATE (Xs, float, NAVE);
  ALLOCATE (Ys, float, 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);
  
  /* check on photcode */
  Nsec = GetPhotcodeNsec (thiscode[0].code);

  /* project onto rectilinear grid with 1 arcsec pixels, sort by X */
  /* reference for coords is this image */
  tcoords = image[0].coords;
  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");
  
  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 (&Xs[i], &Ys[i], catalog[0].average[i].R, catalog[0].average[i].D, &image[0].coords);
    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 = NSIGMA * 0.02 * image[0].cerror;
  } else {
    RADIUS = DEFAULT_RADIUS;
  }
  RADIUS2 = RADIUS*RADIUS;
  /* correct instrumental mags for exposure time */
  if (image[0].exptime > 0) {
    MTIME = 2500*log10(image[0].exptime);
  } else {
    MTIME = 0;
  }

  /** 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 ++;
	/** insert star in measurement list */
	n = N2[j];
	N = N1[i];
	m = catalog[0].average[n].offset;  /* first measurement of this star */
	for (k = 0; k < catalog[0].average[n].Nm - 1; k++)
	  m = next[m];
	next[Nmeas] = next[m];
	next[m] = Nmeas;
	if (next[Nmeas] == -1) { /* last just was moved */
	  last = Nmeas;
	}
	
	if (CALIBRATE) {
	  found = FALSE;
	  found0 = FALSE;
	  found1 = FALSE;
	  m = catalog[0].average[n].offset;  /* first measurement of this star */
	  for (k = 0; !found && (k < catalog[0].average[n].Nm); k++) {
	    if (catalog[0].measure[m+k].source == CalReference) { found0 = TRUE; CalM = catalog[0].measure[m+k].M; }
	    if (catalog[0].measure[m+k].source == CalColor)     { found1 = TRUE; CalC = catalog[0].measure[m+k].M; }
	    if (found0 && found1) {
	      found = TRUE;
	      SaveCalibration (MIN (NO_MAG, 1000*stars[N].M + MTIME), 
			       MIN (255, stars[N].dM),
			       CalM, CalC,
			       image[0].secz, N,
			       stars[N].R, stars[N].D, Xs[n], Ys[n]);
	    }
	  }
	}
	      
	/** 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      = MIN (NO_MAG, 1000*stars[N].M + MTIME);
	catalog[0].measure[Nmeas].Mcal   = image[0].Mcal;
	catalog[0].measure[Nmeas].dM     = MIN (255, stars[N].dM);  /* error in input files stored in thousandths of mag */
	catalog[0].measure[Nmeas].t      = image[0].tzero + 1e-4*stars[N].Y*image[0].trate;  /* trate is in 0.1 msec / row */
	catalog[0].measure[Nmeas].averef = n;
	catalog[0].measure[Nmeas].source = image[0].source;  /* photometry source */
	catalog[0].measure[Nmeas].dophot = stars[N].dophot;  
	catalog[0].measure[Nmeas].flags  = 0;
	catalog[0].measure[Nmeas].dt     = MTIME;

	catalog[0].measure[Nmeas].Mgal     = MIN (NO_MAG, 1000*stars[N].Mgal + MTIME);
	catalog[0].measure[Nmeas].airmass  = image[0].secz;
	catalog[0].measure[Nmeas].FWx      = MIN (NO_MAG, 100*stars[N].fx);
	catalog[0].measure[Nmeas].fwy      = MIN (255, 100*(stars[N].fy / stars[N].fx));
	catalog[0].measure[Nmeas].theta    = MIN (255, (255/360)*stars[N].df);
	/* refers to same number as first measurement */
	
	/* check for entries in the secfilt lists */
	Mrel = catalog[0].measure[Nmeas].M + catalog[0].measure[Nmeas].airmass * thiscode[0].K;
	Mval = (Nsec == -1) ? &catalog[0].average[n].M : &catalog[0].secfilt[n*Nsecfilt+Nsec].M;
	if (*Mval == NO_MAG) *Mval = Mrel;

	/*** 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++;
  }

  /* add reference for undetected catalog stars */
  for (j = 0; (j < Nave) && !SKIP_MISSED; j++) {
    n = N2[j];
    if (catalog[0].found[n] < 0) { 
      /* should the catalog star be on this image? project into image coords */
      if (in_image (catalog[0].average[n].R, catalog[0].average[n].D, image)) {
	/* find last missing meas */
	if (catalog[0].average[n].Nn < 1) { /* no previous missing obs */
	  catalog[0].average[n].missing = Nmiss;
	  next_miss[last_miss] = Nmiss;
	  next_miss[Nmiss] = -1;
	  last_miss = Nmiss;
	} else {
	  m = catalog[0].average[n].missing;
	  for (k = 0; k < catalog[0].average[n].Nn - 1; k++)
	    m = next_miss[m];
	  next_miss[Nmiss] = next_miss[m];
	  next_miss[m] = Nmiss;
	  if (next_miss[Nmiss] == -1) { /* last just was moved */
	    last_miss = Nmiss;
	  }
	}
	RD_to_XY (&X, &Y, catalog[0].average[n].R, catalog[0].average[n].D, &image[0].coords);	  
	catalog[0].missing[Nmiss].t  = image[0].tzero + 1e-4*Y*image[0].trate;  /* trate is in 0.1 msec / row */
	catalog[0].average[n].Nn ++;
	Nmiss ++;
	if (Nmiss == NMISS) {
	  NMISS = Nmiss + 1000;
	  REALLOCATE (next_miss, int, NMISS);
	  REALLOCATE (catalog[0].missing, Missing, NMISS);
	}
      }
    }
  }

  /* incorporate unmatched image stars, if this star is in field of this catalog */
  for (i = 0; i < Nstars; i++) {
    N = N1[i];
    if ((stars[N].found < 0) && IN_CATALOG (stars[N].R, stars[N].D)) {
      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;
      }
      
      Mrel = catalog[0].measure[Nmeas].M + catalog[0].measure[Nmeas].airmass * thiscode[0].K;
      Mval = (Nsec == -1) ? &catalog[0].average[Nave].M : &catalog[0].secfilt[Nave*Nsecfilt+Nsec].M;
      if (*Mval == NO_MAG) *Mval = Mrel;

      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        = MIN (NO_MAG, 1000.0*stars[N].M + MTIME);
      catalog[0].measure[Nmeas].Mcal     = image[0].Mcal;
      catalog[0].measure[Nmeas].dM       = MIN (255, stars[N].dM);
      catalog[0].measure[Nmeas].t        = image[0].tzero + 1e-4*stars[N].Y*image[0].trate; /* trate is in 0.1 msec / row */
      catalog[0].measure[Nmeas].averef   = Nave;
      catalog[0].measure[Nmeas].source   = image[0].source;  /* photometry source */
      catalog[0].measure[Nmeas].dophot   = stars[N].dophot;  
      catalog[0].measure[Nmeas].flags    = 0;
      catalog[0].measure[Nmeas].dt       = MTIME;

      catalog[0].measure[Nmeas].Mgal     = MIN (NO_MAG, 1000.0*stars[N].Mgal + MTIME);
      catalog[0].measure[Nmeas].airmass  = image[0].secz;
      catalog[0].measure[Nmeas].FWx      = MIN (0x7fff, 100*stars[N].fx);
      catalog[0].measure[Nmeas].fwy      = MIN (255, 100*(stars[N].fy / stars[N].fx));
      catalog[0].measure[Nmeas].theta    = MIN (255, (255/360)*stars[N].df);
      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);
      }
      /** now add references from all previous non-detection observations of this spot on the sky */
      for (j = 0; (j < Nimage) && !SKIP_MISSED; 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);
	  }
  
	}
      }
      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;
    }
  }

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

int in_image (double r, double d, Image *image) {

  double X, Y;

  RD_to_XY (&X, &Y, r, d, &image[0].coords);
  if (((X) >= 0) && ((X) < image[0].NX) && ((Y) >= 0) && ((Y) < image[0].NY))
    return (TRUE);
  else 
    return (FALSE);

}

/* 
   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

       */

     
# ifdef OLDJUNKHERE

 
  /* fix measurement list order */
  for (i = 0; i < Nmeas - 1; i++) {
    if (next[i] != i + 1) {
      n0 = next[i];
      n1 = next[i+1];
      n2 = next[n0];
      t = i;
      while ((t != -1) && (next[t] != i+1))
	t = next[t];
      if ((t == -1) || (n2 == i+1)) {
	t = n2;
      }
      next[i] = i+1;
      next[i+1] = n2;
      next[n0] = n1;
      next[t] = n0;
      tmpmeasure = catalog[0].measure[n0];
      catalog[0].measure[n0] = catalog[0].measure[i+1];
      catalog[0].measure[i+1] = tmpmeasure;
    }
  }
 
  /* fix offset references */
  n = (0x00ffffff & catalog[0].measure[0].average);
  catalog[0].average[n].offset = 0;
  for (i = 1; i < Nmeas; i++) {
    if (n != (0x00ffffff & catalog[0].measure[i].average)) {
      n = (0x00ffffff & catalog[0].measure[i].average);
      catalog[0].average[n].offset = i;
    }
  }

  /* fix missing list order */
  for (i = 0; i < Nmiss - 1; i++) {
    if (next_miss[i] != i + 1) {
      n0 = next_miss[i];
      n1 = next_miss[i+1];
      n2 = next_miss[n0];
      t = i;
      while ((t != -1) && (next_miss[t] != i+1))
	t = next_miss[t];
      if ((t == -1) || (n2 == i+1)) {
	t = n2;
      }
      next_miss[i] = i+1;
      next_miss[i+1] = n2;
      next_miss[n0] = n1;
      next_miss[t] = n0;
      tmpmissing = catalog[0].missing[n0];
      catalog[0].missing[n0] = catalog[0].missing[i+1];
      catalog[0].missing[i+1] = tmpmissing;
    }
  }
 
  /* fix offset references */
  n = (0x00ffffff & catalog[0].measure[0].average);
  catalog[0].average[n].offset = 0;
  for (i = 1; i < Nmeas; i++) {
    if (n != (0x00ffffff & catalog[0].measure[i].average)) {
      n = (0x00ffffff & catalog[0].measure[i].average);
      catalog[0].average[n].offset = i;
    }
  }

# endif
