Changeset 38143

Timestamp:

Apr 17, 2015, 12:37:52 PM (11 years ago)

Author:

eugene

Message:

apply zero point calibration to R5,R6,R7 data and set mean in Jy

Location:

branches/eam_branches/ipp-20150405/Ohana/src/dvolens

Files:

• Makefile (modified) (2 diffs)
• include/dvolens.h (modified) (2 diffs)
• src/myIndex.c (added)
• src/update_objects_catalog.c (modified) (7 diffs)

branches/eam_branches/ipp-20150405/Ohana/src/dvolens/Makefile

@@ -29 +29 @@
 $(SRC)/help.$(ARCH).o            \
 $(SRC)/initialize.$(ARCH).o      \
+$(SRC)/myIndex.$(ARCH).o         \
 $(SRC)/update_objects.$(ARCH).o  \
 $(SRC)/update_objects_catalog.$(ARCH).o  \
@@ -44 +45 @@
 $(SRC)/help.$(ARCH).o            \
 $(SRC)/initialize.$(ARCH).o      \
+$(SRC)/myIndex.$(ARCH).o         \
 $(SRC)/update_objects.$(ARCH).o  \
 $(SRC)/update_objects_catalog.$(ARCH).o  \

branches/eam_branches/ipp-20150405/Ohana/src/dvolens/include/dvolens.h

@@ -4 +4 @@
 # 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 +75 @@
 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);

branches/eam_branches/ipp-20150405/Ohana/src/dvolens/src/update_objects_catalog.c

@@ -1 +1 @@
 # include "dvolens.h"
+# define SCALE 0.001
+float MagToFlux (float Mag); // in libdvo, but not exposed?
 
 int update_objects_catalog (Catalog *catalog) {
@@ -26 +28 @@
   // will either have 0 or Nsecfilt lensobj entries
   int Nlensobj = 0;
+
+  myIndexType *measureIndex = myIndexAlloc();
 
   for (i = 0; i < catalog->Naverage; i++) {
@@ -44 +48 @@
     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?)
     for (Lj = 0; Lj < average->Nlensing; Lj++) {
@@ -52 +73 @@
 
       // pointer from lensing entry to corresponding measure entry (keep updated?)
-      int Mseq = lensing->measureSeq;
-      myAssert (Mseq < average->Nmeasure, "oops");
-
-      Measure *measure = &catalog->measure[Moff + Mseq];
+      // 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");
 
@@ -88 +114 @@
       // 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;
+      if ((measure->dbFlags & ID_MEAS_WARP_USED) == 0) continue;
 
       Mxx_obj[Nsec] += measure->Mxx;
@@ -118 +144 @@
       lensobj[Nsec]. E2_sh_psf += lensing-> E2_sh_psf;
 
-      myAbort ("NOTE: relphot needs to set measure->Mcal -- is this happening or not??");
-
-      int Np = photcodes[0].hashcode[measure[0].photcode];
-      if (Np == -1) return (NAN);
-      PhotCode *code = &photcodes[0].code[Np];
-      
-      float Mdef = code[0].K*(measure[0].airmass - 1.000) + SCALE*code[0].C;
-      float Moff = Mdef + 8.9 - measure->Mcal;
-      float Foff = 3630.8 * MagToFlux(Moff);
+      // 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
-      // Foff * lensing->F_ApR5 is in Jy
-      // XXX this is unweighted -- apply a weighting?
-
-      lensobj[Nsec]. F_ApR5 += Foff * lensing-> F_ApR5;
-      lensobj[Nsec].dF_ApR5 += SQ(Foff * lensing->dF_ApR5);
-      lensobj[Nsec].sF_ApR5 += SQ(Foff * lensing->sF_ApR5);
-      lensobj[Nsec].fF_ApR5 += Foff * lensing->fF_ApR5;
-
-      lensobj[Nsec]. F_ApR6 += Foff * lensing-> F_ApR6;
-      lensobj[Nsec].dF_ApR6 += SQ(Foff * lensing->dF_ApR6);
-      lensobj[Nsec].sF_ApR6 += SQ(Foff * lensing->sF_ApR6);
-      lensobj[Nsec].fF_ApR6 += Foff * lensing->fF_ApR6;
-
-      lensobj[Nsec]. F_ApR7 += Foff * lensing-> F_ApR7;
-      lensobj[Nsec].dF_ApR7 += SQ(Foff * lensing->dF_ApR7);
-      lensobj[Nsec].sF_ApR7 += SQ(Foff * lensing->sF_ApR7);
-      lensobj[Nsec].fF_ApR7 += Foff * lensing->fF_ApR7;
+      // 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   ++;
@@ -220 +239 @@
   free (Myy_obj);
 
+  myIndexFree (measureIndex);
+
   catalog->Nlensobj = Nlensobj;
   catalog->Nlensobj_disk = 0;