Index: trunk/Ohana/src/relastro/Makefile
===================================================================
--- trunk/Ohana/src/relastro/Makefile	(revision 28043)
+++ trunk/Ohana/src/relastro/Makefile	(revision 28184)
@@ -56,6 +56,8 @@
 $(SRC)/save_catalogs.$(ARCH).o       \
 $(SRC)/write_coords.$(ARCH).o        \
-$(SRC)/CoordOps.$(ARCH).o        \
-$(SRC)/FixProblemImages.$(ARCH).o        \
+$(SRC)/CoordOps.$(ARCH).o            \
+$(SRC)/FixProblemImages.$(ARCH).o    \
+$(SRC)/high_speed_catalogs.$(ARCH).o  \
+$(SRC)/high_speed_objects.$(ARCH).o  \
 $(SRC)/relastroVisual.$(ARCH).o
 
Index: trunk/Ohana/src/relastro/doc/high-speed.txt
===================================================================
--- trunk/Ohana/src/relastro/doc/high-speed.txt	(revision 28184)
+++ trunk/Ohana/src/relastro/doc/high-speed.txt	(revision 28184)
@@ -0,0 +1,7 @@
+
+Niall Deacon's 2MASS / PS1 high-speed search uses the following set of criteria (roughly):
+
+2 sets of detections:
+
+A : 2MASS only, J significant, H & K anything, set of phot_flags?
+B : PS1 only, 2 y-band, 2 z-band, S/N > 5, not a galaxy 
Index: trunk/Ohana/src/relastro/include/relastro.h
===================================================================
--- trunk/Ohana/src/relastro/include/relastro.h	(revision 28043)
+++ trunk/Ohana/src/relastro/include/relastro.h	(revision 28184)
@@ -14,5 +14,5 @@
 typedef enum {FIT_NONE, FIT_AVERAGE, FIT_PM_ONLY, FIT_PAR_ONLY, FIT_PM_AND_PAR} FitMode;
 
-typedef enum {TARGET_NONE, TARGET_OBJECTS, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS} FitTarget;
+typedef enum {TARGET_NONE, TARGET_OBJECTS, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS, TARGET_HIGH_SPEED} FitTarget;
 
 typedef struct {
@@ -89,4 +89,6 @@
 int SRC_MEAS_TOOFEW; //catalog objects wich fewer detections then this are ignored
 double MIN_ERROR;
+
+double RADIUS; // match radius for high-speed objects
 
 int    VERBOSE;
@@ -112,4 +114,8 @@
 int           NphotcodesKeep,      NphotcodesSkip;
 PhotCode     **photcodesKeep,     **photcodesSkip;
+
+char          *PHOTCODE_A_LIST,  *PHOTCODE_B_LIST;
+int           NphotcodesGroupA,  NphotcodesGroupB;
+PhotCode     **photcodesGroupA, **photcodesGroupB;
 
 int AreaSelect;
@@ -321,2 +327,6 @@
 int saveCoords (Coords *coords, off_t N);
 void resetImageRaw (Catalog *catalog, int Ncatalog, off_t im);
+
+int high_speed_catalogs ();
+int high_speed_objects (SkyRegion *region, Catalog *catalog);
+int MeasMatchesPhotcode(Measure *measure, PhotCode **photcodeSet, int Nset);
Index: trunk/Ohana/src/relastro/src/args.c
===================================================================
--- trunk/Ohana/src/relastro/src/args.c	(revision 28043)
+++ trunk/Ohana/src/relastro/src/args.c	(revision 28184)
@@ -28,4 +28,17 @@
     }
   }
