# include "relphot.h"

# if (0)
# define TEST_OBJ_ID 0x0000000e
# define TEST_CAT_ID 0x000076ee
# else
# define TEST_OBJ_ID 0
# define TEST_CAT_ID 0
# endif

int magStatsByRankingClipped (StatDataSet *dataset, StatType *stats);
int magStatsByRankingIRLS (StatDataSet *dataset, StatType *stats);
int markMeasureByRanking (StatDataSet *dataset, Measure *measure, int minrank, DVOMeasureFlags keepflag, DVOMeasureFlags maskflag);
void GetPhotFlagStats (uint32_t *photFlagUpper, uint32_t *photFlagLower, uint32_t *photflag_list, int Nphotflag);
void sort_entry_by_offset (double *offset, int *entry, int N);
void sort_StatDataSet (StatDataSet *dataset);
int liststats_irls (StatDataSet *dataset, int Npoints, StatType *stats);

# define MAG_STATS_BY_RANKING magStatsByRankingIRLS

# define UBERCAL_WEIGHT 100.0

# define SKIP_THIS_MEAS(REASON) {			\
    results->REASON ++;					\
    continue; }

# define SKIP_THIS_MEAS_STACK(REASON) {			\
    continue; }

# define CHECK_VALID_MAG(MAG,D_MAG) (isfinite(MAG) && isfinite(D_MAG) && (MAG > -5.0) && (MAG < 30.0))
# define CHECK_VALID_FLUX(FLUX,D_FLUX) (isfinite(FLUX) && isfinite(D_FLUX))

