# include "relastro.h"

# define NEXT_I { i++; continue; }
# define NEXT_J { j++; continue; }

// XXX these were used to spot-check specific objects
# define REF_RA1 192.455781 
# define REF_DEC1 2.563102 
# define REF_RA2 192.444302113 
# define REF_DEC2 2.5786559112 
int high_speed_objects (SkyRegion *region, Catalog *catalog) {

  // XXX seems silly to require region here just to find the catalog bounds / center

  off_t i, j, m, J, ni, nj, *N1, Nslow, NgroupA, NgroupB;
  int *slowMoving, *groupA, *groupB, status, foundA, foundB, Nmatch, valid;
  double *X1, *Y1;
  double dX, dY, dR, RADIUS2, yJ;
  Coords tcoords;

  int zcode, zNsec, ycode, yNsec, jcode, jNsec, hcode, hNsec, kcode, kNsec, USNO_R, USNO_N, Nsecfilt;

  Nsecfilt = GetPhotcodeNsecfilt();

  ycode = GetPhotcodeCodebyName("y");
  yNsec = GetPhotcodeNsec(ycode);
  
  zcode = GetPhotcodeCodebyName("z");
  zNsec = GetPhotcodeNsec(zcode);
  
  jcode = GetPhotcodeCodebyName("J");
  jNsec = GetPhotcodeNsec(jcode);
  
  hcode = GetPhotcodeCodebyName("H");
  hNsec = GetPhotcodeNsec(hcode);
  
  kcode = GetPhotcodeCodebyName("K");
  kNsec = GetPhotcodeNsec(kcode);
  
  USNO_R = GetPhotcodeCodebyName("USNO_RED");
  USNO_N = GetPhotcodeCodebyName("USNO_N");

  // high-speed between different surveys (easier case):
  // we have two sets of photcodes (A) and (B) and are looking for objects
  // with detections in only (A) and separately only (B).

  // we need at least 2 objects if we are going to match anything...
  if (catalog[0].Naverage < 2) return (TRUE);

  // mask with which to mark objects to be ignored
  ALLOCATE (slowMoving, int, catalog[0].Naverage);
  memset (slowMoving, 0, catalog[0].Naverage*sizeof(int));

  // record to which photcode group the object belongs:
  NgroupA = 0;
  ALLOCATE (groupA, int, catalog[0].Naverage);
  memset (groupA, 0, catalog[0].Naverage*sizeof(int));

  NgroupB = 0;
  ALLOCATE (groupB, int, catalog[0].Naverage);
  memset (groupB, 0, catalog[0].Naverage*sizeof(int));

  fprintf (stderr, "checking %lld objects\n", (long long) catalog[0].Naverage);
  Nslow = 0;

  fprintf (stdout, "#      RA_A  DEC_A              RA_B  DEC_B      :     pmx     pmy      dR  :    z      dz       y      dy       J      dJ       H      dH       K      dK\n");
  //................270.2346670  20.2471540  270.2170434  20.2717396 :  -59.11   88.78   106.66 :   0.000  0.000   19.582  0.112   16.041  0.110   15.388  0.098   14.858  0.001

  // mark (exclude) objects with both sets of target photcodes
  for (i = 0; i < catalog[0].Naverage; i++) {

    // if (i % 100 == 0) fprintf (stderr, ".");

    // radius = hypot((REF_DEC1 - catalog[0].average[i].D), (REF_RA1 - catalog[0].average[i].R));
    // if (radius < 0.0001) {
    //   fprintf (stderr, "found it\n");
    // }
    // radius = hypot((REF_DEC2 - catalog[0].average[i].D), (REF_RA2 - catalog[0].average[i].R));
    // if (radius < 0.0001) {
    //   fprintf (stderr, "found it\n");
    // }

    // do any of the measures for this object match group A?
    m = catalog[0].average[i].measureOffset;
    foundA = FALSE;
    for (j = 0; !foundA && (j < catalog[0].average[i].Nmeasure); j++, m++) {
      if (MeasMatchesPhotcode(&catalog[0].measure[m], photcodesGroupA, NphotcodesGroupA)) {
	foundA = TRUE;
      }
    }

    // do any of the measures for this object match group B?
    m = catalog[0].average[i].measureOffset;
    foundB = FALSE;
    for (j = 0; !foundB && (j < catalog[0].average[i].Nmeasure); j++, m++) {
      if (MeasMatchesPhotcode(&catalog[0].measure[m], photcodesGroupB, NphotcodesGroupB)) {
	foundB = TRUE;
      }
    }

    if (foundA && foundB) {
      slowMoving[i] = TRUE;
      Nslow ++;
      continue;
    }

    // additional constraints:
    // * group A : require 

    // 2MASS detections:
    if (foundA && !foundB) {
      // average-based tests:
      valid = TRUE;
      valid &= ((catalog[0].average[i].flags & 0x40000) == 0);
      valid &= ((catalog[0].average[i].flags & 0x80000)  > 0);
      valid &= (catalog[0].secfilt[i*Nsecfilt + yNsec].M < 1.0);
      valid &= (catalog[0].secfilt[i*Nsecfilt + zNsec].M < 1.0);
      // XXX color restrictions are applied in the pair-wise matching below
      
      // measure-based tests:
      m = catalog[0].average[i].measureOffset;
      for (j = 0; valid && (j < catalog[0].average[i].Nmeasure); j++, m++) {
	if (catalog[0].measure[m].photcode == USNO_R) {
	  valid &= ((catalog[0].measure[m].M > 20.0) || isnan(catalog[0].measure[m].M));
	}
	if (catalog[0].measure[m].photcode == USNO_N) {
	  valid &= ((catalog[0].measure[m].M > 18.0) || isnan(catalog[0].measure[m].M));
	}
      }

      // ((objflags & 524288) == 524288) && 
      // ((objflags & 262144) == 0) && 
      // (z<1.0) && 
      // (y<1.0) && 
      // ((yp[$i]-J)>1.25) && ((yp[$i]-J)<5.0) && 
      // ((USNO_RED<1.0)||(USNO_RED>20.0)) && 
      // ((USNO_N>18.0)||(USNO_N<1.0))

      if (valid) {
	groupA[i] = TRUE;
	NgroupA ++;
	continue;
      } else {
	// fprintf (stderr, "skip %lld (group A)\n", (long long) i);
      }
    }

    // PS1 detections:
    if (foundB && !foundA) {

      // average-based tests:
      valid = TRUE;
      valid &= ((catalog[0].average[i].flags & 0x10000) == 0);
      valid &= ((catalog[0].average[i].flags & 0x20000)  > 0);
      valid &= (catalog[0].secfilt[i*Nsecfilt + jNsec].M < 1.0);
      valid &= (catalog[0].secfilt[i*Nsecfilt + yNsec].Nused > 1);
      valid &= (catalog[0].secfilt[i*Nsecfilt + yNsec].dM < 0.2);

      if ((catalog[0].secfilt[i*Nsecfilt + zNsec].M < 1.0) || (catalog[0].secfilt[i*Nsecfilt + zNsec].Nused == 1)) {
	valid &= TRUE;
      } else {
	valid &= ((catalog[0].secfilt[i*Nsecfilt + zNsec].M - catalog[0].secfilt[i*Nsecfilt + yNsec].M) > 0.2);
      }

      // measure-based tests:
      m = catalog[0].average[i].measureOffset;
      for (j = 0; valid && (j < catalog[0].average[i].Nmeasure); j++, m++) {
	if (catalog[0].measure[m].photcode == USNO_R) {
	  valid &= ((catalog[0].measure[m].M > 20.0) || isnan(catalog[0].measure[m].M));
	}
	if (catalog[0].measure[m].photcode == USNO_N) {
	  valid &= ((catalog[0].measure[m].M > 18.5) || isnan(catalog[0].measure[m].M));
	}
      }

      // (abs(glat)>15.0) && 
      // ((objflags & 65536) == 0) && 
      // ((objflags & 131072) == 131072) &&
      // (J<1.0) && 
      // (y:nphot>1) && 
      // (y:err<0.2) && 
      // ((USNO_RED>20.0)||(USNO_RED<1.0)) && 
      // ((USNO_N>18.5)||(USNO_N<1.0)) && 
      // ((z<1.0)||(z:nphot == 1)||((z:nphot>1)&&((z-y)>0.2)))

      if (valid) {
	groupB[i] = TRUE;
	NgroupB ++;
	continue;
      } else {
	// fprintf (stderr, "skip %lld (group B)\n", (long long) i);
      }
    }

    // this object does not have a detection from either groupA or groupB -- skip as if slow
    slowMoving[i] = TRUE;
    Nslow ++;
  }
  fprintf (stderr, "%lld of %lld are slow or invalid; %lld in group A, %lld in group B\n", (long long) Nslow, (long long) catalog[0].Naverage, (long long) NgroupA, (long long) NgroupB);
    
  // double loop over unmarked objects (sorted in RA / X)
  // record all pairs within the desired match distance

  // create tmp local positional index
  ALLOCATE (X1, double, catalog[0].Naverage);
  ALLOCATE (Y1, double, catalog[0].Naverage);
  ALLOCATE (N1, off_t,  catalog[0].Naverage);

  // define a local projection
  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 */
  for (i = 0; i < catalog[0].Naverage; i++) {
    status = RD_to_XY (&X1[i], &Y1[i], catalog[0].average[i].R, catalog[0].average[i].D, &tcoords);
    N1[i] = i;
    assert (status);
  }
  sort_coords_index (X1, Y1, N1, catalog[0].Naverage);

  Nmatch = 0;
  RADIUS2 = SQ(RADIUS);

  // mark (exclude) objects with both sets of target photcodes
  for (i = j = 0; (i < catalog[0].Naverage) && (j < catalog[0].Naverage);) {

    ni = N1[i];
    nj = N1[j];

    // if (ni == 131030) {
    //   fprintf (stderr, "got 2mass\n");
    // }
    // if (nj == 53581) {
    //   fprintf (stderr, "got ps1\n");
    // }

    if (slowMoving[ni]) NEXT_I;
    if (slowMoving[nj]) NEXT_J;

    // i => groupA, j => groupB
    if (!groupA[ni]) NEXT_I;
    if (!groupB[nj]) NEXT_J;

    if (!finite(X1[i]) || !finite(Y1[i])) NEXT_I;
    if (!finite(X1[j]) || !finite(Y1[j])) NEXT_J;

    // look for pairs that are within the maximum separation
    dX = X1[i] - X1[j];

    if (dX <= -1.02*RADIUS) NEXT_I; // negative dX: i is too small
    if (dX >= +1.02*RADIUS) NEXT_J; // positive dX: j is too small

    // within match range; look for valid matches
    for (J = j; (dX > -1.02*RADIUS) && (J < catalog[0].Naverage); J++) {
      if (J == i) continue;  // avoid auto-matches

      dX = X1[i] - X1[J];

      nj = N1[J];
      if (!groupB[nj]) continue;

      dY = Y1[i] - Y1[J];
      dR = dX*dX + dY*dY;
      if (dR > RADIUS2) continue;

      // y_groupA - J_groupB 
      yJ = catalog[0].secfilt[nj*Nsecfilt + yNsec].M - catalog[0].secfilt[ni*Nsecfilt + jNsec].M;
      if (yJ < 1.25) continue;
      if (yJ > 5.0) continue;

      /*** a match is found (just print it for the moment) ***/
      Nmatch ++;

      // for the moment, just write the matches to stdout:
      float pmx = -dX; // proper-motion displacement in ra  direction (B - A) [arcsec]
      float pmy = -dY; // proper-motion displacement in dec direction (B - A) [arcsec]
      fprintf (stdout, "%11.7f %11.7f  %11.7f %11.7f : %7.2f %7.2f  %7.2f : %7.3f %6.3f  %7.3f %6.3f  %7.3f %6.3f  %7.3f %6.3f  %7.3f %6.3f\n", 
	       catalog[0].average[ni].R, catalog[0].average[ni].D,
	       catalog[0].average[nj].R, catalog[0].average[nj].D,
	       pmx, pmy, sqrt(dR),
	       catalog[0].secfilt[nj*Nsecfilt + zNsec].M,
	       catalog[0].secfilt[nj*Nsecfilt + zNsec].dM,
	       catalog[0].secfilt[nj*Nsecfilt + yNsec].M, 
	       catalog[0].secfilt[nj*Nsecfilt + yNsec].dM,
	       catalog[0].secfilt[ni*Nsecfilt + jNsec].M,
	       catalog[0].secfilt[ni*Nsecfilt + jNsec].dM,
	       catalog[0].secfilt[ni*Nsecfilt + hNsec].M,
	       catalog[0].secfilt[ni*Nsecfilt + hNsec].dM,
	       catalog[0].secfilt[ni*Nsecfilt + kNsec].M,
	       catalog[0].secfilt[ni*Nsecfilt + kNsec].dM);
    }
    i++;
  }
  fprintf (stderr, "found %d matches\n", Nmatch);
  
  free (slowMoving);
  free (groupA);
  free (groupB);
  free (X1);
  free (Y1);
  free (N1);

  return TRUE;
}

// return TRUE if any measure->photcodes match the photcodeSet
int MeasMatchesPhotcode(Measure *measure, PhotCode **photcodeSet, int Nset) {

  int found, k;

  assert (Nset > 0);

  if (!finite(measure[0].dR) || !finite(measure[0].dD)) return FALSE;
  if (!finite(measure[0].M)) return FALSE;
  
  /* select measurements by photcode, or equiv photcode, if specified */
  if (Nset > 0) {
    found = FALSE;
    for (k = 0; (k < Nset) && !found; k++) {
      if (photcodeSet[k][0].code == measure[0].photcode) found = TRUE;
      if (photcodeSet[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
    }
    if (!found) return FALSE;
  }
  
  return TRUE;
}