+  if ((N = get_argument (argc, argv, "-high-speed"))) {
+    // XXX include a parallax / no-parallax option
+    if (N >= argc - 3) usage();
+    FIT_TARGET = TARGET_HIGH_SPEED;
+    remove_argument (N, &argc, argv);
+    PHOTCODE_A_LIST = strcreate(argv[N]);
+    remove_argument (N, &argc, argv);
+    PHOTCODE_B_LIST = strcreate(argv[N]);
+    remove_argument (N, &argc, argv);
+    RADIUS = atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
   if ((N = get_argument (argc, argv, "-update-simple"))) {
     remove_argument (N, &argc, argv);
@@ -271,6 +284,6 @@
 void usage () {
   fprintf (stderr, "ERROR: USAGE: relastro -region RA RA DEC DEC\n");
-  fprintf (stderr, "       OR:    relastro -catalog (ra) (dec)\n");
-  fprintf (stderr, "  working options: \n");
+  fprintf (stderr, "       OR:    relastro -catalog (ra) (dec)\n\n");
+  fprintf (stderr, "  specify one of the following modes: \n");
   fprintf (stderr, "  -update-objects\n");
   fprintf (stderr, "    -pm\n");
@@ -280,4 +293,6 @@
   fprintf (stderr, "  -update-chips\n");
   fprintf (stderr, "  -update-mosaics\n");
+  fprintf (stderr, "  -high-speed (code[,code,code]) (code[,code,code]) (radius)\n\n");
+  fprintf (stderr, "  additional options: \n");
   fprintf (stderr, "  -time (start)(stop)\n");
   fprintf (stderr, "  +photcode (code)[,code,code...]\n");
Index: trunk/Ohana/src/relastro/src/high_speed_catalogs.c
===================================================================
--- trunk/Ohana/src/relastro/src/high_speed_catalogs.c	(revision 28184)
+++ trunk/Ohana/src/relastro/src/high_speed_catalogs.c	(revision 28184)
@@ -0,0 +1,55 @@
+# include "relastro.h"
+
+int high_speed_catalogs () {
+
+  int i;
+
+  SkyTable *sky = NULL;
+  SkyList *skylist = NULL;
+  Catalog catalog;
+
+  // load the current sky table (layout of all SkyRegions) 
+  sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE);
+  SkyTableSetFilenames (sky, CATDIR, "cpt");
+  
+  // determine the populated SkyRegions overlapping the requested area (default depth)
+  if (UserCatalog) {
+    skylist = SkyRegionByPoint (sky, -1, UserCatalogRA, UserCatalogDEC);
+  } else {
+    skylist = SkyListByPatch (sky, -1, &UserPatch);
+  }
+
+  // load data from each region file, only use bright stars
+  for (i = 0; i < skylist[0].Nregions; i++) {
+
+    // set up the basic catalog info
+    catalog.filename  = skylist[0].filename[i];
+    catalog.catformat = dvo_catalog_catformat (CATFORMAT);    // set the default catformat from config data
+    catalog.catmode   = dvo_catalog_catmode (CATMODE);        // set the default catmode from config data
+    catalog.catflags  = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
+    catalog.Nsecfilt  = GetPhotcodeNsecfilt ();
+
+    if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) {
+      fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog.filename);
+      exit (1);
+    }
+    if (!catalog.Naves_disk) {
+      if (VERBOSE) fprintf (stderr, "no data in %s, skipping\n", catalog.filename);
+      dvo_catalog_unlock (&catalog);
+      dvo_catalog_free (&catalog);
+      continue;
+    }
+
+    high_speed_objects (skylist[0].regions[i], &catalog);
+
+    if (!UPDATE) {
+      dvo_catalog_unlock (&catalog);
+      dvo_catalog_free (&catalog);
+      continue;
+    }
+    
+    save_catalogs (&catalog, 1);
+  }
+  
+  return (TRUE);
+}
Index: trunk/Ohana/src/relastro/src/high_speed_objects.c
===================================================================
--- trunk/Ohana/src/relastro/src/high_speed_objects.c	(revision 28184)
+++ trunk/Ohana/src/relastro/src/high_speed_objects.c	(revision 28184)
@@ -0,0 +1,191 @@
+# include "relastro.h"
+
+# define NEXT_I { i++; continue; }
+# define NEXT_J { j++; continue; }
+
+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, ni, nj, Nj, *N1, Nslow, NgroupA, NgroupB;
+  int *slowMoving, *groupA, *groupB, status, foundA, foundB, Nmatch;
+  double *X1, *Y1;
+  double dX, dY, dR, RADIUS2;
+  Coords tcoords;
+
+  // 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;
+
+  // mark (exclude) objects with both sets of target photcodes
+  for (i = 0; i < catalog[0].Naverage; i++) {
+
+    if (i % 100 == 0) fprintf (stderr, ".");
+
+    // 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 
+
+    if (foundA && !foundB) {
+      groupA[i] = TRUE;
+      NgroupA ++;
+      continue;
+    }
+
+    if (!foundA && foundB) {
+      groupB[i] = TRUE;
+      NgroupB ++;
+      continue;
+    }
+
+    // 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; %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 (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];
+      dY = Y1[i] - Y1[J];
+      dR = dX*dX + dY*dY;
+      if (dR > RADIUS2) continue;
+
+      /*** a match is found (just print it for the moment) ***/
+      Nmatch ++;
+      Ni = N1[i];
+      Nj = N1[J];
+
+      if (0) fprintf (stderr, "match %d = %lld %lld : %f %f : %f %f : %f\n", 
+	       Nmatch, (long long) i, (long long) J,
+	       catalog[0].average[Ni].R, catalog[0].average[Ni].D,
+	       catalog[0].average[Nj].R, catalog[0].average[Nj].D,
+	       sqrt(dR));
+    }
+    i++;
+  }
+  fprintf (stderr, "found %d matches\n", Nmatch);
+  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;
+}
Index: trunk/Ohana/src/relastro/src/initialize.c
===================================================================
--- trunk/Ohana/src/relastro/src/initialize.c	(revision 28043)
+++ trunk/Ohana/src/relastro/src/initialize.c	(revision 28184)
@@ -56,4 +56,46 @@
   }
 