// set A to B if B > A or if A is NAN, as long as B is not NAN
# define MAX_NOT_NAN(A,B) { if (isfinite(B) && (!isfinite(A) || (B > A))) A = B; }

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 setMrel_catalog_alt (Catalog *catalog, int Nc, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {

  off_t j;

  // the same STATMODE is used for chip and warp averages
  liststats_setmode (&results->psfstats,  STATMODE);
  liststats_setmode (&results->apstats,   STATMODE);
  liststats_setmode (&results->kronstats, STATMODE);

  SetMrelInfoReset (results); // reset the countesrs

  for (j = 0; j < catalog[Nc].Naverage; j++) {

    if (STAGES & STAGE_CHIP) {
      setMrelAverageExposure (&catalog[Nc], Nc, j, Nsecfilt, isSetMrelFinal, flatcorr, results);
    }

    // only apply Stack operation on setMrelFinal in first pass 
    if (isSetMrelFinal && (STAGES & STAGE_STACK) && !IS_DIFF_DB) {
      setMrelAverageStack (&catalog[Nc], Nc, j, Nsecfilt, flatcorr);
    }

    // only measure force-warp mean values if issetMrelFinal (make it optional?)
    if (isSetMrelFinal && (STAGES & STAGE_WARP)) {
      setMrelAverageForcedWarp (&catalog[Nc], Nc, j, Nsecfilt, flatcorr, results);
    }
  }
  return (TRUE);
}

// 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.

// 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

// set mean of chip measurements (selected by photcode range for now):
int setMrelAverageExposure (Catalog *catalog, int cat, off_t ave, int Nsecfilt, int isSetMrelFinal, FlatCorrectionTable *flatcorr, SetMrelInfo *results) {

  int Nsec;
  off_t k;

  // we are guaranteed to have averageT, but not average
  AverageTiny *averageT = &catalog[0].averageT[ave];

  off_t measureOffset = averageT->measureOffset;
  MeasureTiny *measureT = &catalog[0].measureT[measureOffset];

  // we are NOT guaranteed to have average, measure, secfilt
  Average *average  = catalog[0].average     ? &catalog[0].average[ave]               : NULL;
  Measure *measure  = catalog[0].measure     ? &catalog[0].measure[measureOffset]     : NULL;
  SecFilt *secfilt  = catalog[0].secfilt     ? &catalog[0].secfilt[ave*Nsecfilt]      : NULL;
  char *measureRank = catalog[0].measureRank ? &catalog[0].measureRank[measureOffset] : NULL;

  // we are measuring means for 3 types of mags: psf, ap, kron.

  // in the final assignment, set the mean mags even if only 1 measurements exists
  int Nminmeas = isSetMrelFinal ? 1 : STAR_TOOFEW + 1;

  // isSetMrelFinal : in the final pass, we set psf, kron, ap mag values and extra stats; other passes only do psf mags

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

  SetMrelInfoResetObject (results); // reset the per-object arrays

  int NextPS1 = 0;
  int NpsfPS1 = 0;
  int GoodPS1 = FALSE;
  int Good2MASS = FALSE;
  int Galaxy2MASS = FALSE;
  int haveTYCHO = FALSE;
  int haveHSC = FALSE;
  int haveCFH = FALSE;

  float stargalmax = 0.0;

  // assign measurements to the photcode lists 
  for (k = 0; k < averageT[0].Nmeasure; k++) {

    // these bits should not be set unless we use them in this pass
    // (note that we can only un-set them in the final pass when we actually have measure, not just measureT
    if (measure) measure[k].dbFlags &= ~(ID_MEAS_PHOTOM_PSF | ID_MEAS_PHOTOM_APER | ID_MEAS_PHOTOM_KRON);

    // skip measurements that do not have a valid photcode (raise exception?)
    PhotCode *code = GetPhotcodebyCode (measureT[k].photcode);
    myAssert (code, "invalid photcode??");

    // SKIP gpc1 stack data (hard-wired photcodes)
    if (isGPC1stack(measureT[k].photcode)) continue;
    
    // SKIP gpc1 forced-warp data (hard-wired photcodes)
    if (isGPC1warp(measureT[k].photcode)) continue;
    
    if (isTYCHO(measureT[k].photcode)) { haveTYCHO = TRUE; }
    if (isHSCchip(measureT[k].photcode)) { haveHSC = TRUE; }
    if (isCFHchip(measureT[k].photcode)) { haveCFH = TRUE; }

    if (is2MASS(measureT[k].photcode)) {
      if (measureT[k].photFlags & 0x00c00000) {
	Galaxy2MASS = TRUE; // per object value
      }
      if (measureT[k].photFlags & 0x00000007) {
	Good2MASS = TRUE; // per object value 
      }
    }

    // assign the Nsec value so we can assign to the right lists
    int Nsec = GetPhotcodeNsec (code->equiv);
    if (Nsec < 0) continue; // skip measurements which do not have an equiv average photcode

    results->Nmeas[Nsec] ++;

    // various steps only relevant to the final pass 
    if (isSetMrelFinal) {
      MAX_NOT_NAN (results->psfQfMax[Nsec],     measure[k].psfQF);
      MAX_NOT_NAN (results->psfQfPerfMax[Nsec], measure[k].psfQFperf);

      // are we a PS1 exposure photcode? (hard-wired photcodes)
      if (isGPC1chip(measure[k].photcode)) { 
	results->NexpPS1[Nsec] ++; 
	results->havePS1[Nsec] = TRUE; 
	
	// only count psfQF > 0.85 here
	if (isfinite(measure[k].extNsigma) && isfinite(measure[k].psfQF) && (measure[k].psfQF > 0.85)) {
	  stargalmax = MAX (stargalmax, measure[k].extNsigma);
	  results->stargal_list[results->Nstargal] = measure[k].extNsigma;
	  results->Nstargal ++;
	}

	// for all GPC1 measurements, track the photflags
	results->photflag_list[results->Nphotflags] = measure[k].photFlags;
	results->Nphotflags ++;
      }

      // force the use of SYN even if we have PS1 mags?
      if (isGPC1synth(measure[k].photcode)) {
	results->haveSYN[Nsec] = TRUE;
	results->measSYN[Nsec] = k;
	measureT[k].dbFlags |= ID_MEAS_SYNTH_MAG; // redundant with photcode
	// note that synthetic mags have the real mags in measure.M (no zero point offset)
	if (measure[k].M < results->minSYN[Nsec]) results->needSYN[Nsec] = TRUE;
	continue;
      }
    }

    // ** Choose the calibration (depends on the mode : do I have an image reference or not?) 
    float Mcal = 0.0, Mmos = 0.0, Mgrid = 0.0, Mflat = 0.0;
    off_t meas = measureOffset + k;
    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 {
      // getMcal_alt returns image[].Mcal modified by flatcorr(image.photom_map_id,x,y)
      // NOTE: getMcal_alt does not include measure.Mflat
      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);
      Mflat = isnan (measureT[k].Mflat) ? 0.0 : measureT[k].Mflat;
    }

    int myUbercalDist = getUbercalDist(meas, cat);
    results->minUbercalDist[Nsec] = MIN(results->minUbercalDist[Nsec], myUbercalDist);

    int isUbercal = (measureT[k].dbFlags & ID_MEAS_PHOTOM_UBERCAL);

    if (isUbercal) results->haveUbercal[Nsec] = TRUE; // haveUbercal is set per secfilt, isUbercal is per measure XXX define this array

    int useUbercalWeight = isUbercal || (refPhotcode && (code->code == refPhotcode->code));

    float Map  = NAN;
    float dMap = NAN;
    if (isSetMrelFinal) {
      Map = PhotCat (&measure[k], MAG_CLASS_APER);

      // NOTE: for PV3, apFluxErr is broken (see pmSourcePhotometry.c:245).  we are going to
      // use PSF flux error instead here:
      // dMap = MAX (hypot(measure[k].dMap,   code->photomErrSys), MIN_ERROR);
      dMap = MAX (hypot(measure[k].dM,   code->photomErrSys), MIN_ERROR);

      if (CHECK_VALID_MAG(Map, dMap)) {
	int Nap = results->aperData[Nsec].Nlist;
	results->aperData[Nsec].flxlist[Nap] = Map - Mflat - Mcal - Mmos - Mgrid; // this is consistent with PhotRel
	results->aperData[Nsec].errlist[Nap] = dMap;
	results->aperData[Nsec].wgtlist[Nap] = useUbercalWeight ? UBERCAL_WEIGHT : 1.0;
	results->aperData[Nsec].ranking[Nap] = measureRank[k];
	results->aperData[Nsec].measSeq[Nap] = k;
	results->aperData[Nsec].msklist[Nap] = 0;
	results->aperData[Nsec].Nlist ++;
      }
    }

    float Mkron, dMkron;
    if (isSetMrelFinal) {
      Mkron = PhotCat (&measure[k], MAG_CLASS_KRON);
      dMkron = MAX (hypot(measure[k].dMkron, code->photomErrSys), MIN_ERROR);
      if (CHECK_VALID_MAG(Mkron, dMkron)) {
	int Nkron = results->kronData[Nsec].Nlist;
	results->kronData[Nsec].flxlist[Nkron] = Mkron - Mflat - Mcal - Mmos - Mgrid; // this is consistent with PhotRel
	results->kronData[Nsec].errlist[Nkron] = dMkron;
	results->kronData[Nsec].wgtlist[Nkron] = isUbercal ? UBERCAL_WEIGHT : 1.0;
	results->kronData[Nsec].ranking[Nkron] = measureRank[k];
	results->kronData[Nsec].measSeq[Nkron] = k;
	results->kronData[Nsec].msklist[Nkron] = 0;
	results->kronData[Nsec].Nlist ++;
      }
    }

    float Mpsf = PhotSysTiny (&measureT[k], &averageT[0], &secfilt[0], MAG_CLASS_PSF);
    float dMpsf  = MAX (hypot(measureT[k].dM, code->photomErrSys), MIN_ERROR);
    if (CHECK_VALID_MAG(Mpsf, dMpsf)) {
      int Npsf = results->psfData[Nsec].Nlist;
      results->psfData[Nsec].flxlist[Npsf] = Mpsf - Mflat - Mcal - Mmos - Mgrid; // this is consistent with PhotRel
      results->psfData[Nsec].errlist[Npsf] = dMpsf;
      results->psfData[Nsec].wgtlist[Npsf] = isUbercal ? UBERCAL_WEIGHT : 1.0;
      results->psfData[Nsec].ranking[Npsf] = measureRank ? measureRank[k] : 0;
      results->psfData[Nsec].measSeq[Npsf] = k;
      results->psfData[Nsec].msklist[Npsf] = 0;
      results->psfData[Nsec].Nlist ++;
    }

    // count the extended detections (all PS1 bands)
    if (isSetMrelFinal && isGPC1chip(measureT[k].photcode) && !isnan(Map) && !isnan(Mpsf)) {
      float dMagAp = Mpsf - Map;
      float SigmaAp = hypot(0.1, 2.5*dMpsf);
      // 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) {
	results->Next[Nsec] ++;
	NextPS1 ++;
      } else {
	NpsfPS1 ++;
      }
    }
  }

  SynthZeroPoints *synthzpts = SynthZeroPointsGet();

  // find max values across filters
  float psfQfMax     = 0.0;
  float psfQfPerfMax = 0.0;

  // now calculate the mean stats for the Nsec bands.
  for (Nsec = 0; Nsec < Nsecfilt; Nsec++) {

    // if we detected this object in PS1, or do not request -preserve-ps1, keep the mean photometry values
    if (!PRESERVE_PS1 || !(secfilt[Nsec].flags & ID_SECF_HAS_PS1)) {
      dvo_secfilt_init (&secfilt[Nsec], SECFILT_RESET_CHIP); // this does not reset astrometry or STACK bits
    }

    if (haveTYCHO) {
      secfilt[Nsec].flags |= ID_SECF_HAS_TYCHO;
    }
    if (haveHSC) {
      secfilt[Nsec].flags |= ID_SECF_HAS_HSC;
    }
    if (haveCFH) {
      secfilt[Nsec].flags |= ID_SECF_HAS_CFH;
    }

    if (PRESERVE_PS1 && (secfilt[Nsec].flags & ID_SECF_HAS_PS1)) continue;
    // if -preserve-ps1 is set and this object has PS1 data, skip the rest of the steps:

    // XXX hardwired grizy = (01234) JHK = (567) w = (8)
    if ((Nsec < 5) || (Nsec == 8)) {
      secfilt[Nsec].Ncode = results->NexpPS1[Nsec]; 
    } else {
      secfilt[Nsec].Ncode = results->Nmeas[Nsec]; // 2MASS data if it exists
    }

    // use the synthetic magnitude for mean photometry (apply zpt calibration if no tycho
    // photometry exists)
    if (results->needSYN[Nsec] || (!results->havePS1[Nsec] && results->haveSYN[Nsec])) {
      if (!isSetMrelFinal) continue;

      // use the single SYNTH value instead of the other mags here
      float Mpsf = measure[results->measSYN[Nsec]].M;
      float ZP = 0.0;

      // if we have loaded the synth zpt correction, and do NOT have a tycho measurement, apply the correction
      if (!haveTYCHO && synthzpts) {
	// need to look up the (X,Y) coords from this (R,D) location
	double X, Y;
	double R = ohana_normalize_angle_to_midpoint (measure->R, 0.0); // XXX 0.0 or 180.0 for center?
	RD_to_XY (&X, &Y, R, measure->D, &synthzpts->coords);
	if (X < 0) return FALSE;
	if (Y < 0) return FALSE;
	if (X >= synthzpts->Nx) return FALSE;
	if (Y >= synthzpts->Ny) return FALSE;
	
	int Xpix = X;
	int Ypix = Y;
	int Npix = Xpix + Ypix*synthzpts->Nx;
	
	float *value = (float *) synthzpts->matrix[Nsec].buffer;
	ZP = !isnan(value[Npix]) ? value[Npix] : 0.0;
	if (isfinite(value[Npix])) {
	  secfilt[Nsec].flags |= ID_SECF_FIX_SYNTH;
	}
      }
      secfilt[Nsec].M      = Mpsf + ZP;
      secfilt[Nsec].dM     = 0.6;
      secfilt[Nsec].Mchisq = 0.0;
      secfilt[Nsec].flags |= ID_SECF_USE_SYNTH;
      continue;
    } 

    // if too few valid measurements meet the minimum criteria, go to the next entry
    StatType *psfstats = &results->psfstats;
    int Nranking = MAG_STATS_BY_RANKING (&results->psfData[Nsec], psfstats);
    if (Nranking < Nminmeas) { 
      secfilt[Nsec].flags |= ID_OBJ_FEW;
    } else {
      secfilt[Nsec].M  	 = psfstats->mean;
      secfilt[Nsec].dM 	 = psfstats->error;
      secfilt[Nsec].Mchisq = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN;
    }
    int minRankPSF = (Nranking > 0) ? results->psfData[Nsec].ranking[0] : 10;

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

    if (isSetMrelFinal) {
      if ((average[0].objID == TEST_OBJ_ID) && (average[0].catID == TEST_CAT_ID)) {
	fprintf (stderr, "test obj\n");
      }

      // mark the measurements matching this ranking 
      markMeasureByRanking (&results->psfData[Nsec], measure, minRankPSF, ID_MEAS_PHOTOM_PSF, ID_MEAS_MASKED_PSF);

      if (Nranking) {
	secfilt[Nsec].Mstdev = psfstats->sigma;
	secfilt[Nsec].Nused  = psfstats->Nmeas;
	secfilt[Nsec].Mmax   = psfstats->max;
	secfilt[Nsec].Mmin   = psfstats->min;
      }
      secfilt[Nsec].psfQfMax     = results->psfQfMax[Nsec];
      secfilt[Nsec].psfQfPerfMax = results->psfQfPerfMax[Nsec];

      MAX_NOT_NAN (psfQfMax,     secfilt[Nsec].psfQfMax);
      MAX_NOT_NAN (psfQfPerfMax, secfilt[Nsec].psfQfPerfMax);

      StatType *apstats = &results->apstats;
      Nranking = MAG_STATS_BY_RANKING (&results->aperData[Nsec], apstats);
      if (Nranking) {
	secfilt[Nsec].Map     = apstats->mean; 
	secfilt[Nsec].dMap    = apstats->error; 
	secfilt[Nsec].sMap    = apstats->sigma; 
	secfilt[Nsec].NusedAp = Nranking;
      }
      int minRankAper = (Nranking > 0) ? results->aperData[Nsec].ranking[0] : 10;
      markMeasureByRanking (&results->aperData[Nsec], measure, minRankAper, ID_MEAS_PHOTOM_APER, ID_MEAS_MASKED_APER);

      StatType *kronstats = &results->kronstats;
      Nranking = MAG_STATS_BY_RANKING (&results->kronData[Nsec], kronstats);
      if (Nranking) {
	secfilt[Nsec].Mkron     = kronstats->mean; 
	secfilt[Nsec].dMkron    = kronstats->error; 
	secfilt[Nsec].sMkron    = kronstats->sigma; 
	secfilt[Nsec].NusedKron = Nranking;
      }
      int minRankKron = (Nranking > 0) ? results->kronData[Nsec].ranking[0] : 10;
      markMeasureByRanking (&results->kronData[Nsec], measure, minRankKron, ID_MEAS_PHOTOM_KRON, ID_MEAS_MASKED_KRON);

      // does this object appear extended in > 50% of measurements?
      if ((results->Next[Nsec] > 0) && (results->Next[Nsec] > 0.5*psfstats->Nmeas)) {
	secfilt[Nsec].flags |= ID_SECF_OBJ_EXT;
      }

      switch (minRankPSF) {
	case 0:
	  secfilt[Nsec].flags |= ID_SECF_RANK_0;
	  if (results->havePS1[Nsec]) GoodPS1 = TRUE;
	  break;
	case 1:
	  secfilt[Nsec].flags |= ID_SECF_RANK_1;
	  if (results->havePS1[Nsec]) GoodPS1 = TRUE;
	  break;
	case 2:
	  secfilt[Nsec].flags |= ID_SECF_RANK_2;
	  if (results->havePS1[Nsec]) GoodPS1 = TRUE;
	  break;
	case 3:
	  secfilt[Nsec].flags |= ID_SECF_RANK_3;
	  break;
	case 4:
	  secfilt[Nsec].flags |= ID_SECF_RANK_4;
	  break;
	default:
	  break;
      }
      if (results->havePS1[Nsec]) {
	secfilt[Nsec].flags |= ID_SECF_HAS_PS1;
      }	
      if (results->haveUbercal[Nsec]) {
	secfilt[Nsec].flags |= ID_SECF_USE_UBERCAL;
      }
    }
  }

  if (isSetMrelFinal) {
    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;
    }

    GetPhotFlagStats (&average[0].photFlagsUpper, &average[0].photFlagsLower, results->photflag_list, results->Nphotflags);

    // calculate the median of the stargal (ExtNsigma) values
    dsort(results->stargal_list, results->Nstargal);
    double stargal_median;
    if (results->Nstargal % 2) {
      int Ncenter = results->Nstargal / 2;
      stargal_median = results->stargal_list[Ncenter];
    } else {
      int Ncenter = results->Nstargal / 2 - 1;
      stargal_median = 0.5*(results->stargal_list[Ncenter] + results->stargal_list[Ncenter + 1]);
    }
    average[0].stargal   = stargal_median;

    average[0].psfQF     = psfQfMax;
    average[0].psfQFperf = psfQfPerfMax;
  }
  return (TRUE);
}

