Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/include/relastro.h
===================================================================
--- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/include/relastro.h	(revision 33569)
+++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/include/relastro.h	(revision 33570)
@@ -4,4 +4,5 @@
 # include <signal.h>
 # include <assert.h>
+# include <pthread.h>
 
 // choose off_t or int depending on full-scale relphot analysis resources
@@ -159,4 +160,5 @@
 int    RESET;
 int    NLOOP;
+int    NTHREADS;
 int    UPDATE;
 int    PLOTSTUFF;
Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateChips.c
===================================================================
--- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateChips.c	(revision 33569)
+++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/UpdateChips.c	(revision 33570)
@@ -1,13 +1,38 @@
 # include "relastro.h"
+
 int plotChipFits (double *Ro, double *Do, char *mode, int Nimage);
 int saveCenter (Image *image, double *Ro, double *Do, int im);
-
+off_t getNextImageForThread ();
+int UpdateChips_threaded (Catalog *catalog, int Ncatalog);
+void *UpdateChips_worker (void *data);
+
+enum {THREAD_RUN, THREAD_DONE};
+
+typedef struct {
+  int entry;
+  int state;
+  double *Ro;
+  double *Do;
+  char *mode;
+  off_t Nskip;
+  off_t Nmosaic;
+  off_t NnewFit;
+  off_t NoldFit;
+  Catalog *catalog;
+  int Ncatalog;
+} ThreadInfo;
+
+static Image *image = NULL;
+static off_t Nimage = 0;
+static off_t nextImage = 0;
+
+// update astrometry of all chips relative to the average positions
+// if NTHREADS is non-zero, call the threaded version of this function
 int UpdateChips (Catalog *catalog, int Ncatalog) {
 
-  int Nskip, Nmosaic, NnewFit, NoldFit;
+  off_t Nskip, Nmosaic, NnewFit, NoldFit;
 
   /* we can measure new image parameters for each non-mosaic chip independently */
-  off_t i, Nimage, Nraw, Nref, nFitAstr;
-  Image *image;
+  off_t i, Nraw, Nref, nFitAstr;
   StarData *raw, *ref;
   Coords *oldCoords;
@@ -15,4 +40,9 @@
   double *Ro, *Do;
   char *mode;
+
+  if (NTHREADS) {
+    UpdateChips_threaded (catalog, Ncatalog);
+    return TRUE;
+  }
 
   Nskip = Nmosaic = NnewFit = NoldFit = 0;
@@ -108,9 +138,215 @@
 
   plotChipFits (Ro, Do, mode, Nimage);
-
-  fprintf (stderr, "UpdateChips: %d fitted, %d keep old, %d skipped, %d mosaic (skipped)\n", NnewFit, NoldFit, Nskip, Nmosaic);
+  free (Ro);
+  free (Do);
+  free (mode);
+
+  fprintf (stderr, "UpdateChips: %d fitted, %d keep old, %d skipped, %d mosaic (skipped)\n", 
+	   (int) NnewFit, (int) NoldFit, (int) Nskip, (int) Nmosaic);
   return (TRUE);
 }
 
