# include "relphot.h"
# define UBERCAL_WEIGHT 100.0

# define SKIP_THIS_MEAS(REASON) {				\
    measureT[k].dbFlags |= ID_MEAS_SKIP_PHOTOM;	\
    if (measure) {					\
      measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM;		\
    }								\
    results->REASON ++;						\
    continue; }

# define SKIP_THIS_MEAS_STACK(REASON) {			\
    measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM;		\
    continue; }

static float MagToFlux (float Mag) {
  float Flux = pow(10.0, -0.4*(Mag));
  return (Flux);
}

int print_measure_set_alt (Average *average, SecFilt *secfilt, Measure *measure) {

  int Ns;
  off_t k;

  int Nsecfilt = GetPhotcodeNsecfilt ();

  for (k = 0; k < average[0].Nmeasure; k++) {
    fprintf (stderr, "meas: %08x\n", measure[k].dbFlags);
  }

  for (Ns = 0; Ns < Nsecfilt; Ns++) {
    fprintf (stderr, "secf: %08x\n", secfilt[Ns].flags);
  }
  return (TRUE);
}

int setMrelAverageForcedWarp (off_t measureOffset, int cat, int pass, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t *found);

int setMrel_catalog_alt (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {

  off_t j;

  liststats_setmode (&results->psfstats, STATMODE);
  liststats_setmode (&results->apstats, STATMODE);
  liststats_setmode (&results->kronstats, STATMODE);

  SetMrelInfoInit (results, FALSE); // do not allocate list,dlist,wlist arrays

  int isSetMrelFinal = (pass >= 0);

  char *primaryCell = NULL;
  if (isSetMrelFinal) {
    ALLOCATE (primaryCell, char, DVO_MAX_PATH);
  }

  for (j = 0; j < catalog[Nc].Naverage; j++) {
    // we are guaranteed to have averageT, but not average
    off_t m = catalog[Nc].averageT[j].measureOffset;
    Average *average = catalog[Nc].average ? &catalog[Nc].average[j] : NULL;
    Measure *measure = catalog[Nc].measure ? &catalog[Nc].measure[m] : NULL;
    setMrelAverageExposure (m, Nc, pass, flatcorr, results, average, &catalog[Nc].averageT[j], &catalog[Nc].secfilt[j*Nsecfilt], measure, &catalog[Nc].measureT[m], &catalog[Nc].found_t[Nsecfilt*j]);

    // only apply Stack operation on setMrelFinal in first pass 
    if (isSetMrelFinal && (pass == 0)) {
      setMrelAverageStack (m, Nc, flatcorr, results, &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], &catalog[Nc].measure[m], primaryCell);
      setGlobalObjStats (&catalog[Nc].average[j], &catalog[Nc].measure[m]);
    }

    // only measure force-warp mean values if issetMrelFinal (make it optional?)
    if (isSetMrelFinal) {
      setMrelAverageForcedWarp (m, Nc, pass, results, &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], &catalog[Nc].measure[m], &catalog[Nc].foundWarp_t[Nsecfilt*j]);
      setGlobalObjStats (&catalog[Nc].average[j], &catalog[Nc].measure[m]);
    }
  }
  if (primaryCell) free (primaryCell);
  return (TRUE);
}