// only apply Stack operation on setMrelFinal in first pass 
// this function has 3 goals, not to be confused:
// 1) find and mark the PRIMARY detections (regardless of the quality of the detection)
// 2) select the BEST detections per filter (regardless of PRIMARY)
// 3) apply the zero point and AB->Jy transformations
int setMrelAverageStack (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr) {

  // we are guaranteed to have average, measure, secfilt
  Average *average = &catalog[0].average[ave];
  off_t measureOffset = average->measureOffset;
  int Nmeasure = average->Nmeasure;

  Measure *measure = &catalog[0].measure[measureOffset];
  SecFilt *secfilt = &catalog[0].secfilt[ave*Nsecfilt];

  off_t k;

  float Mcal = 0, Mmos = 0, Mgrid = 0, Finst = 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);

  average[0].tessID       = tessID;
  average[0].projectionID = projectID;
  average[0].skycellID    = skycellID;

  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;

    float psfQFbest = 0.0;

    off_t stackBestMeasure = -1;
    off_t stackPrimaryMeasure = -1;

    int isBad = FALSE;
    int isSuspect = FALSE;

    int Nstack = 0; // number for this photcode
    int NstackDet = 0; // number for this photcode

    secfilt[Nsec].stackBestOff = -1;
    secfilt[Nsec].stackPrmryOff = -1;

    off_t meas = measureOffset;
    for (k = 0; k < Nmeasure; k++, meas++) {

      // only examine gpc1 stack data
      if (!isGPC1stack(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; }

      Nstack ++;
      if ((measure[k].photFlags2 & 0x00000004) == 0) NstackDet ++;

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

      haveStack = TRUE;
      haveStackObject = TRUE;
      
      // ** find the PRIMARY measurement

      // if we request the primary (USE_TREE_FOR_PRIMARY), this is true if the measurement is from the 
      // primary skycell for this position
      if (MatchImageSkycellID (meas, cat, tessID, projectID, skycellID)) {
	stackPrimaryMeasure = k;
	measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY;
	secfilt[Nsec].stackPrmryOff = meas;
	myAssert (secfilt[Nsec].stackPrmryOff <= catalog[0].Nmeasure, "stackPrmryOff out of range");
      }

      // ** now choose the BEST measurements (may also be PRIMARY)

      // ensure that we at least have a single best measure 
      if (stackBestMeasure == -1) stackBestMeasure = k;

      // choose the best psfQFperf value for the BEST measurement
      if (isfinite(measure[k].psfQFperf) && (measure[k].psfQFperf > psfQFbest)) {
	psfQFbest = measure[k].psfQFperf;
	stackBestMeasure = k;
      }
      // ... UNLESS psfQFperf > 0.98 for the primary, in which case just use the primary.
      if ((measure[k].dbFlags & ID_MEAS_STACK_PRIMARY) && isfinite(measure[k].psfQFperf) && (measure[k].psfQFperf > 0.98)) {
	psfQFbest = 1000; // force this to be the best entry
	stackBestMeasure = k;
      }

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

      // 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

      // match measurement to its image (this is just a check, right?)
      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 = measure[k].Mcal; // check that this is zero for loaded REF value
      } else {
	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);
      }

      // NOTE: negative and insignificant flux values are allowed, but not NAN flux values
      Finst = PhotFluxInst (&measure[k], MAG_CLASS_PSF);
      if (isnan(Finst)) SKIP_THIS_MEAS_STACK(Ninst);

      // 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 = (stackBestMeasure >= 0) ? stackBestMeasure : stackPrimaryMeasure;
    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;

    // NOTE: for PV3, apFluxErr is broken (see pmSourcePhotometry.c:245).  we are going to
    // use PSF flux error instead here:
    // secfilt[Nsec].dFapStk   = zpFactor * measure[k].dFluxAp;
    secfilt[Nsec].dFapStk   = zpFactor * measure[k].dFluxPSF;

    // Jy to AB mags
    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].stackBestOff = k + measureOffset;
    myAssert (secfilt[Nsec].stackBestOff <= catalog[0].Nmeasure, "stackBestOff out of range");

    secfilt[Nsec].Nstack    = Nstack;
    secfilt[Nsec].NstackDet = NstackDet;

    // this is the measurement used by secfilt[]
    measure[k].dbFlags |= ID_MEAS_STACK_PHOT_SRC;
    if (k == stackPrimaryMeasure) {
      secfilt[Nsec].flags |= ID_SECF_STACK_PRIMARY;
    }

    // stack measurement used for 'best' was a detection (not forced from the other bands)
    if ((measure[k].photFlags2 & 0x00000004) == 0) {
      secfilt[Nsec].flags |= ID_SECF_STACK_BESTDET;
    }
    // stack measurement used for 'primary' was a detection (not forced from the other bands)
    if ((stackPrimaryMeasure >= 0) && ((measure[stackPrimaryMeasure].photFlags2 & 0x00000004) == 0)) {
      secfilt[Nsec].flags |= ID_SECF_STACK_PRIMDET;
    }

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

  // in PSPS, there will be a stackObject row with the primary values.  if all filter primary values which
  // exist are also best values, then this row will be the best entry

  // the secfilt[Ns] entries are defined to be the Best values, so

  // conversely, if any secfilt[Ns] values are not PRIMARY, then the PRIMARY row is not the BEST row

  int PrimaryIsBest = TRUE;
  for (Ns = 0; Ns < Nphotcodes; Ns++) {
    if (!(secfilt[Ns].flags & ID_SECF_HAS_PS1_STACK)) continue; // no stack detection in PS1, nothing is best
    if (secfilt[Ns].flags & ID_SECF_STACK_PRIMARY) continue;    // primary stack detection is best
    PrimaryIsBest = FALSE;
  }

  if (PrimaryIsBest) {
    average[0].flags |= ID_OBJ_BEST_STACK;
  }

  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);
}