+int UpdateChips_threaded (Catalog *catalog, int Ncatalog) {
+
+  int i;
+  off_t Nskip, Nmosaic, NnewFit, NoldFit;
+  double *Ro, *Do;
+  char *mode;
+
+  image = getimages (&Nimage, NULL);
+
+  // save fit results for summary plot
+  ALLOCATE (Ro, double, Nimage);
+  ALLOCATE (Do, double, Nimage);
+  ALLOCATE (mode, char, Nimage);
+
+  pthread_attr_t attr;
+  pthread_attr_init (&attr);
+  pthread_attr_setdetachstate (&attr, PTHREAD_CREATE_DETACHED);
+  
+  pthread_t *threads;
+  ALLOCATE (threads, pthread_t, NTHREADS);
+
+  ThreadInfo *threadinfo;
+  ALLOCATE (threadinfo, ThreadInfo, NTHREADS);
+
+  // launch N worker threads
+  for (i = 0; i < NTHREADS; i++) {
+    threadinfo[i].entry = i;
+    threadinfo[i].state = THREAD_RUN;
+    threadinfo[i].Ro = Ro;
+    threadinfo[i].Do = Do;
+    threadinfo[i].mode = mode;
+    threadinfo[i].Nskip = 0;
+    threadinfo[i].Nmosaic = 0;
+    threadinfo[i].NnewFit = 0;
+    threadinfo[i].NoldFit = 0;
+    threadinfo[i].catalog  =  catalog;
+    threadinfo[i].Ncatalog = Ncatalog;
+    pthread_create (&threads[i], NULL, UpdateChips_worker, &threadinfo[i]);
+  }
+  pthread_attr_destroy (&attr);
+
+  // wait until all threads have finished
+  while (1) {
+    int allDone = TRUE;
+    for (i = 0; i < NTHREADS; i++) {
+      if (threadinfo[i].state == THREAD_RUN) allDone = FALSE;
+    }
+    if (allDone) {
+      break;
+    }
+    usleep (500000);
+  }
+
+  // all threads are done, free the threads array and grab the info
+  free (threads);
+  
+  plotChipFits (Ro, Do, mode, Nimage);
+  free (Ro);
+  free (Do);
+  free (mode);
+
+  // report stats & summary from the threads
+  Nskip = Nmosaic = NnewFit = NoldFit = 0;
+  for (i = 0; i < NTHREADS; i++) {
+    fprintf (stderr, "UpdateChips thread %d : %d fitted, %d keep old, %d skipped, %d mosaic (skipped)\n", 
+	     i, (int) threadinfo[i].NnewFit, (int) threadinfo[i].NoldFit, (int) threadinfo[i].Nskip, (int) threadinfo[i].Nmosaic);
+
+    NnewFit += threadinfo[i].NnewFit;
+    NoldFit += threadinfo[i].NoldFit;
+    Nskip   += threadinfo[i].Nskip;
+    Nmosaic += threadinfo[i].Nmosaic;
+  }
+  free (threadinfo);
+
+  fprintf (stderr, "UpdateChips: %d fitted, %d keep old, %d skipped, %d mosaic (skipped)\n", 
+	   (int) NnewFit, (int) NoldFit, (int) Nskip, (int) Nmosaic);
+  return (TRUE);
+}
+
+void *UpdateChips_worker (void *data) {
+
+  /* we can measure new image parameters for each non-mosaic chip independently */
+  off_t Nraw, Nref, nFitAstr;
+  StarData *raw, *ref;
+  Coords *oldCoords;
+  float dXpixSys, dYpixSys;
+
+  ThreadInfo *threadinfo = data;
+
+  while (1) {
+
+    off_t i = getNextImageForThread();
+    if (i == -1) {
+      threadinfo->state = THREAD_DONE;
+      return NULL;
+    }
+
+    /* skip all except WRP images */
+    if (strcmp(&image[i].coords.ctype[4], "-WRP")) {
+      threadinfo->Nmosaic ++;
+      threadinfo->mode[i] = 0;
+      continue;
+    }
+
+    /* convert measure coordinates to raw entries */
+    raw = getImageRaw (threadinfo->catalog, threadinfo->Ncatalog, i, &Nraw, MODE_MOSAIC);
+    if (!raw) {
+      threadinfo->Nskip ++;
+      threadinfo->mode[i] = 0;
+      continue;
+    }
+
+    /* convert average coordinates to ref entries */
+    ref = getImageRef (threadinfo->catalog, threadinfo->Ncatalog, i, &Nref, MODE_MOSAIC);
+    if (!ref) {
+      threadinfo->Nskip ++;
+      threadinfo->mode[i] = 0;
+      free (raw);
+      continue;
+    }
+
+    // note that Nraw & Nref must be equal: if not, we made a programming error in one of these two functions.
+    assert (Nraw == Nref);
+
+    // save these in case of failure
+    saveCoords (&image[i].coords, i);
+    dXpixSys = image[i].dXpixSys;
+    dYpixSys = image[i].dYpixSys;
+    nFitAstr = image[i].nFitAstrom;
+
+    // FitChip does iterative, clipped fitting
+    // fprintf (stderr, "image "OFF_T_FMT" : Nstars: "OFF_T_FMT"\n",  i,  Nraw);
+    if (!FitChip (raw, ref, Nraw, &image[i])) {
+      if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT" of %d\n", image[i].name,  i,  Nraw, image[i].nstar);
+
+      // restore status quo ante
+      oldCoords = getCoords (i);
+      memcpy (&image[i].coords, oldCoords, sizeof(Coords));
+      image[i].dXpixSys = dXpixSys;
+      image[i].dYpixSys = dYpixSys; 
+      image[i].nFitAstrom = nFitAstr;
+
+      saveCenter (image, &threadinfo->Ro[i], &threadinfo->Do[i], i);
+      threadinfo->mode[i] = 1;
+      threadinfo->NoldFit ++;
+      free (raw);
+      free (ref);
+      continue;
+    }
+
+    if (!checkStarMap (i)) {
+      if (VERBOSE) fprintf (stderr, "fit diverges too much for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name,  i,  Nraw);
+      // restore status quo ante
+      oldCoords = getCoords (i);
+      memcpy (&image[i].coords, oldCoords, sizeof(Coords));
+      image[i].dXpixSys = dXpixSys;
+      image[i].dYpixSys = dYpixSys; 
+      image[i].nFitAstrom = nFitAstr;
+
+      saveCenter (image, &threadinfo->Ro[i], &threadinfo->Do[i], i);
+      threadinfo->mode[i] = 2;
+      image[i].flags |= ID_IMAGE_ASTROM_POOR;
+      threadinfo->NoldFit ++;
+      free (raw);
+      free (ref);
+      continue;
+    } 
+
+    saveCenter (image, &threadinfo->Ro[i], &threadinfo->Do[i], i);
+    threadinfo->mode[i] = 3;
+    threadinfo->NnewFit ++;
+    free (raw);
+    free (ref);
+  }
+  
+  // we should never reach here...
+  return NULL;
+}
+
+// mutex to lock UpdateChips_worker operations 
+static pthread_mutex_t UpdateChips_mutex = PTHREAD_MUTEX_INITIALIZER;
+
+// we have an array of chips (image, Nimage).  we need to hand out images one at a time to
+// the worker threads as they need
+off_t getNextImageForThread () {
+
+  pthread_mutex_lock (&UpdateChips_mutex);
+  if (nextImage >= Nimage) {
+    pthread_mutex_unlock (&UpdateChips_mutex);
+    return (-1);
+  }
+  off_t thisImage = nextImage;
+  nextImage ++;
+
+  pthread_mutex_unlock (&UpdateChips_mutex);
+  return (thisImage);
+}
+
+// This function uses the mosaic array in MosaicOps.c to associate PHU + CHIP this is
+// thread safe : the lookup (getMosaicForImage) is selecting an element of a previously
+// allocated and assigned static array.  This function does NOT use the thread-unsafe (and
+// somewhat slower) functions in libdvo/src/mosaic_astrom.c to associate CHIP and PHU
 int saveCenter (Image *image, double *Ro, double *Do, int im) {
 
@@ -139,9 +375,9 @@
     // note that for a Simple image, L,M = P,Q
     XY_to_LM (&L, &M, X, Y, imcoords);
-    LM_to_RD (&R, &D, L, M, imcoords);
+    LM_to_RD (&R, &D, L, M, imcoords); // because of the block above, ctype is not -WRP 
   } else {
     XY_to_LM (&L, &M, X, Y, imcoords);
     XY_to_LM (&P, &Q, L, M, moscoords);
-    LM_to_RD (&R, &D, P, Q, moscoords);
+    LM_to_RD (&R, &D, P, Q, moscoords); // moscoords.ctype is -DIS
   }
 
Index: /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/args.c
===================================================================
--- /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/args.c	(revision 33569)
+++ /branches/eam_branches/ipp-20111122/Ohana/src/relastro/src/args.c	(revision 33570)
@@ -360,4 +360,11 @@
     remove_argument (N, &argc, argv);
     NLOOP = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  NTHREADS = 0;
+  if ((N = get_argument (argc, argv, "-nthreads"))) {
+    remove_argument (N, &argc, argv);
+    NTHREADS = atof (argv[N]);
     remove_argument (N, &argc, argv);
   }