+  NphotcodesGroupA = 0;
+  photcodesGroupA = NULL;
+  if (PHOTCODE_A_LIST != NULL) {
+    NPHOTCODES = 10;
+    ALLOCATE (photcodesGroupA, PhotCode *, NPHOTCODES);
+
+    /* parse the comma-separated list of photcodesGroupA */
+    list = PHOTCODE_A_LIST;
+    while ((codename = strtok_r (list, ",", &ptr)) != NULL) {
+      list = NULL; // pass NULL on successive strtok_r calls
+      fprintf (stderr, "PHOTCODE_A_LIST: %s\n", PHOTCODE_A_LIST);
+      fprintf (stderr, "codename: %s\n", codename);
+      if ((photcodesGroupA[NphotcodesGroupA] = GetPhotcodebyName (codename)) == NULL) {
+        fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename);
+        exit (1);
+      }
+      NphotcodesGroupA ++;
+      CHECK_REALLOCATE (photcodesGroupA, PhotCode *, NPHOTCODES, NphotcodesGroupA, 10);
+    }
+  }
+
+  NphotcodesGroupB = 0;
+  photcodesGroupB = NULL;
+  if (PHOTCODE_B_LIST != NULL) {
+    NPHOTCODES = 10;
+    ALLOCATE (photcodesGroupB, PhotCode *, NPHOTCODES);
+
+    /* parse the comma-separated list of photcodesGroupB */
+    list = PHOTCODE_B_LIST;
+    while ((codename = strtok_r (list, ",", &ptr)) != NULL) {
+      list = NULL; // pass NULL on successive strtok_r calls
+      fprintf (stderr, "PHOTCODE_B_LIST: %s\n", PHOTCODE_B_LIST);
+      fprintf (stderr, "codename: %s\n", codename);
+      if ((photcodesGroupB[NphotcodesGroupB] = GetPhotcodebyName (codename)) == NULL) {
+        fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename);
+        exit (1);
+      }
+      NphotcodesGroupB ++;
+      CHECK_REALLOCATE (photcodesGroupB, PhotCode *, NPHOTCODES, NphotcodesGroupB, 10);
+    }
+  }
+
   initstats (STATMODE);
 
Index: trunk/Ohana/src/relastro/src/relastro.c
===================================================================
--- trunk/Ohana/src/relastro/src/relastro.c	(revision 28043)
+++ trunk/Ohana/src/relastro/src/relastro.c	(revision 28184)
@@ -23,4 +23,10 @@
   if (FIT_TARGET == TARGET_OBJECTS) {
     relastro_objects ();
+    exit (0);
+  }
+
+  /* the object analysis is a separate process iterating over catalogs */
+  if (FIT_TARGET == TARGET_HIGH_SPEED) {
+    high_speed_catalogs ();
     exit (0);
   }