// 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
// analysis is done on flux, not mags (as the faintest objects will be nearly insignificant)
int setMrelAverageForcedWarp (Catalog *catalog, int cat, off_t ave, int Nsecfilt, FlatCorrectionTable *flatcorr, SetMrelInfo *results) {
  OHANA_UNUSED_PARAM(flatcorr);

  // we are guaranteed to have average, measure, secfilt
  Average *average    = &catalog[0].average[ave];
  off_t measureOffset = average->measureOffset;

  Measure *measure    = &catalog[0].measure[measureOffset];
  SecFilt *secfilt    = &catalog[0].secfilt[ave*Nsecfilt];
  char *measureRank   = &catalog[0].measureRank[measureOffset];

  off_t k, Nsec;
  float Mcal = 0;

  SetMrelInfoResetObject (results); // reset the per-object arrays

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

  // for a diff db, we just want primary, if available or not. for warp we want to use the
  // skycell which provided the stack best entry
  int tessID, projID, skycellID;
  if (IS_DIFF_DB) {
    get_tess_ids(&tessID, &projID, &skycellID, average[0].R, average[0].D);
  }
  for (Nsec = 0; Nsec < Nsecfilt; Nsec++) {
    if (!IS_DIFF_DB) {
      myAssert (secfilt[Nsec].stackBestOff <= catalog->Nmeasure, "stackBestOff out of range");
      if (!FindImageSkycellID (secfilt[Nsec].stackBestOff, cat, &tessID, &projID, &skycellID)) continue;
    }
    results->tessID[Nsec]    = tessID;
    results->projID[Nsec]    = projID;
    results->skycellID[Nsec] = skycellID;
  }

  average->NwarpOK = 0; // count the number of warps with psfQF > 0.0

  // assign measurements to the photcode lists 
  for (k = 0; k < average->Nmeasure; k++) {
    off_t meas = measureOffset + k;

    measure[k].dbFlags &= ~ID_MEAS_WARP_USED; // this should not be set unless we use it in this pass

    // skip measurements that do not have a valid photcode (raise exception?)
    PhotCode *code = GetPhotcodebyCode (measure[k].photcode);
    myAssert (code, "invalid photcode??");

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

    // assign the Nsec value so we can assign to the right lists
    int Nsec = GetPhotcodeNsec (code->equiv);
    if (Nsec < 0) continue; // skip measurements which do not have an equiv average photcode

    results->Nmeas[Nsec] ++;
    if (isfinite(measure[k].psfQFperf) && (measure[k].psfQFperf > 0.85)) results->NmeasGood[Nsec] ++;
    if (isfinite(measure[k].psfQF)     && (measure[k].psfQF     > 0.00)) average->NwarpOK ++;

    // use primary skycell for DIFF, stack skycell for WARP
    if (!MatchImageSkycellID (meas, cat, results->tessID[Nsec], results->projID[Nsec], results->skycellID[Nsec])) {
      continue;
    }

    // assign the max psfQF values for each Nsec value
    if (IS_DIFF_DB) {
      measure[k].dbFlags |= ID_MEAS_STACK_PRIMARY;
      MAX_NOT_NAN (results->psfQfMax[Nsec],     measure[k].psfQF);
      MAX_NOT_NAN (results->psfQfPerfMax[Nsec], measure[k].psfQFperf);
    }

    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);

    // 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);
    float dFpsf = PhotFluxCatErr (&measure[k], MAG_CLASS_PSF);
    if (CHECK_VALID_FLUX(Fpsf, dFpsf)) {
      dFpsf = MAX (hypot(dFpsf, code->photomErrSys*Fpsf), MIN_ERROR*Fpsf); // bump up the error by a systematic floor
      int Npsf = results->psfData[Nsec].Nlist;
      results->psfData[Nsec].flxlist[Npsf] = Fpsf * Fcal;
      results->psfData[Nsec].errlist[Npsf] = dFpsf * Fcal;
      results->psfData[Nsec].wgtlist[Npsf] = 1.0;
      results->psfData[Nsec].ranking[Npsf] = measureRank ? measureRank[k] : 0;
      results->psfData[Nsec].measSeq[Npsf] = k;
      results->psfData[Nsec].msklist[Npsf] = 0;
      results->psfData[Nsec].Nlist ++;
    }     

    float Fap = PhotFluxCat (&measure[k], MAG_CLASS_APER);

    // NOTE: for PV3, apFluxErr is broken (see pmSourcePhotometry.c:245).  we are going to
    // use PSF flux error instead here:
    // float dFap = PhotFluxCatErr (&measure[k], MAG_CLASS_APER);
    float dFap = PhotFluxCatErr (&measure[k], MAG_CLASS_PSF);

    if (CHECK_VALID_FLUX(Fap, dFap)) {
      dFap = MAX (hypot(dFap, code->photomErrSys*Fap), MIN_ERROR*Fap); // bump up the error by a systematic floor
      int Naper = results->aperData[Nsec].Nlist;
      results->aperData[Nsec].flxlist[Naper] = Fap * Fcal;
      results->aperData[Nsec].errlist[Naper] = dFap * Fcal;
      results->aperData[Nsec].wgtlist[Naper] = 1.0;
      results->aperData[Nsec].ranking[Naper] = measureRank ? measureRank[k] : 0;
      results->aperData[Nsec].measSeq[Naper] = k;
      results->aperData[Nsec].msklist[Naper] = 0;
      results->aperData[Nsec].Nlist ++;
    }

    float Fkron = PhotFluxCat (&measure[k], MAG_CLASS_KRON);
    float dFkron = PhotFluxCatErr (&measure[k], MAG_CLASS_KRON);
    if (CHECK_VALID_FLUX(Fkron, dFkron)) {
      dFkron = MAX (hypot(dFkron, code->photomErrSys*Fkron), MIN_ERROR*Fkron); // bump up the error by a systematic floor
      int Nkron = results->kronData[Nsec].Nlist;
      results->kronData[Nsec].flxlist[Nkron] = Fkron * Fcal;
      results->kronData[Nsec].errlist[Nkron] = dFkron * Fcal;
      results->kronData[Nsec].wgtlist[Nkron] = 1.0;
      results->kronData[Nsec].ranking[Nkron] = measureRank ? measureRank[k] : 0;
      results->kronData[Nsec].measSeq[Nkron] = k;
      results->kronData[Nsec].msklist[Nkron] = 0;
      results->kronData[Nsec].Nlist ++;
    }
  }

  // find max values across filters
  float psfQfMax     = 0.0;
  float psfQfPerfMax = 0.0;

  // now calculate the mean stats for the Nsec bands.
  for (Nsec = 0; Nsec < Nsecfilt; Nsec++) {
    dvo_secfilt_init (&secfilt[Nsec], SECFILT_RESET_WARP);

    if (IS_DIFF_DB) {
      // for non DIFF_DB, these are set in setMrelAverageExposure
      secfilt[Nsec].psfQfMax     = results->psfQfMax[Nsec];
      secfilt[Nsec].psfQfPerfMax = results->psfQfPerfMax[Nsec];

      MAX_NOT_NAN (psfQfMax,     secfilt[Nsec].psfQfMax);
      MAX_NOT_NAN (psfQfPerfMax, secfilt[Nsec].psfQfPerfMax);
    }

    int Nranking, minRank;

    // if too few valid measurements meet the minimum criteria, go to the next entry
    StatType *psfstats = &results->psfstats;
    Nranking = MAG_STATS_BY_RANKING (&results->psfData[Nsec], psfstats);
    if (Nranking) {
      secfilt[Nsec].FpsfWrp  = psfstats->mean;
      secfilt[Nsec].dFpsfWrp = psfstats->error;
      secfilt[Nsec].sFpsfWrp = (psfstats->Nmeas > 1) ? psfstats->chisq : NAN;
      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
    }
    minRank = (Nranking > 0) ? results->psfData[Nsec].ranking[0] : 10;
    markMeasureByRanking (&results->psfData[Nsec], measure, minRank, ID_MEAS_WARP_USED | ID_MEAS_PHOTOM_PSF, ID_MEAS_MASKED_PSF);

    // if too few valid measurements meet the minimum criteria, go to the next entry
    StatType *apstats = &results->apstats;
    Nranking = MAG_STATS_BY_RANKING (&results->aperData[Nsec], apstats);
    if (Nranking) {
      secfilt[Nsec].FapWrp     = apstats->mean;
      secfilt[Nsec].dFapWrp    = apstats->error;
      secfilt[Nsec].sFapWrp    = (apstats->Nmeas > 1) ? apstats->chisq : NAN;
      secfilt[Nsec].NusedApWrp = apstats->Nmeas;
      secfilt[Nsec].MapWrp     = isnan(secfilt[Nsec].FapWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FapWrp); // 8.9 since flux is in Jy
    }
    minRank = (Nranking > 0) ? results->aperData[Nsec].ranking[0] : 10;
    markMeasureByRanking (&results->aperData[Nsec], measure, minRank, ID_MEAS_WARP_USED | ID_MEAS_PHOTOM_APER, ID_MEAS_MASKED_APER);

    // if too few valid measurements meet the minimum criteria, go to the next entry
    StatType *kronstats = &results->kronstats;
    Nranking = MAG_STATS_BY_RANKING (&results->kronData[Nsec], kronstats);
    if (Nranking) {
      secfilt[Nsec].FkronWrp     = kronstats->mean;
      secfilt[Nsec].dFkronWrp    = kronstats->error;
      secfilt[Nsec].sFkronWrp    = (kronstats->Nmeas > 1) ? kronstats->chisq : NAN;
      secfilt[Nsec].NusedKronWrp = kronstats->Nmeas;
      secfilt[Nsec].MkronWrp     = isnan(secfilt[Nsec].FkronWrp) ? NAN : 8.9 - 2.5*log10(secfilt[Nsec].FkronWrp); // 8.9 since flux is in Jy
    }
    minRank = (Nranking > 0) ? results->kronData[Nsec].ranking[0] : 10;
    markMeasureByRanking (&results->kronData[Nsec], measure, minRank, ID_MEAS_WARP_USED | ID_MEAS_PHOTOM_KRON, ID_MEAS_MASKED_KRON);

    secfilt[Nsec].Nwarp     = results->Nmeas[Nsec];
    secfilt[Nsec].NwarpGood = results->NmeasGood[Nsec];
  }
  return (TRUE);
}

