Index: branches/eam_branches/ipp-20150405/Ohana/src/dvolens/Makefile =================================================================== --- branches/eam_branches/ipp-20150405/Ohana/src/dvolens/Makefile (revision 38142) +++ branches/eam_branches/ipp-20150405/Ohana/src/dvolens/Makefile (revision 38143) @@ -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: branches/eam_branches/ipp-20150405/Ohana/src/dvolens/include/dvolens.h =================================================================== --- branches/eam_branches/ipp-20150405/Ohana/src/dvolens/include/dvolens.h (revision 38142) +++ branches/eam_branches/ipp-20150405/Ohana/src/dvolens/include/dvolens.h (revision 38143) @@ -4,4 +4,16 @@ # include # include + +# 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: branches/eam_branches/ipp-20150405/Ohana/src/dvolens/src/myIndex.c =================================================================== --- branches/eam_branches/ipp-20150405/Ohana/src/dvolens/src/myIndex.c (revision 38143) +++ branches/eam_branches/ipp-20150405/Ohana/src/dvolens/src/myIndex.c (revision 38143) @@ -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: branches/eam_branches/ipp-20150405/Ohana/src/dvolens/src/update_objects_catalog.c =================================================================== --- branches/eam_branches/ipp-20150405/Ohana/src/dvolens/src/update_objects_catalog.c (revision 38142) +++ branches/eam_branches/ipp-20150405/Ohana/src/dvolens/src/update_objects_catalog.c (revision 38143) @@ -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) { @@ -26,4 +28,6 @@ // will either have 0 or Nsecfilt lensobj entries int Nlensobj = 0; + + myIndexType *measureIndex = myIndexAlloc(); for (i = 0; i < catalog->Naverage; i++) { @@ -44,4 +48,21 @@ 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,8 +73,13 @@ // 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,5 +114,5 @@ // 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,32 +144,25 @@ 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,4 +239,6 @@ free (Myy_obj); + myIndexFree (measureIndex); + catalog->Nlensobj = Nlensobj; catalog->Nlensobj_disk = 0;