# include "photdbc.h"

void join_stars (Catalog *catalog) {

  int i, j, k, m, M, Ni, Nj, first_j;
  int Naves, Nmeas, Ncurr;
  int Naverage, Nmeasure, *found, *index;
  float *X, *Y, dX, dY, dR, RADIUS2;
  double Sr, Sd, Rmid, Dmid;

  Average *naverage, *average;
  Measure *nmeasure, *measure;
  SecFilt *nsecfilt, *secfilt;
  Mpointer *mpointer;
  Coords tcoords;

  if (VERBOSE) fprintf (stderr, "joining overlapping stars\n");

  average = catalog[0].average;
  measure = catalog[0].measure;
  secfilt = catalog[0].secfilt;

  Naverage = catalog[0].Naverage;
  Nmeasure = catalog[0].Nmeasure;
  RADIUS2 = JOIN_RADIUS*JOIN_RADIUS;

  /* reference for coords is this region */
  Rmid = Dmid = 0;
  for (i = 0; i < Naverage; i++) {
    Rmid += average[i].R_PS;
    Dmid += average[i].D_PS;
  }
  Rmid /= Naverage;
  Dmid /= Naverage;

  /* coordinate system for projection */
  tcoords.crval1 = Rmid;
  tcoords.crval2 = Dmid;
  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;
  tcoords.Npolyterms = 0;
  strcpy (tcoords.ctype, "RA---TAN");
  
  /* project & sort coordinates in local linear frame */
  ALLOCATE (X, float, Naverage);
  ALLOCATE (Y, float, Naverage);
  ALLOCATE (index, int, Naverage);
  for (i = 0; i < Naverage; i++) {
    index[i] = i;
    fRD_to_XY (&X[i], &Y[i], average[i].R_PS, average[i].D_PS, &tcoords);
  }
  sort_lists (X, Y, index, Naverage);
  
  /* flags to mark if a star has been handled */
  ALLOCATE (found, int, Naverage);
  for (i = 0; i < Naverage; i++) found[i] = FALSE;

  ALLOCATE (naverage, Average, Naverage);
  ALLOCATE (mpointer, Mpointer, Nmeasure);
  for (i = 0; i < Nmeasure; i++) mpointer[i].averef = -1;

  Naves =  0; // counter for new averages
  Nmeas =  0; // counter for new measures
  for (i = j = 0; (i < Naverage) && (j < Naverage); ) {
    
    Ni = index[i];
    Nj = index[j];

    /* a new star, add it to naverage[] */
    if (!found[i]) { 	
      naverage[Naves]        = average[Ni];
      naverage[Naves].offset = Nmeas;
      for (k = 0; k < average[Ni].Nm; k++) {
	m = average[Ni].offset + k;
	mpointer[Nmeas].measure = m;
	mpointer[Nmeas].averef  = Naves;
	mpointer[Nmeas].R       = average[Ni].R_PS - measure[m].dR_PS / 3600.0;
	mpointer[Nmeas].D       = average[Ni].D_PS - measure[m].dD_PS / 3600.0;
	Nmeas ++;
      }
      Ncurr = Naves;
      found[i] = TRUE;
      Naves ++;
    }

    if (found[j]) { j++; continue; }  // don't duplicate
    if (j == i) {   j++; continue; }  // don't auto-correlate

    dX = X[i] - X[j];
    if (dX <= -JOIN_RADIUS) { /* X[j] is too large */
      while (found[i] && (i < Naverage)) i++;
      continue;
    }
    if (dX >= +JOIN_RADIUS) { /* X[i] is too large */
      j++;
      continue;
    }

    first_j = j;
    for (; (dX > -2*JOIN_RADIUS) && (j < Naverage); j++) {
      Nj = index[j];
      if (found[j]) continue;
      dX = X[i] - X[j];
      dY = Y[i] - Y[j];
      dR = dX*dX + dY*dY;
      if (dR < RADIUS2) {  /* matched star, join to first */

	/* define pointers for new measures */
	for (k = 0; k < average[Nj].Nm; k++) {
	  m = average[Nj].offset + k;
	  mpointer[Nmeas].measure = m;
	  mpointer[Nmeas].averef  = Ncurr;
	  mpointer[Nmeas].R       = average[Nj].R_PS - measure[m].dR_PS / 3600.0;
	  mpointer[Nmeas].D       = average[Nj].D_PS - measure[m].dD_PS / 3600.0;
	  Nmeas ++;
	}
	naverage[Ncurr].Nm += average[Nj].Nm;
	found[j] = TRUE;
	
	/* recalculate naverage[Ncurr].RA,DEC */
	/* this must be done here to keep the average position consistent
	   for the next star found */
	for (Sr = Sd = k = 0; k < naverage[Ncurr].Nm; k++) {
	  m = naverage[Ncurr].offset + k;
	  Sr += mpointer[m].R;
	  Sd += mpointer[m].D;
	}
	Sr = Sr / naverage[Ncurr].Nm;
	Sd = Sd / naverage[Ncurr].Nm;
	naverage[Ncurr].R_PS = Sr;
	naverage[Ncurr].D_PS = Sd;

	/* update original measurement offsets */
	for (k = 0; k < naverage[Ncurr].Nm; k++) {
	  m = naverage[Ncurr].offset + k;
	  M = mpointer[m].measure;
	  measure[M].dR_PS = 3600.0*(Sr - mpointer[m].R);
	  measure[M].dD_PS = 3600.0*(Sd - mpointer[m].D);
	}

	/* update current reference star position */
	fRD_to_XY (&X[i], &Y[i], Sr, Sd, &tcoords);
      }
    }
    while (found[i] && (i < Naverage)) i++;
    j = first_j;
  }
  if (Nmeas != Nmeasure) {
    fprintf (stderr, "failure to match measures\n");
  }
  
  /* create a new Measure table in the appropriate sequence */
  ALLOCATE (nmeasure, Measure, Nmeasure);
  for (i = 0; i < Nmeasure; i++) {
    nmeasure[i]        = measure[mpointer[i].measure];
    nmeasure[i].averef = mpointer[i].averef;
  }

  /* create an empty SecFilt table */
  ALLOCATE (nsecfilt, SecFilt, Naves*catalog[0].Nsecfilt);
  bzero (nsecfilt, Naves*catalog[0].Nsecfilt*sizeof(SecFilt));

  free (catalog[0].average);
  free (catalog[0].measure);
  free (catalog[0].secfilt);

  catalog[0].average = naverage;
  catalog[0].measure = nmeasure;
  catalog[0].secfilt = nsecfilt;
  
  catalog[0].Naverage = Naves;
}

/* notes:
   merge averages of stars within RADIUS to single stars 
   updates the average RA and DEC fields
   creates a new Average table 'naverage' 
   updates the Measure table to the new sequence
   generates a new SecFilt table with empty values 
   does NOT update magnitudes 
*/   
   