int magStatsByRanking (StatDataSet *dataset, StatType *stats) {

  liststats_init (stats);

  if (dataset->Nlist == 0) return 0;

  // we have a list of measurements for this Nsec value, along with errors, weights, and their rank
  // first, sort by the rank (increasing)
  sort_StatDataSet (dataset);

  // find the values with the same minimum rank:
  int Nranking = 0;
  int minRank = dataset->ranking[0]; // MIN (5, result->psfData[Nsec].ranking[0]);  -- only allow rank 5 or <

  int i;
  for (i = 0; (i < dataset->Nlist) && (dataset->ranking[i] == minRank); i++, Nranking++);

  liststats (dataset->flxlist, dataset->errlist, dataset->wgtlist, Nranking, stats);
  return (Nranking);
}

int magStatsByRankingIRLS (StatDataSet *dataset, StatType *stats) {

  liststats_init (stats);

  if (dataset->Nlist == 0) return 0;

  // we have a list of measurements for this Nsec value, along with errors, weights, and their rank
  // first, sort by the rank (increasing)
  sort_StatDataSet (dataset);

  // find the values with the same minimum rank:
  int Nranking = 0;
  int minRank = dataset->ranking[0]; // MIN (5, result->psfData[Nsec].ranking[0]);  -- only allow rank 5 or <

  int i;
  for (i = 0; (i < dataset->Nlist) && (dataset->ranking[i] == minRank); i++, Nranking++);

  liststats_irls (dataset, Nranking, stats);
  return (Nranking);
}

