# include "dvomerge.h"
# define PSPS_ID TRUE

# define IN_REGION(R,D) ( \
((D) >= region[0].Dmin) && ((D) < region[0].Dmax) && \
((R) >= region[0].Rmin)  && ((R) < region[0].Rmax))

// merge the input data into the output catalog
// input entries are matched by position to the new objects
int merge_catalogs_old (SkyRegion *region, Catalog *output, Catalog *input, double RADIUS) {

  off_t i, j, Nin, offset, J, Jmin, status, Nstars;
  double RADIUS2, Rmin, Rin, Din;
  double *X1, *Y1, *X2, *Y2;
  double dX, dY, dR;
  off_t *N1, *N2,  *next_meas;
  off_t Nave, NAVE, Nmeas, NMEAS, Nmatch;
  int Nsecfilt;
  unsigned int objID, catID;
  Coords tcoords;

  Nsecfilt = output[0].Nsecfilt;

  /** allocate local arrays (stars) **/
  ALLOCATE (X1, double, input[0].Naverage);
  ALLOCATE (Y1, double, input[0].Naverage);
  ALLOCATE (N1, off_t,  input[0].Naverage);
  ALLOCATE (input[0].found, off_t, input[0].Naverage);

  /** allocate local arrays (catalog) **/
  NAVE = Nave = output[0].Naverage;
  ALLOCATE (X2, double, NAVE);
  ALLOCATE (Y2, double, NAVE);
  ALLOCATE (N2, off_t,    NAVE);
  ALLOCATE (output[0].found, off_t, NAVE);
  /* for secfilt j and star i, secfilt[i*Nsecfilt+j] */

  /* internal counters */
  Nmatch = 0;
  NMEAS = Nmeas = output[0].Nmeasure;

  // current max obj ID for this catalog
  objID = output[0].objID;
  catID = output[0].catID;

  /* project onto rectilinear grid with 1 arcsec pixels. the choice of ZEA projection has the
   * advantage that every point in R,D has a mapping to a unique X,Y.  However, note that not all
   * possible X,Y points map back to R,D and the local plate scale changes substantially far from
   * the projection pole.  a better mapping might be ARC, not yet implemented (see
   * coordops.update.c).  We use the center of the region (catalog) for crval1,2. 
   */
  tcoords.crval1 = 0.5*(region[0].Rmin + region[0].Rmax);
  if (region[0].Dmax < 90) {
    tcoords.crval2 = 0.5*(region[0].Dmin + region[0].Dmax);
  } else {
    tcoords.crval2 = 90.0;
  }
  tcoords.crpix1 = 0;
  tcoords.crpix2 = 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 = 1;
  strcpy (tcoords.ctype, "RA---ARC");

  /* build spatial index (RA sort) referencing input array sequence */
  Nstars = 0;
  for (i = 0; i < input[0].Naverage; i++) {
    status = RD_to_XY (&X1[Nstars], &Y1[Nstars], input[0].average[i].R, input[0].average[i].D, &tcoords);
    if (!status) continue;
    N1[Nstars] = i;
    Nstars ++;
    input[0].found[i] = -1;
  }
  if (Nstars < 1) {
    if (VERBOSE) fprintf (stderr, "skipping %s, no overlapping stars\n", output[0].filename);
    free (output[0].found);
    free (input[0].found);
    free (X1);
    free (Y1);
    free (N1);
    free (X2);
    free (Y2);
    free (N2);
    return (0);
  }
  if (Nstars > 1) sort_coords_index (X1, Y1, N1, Nstars);

  /* build spatial index (RA sort) */
  for (i = 0; i < Nave; i++) {
    RD_to_XY (&X2[i], &Y2[i], output[0].average[i].R, output[0].average[i].D, &tcoords);
    output[0].found[i] = -1;
    N2[i] = i;
  }
  if (Nave > 1) sort_coords_index (X2, Y2, N2, Nave);

  /* set up pointers for linked list of measure */
  if (output[0].sorted && (output[0].Nmeasure == output[0].Nmeas_disk)) {
    // this version is only valid if we have done a full catalog load, and if the catalog
    // is sorted while processed
    next_meas = init_measure_links (output[0].average, Nave, output[0].measure, Nmeas);
  } else {
    next_meas = build_measure_links (output[0].average, Nave, output[0].measure, Nmeas);
  }    

  /* choose a radius for matches */
  RADIUS2 = RADIUS*RADIUS;

  /** find matched stars **/
  for (i = j = 0; (i < Nstars) && (j < Nave); ) {
    if (!finite(X1[i]) || !finite(Y1[i])) { 
      i++; 
      continue;
    }
    if (!finite(X2[j]) || !finite(Y2[j])) { 
      j++; 
      continue;
    }
    
    /* negative dX: j is too large */
    dX = X1[i] - X2[j];
    if (dX <= -1.02*RADIUS) {
      i++;
      continue;
    }
    /* positive dX, i is too large */
    if (dX >= 1.02*RADIUS) {
      j++;
      continue;
    }

    Jmin = -1;
    Rmin = RADIUS2;
    for (J = j; (dX > -1.02*RADIUS) && (J < Nave); J++) {
      /* find closest match for this detection */
      dX = X1[i] - X2[J];
      dY = Y1[i] - Y2[J];
      dR = dX*dX + dY*dY;
      if (dR > RADIUS2) continue;
      if (dR < Rmin) {
	Rmin = dR;
	Jmin  = J;
      }
    }

    /* no match, try next detection */ 
    if (Jmin == -1) {
      i++;
      continue;
    }

    /*** a match is found, add to average, measure ***/
    Nmatch ++;
    off_t n = N2[Jmin];
    off_t N = N1[i];

    /* make sure there is space for next Nmeasure entries */
    if (Nmeas + input[0].average[N].Nmeasure >= NMEAS) {
      NMEAS = Nmeas + input[0].average[N].Nmeasure + 1000;
      REALLOCATE (next_meas, off_t, NMEAS);
      REALLOCATE (output[0].measure, Measure, NMEAS);
    }

    /** add ALL measurements for this input average object **/
    for (Nin = 0; Nin < input[0].average[N].Nmeasure; Nin ++) {
      offset = input[0].average[N].measureOffset + Nin;

      /* add to end of measurement list */
      add_meas_link (&output[0].average[n], next_meas, Nmeas, NMEAS);
	
      // set the new measurements
      output[0].measure[Nmeas] = input[0].measure[offset];

      Rin = input[0].average[N].R - input[0].measure[offset].dR / 3600.0;
      Din = input[0].average[N].D - input[0].measure[offset].dD / 3600.0;

      /** dR,dD now represent arcsec **/
      output[0].measure[Nmeas].dR       = 3600.0*(output[0].average[n].R - Rin);  // XXX update these based on choice of astromety
      output[0].measure[Nmeas].dD       = 3600.0*(output[0].average[n].D - Din);  // XXX update these based on choice of astromety
      output[0].measure[Nmeas].dbFlags  = 0;  // XXX why reset these?
      output[0].measure[Nmeas].averef   = n;
      output[0].measure[Nmeas].objID    = output[0].average[n].objID;
      output[0].measure[Nmeas].catID    = output[0].catID;

      // fprintf (stderr, "Nave : %lld, Nmeas : %lld, dR: %f, dD: %f, catID: %d\n", (long long) n, (long long) Nmeas, output[0].measure[Nmeas].dR, output[0].measure[Nmeas].dD, output[0].measure[i].catID);

      // rationalize dR
      if (output[0].measure[Nmeas].dR > +180.0*3600.0) {
	// average on high end of boundary, move star up
	Rin += 360.0;
	output[0].measure[Nmeas].dR = 3600.0*(output[0].average[n].R - Rin);
      }
      if (output[0].measure[Nmeas].dR < -180.0*3600.0) {
	// average on low end of boundary, move star down
	Rin -= 360.0;
	output[0].measure[Nmeas].dR = 3600.0*(output[0].average[n].R - Rin);
      }
      if (fabs(output[0].measure[Nmeas].dR) > 10*RADIUS) {
	fprintf (stderr, "error: %10.6f,%10.6f vs %10.6f,%10.6f (%f,%f vs %f,%f)\n", 
		 output[0].average[n].R, output[0].average[n].D, Rin, Din,
		 X1[i], X2[Jmin], Y1[i], Y2[Jmin]);
	// XXX abort on this? -- this is a bad failure...
      }
      input[0].found[N] = Nmeas;
      output[0].average[n].Nmeasure ++;
      Nmeas ++;
    }

    /* Nm is updated, but not written out in -update mode (for existing entries)
       Nm is recalculated in build_meas_links if loaded table is not sorted */
    output[0].found[n] = Nmeas;
    i++;
  }

# if (0)
  fprintf (stderr, "--- 1 ---\n");
  for (i = 0; i < Nmeas; i++) {
    fprintf (stderr, "Nave : %d, Nmeas : %lld, dR: %f, dD: %f, catID: %d\n", output[0].measure[i].averef, (long long) i, output[0].measure[i].dR, output[0].measure[i].dD, output[0].measure[i].catID);
  }
# endif

  /** incorporate unmatched image stars, if this star is in field of this catalog **/
  /* these new entries are all written out in UPDATE mode */ 
  for (i = 0; i < Nstars; i++) {
    off_t N = N1[i];

    /* make sure there is space for next entry */
    if (Nmeas + input[0].average[N].Nmeasure >= NMEAS) {
      NMEAS = Nmeas + input[0].average[N].Nmeasure + 1000;
      REALLOCATE (next_meas, off_t, NMEAS);
      REALLOCATE (output[0].measure, Measure, NMEAS);
    }
    if (Nave >= NAVE) {
      NAVE = Nave + 1000;
      REALLOCATE (output[0].average, Average, NAVE);
      REALLOCATE (output[0].secfilt, SecFilt, NAVE*output[0].Nsecfilt);
    }

    if (input[0].found[N] >= 0) continue;
    if (!IN_REGION (input[0].average[N].R, input[0].average[N].D)) continue;

    output[0].average[Nave].R         	   = input[0].average[N].R;
    output[0].average[Nave].D         	   = input[0].average[N].D;
    output[0].average[Nave].dR        	   = 0;
    output[0].average[Nave].dD        	   = 0;

    output[0].average[Nave].Nmeasure  	   = 1;
    output[0].average[Nave].Nmissing  	   = 0;
    output[0].average[Nave].Nextend        = 0;

    output[0].average[Nave].measureOffset  = Nmeas;
    output[0].average[Nave].missingOffset  = -1;
    output[0].average[Nave].extendOffset   = -1;

    output[0].average[Nave].uR        	   = 0;
    output[0].average[Nave].uD        	   = 0;
    output[0].average[Nave].duR       	   = 0;
    output[0].average[Nave].duD       	   = 0;
    output[0].average[Nave].P         	   = 0;
    output[0].average[Nave].dP        	   = 0;

    output[0].average[Nave].Xp        	   = 0;
    output[0].average[Nave].ChiSqAve   	   = 0.0;
    output[0].average[Nave].ChiSqPM   	   = 0.0;
    output[0].average[Nave].ChiSqPar   	   = 0.0;
    output[0].average[Nave].Tmean   	   = 0;
    output[0].average[Nave].Trange   	   = 0;
    output[0].average[Nave].Npos    	   = 0;

    output[0].average[Nave].objID     	   = objID; // we create objID values in the context of the output db
    output[0].average[Nave].catID     	   = catID; // we create catID values in the context of the output db
    output[0].average[Nave].flags          = 0;  // XXX why reset these?
    if (PSPS_ID) {
      output[0].average[Nave].extID = CreatePSPSObjectID(output[0].average[Nave].R, output[0].average[Nave].D);
    } else {
      output[0].average[Nave].extID = 0;
    }

    objID ++;

    for (j = 0; j < Nsecfilt; j++) {
      output[0].secfilt[Nave*Nsecfilt+j].M  = NAN;
      output[0].secfilt[Nave*Nsecfilt+j].dM = NAN;
      output[0].secfilt[Nave*Nsecfilt+j].Xm = NAN_S_SHORT;
      output[0].secfilt[Nave*Nsecfilt+j].M_20 	= NAN_S_SHORT;
      output[0].secfilt[Nave*Nsecfilt+j].M_80 	= NAN_S_SHORT;
      output[0].secfilt[Nave*Nsecfilt+j].Ncode = 0;
      output[0].secfilt[Nave*Nsecfilt+j].Nused = 0;
    }

    /** add measurements for this input average object **/
    for (Nin = 0; Nin < input[0].average[N].Nmeasure; Nin ++) {
      offset = input[0].average[N].measureOffset + Nin;

      // supply the measurments from this detection
      output[0].measure[Nmeas]           = input[0].measure[offset];

      // the following measure elements cannot be set until here:
      output[0].measure[Nmeas].dbFlags  = 0;
      output[0].measure[Nmeas].averef   = Nave;
      output[0].measure[Nmeas].objID    = output[0].average[Nave].objID;
      output[0].measure[Nmeas].catID    = output[0].catID;

      /* next[Nmeas] should always be -1 in this context (it is always the only
	 measurement for the star) */
      input[0].found[N] = Nmeas;
      next_meas[Nmeas] = -1;
      Nmeas ++;
    }
    int Ngroup = input[0].average[N].Nmeasure;
    for (j = 0; j < Ngroup - 1; j++) {
	next_meas[Nmeas - Ngroup + j] = Nmeas - Ngroup + j + 1;
    }
    Nave ++;
  }
      
  REALLOCATE (output[0].average, Average, Nave);
  REALLOCATE (output[0].measure, Measure, Nmeas);
 
# if (0)
  fprintf (stderr, "--- 2 ---\n");
  for (i = 0; i < Nmeas; i++) {
    fprintf (stderr, "Nave : %d, Nmeas : %lld, dR: %f, dD: %f, catID: %d\n", output[0].measure[i].averef, (long long) i, output[0].measure[i].dR, output[0].measure[i].dD, output[0].measure[i].catID);
  }
# endif

# define NOSORT FALSE
  if (NOSORT) {
    output[0].sorted = FALSE;
  } else {
    output[0].sorted = TRUE;
    output[0].measure = sort_measure (output[0].average, Nave, output[0].measure, Nmeas, next_meas);
  }

# if (0)
  fprintf (stderr, "--- 3 ---\n");
  for (i = 0; i < Nmeas; i++) {
    fprintf (stderr, "Nave : %d, Nmeas : %lld, dR: %f, dD: %f, catID: %d\n", output[0].measure[i].averef, (long long) i, output[0].measure[i].dR, output[0].measure[i].dD, output[0].measure[i].catID);
  }
# endif

  /* note stars which have been found in this catalog */
  for (i = 0; i < input[0].Naverage; i++) {
    if (input[0].found[i] > -1) {
      input[0].found[i] = -2;
    } else {
      input[0].found[i] = -3;
    }
  }

  /* check if the catalog has changed?  if no change, no need to write */
  output[0].objID    = objID; // new max value, save on catalog close
  output[0].Naverage = Nave;
  output[0].Nmeasure = Nmeas;
  output[0].Nsecf_mem = Nave*Nsecfilt;
  if (VERBOSE) fprintf (stderr, "Nstars, Nave, Nmeas: %lld %lld %lld, (%lld matches)\n", (long long) Nstars, (long long) Nave, (long long) Nmeas, (long long) Nmatch);

  free (next_meas);

  free (output[0].found);
  free (X2);
  free (Y2);
  free (N2);
  free (input[0].found);
  free (X1);
  free (Y1);
  free (N1);
  return (Nmatch);
}

/* 
   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
   
   output[0].found[Ncat] but stars[Nstar].found
   
*/