// set mean of chip measurements (selected by photcode range for now):
int setMrelAverageExposure (off_t measureOffset, int cat, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, AverageTiny *averageT, SecFilt *secfilt, Measure *measure, MeasureTiny *measureT, off_t *found) {

  off_t k;
  float Msys = 0, Mcal= 0, Mmos = 0, Mgrid = 0;

  // we are measuring means for 3 types of mags: psf, ap, kron.  I am using the psf mag
  // error for the ap mags, but krons have their own errors.  it is an open question if I
  // should be doing weighted or unweighted fits (this is a user option)
  double *Mpsflist  = results->Mpsflist;
  double *dpsflist  = results->dpsflist;
  double *wpsflist  = results->wpsflist;

  double *Maplist   = results->Maplist;
  double *daplist   = results->daplist;
  double *waplist   = results->waplist;

  double *Mkronlist = results->Mkronlist;
  double *dkronlist = results->dkronlist;
  double *wkronlist = results->wkronlist;

  StatType *psfstats  = &results->psfstats;
  StatType *apstats   = &results->apstats;
  StatType *kronstats = &results->kronstats;

  int isSetMrelFinal = (pass >= 0);

  // option for a test print
  if (FALSE && (average[0].objID == 0x7146) && (average[0].catID == 0x49d8)) {
    fprintf (stderr, "test obj\n");
    print_measure_set_alt (average, secfilt, measure);
  }

  int GoodPS1 = FALSE;
  int Good2MASS = FALSE;
  int Galaxy2MASS = FALSE;

  int NextPS1 = 0;
  int NpsfPS1 = 0;

  int Ns;
  for (Ns = 0; Ns < Nphotcodes; Ns++) {

    int thisCode = photcodes[Ns][0].code;
    int Nsec = GetPhotcodeNsec(thisCode);

    /* calculate the average mag in this SEC photcode for a single star */

    /* star/photcodes already calibrated */
    if ( isSetMrelFinal && found[Nsec]) continue;  
      
    // skip bad stars
    if (!isSetMrelFinal && (secfilt[Nsec].flags & STAR_BAD)) continue;

    int NexpPS1 = 0;
    int Ncode = 0;
    int Next = 0;
    int haveSynth = FALSE;
    int havePS1   = FALSE;

    int forceSynth = FALSE;
    int forceSynthEntry = -1;

    int haveUbercal = FALSE;

    int minUbercalDist = 1000;
    
    off_t meas = measureOffset;

    int Nap = 0;
    int Npsf = 0;
    int Nkron = 0;
    for (k = 0; k < averageT[0].Nmeasure; k++, meas++) {

      // skip measurements that do not match the current photcode
      PhotCode *code = GetPhotcodebyCode (measureT[k].photcode);
      if (!code) continue;
      if (code->equiv != thisCode) { continue; }
      Ncode ++;

      // are we a PS1 exposure photcode?
      if (isGPC1chip(measureT[k].photcode)) NexpPS1 ++;

      // SKIP gpc1 stack data
      if (isGPC1stack(measureT[k].photcode)) continue;

      // SKIP gpc1 forced-warp data
      if (isGPC1warp(measureT[k].photcode)) continue;

      if (measureT[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad); 

      if (getImageEntry (meas, cat) < 0) {
	// measurements without an image are either external reference photometry or
	// data for which the associated image has not been loaded (probably because of
	// overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
	Mmos = Mgrid = 0;
	Mcal = measureT[k].Mcal; // check that this is zero for loaded REF value
      } else {
	Mcal  = getMcal_alt (meas, cat, flatcorr, measureT[k].Xccd, measureT[k].Yccd);
	if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
	Mmos  = getMmos  (meas, cat);
	if (isnan(Mmos))  SKIP_THIS_MEAS(Nmos);
	Mgrid = getMgrid (meas, cat);
	if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid);
      }

      // skip some absurd values NAN, < 0.0, > 30.0
      Msys = PhotSysTiny (&measureT[k], &averageT[0], &secfilt[0], MAG_CLASS_PSF);
      if (isnan(Msys)) SKIP_THIS_MEAS(Nsys);
      if (Msys <  0.0) SKIP_THIS_MEAS(Nsys);
      if (Msys > 30.0) SKIP_THIS_MEAS(Nsys);

      int myUbercalDist = getUbercalDist(meas, cat);
      minUbercalDist = MIN(minUbercalDist, myUbercalDist);

      int isUbercal = (measureT[k].dbFlags & ID_MEAS_PHOTOM_UBERCAL);
      float dMpsf = MAX (hypot(measureT[k].dM, code->photomErrSys), MIN_ERROR);

      if (isUbercal) haveUbercal = TRUE; // haveUbercal is set per secfilt, isUbercal is per measure

      if (isSetMrelFinal) {
	// special options for PS1 data
	if ((measure[k].photcode >= 10000) && (measure[k].photcode <= 10500)) {
	  // count the extended detections
	  if (!isnan(measure[k].Map)) {
	    float dMagAp = measure[k].M - measure[k].Map;
	    float SigmaAp = hypot(0.1, 2.5*measure[k].dM);
	    // XXX this is still quite ad hoc, but at least it:
	    // (a) converges to 0.1 mag offset at the bright end
	    // (b) converges to 0.5 mag offset at the faint end (dM = 0.2)
	    if (dMagAp > SigmaAp) {
	      Next ++;
	      NextPS1 ++;
	    } else {
	      NpsfPS1 ++;
	    }
	  }
	  havePS1 = TRUE;
	}

	// count extended detections for 2MASS (XXX NOTE hardwired photcodes 2011, 2012, 2013)
	if ((measure[k].photcode >= 2011) && (measure[k].photcode <= 2013)) {
	  if (measure[k].photFlags & 0x00c00000) {
	    Next ++;
	    Galaxy2MASS = TRUE;
	  }
	  if (pass == 0) {
	    if (measure[k].photFlags & 0x00000007) {
	      Good2MASS = TRUE;
	    } else {
	      // detections without one of these bits should only be used in PASS_1
	      SKIP_THIS_MEAS(Nbad);
	    }
	  }
	}

	// Blindly accepth the SYNTH mags if we are above saturation, otherwise, 
	// ignore SYNTH photcodes until PASS == 4 (where we also accept saturated stars)
	if ((measure[k].photcode >= 3001) && (measure[k].photcode <= 3005)) {
	  // something of a hack: force object to use synth values if synth mags >>
	  // saturation (3pi instrumental mags < -15)
	  float MaxMagForceSynth = NAN;
	  switch (measure[k].photcode) {
	    case 3001:
	      MaxMagForceSynth = 13.64;
	      break;
	    case 3002:
	      MaxMagForceSynth = 13.76;
	      break;
	    case 3003:
	      MaxMagForceSynth = 13.74;
	      break;
	    case 3004:
	      MaxMagForceSynth = 12.94;
	      break;
	    case 3005:
	      MaxMagForceSynth = 12.01;
	      break;
	  }
	  if (measureT[k].M < MaxMagForceSynth) {
	    forceSynth = TRUE;
	    forceSynthEntry = Npsf;
	  } else {
	    if (pass < 4) {
	      SKIP_THIS_MEAS(Nbad);
	    }
	    haveSynth = TRUE;
	  }
	}

	// Map (Maplist) and Mkron (Mkronlist,dkronlist) are used to calculate mean mags per filter
	float Map = PhotCat (&measure[k], MAG_CLASS_APER);
	if (!isnan(Map)) {
	    Maplist[Nap] = Map - Mcal - Mmos - Mgrid;
	    daplist[Nap] = dMpsf;// XXX check on this...
	    waplist[Nap] = isUbercal ? UBERCAL_WEIGHT : 1.0;
	    Nap ++;
	}

	float Mkron = PhotCat (&measure[k], MAG_CLASS_KRON);
	if (!isnan(Mkron)) {
	  Mkronlist[Nkron] = Mkron - Mcal - Mmos - Mgrid;
	  dkronlist[Nkron] = measure[k].dMkron;
	  wkronlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0;
	  Nkron ++;
	}
      } // if (isSetMrelFinal)

      // dlist gives the error per measurement, wlist gives the weight
      // we can modify the error and weight in a few ways:
      // 1) MIN_ERROR guarantees a floor
      // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode
      // 3) UBERCAL measurements can have their weight increased by a big factor to help tie down the averages
      // 4) some reference photcode of some kind can be specified as fixed and have a high weight
      Mpsflist[Npsf] = Msys - Mcal - Mmos - Mgrid;
      dpsflist[Npsf] = dMpsf;
      wpsflist[Npsf] = isUbercal ? UBERCAL_WEIGHT : 1.0;

      // NOTE: 
      // Msys is measure[i].M + zp corrections
      // Mcal is image[j].Mcal
      // Mmos and Mgrid are offsets for mosaic and grid

      // tie down reference photometry if the -refcode (code) option is selected
      // eg, -refcode g_SDSS
      // this probably makes no sense in the context of multifilter analysis
      // XXX probably need to use the photocde table to assign reference mag weights.
      if (refPhotcode) {
	if (code->code == refPhotcode->code) {
	  wpsflist[Npsf] = UBERCAL_WEIGHT;
	}
      }
      Npsf ++;
    }

    int Nminmeas = isSetMrelFinal ? 1 : STAR_TOOFEW + 1;

    // XXX : ugh : another hard-wired photcode entry...
    if (isSetMrelFinal && (pass == 0)) {
      if ((thisCode < 6) || (thisCode == 9)) {
	secfilt[Nsec].Ncode = NexpPS1; 
      } else {
	secfilt[Nsec].Ncode = Ncode; // 2MASS data if it exists
      }
    }

    // when performing the grid analysis, STAR_TOOFEW should be set to 1;
    if (Npsf < Nminmeas) { /* too few measurements */
      // fprintf (f, "%10.6f %10.6f %d %d %d\n", averageT[0].R, averageT[0].D, measureT[0].imageID, Npsf, STAR_TOOFEW); 
      secfilt[Nsec].flags |= ID_STAR_FEW;
      if (Ncode == 0) {
	results->Ncode ++;
      } else {
	results->Nfew ++;
      }
      continue;
    } else {
      secfilt[Nsec].flags &= ~ID_STAR_FEW;
    }	

    if (forceSynth) {
      // use the single SYNTH value instead of the other mags here
      myAssert ((forceSynthEntry < Npsf) && (forceSynthEntry >= 0), "programming error");
      Mpsflist[0] = Mpsflist[forceSynthEntry];
      dpsflist[0] = dpsflist[forceSynthEntry];
      wpsflist[0] = wpsflist[forceSynthEntry];
      Npsf = 1;
    }
    liststats (Mpsflist, dpsflist, wpsflist, Npsf, psfstats);

    secfilt[Nsec].M  	 = psfstats->mean;
    secfilt[Nsec].dM 	 = psfstats->error;
    secfilt[Nsec].Mchisq = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN;

    // when running -averages, we have no information about the images, so we cannot set this
    if (minUbercalDist > -1) {
      secfilt[Nsec].ubercalDist = minUbercalDist;
    }

    if (isSetMrelFinal) {
      found[Nsec] = TRUE;

      secfilt[Nsec].Mstdev = psfstats->sigma; // Mstdev is in millimags (not enough space for more precision)
      // secfilt[Nsec].Ncode = Ncode;
      secfilt[Nsec].Nused = psfstats->Nmeas;

      secfilt[Nsec].Mmax = psfstats->max;
      secfilt[Nsec].Mmin = psfstats->min;

      // NOTE : use the modified weight for apmags as well as psf mags
      liststats (Maplist, daplist, waplist, Nap, apstats);
      secfilt[Nsec].Map  = Nap > 0 ? apstats->mean : NAN; 
      secfilt[Nsec].dMap  = Nap > 0 ? apstats->error : NAN; 
      secfilt[Nsec].sMap  = Nap > 0 ? apstats->sigma : NAN; 
      secfilt[Nsec].NusedAp  = Nap;

      liststats (Mkronlist, dkronlist, wkronlist, Nkron, kronstats);
      secfilt[Nsec].Mkron  = Nkron > 0 ? kronstats->mean  : NAN; 
      secfilt[Nsec].dMkron = Nkron > 0 ? kronstats->error : NAN; 
      secfilt[Nsec].sMkron = Nkron > 0 ? kronstats->sigma : NAN; 
      secfilt[Nsec].NusedKron  = Nkron;

      // NOTE: for 2MASS measurements, Next should be 1, as should N
      if ((Next > 0) && (Next > 0.5*Npsf)) {
	secfilt[Nsec].flags |= ID_SECF_OBJ_EXT;
      }

      switch (pass) {
	case 0:
	  secfilt[Nsec].flags |= ID_PHOTOM_PASS_0;
	  if (havePS1) GoodPS1 = TRUE;
	  break;
	case 1:
	  secfilt[Nsec].flags |= ID_PHOTOM_PASS_1;
	  if (havePS1) GoodPS1 = TRUE;
	  break;
	case 2:
	  secfilt[Nsec].flags |= ID_PHOTOM_PASS_2;
	  if (havePS1) GoodPS1 = TRUE;
	  break;
	case 3:
	  secfilt[Nsec].flags |= ID_PHOTOM_PASS_3;
	  break;
	case 4:
	  secfilt[Nsec].flags |= ID_PHOTOM_PASS_4;
	  break;
      }
      if (haveSynth) {
	secfilt[Nsec].flags |= ID_SECF_USE_SYNTH;
      }	
      if (havePS1) {
	secfilt[Nsec].flags |= ID_SECF_HAS_PS1;
      }	
      if (haveUbercal) {
	secfilt[Nsec].flags |= ID_SECF_USE_UBERCAL;
      }	
    }
  }

  if (isSetMrelFinal) {
    if (pass == 0) {
      DVOAverageFlags flagBits = ID_OBJ_EXT | ID_OBJ_EXT_ALT | ID_OBJ_GOOD | ID_OBJ_GOOD_ALT;
      // we attempt to set a few flags here; reset those bits before trying:
      average[0].flags &= ~flagBits;
    }

    if (NextPS1 && (NextPS1 > NpsfPS1)) {
      average[0].flags |= ID_OBJ_EXT;
    }
    if (GoodPS1) {
      average[0].flags |= ID_OBJ_GOOD;
    }
    if (Galaxy2MASS) {
      average[0].flags |= ID_OBJ_EXT_ALT;
    }
    if (Good2MASS) {
      average[0].flags |= ID_OBJ_GOOD_ALT;
    }
  }
  return (TRUE);
}