// outlier warp measurements are driving bad mean values.  
int magStatsByRankingClipped (StatDataSet *dataset, StatType *stats) {

  liststats_init (stats);

  if (dataset->Nlist == 0) return 0;

  // we have a list of measurements for this Nsec value, along with errors, weights, and their rank
  // first, sort by the rank (increasing)
  sort_StatDataSet (dataset);

  // find the values with the same minimum rank:
  int Nranking = 0;
  int minRank = dataset->ranking[0]; // MIN (5, result->psfData[Nsec].ranking[0]);  -- only allow rank 5 or <

  int i;
  for (i = 0; (i < dataset->Nlist) && (dataset->ranking[i] == minRank); i++, Nranking++);

  // we want to measure the weighted average, but first we want to clip points which are > 3 sigma from the median
  // since we are only allowed to clip at most 25% of the points, we need to have at least 4 points to do this clipping
  if (Nranking > 3) {

    double *fluxes = NULL;
    double *offset = NULL;
    int    *entry  = NULL;

    // now find the distances
    ALLOCATE (entry,  int,    Nranking);
    ALLOCATE (offset, double, Nranking);
    ALLOCATE (fluxes, double, Nranking);

    // first, measure the median of the flxlist:
    for (i = 0; i < Nranking; i++) {
      fluxes[i] = dataset->flxlist[i];
    }
    dsort (fluxes, Nranking);
    int Nmidpoint = (int)(0.5*Nranking); // if Nranking is odd (e.g., 5), Nmidpoint is central bin (e.g., 2); if Nranking is even (e.g., 6), Nmidpoint is central bin + 1 (e.g., 3);
    double medianFlux = (Nranking % 2) ? fluxes[Nmidpoint] : 0.5*(fluxes[Nmidpoint] + fluxes[Nmidpoint - 1]);

    for (i = 0; i < Nranking; i++) {
      offset[i] = fabs((dataset->flxlist[i] - medianFlux) / dataset->errlist[i]);
      entry[i] = i;
    }

    sort_entry_by_offset (offset, entry, Nranking);

    // we are only going to clip at most 25% of the points.  Start at the 75% point and check for outliers
    int Nreject = 0;
    for (i = 0.75*Nranking; i < Nranking; i++) {
      if (offset[i] > 4.0) {
	// this is an outlier -- add 10 to the ranking
	dataset->ranking[entry[i]] += 10;
	Nreject ++;
      }
    }
    Nranking -= Nreject; // the excluded points are now removed from the count

    // re-sort by the rank so liststats below only operates on the unrejected points
    sort_StatDataSet (dataset);

    free (fluxes);
    free (offset);
    free (entry);
  }

  // get stats on only the first Nranking entries
  liststats (dataset->flxlist, dataset->errlist, dataset->wgtlist, Nranking, stats);

  return (Nranking);
}

