# include "dvolens.h"

int update_objects_catalog (Catalog *catalog) {

  off_t i, j, Lj, Mj;

  if (VERBOSE2) fprintf (stderr, "lensobj catalog %d : "OFF_T_FMT" ave, "OFF_T_FMT" meas, "OFF_T_FMT" lensing\n", 
			 catalog->catID, catalog->Naverage, catalog->Nmeasure, catalog->Nlensing);

  // I am making this measurement per filter
  int Nsecfilt = GetPhotcodeNsecfilt ();
  
  // I think we have already allocated lensobj
  REALLOCATE (catalog->lensobj, Lensobj, catalog->Naverage * Nsecfilt);

  // Nlensobj tracks how many lensobj entries have I actually generated so far for this
  // catalog This is not necessarily 1-to-1 with secfilt entries, but any average entry
  // will either have 0 or Nsecfilt lensobj entries
  int Nlensobj = 0;

  for (i = 0; i < catalog->Naverage; i++) {
    Average *average = &catalog->average[i];
    if (average->Nlensing == 0) continue;

    // NOTE Nlensobj increments by Nsecfilt each time
    Lensobj *lensobj = &catalog->lensobj[Nlensobj];
    for (j = 0; j < Nsecfilt; j++) {
      dvo_lensobj_init (&lensobj[j], TRUE); // init accumulated values to 0
    }

    // reset the average values

    off_t Loff = average->lensingOffset;
    off_t Moff = average->measureOffset;

    // 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++) {
      
      // find the corresponding measure
      Lensing *lensing = &catalog->lensing[Loff + Lj];
      if (!isfinite(lensing->X11_sm_obj)) continue;

      Measure *measure = NULL;
      int found = FALSE;
      for (Mj = 0; !found && (Mj < average->Nmeasure); Mj++) {
	measure = &catalog->measure[Moff + Mj];
	if (measure->imageID != lensing->imageID) continue;
	found = TRUE;
      }
      if (!found) {
	fprintf (stderr, "error: cannot match measurement with lensing parameter\n");
	fprintf (stderr, "objID: %d, catID: %d, detID: %d, imageID: %d, averef: %d\n",
		 lensing->objID, lensing->catID, lensing->detID, lensing->imageID, lensing->averef); 
	abort();
      }
      
      // skip measurements that do not match the current photcode
      PhotCode *code = GetPhotcodebyCode (measure->photcode);
      myAssert (code, "missing photcode?");
      myAssert (code->equiv > -1, "photcode not equivalent to secfilt?");
      
      int Nsec = GetPhotcodeNsec (code->equiv);
      myAssert (Nsec > -1, "cannot find Nsec?");
	
      // if (measure->dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad); 
      if (measure->psfQF < 0.85) continue; 
      if (measure->psfQFperf < 0.85) continue; 
      

      lensobj[Nsec].X11_sm_obj += lensing->X11_sm_obj;
      lensobj[Nsec].X12_sm_obj += lensing->X12_sm_obj;
      lensobj[Nsec].X22_sm_obj += lensing->X22_sm_obj;
      lensobj[Nsec]. E1_sm_obj += lensing-> E1_sm_obj;
      lensobj[Nsec]. E2_sm_obj += lensing-> E2_sm_obj;

      lensobj[Nsec].X11_sh_obj += lensing->X11_sh_obj;
      lensobj[Nsec].X12_sh_obj += lensing->X12_sh_obj;
      lensobj[Nsec].X22_sh_obj += lensing->X22_sh_obj;
      lensobj[Nsec]. E1_sh_obj += lensing-> E1_sh_obj;
      lensobj[Nsec]. E2_sh_obj += lensing-> E2_sh_obj;

      lensobj[Nsec].X11_sm_psf += lensing->X11_sm_psf;
      lensobj[Nsec].X12_sm_psf += lensing->X12_sm_psf;
      lensobj[Nsec].X22_sm_psf += lensing->X22_sm_psf;
      lensobj[Nsec]. E1_sm_psf += lensing-> E1_sm_psf;
      lensobj[Nsec]. E2_sm_psf += lensing-> E2_sm_psf;

      lensobj[Nsec].X11_sh_psf += lensing->X11_sh_psf;
      lensobj[Nsec].X12_sh_psf += lensing->X12_sh_psf;
      lensobj[Nsec].X22_sh_psf += lensing->X22_sh_psf;
      lensobj[Nsec]. E1_sh_psf += lensing-> E1_sh_psf;
      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;

      lensobj[Nsec].Nmeas   ++;
      lensobj[Nsec].photcode = code->equiv;
      lensobj[Nsec].objID = lensing->objID;
      lensobj[Nsec].catID = lensing->catID;
    }
    
    for (j = 0; j < Nsecfilt; j++) {
      if (!lensobj[j].Nmeas) {
	dvo_lensobj_init (&lensobj[j], FALSE);
	continue;
      }
      float Nmeas = lensobj[j].Nmeas;
      lensobj[j].X11_sm_obj /= Nmeas;
      lensobj[j].X12_sm_obj /= Nmeas;
      lensobj[j].X22_sm_obj /= Nmeas;
      lensobj[j]. E1_sm_obj /= Nmeas;
      lensobj[j]. E2_sm_obj /= Nmeas;
      lensobj[j].X11_sh_obj /= Nmeas;
      lensobj[j].X12_sh_obj /= Nmeas;
      lensobj[j].X22_sh_obj /= Nmeas;
      lensobj[j]. E1_sh_obj /= Nmeas;
      lensobj[j]. E2_sh_obj /= Nmeas;
      lensobj[j].X11_sm_psf /= Nmeas;
      lensobj[j].X12_sm_psf /= Nmeas;
      lensobj[j].X22_sm_psf /= Nmeas;
      lensobj[j]. E1_sm_psf /= Nmeas;
      lensobj[j]. E2_sm_psf /= Nmeas;
      lensobj[j].X11_sh_psf /= Nmeas;
      lensobj[j].X12_sh_psf /= Nmeas;
      lensobj[j].X22_sh_psf /= Nmeas;
      lensobj[j]. E1_sh_psf /= Nmeas;
      lensobj[j]. E2_sh_psf /= Nmeas;

      lensobj[j]. F_ApR5 /= Nmeas;
      lensobj[j].dF_ApR5 /= Nmeas;
      lensobj[j].sF_ApR5 /= Nmeas;
      lensobj[j].fF_ApR5 /= Nmeas;
      lensobj[j]. F_ApR6 /= Nmeas;
      lensobj[j].dF_ApR6 /= Nmeas;
      lensobj[j].sF_ApR6 /= Nmeas;
      lensobj[j].fF_ApR6 /= 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);

      // somewhere in here I should calculate the values of E1, E2, and gamma
    }

    average->Nlensobj = Nsecfilt;
    average->lensobjOffset = Nlensobj;
    Nlensobj += Nsecfilt;
    // XVERB |= (catalog->average[j].objID == OBJ_ID_SRC) && (catalog->average[j].catID == CAT_ID_SRC);
    // XVERB |= (catalog->average[j].objID == OBJ_ID_DST) && (catalog->average[j].catID == CAT_ID_DST);
  }

  catalog->Nlensobj = Nlensobj;
  catalog->Nlensobj_disk = 0;
  catalog->Nlensobj_off  = 0;
  REALLOCATE (catalog->lensobj, Lensobj, Nlensobj);

  return (TRUE);
}