// only apply Stack operation on setMrelFinal in first pass 
int setMrelAverageStack (off_t measureOffset, int cat, FlatCorrectionTable *flatcorr, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, char *primaryCell) {

  off_t k, ID;

  float Msys = 0, Mcal= 0, Mmos = 0, Mgrid = 0;

  // set the primary projection cell and skycell for this coordinate
  int tessID    = -1;
  int projectID = -1;
  int skycellID = -1;
  get_tess_ids(&tessID, &projectID, &skycellID, average[0].R, average[0].D);

  int NstackGood = 0;
  int NstackSuspect = 0;
  int haveStackObject = FALSE;

  int Ns;
  for (Ns = 0; Ns < Nphotcodes; Ns++) {

    int thisCode = photcodes[Ns][0].code;
    int Nsec = GetPhotcodeNsec(thisCode);

    /* calculate the average mag in this SEC photcode for a single star */

    int haveStack = FALSE;

    // need to find the measurement closest to the center of its skycell, as well as the
    // closest for the subset of primary projection cells

    float stackCenterOffsetMin = 1e9;
    off_t stackCenterMeasureMin = -1;

    float stackPrimaryOffsetMin = 1e9;
    off_t stackPrimaryMeasureMin = -1;

    int isBad = FALSE;
    int isSuspect = FALSE;

    off_t meas = measureOffset;
    for (k = 0; k < average[0].Nmeasure; k++, meas++) {

      // skip measurements that do not match the current photcode
      PhotCode *code = GetPhotcodebyCode (measure[k].photcode);
      if (!code) continue;
      if (code->equiv != thisCode) { continue; }

      // clear this bit for all measurements
      measure[k].dbFlags &= ~ID_MEAS_STACK_PRIMARY;

      // only examine gpc1 stack data
      if (!isGPC1stack(measure[k].photcode)) continue;

      haveStack = TRUE;
      haveStackObject = TRUE;
      
      if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS_STACK(Nbad); 

      if (getImageEntry (meas, cat) < 0) continue;
      // measurements without an image are either external reference photometry or
      // data for which the associated image has not been loaded (probably because of
      // overlaps).  we only want measurements associated with stack images in this loop

      Mcal  = getMcal_alt  (meas, cat, flatcorr, measure[k].Xccd, measure[k].Yccd);
      if (isnan(Mcal))  SKIP_THIS_MEAS_STACK(Ncal);
      Mmos  = getMmos  (meas, cat);
      if (isnan(Mmos))  SKIP_THIS_MEAS_STACK(Nmos);
      Mgrid = getMgrid (meas, cat);
      if (isnan(Mgrid)) SKIP_THIS_MEAS_STACK(Ngrid);

      Msys = PhotSys (&measure[k], &average[0], &secfilt[0], MAG_CLASS_PSF);
      if (isnan(Msys)) SKIP_THIS_MEAS_STACK(Nsys);

      unsigned int stackImageID;

      // which stack image should we use for the mean value?
      // if we request the primary (USE_TREE_FOR_PRIMARY), then find the min distances for data from the primary cell
      if (MatchImageSkycellID (meas, cat, tessID, projectID, skycellID)) {
	float stackPrimaryOffset = getCenterOffset (meas, cat, &measure[k], &stackImageID);
	if (stackPrimaryOffset < stackPrimaryOffsetMin) {
	  stackPrimaryOffsetMin = stackPrimaryOffset;
	  // stackPrimaryIDmin = stackImageID;
	  stackPrimaryMeasureMin = k;
	}
      }

      // get the center distance for the generic case:
      float stackCenterOffset = getCenterOffset (meas, cat, &measure[k], &stackImageID);
      if (stackCenterOffset < stackCenterOffsetMin) {
	stackCenterOffsetMin = stackCenterOffset;
	// stackCenterIDmin = stackImageID;
	stackCenterMeasureMin = k;
      }

      // data quality assessment
      isBad |= (measure[k].photFlags & code->photomBadMask);
      isBad |= (measure[k].psfQF < 0.85);
      isBad |= isnan(measure[k].psfQF);
      isBad |= measure[k].dM > 0.2; // S/N < 5.0

      isSuspect |= (measure[k].photFlags & code->photomPoorMask);
      isSuspect |= (measure[k].psfQFperf < 0.85);
    }

    if (!haveStack) continue;

    if (!isSuspect && !isBad) {
      NstackGood ++;
    }
    if (isSuspect && !isBad) {
      NstackSuspect ++;
    }

    // measurements which are bad will not have a valid stack entry and are skipped
    k = (stackPrimaryMeasureMin >= 0) ? stackPrimaryMeasureMin : stackCenterMeasureMin;
    if (k < 0) continue;


    // we are now populating stackDetectID not stack Image ID in secfilt
    // ID = (stackPrimaryMeasureMin >= 0) ? stackPrimaryIDmin      : stackCenterIDmin;
    ID = measure[k].extID; // for the stack, this is the stackDetectID

    // get the zero point for the selected image
    float zp = PhotZeroPoint (&measure[k], &average[0], &secfilt[0]) - (Mcal + Mmos + Mgrid);

    // flux_cgs : erg sec^1 cm^-2 Hz^-1
    // mag_inst : -2.5 log (cts/sec)
    // mag_inst : -2.5 log (flux_inst)
    // flux_inst = ten(-0.4*mag_inst)

    // mag_AB = -2.5 log (flux_cgs) - 48.6 (~by definition) [~Vega flux in V-band]
    // flux_cgs = ten(-0.4*(mag_AB + 48.6))

    // flux_AB : ten(-0.4*mag_AB)

    // flux_cgs = ten(-0.4*48.6) * flux_AB
    // flux_AB  = ten(+0.4*48.6) * flux_cgs
	    
    // flux_Jy : flux_cgs * 10^23

    // flux_AB = ten(+0.4*48.6) * ten(-23) * flux_Jy

    // mag_AB = mag_inst + ZP

    // flux_inst = ten(-0.4*(mag_AB - ZP)) = ten(0.4*ZP) * flux_AB

    // flux_AB = flux_inst * ten(-0.4*ZP)

    // flux_inst * ten(-0.4*ZP) = ten(+0.4*48.6 - 23) * flux_Jy

    // flux_inst = flux_Jy * ten(0.4*ZP + 0.4*48.6 - 23)
    // flux_inst = flux_Jy * ten(0.4*ZP - 3.56)
    // flux_Jy = flux_inst * ten(-0.4*ZP + 3.56)

    // zpFactor to go from instrumental flux to Janskies
    float zpFactor = pow(10.0, -0.4*zp + 3.56);

    // need to put in AB mag factor to get to Janskies (or uJy?)
    secfilt[Nsec].FpsfStk   = zpFactor * measure[k].FluxPSF;  
    secfilt[Nsec].dFpsfStk  = zpFactor * measure[k].dFluxPSF; 
    secfilt[Nsec].FkronStk  = zpFactor * measure[k].FluxKron; 
    secfilt[Nsec].dFkronStk = zpFactor * measure[k].dFluxKron;
    secfilt[Nsec].FapStk    = zpFactor * measure[k].FluxAp;
    secfilt[Nsec].dFapStk   = zpFactor * measure[k].dFluxAp;

    secfilt[Nsec].MpsfStk   = (measure[k].FluxPSF  > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FpsfStk) : NAN;
    secfilt[Nsec].MkronStk  = (measure[k].FluxKron > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FkronStk) : NAN;
    secfilt[Nsec].MapStk    = (measure[k].FluxAp   > 0.0) ? 8.9 - 2.5*log10(secfilt[Nsec].FapStk) : NAN;

    secfilt[Nsec].stackDetectID = ID;

    measure[k].dbFlags |= ID_MEAS_STACK_PHOT_SRC;
    if (stackPrimaryMeasureMin >= 0) {
      measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY;
      secfilt[Nsec].flags |= ID_SECF_STACK_PRIMARY;
    }

    secfilt[Nsec].flags |= ID_SECF_HAS_STACK;
  } // Nsecfilt loop

  if (NstackGood >= 2) {
    average[0].flags |= ID_OBJ_GOOD_STACK;
  } else if (NstackGood + NstackSuspect >= 2) {
    average[0].flags |= ID_OBJ_SUSPECT_STACK;
  } else if (haveStackObject) {
    average[0].flags |= ID_OBJ_BAD_STACK;
  }
  return (TRUE);
}