int markMeasureByRanking (StatDataSet *dataset, Measure *measure, int minrank, DVOMeasureFlags keepflag, DVOMeasureFlags maskflag) {
  int i;

  for (i = 0; i < dataset->Nlist; i++) {
    int k = dataset->measSeq[i];
    if (dataset->ranking[i] > minrank) {
      measure[k].dbFlags &= ~keepflag;
    } else {
      measure[k].dbFlags |= keepflag;
    }
    if (dataset->msklist[i]) measure[k].dbFlags |= maskflag;
  }
  return TRUE;
}

/* sort a coordinate pair (X,Y) and the associated index (S) */
void sort_StatDataSet (StatDataSet *dataset) {
  
# define SWAPFUNC(A,B){ double dtmp; int itmp; \
    dtmp = dataset->flxlist[A]; dataset->flxlist[A] = dataset->flxlist[B]; dataset->flxlist[B] = dtmp;     \
    dtmp = dataset->errlist[A]; dataset->errlist[A] = dataset->errlist[B]; dataset->errlist[B] = dtmp;     \
    dtmp = dataset->wgtlist[A]; dataset->wgtlist[A] = dataset->wgtlist[B]; dataset->wgtlist[B] = dtmp;     \
    itmp = dataset->ranking[A]; dataset->ranking[A] = dataset->ranking[B]; dataset->ranking[B] = itmp;     \
    itmp = dataset->measSeq[A]; dataset->measSeq[A] = dataset->measSeq[B]; dataset->measSeq[B] = itmp;     \
  }
