Index: trunk/Ohana/src/dvolens/Makefile
===================================================================
--- trunk/Ohana/src/dvolens/Makefile	(revision 38042)
+++ trunk/Ohana/src/dvolens/Makefile	(revision 38153)
@@ -29,4 +29,5 @@
 $(SRC)/help.$(ARCH).o		 \
 $(SRC)/initialize.$(ARCH).o	 \
+$(SRC)/myIndex.$(ARCH).o	 \
 $(SRC)/update_objects.$(ARCH).o	 \
 $(SRC)/update_objects_catalog.$(ARCH).o	 \
@@ -44,4 +45,5 @@
 $(SRC)/help.$(ARCH).o		 \
 $(SRC)/initialize.$(ARCH).o	 \
+$(SRC)/myIndex.$(ARCH).o	 \
 $(SRC)/update_objects.$(ARCH).o	 \
 $(SRC)/update_objects_catalog.$(ARCH).o	 \
Index: trunk/Ohana/src/dvolens/include/dvolens.h
===================================================================
--- trunk/Ohana/src/dvolens/include/dvolens.h	(revision 38042)
+++ trunk/Ohana/src/dvolens/include/dvolens.h	(revision 38153)
@@ -4,4 +4,16 @@
 # include <signal.h>
 # include <pthread.h>