# undef SKIP_THIS_MEAS
# define SKIP_THIS_MEAS(REASON) {				\
    measure[k].dbFlags |= ID_MEAS_SKIP_PHOTOM;	\
    results->REASON ++;						\
    continue; }

// set mean of forced-warp measurements (selected by photcode range for now):
// somewhat simplified relative to chip-photometry:
// * no grid, no mosaic, no 2MASS, no SYNTH, no Ubercal, no flatcorr
int setMrelAverageForcedWarp (off_t measureOffset, int cat, int pass, SetMrelInfo *results, Average *average, SecFilt *secfilt, Measure *measure, off_t *found) {

  off_t k;
  float Fsys = 0, Mcal= 0;

  // we are measuring means for 3 types of FLUXes: psf, ap, kron.  I am using the psf mag
  // error for the ap mags, but krons have their own errors.  it is an open question if I
  // should be doing weighted or unweighted fits (this is a user option)
  double *Fpsflist  = results->Mpsflist;
  double *dpsflist  = results->dpsflist;
  double *wpsflist  = results->wpsflist;

  double *Faplist   = results->Maplist;
  double *daplist   = results->daplist;
  double *waplist   = results->waplist;

  double *Fkronlist = results->Mkronlist;
  double *dkronlist = results->dkronlist;
  double *wkronlist = results->wkronlist;

  StatType *psfstats  = &results->psfstats;
  StatType *apstats   = &results->apstats;
  StatType *kronstats = &results->kronstats;

  // option for a test print
  if (FALSE && (average[0].objID == 0x7146) && (average[0].catID == 0x49d8)) {
    fprintf (stderr, "test obj\n");
    print_measure_set_alt (average, secfilt, measure);
  }

  int Ns;
  for (Ns = 0; Ns < Nphotcodes; Ns++) {

    int thisCode = photcodes[Ns][0].code;
    int Nsec = GetPhotcodeNsec(thisCode);

    /* calculate the average mag in this SEC photcode for a single star */

    /* star/photcodes already calibrated */
    if (found[Nsec]) continue;  
      
    off_t meas = measureOffset;

    int Nap = 0;
    int Npsf = 0;
    int Nkron = 0;
    for (k = 0; k < average[0].Nmeasure; k++, meas++) {

      // only examine gpc1 forced-warp data
      if (!isGPC1warp(measure[k].photcode)) continue;

      // skip measurements that do not match the current photcode
      PhotCode *code = GetPhotcodebyCode (measure[k].photcode);
      if (!code) continue;
      if (code->equiv != thisCode) { continue; }

      if (measure[k].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad); 

      if (getImageEntry (meas, cat) < 0) {
	// measurements without an image are either external reference photometry or
	// data for which the associated image has not been loaded (probably because of
	// overlaps).  Msys + measure.Mcal is our best guess of the true magnitude
	Mcal = measure[k].Mcal; // check that this is zero for loaded REF value
      } else {
	// use getMcal not getMcal_alt?
	Mcal  = getMcal_alt (meas, cat, NULL, measure[k].Xccd, measure[k].Yccd);
	// Mcal  = getMcal (meas, cat);
	if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
      }
      float Fcal = MagToFlux(-Mcal);

      // in the calculations below, 
      // ...list gives the error per measurement, wlist gives the weight
      // we can modify the error and weight in a few ways:
      // 1) MIN_ERROR guarantees a floor
      // 2) photomErrSys is added in quadrature as a sytematic error, set per photcode

      // skip some absurd values NAN
      // NOTE : I am using PhotCat not PhotSys for now since GPC1 chip-to-chip color terms
      // are small (and not measured)
      float Fpsf = PhotFluxCat (&measure[k], MAG_CLASS_PSF);
      if (isnan(Fsys)) SKIP_THIS_MEAS(Nsys);
      // if (Msys <  0.0) SKIP_THIS_MEAS(Nsys);
      // if (Msys > 30.0) SKIP_THIS_MEAS(Nsys);

      float dFpsf = PhotFluxCatErr (&measure[k], MAG_CLASS_PSF);
      dFpsf = MAX (dFpsf, MIN_ERROR*Fpsf); // MIN_ERROR is a fractional error
      Fpsflist[Npsf] = Fpsf * Fcal;
      dpsflist[Npsf] = dFpsf * Fcal;
      wpsflist[Npsf] = 1.0;
      Npsf ++;

      float Fap = PhotFluxCat (&measure[k], MAG_CLASS_APER);
      float dFap = PhotFluxCatErr (&measure[k], MAG_CLASS_APER);
      if (!isnan(Fap)) {
	Faplist[Nap] = Fap * Fcal;
	daplist[Nap] = dFap * Fcal;
	waplist[Nap] = 1.0; // drop weight lists?
	Nap ++;
      }

      float Fkron = PhotFluxCat (&measure[k], MAG_CLASS_KRON);
      float dFkron = PhotFluxCatErr (&measure[k], MAG_CLASS_KRON);
      if (!isnan(Fkron)) {
	Fkronlist[Nkron] = Fkron * Fcal;
	dkronlist[Nkron] = dFkron * Fcal;
	wkronlist[Nkron] = 1.0;
	Nkron ++;
      }
    }
    if (Npsf < 1) continue;

    found[Nsec] = TRUE;

    liststats (Fpsflist, dpsflist, wpsflist, Npsf, psfstats);

    secfilt[Nsec].FpsfWrp  = psfstats->mean;
    secfilt[Nsec].dFpsfWrp = psfstats->error;
    secfilt[Nsec].sFpsfWrp = psfstats->sigma; // Mstdev is in millimags (not enough space for more precision)
    secfilt[Nsec].NusedWrp = psfstats->Nmeas;
    secfilt[Nsec].MpsfWrp  = isnan(secfilt[Nsec].FpsfWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FpsfWrp); // 8.9 since flux is in Jy

    // NOTE : use the modified weight for apmags as well as psf mags
    liststats (Faplist, daplist, waplist, Nap, apstats);
    secfilt[Nsec].FapWrp     = Nap > 0 ? apstats->mean : NAN; 
    secfilt[Nsec].dFapWrp    = Nap > 0 ? apstats->error : NAN; 
    secfilt[Nsec].sFapWrp    = Nap > 0 ? apstats->sigma : NAN; 
    secfilt[Nsec].NusedApWrp = Nap;
    secfilt[Nsec].MapWrp  = isnan(secfilt[Nsec].FapWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FapWrp); // 8.9 since flux is in Jy

    liststats (Fkronlist, dkronlist, wkronlist, Nkron, kronstats);
    secfilt[Nsec].FkronWrp     = Nkron > 0 ? kronstats->mean  : NAN; 
    secfilt[Nsec].dFkronWrp    = Nkron > 0 ? kronstats->error : NAN; 
    secfilt[Nsec].sFkronWrp    = Nkron > 0 ? kronstats->sigma : NAN; 
    secfilt[Nsec].NusedKronWrp = Nkron;
    secfilt[Nsec].MkronWrp  = isnan(secfilt[Nsec].FkronWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FkronWrp); // 8.9 since flux is in Jy
  }
  return (TRUE);
}

int setGlobalObjStats (Average *average, Measure *measure) {

  off_t k;

  float psfQFmax     = 0.0;
  float psfQFperfmax = 0.0;
  float stargalmax   = 0.0;

  for (k = 0; k < average[0].Nmeasure; k++) {

    if (isfinite(measure[k].psfQF)) {
      psfQFmax     = MAX (psfQFmax,     measure[k].psfQF);
    }
    if (isfinite(measure[k].psfQFperf)) {
      psfQFperfmax = MAX (psfQFperfmax, measure[k].psfQFperf);
    }

    // only count psfQF > 0.85 here
    if (isfinite(measure[k].extNsigma) && (measure[k].psfQF > 0.85)) {
      stargalmax   = MAX (stargalmax,   measure[k].extNsigma);
    }
  }
  average[0].psfQF     = psfQFmax;
  average[0].psfQFperf = psfQFperfmax;
  average[0].stargal   = stargalmax;

  return (TRUE);
}