# define COMPARE(A,B)(dataset->ranking[A] < dataset->ranking[B])
  
  OHANA_SORT (dataset->Nlist, COMPARE, SWAPFUNC);
  
# undef SWAPFUNC
# undef COMPARE

}

void sort_entry_by_offset (double *offset, int *entry, int N) {

# define SWAPFUNC(A,B){ double dtmp; int itmp; \
    dtmp = offset[A]; offset[A] = offset[B]; offset[B] = dtmp;     \
    itmp = entry[A];  entry[A]  = entry[B];  entry[B]  = itmp;     \
  }
# define COMPARE(A,B)(offset[A] < offset[B])
  
  OHANA_SORT (N, COMPARE, SWAPFUNC);
  
# undef SWAPFUNC
# undef COMPARE
}

void GetPhotFlagStats (uint32_t *photFlagUpper, uint32_t *photFlagLower, uint32_t *photflag_list, int Nphotflag) {

  int i, bit;

  uint32_t outputUpper = 0;
  uint32_t outputLower = 0;

  uint32_t flag = 0x1; // start with the flag on the 0 bit
  for (bit = 0; bit < 32; bit++) {

    int Nraised = 0;
    for (i = 0; i < Nphotflag; i++) {
      if (photflag_list[i] & flag) Nraised ++;
    }

    float fraction = Nraised / (float) Nphotflag;

    if (fraction < 0.25) {
      // do nothing (no bits raised)
    }
    if ((fraction >= 0.25) && (fraction < 0.50)) {
      outputLower |= flag;
    }
    if ((fraction >= 0.50) && (fraction < 0.75)) {
      outputUpper |= flag;
    }
    if (fraction >= 0.75) {
      outputLower |= flag;
      outputUpper |= flag;
    }
    flag = (flag << 1);
  }
  *photFlagUpper = outputUpper;
  *photFlagLower = outputLower;
}
