# include "delstar.h"

find_matches (Catalog *catalog, Image *image, int start, int end) {

  int i, j, k, n, m, N, M, found;
  int *N1, *N2,  *next, *next_miss, *ave_miss, last, last_miss;
  int Nave, NAVE, Nmeas, NMEAS, Nmiss, NMISS, Nmatch, photcode;
  int Nmeasfound, Nsecfilt;
  int this, prev;
  unsigned int flags;
  Measure *tmpmeasure;
  Missing *tmpmissing;
  Coords tcoords;

  /** allocate local arrays **/
  Nave = catalog[0].Naverage;

  Nmeas = catalog[0].Nmeasure;
  ALLOCATE (next, int, MAX(Nmeas,1));
  
  Nmiss = catalog[0].Nmissing;
  ALLOCATE (next_miss, int, MAX(Nmiss,1));
  ALLOCATE (ave_miss, int, MAX(Nmiss,1));
  
  if (VERBOSE) fprintf (stderr, "starting with Nave, Nmeas, Nmiss: %d %d %d\n", Nave, Nmeas, Nmiss);

  /* set up pointers for linked list of measure */
  for (i = 0; i < Nmeas - 1; i++) {
    next[i] = i+1;
  }
  next[i] = -1;
  last = i;
  /* set up pointers for linked list of missing */
  for (i = 0; i < Nmiss - 1; i++) {
    next_miss[i] = i+1;
  }
  next_miss[i] = -1;
  last_miss = i;
  /* set up references for missing to average */
  for (i = 0; i < Nave; i++) {
    for (j = 0; j < catalog[0].average[i].Nn; j++) {
      ave_miss[catalog[0].average[i].missing + j] = i;
    }
  }
  photcode = -1;
  if ((start == 0) && (end == 0)) {
    start = image[0].tzero;
    end   = image[0].tzero + 1e-4*image[0].NY*image[0].trate;  /* trate is in 0.1 msec / row */
    photcode = image[0].source;
  }
  if (VERBOSE) fprintf (stderr, "deleting for range %d to %d (photcode %d)\n", start, end, photcode);
  Nmeasfound = 0;
  Nsecfilt = catalog[0].Nsecfilt;

  /* fprintf (stderr, "fixing the measures...\n"); */
  /** find measure in time range of image **/
  this = prev = -1;
  for (i = 0; (i < Nmeas); i++) {
    if ((i % 10000) == 0) fprintf (stderr, ". ");
    if ((catalog[0].measure[i].t >= start) && (catalog[0].measure[i].t <= end) && (photcode == catalog[0].measure[i].source)) { 
      Nmeasfound ++;
      /* this star is in this image */
      this = next[i];
      next[i] = -2; /* we delete this one */
      if (prev != -1) { next[prev] = this; }
# if (0) 
      /* why is this section disabled? */
      /* fix the list links: connect the previous valid link to the next valid link */
      for (j = i; (j >= 0) && (next[j] == -2); j--); /* find previous entry to fix link */
      if (j >= 0) { /* if j < 0, there is no previous valid link, ignore this step */
	if (next[j] != i) {
	  fprintf (stderr, "error? (1)  this link seems to have been lost\n");
	  fprintf (stderr, "j: %d, next[j]: %d, i: %d\n", j, next[j], i);
	  exit (1);
	}
	/* find next valid entry to fix link */
	for (k = i; (k < Nmeas) && (next[k] == -2); k++);
	if (k < Nmeas)
	  next[j] = k;
	else 
	  next[j] = -1;  /* last link in list gets a -1 */
      }
# endif      
      /*** fix the corresponding average entry ***/
      n = catalog[0].measure[i].averef;
      if (catalog[0].average[n].Nm == 0) { /* this should never happen */
	fprintf (stderr, "error? we deleted one too many objects?\n");
	exit (1);
      }
      catalog[0].average[n].Nm --;
      /* this was only entry in list: will be deleted below.  meanwhile, delete all missing entries*/
      if ((catalog[0].average[n].Nm < 1) && (catalog[0].average[n].Nn > 0)) { 
	m = catalog[0].average[n].missing;
	for (j = 0; j < catalog[0].average[n].Nn; j++) {
	  M = next_miss[m];
	  next_miss[m] = -2;
	  m = M;
	}
	m = catalog[0].average[n].missing;
	/* fix the list links: connect the previous valid link to the next valid link */
	for (j = m; (j >= 0) && (next_miss[j] == -2); j--); /* find previous entry to fix link */
	if (j >= 0) { /* if j < 0, there is no previous valid link, ignore this step */
	  if (next_miss[j] != m) {
	    fprintf (stderr, "error? (2) this link seems to have been lost\n");
	    fprintf (stderr, "j: %d, next_miss[j]: %d, i: %d\n", j, next_miss[j], i);
	    exit (1);
	  }
	  /* find next valid entry to fix link */
	  for (k = m; (k < Nmiss) && (next_miss[k] == -2); k++);
	  if (k < Nmiss)
	    next_miss[j] = k;
	  else 
	    next_miss[j] = -1;  /* last link in list gets a -1 */
	}
      }
      /* this was first entry in list */
      if ((catalog[0].average[n].offset == i) && (catalog[0].average[n].Nm > 0)) { 
	m = catalog[0].average[n].offset;
	/* find next valid entry -- notice lack of error checking... */
	for (j = 0; (j < Nmeas) && (next[m+j] == -2); j++);
	catalog[0].average[n].offset = m + j;
      }
    } else {
      prev = i;
    }
  }
  fprintf (stderr, "found %d meas to remove\n", Nmeasfound);

  /* fprintf (stderr, "fixing the missing...\n"); */
  /** find missing in time range of image **/
  for (i = 0; (i < Nmiss); i++) {
    if ((next_miss[i] != -2) && (catalog[0].missing[i].t >= start) && (catalog[0].missing[i].t <= end)) { 
      /* this star is in this image */

      next_miss[i] = -2; /* we delete this one */
      /* fix the list links: connect the previous valid link to the next valid link */
      for (j = i; (j >= 0) && (next_miss[j] == -2); j--); /* find previous entry to fix link */
      if (j >= 0) { /* if j < 0, there is no previous valid link, ignore this step */
	if (next_miss[j] != i) {
	  fprintf (stderr, "error? (3) this link seems to have been lost\n");
	  fprintf (stderr, "j: %d, next_miss[j]: %d, i: %d\n", j, next_miss[j], i);
	  exit (1);
	}
	/* find next valid entry to fix link */
	for (k = i; (k < Nmiss) && (next_miss[k] == -2); k++);
	if (k < Nmiss)
	  next_miss[j] = k;
	else 
	  next_miss[j] = -1;  /* last link in list gets a -1 */
      }
      
      /* find the corresponding avearge entry 
      found = FALSE;
      for (n = 0; !found && (n < Nave); n++) {
	if ((catalog[0].average[n].missing > 0) && (catalog[0].average[n].missing <= i) && (catalog[0].average[n].missing + catalog[0].average[n].Nn > i)) 
	  found = TRUE;
      }
      n--; */
      /*** fix the corresponding average entry ***/
      n = ave_miss[i];
      if (catalog[0].average[n].Nn == 0) { /* this should never happen */
	fprintf (stderr, "error? we deleted one too many missing?\n");
	exit (1);
      }
      catalog[0].average[n].Nn --;
      /* this was first entry in list */
      if ((catalog[0].average[n].missing == i) && (catalog[0].average[n].Nn > 0)) { 
	m = catalog[0].average[n].missing;
	for (j = 0; (j < Nmiss) && (next_miss[m+j] == -2); j++);
	catalog[0].average[n].missing = m + j;
      }

    }
  }

  /* currently not worked out, but we will need to delete the references to blended image and cat stars */
# if 0
  /*** 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 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;
  }
# endif  

  /* fix Average list: delete entries with Nm == 0 */
  for (i = j = 0; (i < Nave) && (j < Nave); i++, j++) {
    for (; (j < Nave) && (catalog[0].average[j].Nm == 0); j++);
    if ((i != j) && (j < Nave)) {
      catalog[0].average[i] = catalog[0].average[j];
      for (k = 0; k < catalog[0].Nsecfilt; k++) {
	catalog[0].secfilt[i*Nsecfilt + k] = catalog[0].secfilt[j*Nsecfilt + k];
      }
    }    
    if (j == Nave) i--;
  }
  Nave = i;
  REALLOCATE (catalog[0].average, Average, Nave);
  REALLOCATE (catalog[0].secfilt, SecFilt, MAX (1, Nave*Nsecfilt));

  /* 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++) {
      if ((n == -1) || (n == -2)) {
	fprintf (stderr, "error: linked list is confused\n");
	exit (1);
      }
      tmpmeasure[N] = catalog[0].measure[n]; 
      tmpmeasure[N].averef = i;
      n = next[n];
    }
  }
  Nmeas = N;
  free (catalog[0].measure);
  catalog[0].measure = tmpmeasure;
  REALLOCATE (catalog[0].measure, Measure, Nmeas);
    
  /* fprintf (stderr, "fixing the mising order...\n"); */
  /* 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++) {
	if ((n == -1) || (n == -2)) {
	  fprintf (stderr, "error: linked list is confused\n");
	  exit (1);
	}
	tmpmissing[N] = catalog[0].missing[n]; 
	n = next_miss[n];
      }
    }
  }
  Nmiss = N;
  free (catalog[0].missing);
  catalog[0].missing = tmpmissing;
  REALLOCATE (catalog[0].missing, Missing, Nmiss);

  fprintf (stderr, "\n");

  catalog[0].Naverage = Nave;
  catalog[0].Nmeasure = Nmeas;
  catalog[0].Nmissing = Nmiss;
  if (VERBOSE) fprintf (stderr, "  ending with Nave, Nmeas, Nmiss: %d %d %d\n", Nave, Nmeas, Nmiss);

  free (next);
  free (next_miss);
  free (ave_miss);

}