+
+# ifndef MAX_INT
+# define MAX_INT 2147483647
+# endif
+
+typedef struct {
+  int minID;
+  int maxID;
+  int *index;
+  int Nindex;
+  int NINDEX;
+} myIndexType;
 
 typedef enum {
@@ -63,2 +75,9 @@
 int  	      client_logger_message   PROTO((char *format,...));
 
+myIndexType *myIndexAlloc ();
+int myIndexFree (myIndexType *myIndex);
+void myIndexInit (myIndexType *myIndex);
+int myIndexUpdateLimits (myIndexType *myIndex, int value);
+int myIndexSetRange (myIndexType *myIndex);
+int myIndexSetEntry (myIndexType *myIndex, int value, int entry);
+int myIndexGetEntry (myIndexType *myIndex, int value);
Index: trunk/Ohana/src/dvolens/src/myIndex.c
===================================================================
--- trunk/Ohana/src/dvolens/src/myIndex.c	(revision 38153)
+++ trunk/Ohana/src/dvolens/src/myIndex.c	(revision 38153)
@@ -0,0 +1,84 @@
+# include "dvolens.h"
+
+// this is a copy of the file in dvomerge
+
+myIndexType *myIndexAlloc () {
+
+  myIndexType *myIndex;
+  ALLOCATE (myIndex, myIndexType, 1);
+  ALLOCATE (myIndex->index, int, 1);
+
+  // an uninit'ed index is invalid
+  myIndex->minID = 0;
+  myIndex->maxID = 0;
+
+  return myIndex;
+}
+
+int myIndexFree (myIndexType *myIndex) {
+
+  free (myIndex->index);
+  free (myIndex);
+  return TRUE;
+}
+
+void myIndexInit (myIndexType *myIndex) {
+
+  myIndex->minID = MAX_INT;
+  myIndex->maxID = 0;
+
+  myIndex->NINDEX = 1000;
+  myIndex->Nindex = 0;
+
+  return;
+}
+
+// set minID and maxID externally
+int myIndexUpdateLimits (myIndexType *myIndex, int value) {
+
+  myIndex->maxID = MAX(value, myIndex->maxID);
+  myIndex->minID = MIN(value, myIndex->minID);
+
+  return TRUE;
+}
+
+// once minID and maxID are set, this function defines the range and inits
+int myIndexSetRange (myIndexType *myIndex) {
+
+  myIndex->Nindex = myIndex->maxID - myIndex->minID + 1;
+  myIndex->NINDEX = myIndex->Nindex;
+  
+  REALLOCATE (myIndex->index, int, myIndex->NINDEX);
+
+  int i;
+  for (i = 0; i < myIndex->Nindex; i++) {
+    myIndex->index[i] = -1;
+  }
+
+  return TRUE;
+}
+
+// once minID and maxID are set, this function defines the range and inits
+int myIndexSetEntry (myIndexType *myIndex, int value, int entry) {
+
+  int n = value - myIndex->minID;
+
+  myAssert (n >= 0, "impossible!");
+
+  if (myIndex->index[n] > -1) return -1;
+  // myAssert(myIndex->index[n] == -1, "stomping on an earlier index");
+
+  myIndex->index[n] = entry;
+  return n;
+}
+
+// once minID and maxID are set, this function defines the range and inits
+int myIndexGetEntry (myIndexType *myIndex, int value) {
+
+  if (value < myIndex->minID) return -1; // not in this index
+  if (value > myIndex->maxID) return -1; // not in this index
+
+  int n = value - myIndex->minID;
+
+  return myIndex->index[n];
+}
Index: trunk/Ohana/src/dvolens/src/update_objects_catalog.c
===================================================================
--- trunk/Ohana/src/dvolens/src/update_objects_catalog.c	(revision 38042)
+++ trunk/Ohana/src/dvolens/src/update_objects_catalog.c	(revision 38153)
@@ -1,3 +1,5 @@
 # include "dvolens.h"
+# define SCALE 0.001
+float MagToFlux (float Mag); // in libdvo, but not exposed?
 
 int update_objects_catalog (Catalog *catalog) {
@@ -11,4 +13,12 @@
   int Nsecfilt = GetPhotcodeNsecfilt ();
   
+  float *Mxx_obj = NULL;
+  float *Mxy_obj = NULL;
+  float *Myy_obj = NULL;
+
+  ALLOCATE (Mxx_obj, float, Nsecfilt);
+  ALLOCATE (Mxy_obj, float, Nsecfilt);
+  ALLOCATE (Myy_obj, float, Nsecfilt);
+
   // I think we have already allocated lensobj
   REALLOCATE (catalog->lensobj, Lensobj, catalog->Naverage * Nsecfilt);
@@ -19,4 +29,6 @@
   int Nlensobj = 0;
 
+  myIndexType *measureIndex = myIndexAlloc();
+
   for (i = 0; i < catalog->Naverage; i++) {
     Average *average = &catalog->average[i];
@@ -27,10 +39,29 @@
     for (j = 0; j < Nsecfilt; j++) {
       dvo_lensobj_init (&lensobj[j], TRUE); // init accumulated values to 0
+      Mxx_obj[j] = 0.0;
+      Mxy_obj[j] = 0.0;
+      Myy_obj[j] = 0.0;
     }
 
     // reset the average values
-
     off_t Loff = average->lensingOffset;
     off_t Moff = average->measureOffset;
+
+    // I have Nmeasure entries for this object.  I need to match measure to lensing (by imageID)
+    // I could generated a sorted list (imageID, measureSeq) and use bisection to find the desired imageID
+    // I could make an index measureSeq[imageID-imageIDmin]
+    
+    myIndexInit (measureIndex);
+
+    // generate an index for these measure entries (based on Mj and imageID)
+    for (Mj = 0; Mj < average->Nmeasure; Mj++) {
+      myIndexUpdateLimits (measureIndex, catalog->measure[Moff + Mj].imageID);
+    }
+    myIndexSetRange (measureIndex);
+    for (Mj = 0; Mj < average->Nmeasure; Mj++) {
+      // if there are duplicates (multiple measures from the same image), 
+      // myIndexSetEntry will only select one
+      myIndexSetEntry (measureIndex, catalog->measure[Moff + Mj].imageID, Mj);
+    }
 
     // loop over the lensing measurements.  for each one, I need to find the corresponding measurement (make an index in lensing?)
@@ -41,4 +72,17 @@
       if (!isfinite(lensing->X11_sm_obj)) continue;
 
+      // pointer from lensing entry to corresponding measure entry (keep updated?)
+      // XX ALT int Mseq = lensing->measureSeq;
+      // XX ALT myAssert (Mseq < average->Nmeasure, "oops");
+      // XX ALT Measure *measure = &catalog->measure[Moff + Mseq];
+      // XX ALT myAssert (measure->imageID == lensing->imageID, "oops, deux");
+
+      Mj = myIndexGetEntry (measureIndex, lensing->imageID);
+      myAssert (Mj > -1, "missing index");
+
+      Measure *measure = &catalog->measure[Moff + Mj];
+      myAssert (measure->imageID == lensing->imageID, "oops, deux");
+
+# if (0)      
       Measure *measure = NULL;
       int found = FALSE;
@@ -54,4 +98,5 @@
 	abort();
       }
+# endif
       
       // skip measurements that do not match the current photcode
@@ -67,5 +112,12 @@
       if (measure->psfQFperf < 0.85) continue; 
       
-
+      // I should probably only use lensing entries which correspond to warps 
+      // included in the mean warp flux (setMrelCatalog.c)
+      if ((measure->dbFlags & ID_MEAS_WARP_USED) == 0) continue;
+
+      Mxx_obj[Nsec] += measure->Mxx;
+      Mxy_obj[Nsec] += measure->Mxy;
+      Myy_obj[Nsec] += measure->Myy;
+      
       lensobj[Nsec].X11_sm_obj += lensing->X11_sm_obj;
       lensobj[Nsec].X12_sm_obj += lensing->X12_sm_obj;
@@ -92,13 +144,25 @@
       lensobj[Nsec]. E2_sh_psf += lensing-> E2_sh_psf;
 
-      lensobj[Nsec]. F_ApR5 += lensing-> F_ApR5;
-      lensobj[Nsec].dF_ApR5 += SQ(lensing->dF_ApR5);
-      lensobj[Nsec].sF_ApR5 += SQ(lensing->sF_ApR5);
-      lensobj[Nsec].fF_ApR5 += lensing->fF_ApR5;
-
-      lensobj[Nsec]. F_ApR6 += lensing-> F_ApR6;
-      lensobj[Nsec].dF_ApR6 += SQ(lensing->dF_ApR6);
-      lensobj[Nsec].sF_ApR6 += SQ(lensing->sF_ApR6);
-      lensobj[Nsec].fF_ApR6 += lensing->fF_ApR6;
+      // relphot sets measure->Mcal (setMcalOutput.c, called by setMrelFinal.c)
+      float Mcal = code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C - measure->Mcal;
+      float Fcal = 3630.8 * MagToFlux(Mcal);
+
+      // lensing->F_ApR5, etc need to be in units of DN/sec
+      // Fcal * lensing->F_ApR5 is in Jy
+
+      lensobj[Nsec]. F_ApR5 +=    Fcal * lensing-> F_ApR5;
+      lensobj[Nsec].dF_ApR5 += SQ(Fcal * lensing->dF_ApR5);
+      lensobj[Nsec].sF_ApR5 += SQ(Fcal * lensing->sF_ApR5);
+      lensobj[Nsec].fF_ApR5 +=           lensing->fF_ApR5;
+
+      lensobj[Nsec]. F_ApR6 +=    Fcal * lensing-> F_ApR6;
+      lensobj[Nsec].dF_ApR6 += SQ(Fcal * lensing->dF_ApR6);
+      lensobj[Nsec].sF_ApR6 += SQ(Fcal * lensing->sF_ApR6);
+      lensobj[Nsec].fF_ApR6 +=           lensing->fF_ApR6;
+
+      lensobj[Nsec]. F_ApR7 +=    Fcal * lensing-> F_ApR7;
+      lensobj[Nsec].dF_ApR7 += SQ(Fcal * lensing->dF_ApR7);
+      lensobj[Nsec].sF_ApR7 += SQ(Fcal * lensing->sF_ApR7);
+      lensobj[Nsec].fF_ApR7 +=           lensing->fF_ApR7;
 
       lensobj[Nsec].Nmeas   ++;
@@ -139,4 +203,5 @@
       lensobj[j].sF_ApR5 /= Nmeas;
       lensobj[j].fF_ApR5 /= Nmeas;
+
       lensobj[j]. F_ApR6 /= Nmeas;
       lensobj[j].dF_ApR6 /= Nmeas;
@@ -144,10 +209,21 @@
       lensobj[j].fF_ApR6 /= Nmeas;
 
+      lensobj[j]. F_ApR7 /= Nmeas;
+      lensobj[j].dF_ApR7 /= Nmeas;
+      lensobj[j].sF_ApR7 /= Nmeas;
+      lensobj[j].fF_ApR7 /= Nmeas;
+
       lensobj[j].dF_ApR5 = sqrt(lensobj[j].dF_ApR5);
       lensobj[j].sF_ApR5 = sqrt(lensobj[j].sF_ApR5);
       lensobj[j].dF_ApR6 = sqrt(lensobj[j].dF_ApR6);
       lensobj[j].sF_ApR6 = sqrt(lensobj[j].sF_ApR6);
+      lensobj[j].dF_ApR7 = sqrt(lensobj[j].dF_ApR7);
+      lensobj[j].sF_ApR7 = sqrt(lensobj[j].sF_ApR7);
 
       // somewhere in here I should calculate the values of E1, E2, and gamma
+
+      float e0 = Mxx_obj[j] + Myy_obj[j];
+      lensobj[j].E1 = (Mxx_obj[j] - Myy_obj[j]) / e0;
+      lensobj[j].E2 = 2.0 * Mxy_obj[j] / e0;
     }
 
@@ -159,4 +235,10 @@
   }
 
+  free (Mxx_obj);
+  free (Mxy_obj);
+  free (Myy_obj);
+
+  myIndexFree (measureIndex);
+
   catalog->Nlensobj = Nlensobj;
   catalog->Nlensobj_disk = 0;