# include "relphot.h"

static int Nmax;

enum {THREAD_RUN, THREAD_DONE};

typedef struct {
  int entry;
  int state;
  Catalog *catalog;
  int Ncatalog;
  FlatCorrectionTable *flatcorr;
  SetMrelInfo summary;
} ThreadInfo;

// ** internal functions:
void *setMrel_worker (void *data);
int setMrel_threaded (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr);
int setMrel_catalog (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt);
int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure);

// we want to allocate the stats list,dlist arrays only once (or once per thread).
// this function finds the largest array so we can allocate that max size when needed
void initMrel (Catalog *catalog, int Ncatalog) {

  off_t i, j;
  
  Nmax = 0;
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {
      Nmax = MAX (Nmax, catalog[i].averageT[j].Nmeasure);
    }
  }
}  

float getMrel (Catalog *catalog, off_t meas, int cat) {

  int Nsec, Nsecfilt, photcode;
  int ave;
  float value;

  ave = catalog[cat].measureT[meas].averef;
  photcode = catalog[cat].measureT[meas].photcode;

  int ecode = GetPhotcodeEquivCodebyCode (photcode);
  Nsec = GetPhotcodeNsec(ecode);
  Nsecfilt = GetPhotcodeNsecfilt ();

  // is this star OK?
  if (catalog[cat].secfilt[Nsecfilt*ave+Nsec].flags & STAR_BAD) return (NAN);  
 
  value = catalog[cat].secfilt[Nsecfilt*ave+Nsec].M;
  return (value);
}

void SetMrelInfoInit (SetMrelInfo *results, int allocLists) {
  results->Nfew  = 0;
  results->Ncode  = 0;
  results->Nsys  = 0;
  results->Nbad  = 0;
  results->Ncal  = 0;
  results->Nmos  = 0;
  results->Ngrid = 0;
  if (allocLists) {
    ALLOCATE (results->Mpsflist, double, Nmax);
    ALLOCATE (results->dpsflist, double, Nmax);
    ALLOCATE (results->wpsflist, double, Nmax);
  }
}

void SetMrelInfoFree (SetMrelInfo *results) {
  free (results->Mpsflist);
  free (results->dpsflist);
  free (results->wpsflist);
}

void SetMrelInfoAccum (SetMrelInfo *summary, SetMrelInfo *results) {
  summary->Nfew  += results->Nfew ;
  summary->Ncode += results->Ncode ;
  summary->Nsys  += results->Nsys ;
  summary->Nbad  += results->Nbad ;
  summary->Ncal  += results->Ncal ;
  summary->Nmos  += results->Nmos ;
  summary->Ngrid += results->Ngrid;
}

// mutex to lock setMrel_worker operations 
static pthread_mutex_t setMrel_mutex = PTHREAD_MUTEX_INITIALIZER;
static int nextCatalog = 0;

// we have an array of catalogs (catalog, Ncatalog).  we need to hand out catalogs one at a time to
// the worker threads as they need
off_t getNextCatalogForThread (int Ncatalog) {

  pthread_mutex_lock (&setMrel_mutex);
  if (nextCatalog >= Ncatalog) {
    pthread_mutex_unlock (&setMrel_mutex);
    return (-1);
  }
  int thisCatalog = nextCatalog;
  nextCatalog ++;

  pthread_mutex_unlock (&setMrel_mutex);
  return (thisCatalog);
}

// setMrel and setMrelOutput are extremely similar, but have slightly different implications:
// * setMrel uses the internal Tiny structures only
// * setMrelOutput skips stars for which there are too few good measurements 
// * setMrelOutput is meant to be called repeatedly, relaxing the criteria for 'good' on each pass
// * setMrelOutput updates 2MASS average flags
// * setMrelOutput updates average EXT flags (PS1 and 2MASS)

int setMrel (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {

  int i;

  if (NTHREADS) {
    int status = setMrel_threaded (catalog, Ncatalog, flatcorr);
    return status;
  }

  int Nsecfilt = GetPhotcodeNsecfilt ();

  SetMrelInfo summary, results;
  SetMrelInfoInit (&summary, FALSE);
  SetMrelInfoInit (&results, TRUE);

  for (i = 0; i < Ncatalog; i++) {
    // pass == -1 for anything other than the final pass
    switch (SET_MREL_VERSION) {
      case 0:
	setMrel_catalog (catalog, i, -1, flatcorr, &results, Nsecfilt);
	break;
      case 1:
	setMrel_catalog_alt (catalog, i, -1, flatcorr, &results, Nsecfilt);
	break;
      default:
	fprintf (stderr, "invalid setMrel version %d (use 0 or 1)\n", SET_MREL_VERSION);
	exit (2);
    }
    SetMrelInfoAccum (&summary, &results);
  }
  if (VERBOSE2) fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys);

  SetMrelInfoFree (&results);

  return (TRUE);
}

int setMrelOutput (Catalog *catalog, int Ncatalog, int pass, FlatCorrectionTable *flatcorr) {

  int i;

  // fprintf (stderr, "this version of setMrel is invalid for now\n");
  // exit (1);
	
  int Nsecfilt = GetPhotcodeNsecfilt ();

  SetMrelInfo summary, results;
  SetMrelInfoInit (&summary, FALSE);
  SetMrelInfoInit (&results, TRUE); // allocates results->list,dlist,wlist

  ALLOCATE (results.Maplist,  double, Nmax);
  ALLOCATE (results.daplist,  double, Nmax);
  ALLOCATE (results.waplist,  double, Nmax);

  ALLOCATE (results.Mkronlist, double, Nmax);
  ALLOCATE (results.dkronlist, double, Nmax);
  ALLOCATE (results.wkronlist, double, Nmax);

  ALLOCATE (results.psfqf_list, double, Nmax);
  ALLOCATE (results.psfqfperf_list, double, Nmax);
  ALLOCATE (results.stargal_list, double, Nmax);

  for (i = 0; i < Ncatalog; i++) {
    switch (SET_MREL_VERSION) {
      // case 0:
      // setMrel_catalog (catalog, i, pass, flatcorr, &results, Nsecfilt);
      // break;
      case 1:
	setMrel_catalog_alt (catalog, i, pass, flatcorr, &results, Nsecfilt);
	break;
      default:
	fprintf (stderr, "invalid setMrel version %d (use 1)\n", SET_MREL_VERSION);
	exit (2);
    }
    SetMrelInfoAccum (&summary, &results);
  }
  if (VERBOSE2) fprintf (stderr, "%d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys);

  SetMrelInfoFree (&results);
  free (results.Maplist);
  free (results.daplist);
  free (results.waplist);

  free (results.Mkronlist);
  free (results.dkronlist);
  free (results.wkronlist);

  free (results.psfqf_list);
  free (results.psfqfperf_list);
  free (results.stargal_list);
  return (TRUE);
}

int setMrel_threaded (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {

  int i;

  SetMrelInfo summary;
  SetMrelInfoInit (&summary, FALSE);

  pthread_attr_t attr;
  pthread_attr_init (&attr);
  pthread_attr_setdetachstate (&attr, PTHREAD_CREATE_DETACHED);
  
  pthread_t *threads;
  ALLOCATE (threads, pthread_t, NTHREADS);

  ThreadInfo *threadinfo;
  ALLOCATE (threadinfo, ThreadInfo, NTHREADS);

  // each time this function is called, we cycle through the available catalogs.
  // make sure we start at 0
  nextCatalog = 0;;

  // launch N worker threads
  for (i = 0; i < NTHREADS; i++) {
    threadinfo[i].entry = i;
    threadinfo[i].state = THREAD_RUN;
    threadinfo[i].catalog  =  catalog;
    threadinfo[i].Ncatalog = Ncatalog;
    threadinfo[i].flatcorr = flatcorr;
    SetMrelInfoInit (&threadinfo[i].summary, FALSE);
    pthread_create (&threads[i], NULL, setMrel_worker, &threadinfo[i]);
  }
  pthread_attr_destroy (&attr);

  // wait until all threads have finished
  while (1) {
    int allDone = TRUE;
    for (i = 0; i < NTHREADS; i++) {
      if (threadinfo[i].state == THREAD_RUN) allDone = FALSE;
    }
    if (allDone) {
      break;
    }
    usleep (500000);
  }

  // all threads are done, free the threads array and grab the info
  free (threads);
  
  // report stats & summary from the threads
  for (i = 0; i < NTHREADS; i++) {
    if (VERBOSE2) fprintf (stderr, "setMrel thread %d : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", 
	     i, 
	     threadinfo[i].summary.Ncode, 
	     threadinfo[i].summary.Nfew, 
	     threadinfo[i].summary.Nbad, 
	     threadinfo[i].summary.Ncal, 
	     threadinfo[i].summary.Nmos, 
	     threadinfo[i].summary.Ngrid, 
	     threadinfo[i].summary.Nsys);
    SetMrelInfoAccum (&summary, &threadinfo[i].summary);
  }
  if (VERBOSE2) fprintf (stderr, "total : %d stars with no data in photcode, %d stars marked having too few measurements (Nbad: %d, Ncal: %d, Nmos: %d, Ngrid: %d, Nsys: %d)\n", summary.Ncode, summary.Nfew, summary.Nbad, summary.Ncal, summary.Nmos, summary.Ngrid, summary.Nsys);
  free (threadinfo);

  return TRUE;
}

void *setMrel_worker (void *data) {

  ThreadInfo *threadinfo = data;

  int Nsecfilt = GetPhotcodeNsecfilt ();

  SetMrelInfo results;
  SetMrelInfoInit (&results, TRUE); // allocate list, dlist arrays here

  while (1) {

    off_t i = getNextCatalogForThread(threadinfo->Ncatalog);
    if (i == -1) {
      threadinfo->state = THREAD_DONE;
      return NULL;
    }

    Catalog *catalog = threadinfo->catalog;
    FlatCorrectionTable *flatcorr = threadinfo->flatcorr;

    // pass == -1 for anything other than the final pass
    switch (SET_MREL_VERSION) {
      case 0:
	setMrel_catalog (catalog, i, -1, flatcorr, &results, Nsecfilt);
	break;
      case 1:
	setMrel_catalog_alt (catalog, i, -1, flatcorr, &results, Nsecfilt);
	break;
      default:
	fprintf (stderr, "invalid setMrel version %d (use 0 or 1)\n", SET_MREL_VERSION);
	exit (2);
    }
    SetMrelInfoAccum (&threadinfo->summary, &results);
  }

  SetMrelInfoFree (&results);
  return NULL;
}

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


// setMrel_catalog is used in 3 different contexts:

// * during the main relphot iterations:
// ** operations only apply to measureT / averageT
// ** we are applying the analysis to the bright subset catalog
// ** we skip any stars found to be bad (STAR_BAD)

// * during the final pass
// ** operations are applied to all objects
// ** we have special tests for PS1 extended data, 2MASS extended & good data, & synthetic photometry

// * during the stand-alone '-averages' mode (relphot_objects)
// ** no image data is loaded
// ** we only use the provided calibration (measure.Mcal)
// ** we skip outlier rejection of measurements

// set the Mrel values for the specified catalog
// NOTE: here 'catalog' is a pointer to a specific catalog, not the root of the array
int setMrel_catalog (Catalog *catalog, int Nc, int pass, FlatCorrectionTable *flatcorr, SetMrelInfo *results, int Nsecfilt) {

  fprintf (stderr, "??? should you still be using the old version of setMrel_catalog??\n");
  exit (3);

# if (0)

  off_t j, k, m, ID;
  int N;
  float Msys, Mcal= 0, Mmos = 0, Mgrid = 0;

  StatType stats, apstats, kronstats;
  liststats_setmode (&psfstats, STATMODE);
  liststats_setmode (&apstats, STATMODE);
  liststats_setmode (&kronstats, STATMODE);

  double *list      = results->list;
  double *dlist     = results->dlist;
  double *wlist     = results->wlist;
  double *aplist    = results->aplist;
  double *kronlist  = results->kronlist;
  double *dkronlist = results->dkronlist;

  double *psfqf_list     = results->psfqf_list;
  double *psfqfperf_list = results->psfqfperf_list;
  double *stargal_list   = results->stargal_list;

  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++) {
    // option for a test print
    if (FALSE && (catalog[Nc].average[j].objID == 0x7146) && (catalog[Nc].average[j].catID == 0x49d8)) {
      fprintf (stderr, "test obj\n");
      print_measure_set (&catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt], catalog[Nc].measure);
    }

    if (isSetMrelFinal) {
      // set the name of the primary skycell (this is used in a strcmp to match the skycells in stack detections)
      // XXX : this whole function is deprecated
      // BoundaryTreePrimaryCell(primaryCell, catalog[Nc].average[j].R, catalog[Nc].average[j].D);
    }

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

    int NextPS1 = 0;
    int NpsfPS1 = 0;

    int Nmeas = 0; // number of meas used for psfqf, psfqfperf, and stargal

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

      /* star/photcodes already calibrated */
      if ( isSetMrelFinal && catalog[Nc].found_t[Nsecfilt*j+Nsec]) continue;  
      
      // skip bad stars
      if (!isSetMrelFinal && (catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD)) continue;

      int Ncode = 0;
      int Next = 0;
      int haveSynth = FALSE;
      int haveStack = FALSE;
      int havePS1   = 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 forceSynth = FALSE;
      int forceSynthEntry = -1;

      int haveUbercal = FALSE;

      int minUbercalDist = 1000;
    
      int isBad = FALSE;
      int isSuspect = FALSE;

      N = 0;
      m = catalog[Nc].averageT[j].measureOffset;
      for (k = 0; k < catalog[Nc].averageT[j].Nmeasure; k++, m++) {

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

	catalog[Nc].measure[m].dbFlags &= ~ID_MEAS_STACK_PRIMARY;

	if (catalog[Nc].measureT[m].dbFlags & MEAS_BAD) SKIP_THIS_MEAS(Nbad); 

	if (getImageEntry (m, Nc) < 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 = catalog[Nc].measureT[m].Mcal; // check that this is zero for loaded REF value
	} else {
	  Mcal  = getMcal  (m, Nc, flatcorr, catalog);
	  if (isnan(Mcal))  SKIP_THIS_MEAS(Ncal);
	  Mmos  = getMmos  (m, Nc);
	  if (isnan(Mmos))  SKIP_THIS_MEAS(Nmos);
	  Mgrid = getMgrid (m, Nc);
	  if (isnan(Mgrid)) SKIP_THIS_MEAS(Ngrid);
	}

	// skip some absurd values NAN, < 0.0, > 30.0
	Msys = PhotSysTiny (&catalog[Nc].measureT[m], &catalog[Nc].averageT[j], &catalog[Nc].secfilt[j*Nsecfilt], 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(m,Nc);
	minUbercalDist = MIN(minUbercalDist, myUbercalDist);

	if (isSetMrelFinal) {
	  // special options for PS1 data
	  if ((catalog[Nc].measure[m].photcode >= 10000) && (catalog[Nc].measure[m].photcode <= 10500)) {
	    // count the extended detections
	    if (!isnan(catalog[Nc].measure[m].Map)) {
	      float dMagAp = catalog[Nc].measure[m].M - catalog[Nc].measure[m].Map;
	      float SigmaAp = hypot(0.1, 2.5*catalog[Nc].measure[m].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;
	  }

	  // gpc1 stack data
	  if ((catalog[Nc].measure[m].photcode >= 11000) && (catalog[Nc].measure[m].photcode <= 11400)) {
	    haveStack = TRUE;
	    haveStackObject = TRUE;

	    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 (MatchImageName (m, Nc, primaryCell)) {
	      float stackPrimaryOffset = getCenterOffset (m, Nc, &catalog[Nc].measure[m], &stackImageID);
	      if (stackPrimaryOffset < stackPrimaryOffsetMin) {
		stackPrimaryOffsetMin = stackPrimaryOffset;
		// stackPrimaryIDmin = stackImageID;
		stackPrimaryMeasureMin = m;
	      }
	    }

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

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

	    isSuspect |= (catalog[Nc].measure[m].photFlags & code->photomPoorMask);
	    isSuspect |= (catalog[Nc].measure[m].psfQFperf < 0.85);

	    // mean psfQF, psfQFperf, and stargal values are calculated per object (all filters, stacks and exposures)
	    psfqf_list[Nmeas]     = catalog[Nc].measure[m].psfQF;
	    psfqfperf_list[Nmeas] = catalog[Nc].measure[m].psfQFperf;
	    stargal_list[Nmeas]   = catalog[Nc].measure[m].extNsigma;
	    Nmeas++;

	    continue; // do not include stack measurements in the per-exposure statistics
	  }

	  // count extended detections for 2MASS (XXX NOTE hardwired photcodes 2011, 2012, 2013)
	  if ((catalog[Nc].measure[m].photcode >= 2011) && (catalog[Nc].measure[m].photcode <= 2013)) {
	    if (catalog[Nc].measure[m].photFlags & 0x00c00000) {
	      Next ++;
	      Galaxy2MASS = TRUE;
	    }
	    if (pass == 0) {
	      if (catalog[Nc].measure[m].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 ((catalog[Nc].measure[m].photcode >= 3001) && (catalog[Nc].measure[m].photcode <= 3005)) {
	    // something of a hack: force object to use synth values if synth mags >>
	    // saturation (3pi instrumental mags < -15)
	    float MaxMagForceSynth = NAN;
	    switch (catalog[Nc].measure[m].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 (catalog[Nc].measureT[m].M < MaxMagForceSynth) {
	      forceSynth = TRUE;
	      forceSynthEntry = N;
	    } else {
	      if (pass < 4) {
		SKIP_THIS_MEAS(Nbad);
	      }
	      haveSynth = TRUE;
	    }
	  }

	  // Map (aplist) and Mkron (kronlist,dkronlist) are used to calculate mean mags per filter
	  float Map = PhotAper (&catalog[Nc].measure[m]);
	  aplist[N] = Map - Mcal - Mmos - Mgrid;

	  float Mkron = PhotKron (&catalog[Nc].measure[m]);
	  kronlist[N] = Mkron - Mcal - Mmos - Mgrid;
	  dkronlist[N] = catalog[Nc].measure[m].dMkron;

	  // mean psfQF, psfQFperf, and stargal values are calculated per object (all filters, stacks and exposures)
	  psfqf_list[Nmeas]     = catalog[Nc].measure[m].psfQF;
	  psfqfperf_list[Nmeas] = catalog[Nc].measure[m].psfQFperf;
	  stargal_list[Nmeas]   = catalog[Nc].measure[m].extNsigma;
	  Nmeas++;
	}

	// 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
	dlist[N] = MAX (hypot(catalog[Nc].measureT[m].dM, code->photomErrSys), MIN_ERROR);
	wlist[N] = 1.0;
	list[N]  = Msys - Mcal - Mmos - Mgrid;

	// up-weight the ubercal values (or convergence can take a long time...)
	if (catalog[Nc].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
	  haveUbercal = TRUE;
	  wlist[N] = 10.0;
	}

	// 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
	if (refPhotcode) {
	  if (code->code == refPhotcode->code) {
	    wlist[N] = 100.0;
	  }
	}
	N++;
      }

      int Nminmeas = isSetMrelFinal ? 1 : STAR_TOOFEW + 1;

      if (haveStack && (N > 1)) {
	// fprintf (stderr, "multiple stack values for %10.6f %10.6f\n", catalog[Nc].averageT[j].R, catalog[Nc].averageT[j].D);
      }

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

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

      if (forceSynth) {
	// use the single SYNTH value instead of the other mags here
	myAssert ((forceSynthEntry < N) && (forceSynthEntry >= 0), "programming error");
	list[0]  = list[forceSynthEntry];
	dlist[0] = dlist[forceSynthEntry];
	wlist[0] = wlist[forceSynthEntry];
	N = 1;
      }
      liststats (list, dlist, wlist, N, &stats);

      catalog[Nc].secfilt[Nsecfilt*j+Nsec].M      = stats.mean;
      catalog[Nc].secfilt[Nsecfilt*j+Nsec].dM     = stats.error;
      catalog[Nc].secfilt[Nsecfilt*j+Nsec].Mchisq = (stats.Nmeas > 1) ? stats.chisq : NAN;

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

      if (isSetMrelFinal) {
	catalog[Nc].found_t[Nsecfilt*j+Nsec] = TRUE;

	catalog[Nc].secfilt[Nsecfilt*j+Nsec].Mstdev = 1000.0*stats.sigma; // Mstdev is in millimags (not enough space for more precision)
	catalog[Nc].secfilt[Nsecfilt*j+Nsec].Ncode = Ncode;
	catalog[Nc].secfilt[Nsecfilt*j+Nsec].Nused = stats.Nmeas;

	catalog[Nc].secfilt[Nsecfilt*j+Nsec].M_80 = 1000 * stats.Upper80;
	catalog[Nc].secfilt[Nsecfilt*j+Nsec].M_20 = 1000 * stats.Lower20;

	// NOTE : use the modified weight for apmags as well as psf mags
	liststats (aplist, dlist, wlist, N, &apstats);
	catalog[Nc].secfilt[Nsecfilt*j+Nsec].Map  = apstats.mean; 

	liststats (kronlist, dkronlist, wlist, N, &kronstats);
	catalog[Nc].secfilt[Nsecfilt*j+Nsec].Mkron  = kronstats.mean; 
	catalog[Nc].secfilt[Nsecfilt*j+Nsec].dMkron = kronstats.error; 

	if (haveStack) {
	  m  = (stackPrimaryMeasureMin >= 0) ? stackPrimaryMeasureMin : stackCenterMeasureMin;

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

	  // get the zero point for the selected image
	  float zp = PhotZeroPoint (&catalog[Nc].measure[m], &catalog[Nc].average[j], &catalog[Nc].secfilt[j*Nsecfilt]) - (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?)
	  catalog[Nc].secfilt[Nsecfilt*j+Nsec].FluxPSF   = zpFactor * catalog[Nc].measure[m].FluxPSF;  
	  catalog[Nc].secfilt[Nsecfilt*j+Nsec].dFluxPSF  = zpFactor * catalog[Nc].measure[m].dFluxPSF; 
	  catalog[Nc].secfilt[Nsecfilt*j+Nsec].FluxKron  = zpFactor * catalog[Nc].measure[m].FluxKron; 
	  catalog[Nc].secfilt[Nsecfilt*j+Nsec].dFluxKron = zpFactor * catalog[Nc].measure[m].dFluxKron;

	  catalog[Nc].secfilt[Nsecfilt*j+Nsec].stackDetectID = ID;

	  catalog[Nc].measure[m].dbFlags |= ID_MEAS_STACK_PHOT_SRC;
	  if (stackPrimaryMeasureMin >= 0) {
	    catalog[Nc].measure[m].dbFlags |= ID_MEAS_STACK_PRIMARY;
	  }
	}

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

	switch (pass) {
	  case 0:
	    catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_0;
	    if (havePS1) GoodPS1 = TRUE;
	    break;
	  case 1:
	    catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_1;
	    if (havePS1) GoodPS1 = TRUE;
	    break;
	  case 2:
	    catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_2;
	    if (havePS1) GoodPS1 = TRUE;
	    break;
	  case 3:
	    catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_3;
	    break;
	  case 4:
	    catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_PHOTOM_PASS_4;
	    break;
	}
	if (haveSynth) {
	  catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_USE_SYNTH;
	}	
	if (haveStack) {
	  catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_HAS_STACK;
	}	
	if (havePS1) {
	  catalog[Nc].secfilt[Nsecfilt*j+Nsec].flags |= ID_SECF_HAS_PS1;
	}	
	if (haveUbercal) {
	  catalog[Nc].secfilt[Nsecfilt*j+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:
	catalog[Nc].average[j].flags &= ~flagBits;
      }

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

      if (NstackGood >= 2) {
	catalog[Nc].average[j].flags |= ID_OBJ_GOOD_STACK;
      } else if (NstackGood + NstackSuspect >= 2) {
	catalog[Nc].average[j].flags |= ID_OBJ_SUSPECT_STACK;
      } else if (haveStackObject) {
	catalog[Nc].average[j].flags |= ID_OBJ_BAD_STACK;
      }

      // use NAN values for catalog[Nc].average[j].psfQF as signal
      if (!isfinite(catalog[Nc].average[j].psfQF)) {
	// do something here
	float psfQFmax     = 0.0;
	float psfQFperfmax = 0.0;
	float stargalmax   = 0.0;
	int seq;
	for (seq = 0; seq < Nmeas; seq++) {
	  if (isfinite(psfqf_list[seq])) {
	    psfQFmax     = MAX (psfqf_list[seq], psfQFmax);
	  }
	  if (isfinite(psfqfperf_list[seq])) {
	    psfQFperfmax = MAX (psfqfperf_list[seq], psfQFperfmax);
	  }
	  if (isfinite(stargal_list[seq])) {
	    stargalmax   = MAX (stargal_list[seq], stargalmax);
	  }
	}

	catalog[Nc].average[j].psfQF     = psfQFmax;
	catalog[Nc].average[j].psfQFperf = psfQFperfmax;
	catalog[Nc].average[j].stargal   = stargalmax;
      }
    }
  }
  if (primaryCell) free (primaryCell);
  return (TRUE);
# endif
}

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

  off_t k;

  int Nsecfilt = GetPhotcodeNsecfilt ();

  off_t m = average[0].measureOffset;

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

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

/* set measure.Mcal for all measures except ID_MEAS_NOCAL and ID_IMAGE_PHOTOM_NOCAL */
int setMcalOutput (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {

  int i;
  off_t j, k, m;
  float Mcal, Mmos, Mgrid;

  MEAS_BAD = ID_MEAS_NOCAL;
  IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL;

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

      m = catalog[i].averageT[j].measureOffset;
      for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
	if (catalog[i].measureT[m].dbFlags & MEAS_BAD) continue;
	Mcal  = getMcal  (m, i, flatcorr, catalog);
	if (isnan(Mcal)) continue;
	Mmos  = getMmos  (m, i);
	if (isnan(Mmos)) continue;
	Mgrid = getMgrid (m, i);
	if (isnan(Mgrid)) continue;

	// note that measurements for which the image is not selected will not be modified
	// (that is a good thing! -- we keep a prior calibration)

	// set the output calibration
	catalog[i].measure[m].Mcal = Mcal + Mmos + Mgrid;

	if (catalog[i].measureT[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
	  myAssert (isfinite(catalog[i].measure[m].Mcal), "oops, broke an ubercal mag");
	}
      }
    }
  }
  return (TRUE);
}

int markObjects (Catalog *catalog, int Ncatalog) {

  int i, n;
  off_t j;

  // How strongly do I own this object?
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {
      int nOwn = 0;
      int m = catalog[i].averageT[j].measureOffset;
      for (n = 0; n < catalog[i].averageT[j].Nmeasure; n++) {
	if (!catalog[i].measureT[m+n].myDet) continue;
	nOwn ++;
      }
      catalog[i].averageT[j].nOwn = nOwn;
    }
  }
  return TRUE;
}

int dumpObjects (char *filename, Catalog *catalog, int Ncatalog) {

  int i, n;
  off_t j;

  FILE *ftest = fopen (filename, "w");

  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {
      int m = catalog[i].averageT[j].measureOffset;
      for (n = 0; n < catalog[i].averageT[j].Nmeasure; n++) {
	fprintf (ftest, "%08x %08x %10.6f %10.6f  %3d %1d\n", catalog[i].averageT[j].catID, catalog[i].averageT[j].objID, catalog[i].averageT[j].R, catalog[i].averageT[j].D, catalog[i].averageT[j].Nmeasure, catalog[i].measureT[m+n].myDet); 
      }
    }
  }
  fclose (ftest);
  return TRUE;
}

void clean_stars (Catalog *catalog, int Ncatalog) {

  int i, j, Ndel, Nave, Ntot, mark, Ns, Nscat, Nchi, Nnan;
  float dM;
  double MaxScatter, MaxChisq;
  double *xlist, *slist, *dlist;

  StatType stats;
  liststats_setmode (&stats, "MEAN");

  if (VERBOSE) fprintf (stderr, "marking poor stars\n");

  /* find Xm median -> ChiSq lim must be > median */
  for (i = Ntot = 0; i < Ncatalog; i++) {
    Ntot += catalog[i].Naverage; 
  }
  ALLOCATE (xlist, double, Ntot);
  ALLOCATE (slist, double, Ntot);
  ALLOCATE (dlist, double, Ntot);

  int Nsecfilt = GetPhotcodeNsecfilt ();

  // eliminate bad stars using the stats for a single secfilt at a time
  for (Ns = 0; Ns < Nphotcodes; Ns ++) {
    
    int thisCode = photcodes[Ns][0].code;
    int Nsec = GetPhotcodeNsec(thisCode);

    for (i = Ntot = 0; i < Ncatalog; i++) {
      for (j = 0; j < catalog[i].Naverage; j++) {
	if ( catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD ) continue;
	float Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq;
	if (isnan(Mchisq)) continue;
	xlist[Ntot] = Mchisq;
	slist[Ntot] = catalog[i].secfilt[Nsecfilt*j+Nsec].dM;
	dlist[Ntot] = 1;
	Ntot ++;
      }
    }
  
    liststats (xlist, dlist, NULL, Ntot, &stats);
    MaxChisq = MAX (STAR_CHISQ, 2*stats.median);

    liststats (slist, dlist, NULL, Ntot, &stats);
    MaxScatter = MAX (STAR_SCATTER, 2*stats.median);
    fprintf (stderr, "Max Scatter: %f, Max Chisq: %f\n", MaxScatter, MaxChisq);

    Ndel = Nave = Nscat = Nnan = Nchi = 0;
    for (i = 0; i < Ncatalog; i++) {
      for (j = 0; j < catalog[i].Naverage; j++) {
	dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dM;
	float Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq;
	mark = (dM > MaxScatter) || (isnan(Mchisq)) || (Mchisq > MaxChisq);
	if (mark) {
	  catalog[i].secfilt[Nsecfilt*j+Nsec].flags |= ID_STAR_POOR;
	  Ndel ++;
	  if (dM > MaxScatter)   { Nscat ++; }
	  if (isnan(Mchisq))     { Nnan ++; }
	  if (Mchisq > MaxChisq) { Nchi ++; }
	} else {
	  catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_STAR_POOR;
	}
	Nave ++;
      }
    }
    fprintf (stderr, "%d stars marked variable (%d scat, %d nan, %d chi), %d total\n", Ndel, Nscat, Nnan, Nchi, Nave);
  }
  free (xlist);
  free (slist);
  free (dlist);
}

// clean_measures examines the stats for a single star.  It measures the INNER 50% mean
// and sigma, it then re-measures the mean and sigma using all stars with NSIGMA_CLIP (3)
// sigma of the INNER 50% mean.  it then flags any measurements which are more than
// NSIGMA_REJECT (5) sigma of the mean

# define NSIGMA_CLIP 3.0
# define NSIGMA_REJECT 5.0
void clean_measures (Catalog *catalog, int Ncatalog, int final, FlatCorrectionTable *flatcorr) {

  off_t j, k, m, Nmax, Ndel, Nave;
  int i, N, image_bad, TOOFEW;
  off_t *ilist;
  double *tlist, *list, *dlist;
  float Msys, Mcal, Mmos, Mgrid;
  int Ncal, Nmos, Ngrid, Nfew;

  int Nsecfilt = GetPhotcodeNsecfilt ();

  if (VERBOSE) fprintf (stderr, "marking poor measures\n");
  /* Nmeasure is now different, need to reallocate */
  Nmax = 0;
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {
      Nmax = MAX (Nmax, catalog[i].averageT[j].Nmeasure);
    }
  }
  ALLOCATE (list, double, MAX (1, Nmax));
  ALLOCATE (dlist, double, MAX (1, Nmax));
  ALLOCATE (ilist, off_t, Nmax);
  ALLOCATE (tlist, double, Nmax);
  
  /* it makes no sense to mark 3-sigma outliers with <5 measurements */
  TOOFEW = MAX (5, STAR_TOOFEW);

  // stats structures for inner and full stats
  StatType instats, stats;
  liststats_setmode (&instats, "INNER_MEAN");
  liststats_setmode (&stats, "MEAN");

  Ndel = Nave = 0;
  Ncal = Nmos = Ngrid = Nfew = 0;
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {

      int Ns;
      for (Ns = 0; Ns < Nphotcodes; Ns++) {
	
	int thisCode = photcodes[Ns][0].code;
	int Nsec = GetPhotcodeNsec(thisCode);
	
	/* on final processing, skip stars already measured */
	if (final && catalog[i].found_t[Nsecfilt*j + Nsec]) continue;  

	/* skip bad stars to prevent them from becoming good (on inner sample) */
	if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  

	/* accumulate list of valid measurements */
	m = catalog[i].averageT[j].measureOffset;
	N = 0;
	for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
	  // skip measurements that do not match the current photcode
	  int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
	  if (ecode != thisCode) { continue; }

	  // SKIP gpc1 stack data
	  if (isGPC1stack(catalog[i].measureT[m].photcode)) continue;
	  
	  // SKIP gpc1 forced-warp data
	  if (isGPC1warp(catalog[i].measureT[m].photcode)) continue;
	  
	  // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers
	  Mcal  = getMcal  (m, i, flatcorr, catalog);
	  if (isnan(Mcal)) { Ncal ++; continue; }
	  Mmos  = getMmos  (m, i);
	  if (isnan(Mmos)) { Nmos ++; continue; }
	  Mgrid = getMgrid (m, i);
	  if (isnan(Mgrid)) { Ngrid ++; continue; }

	  Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF);
	  list[N] = Msys - Mcal - Mmos - Mgrid;
	  dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR);
	  N++;
	}
	if (N <= TOOFEW) { 
	  Nfew ++; 
	  continue; 
	}

	/* 3-sigma clip based on stats of inner 50% */

	// calculated mean of inner 50%
	liststats (list, dlist, NULL, N, &instats);
	instats.sigma = MAX (MIN_ERROR, instats.sigma); /* if measurements agree too well, sigma -> 0.0 */

	// ignore entries > 3sigma from inner mean
	for (k = m = 0; k < N; k++) {
	  if (fabs (list[k] - instats.median) < NSIGMA_CLIP*instats.sigma) {
	    list[m] = list[k];
	    m++;
	  }
	}
	// recalculate the mean & sigma of the accepted measurements
	liststats (list, dlist, NULL, m, &stats);
	stats.sigma = MAX (MIN_ERROR, stats.sigma);

	/* apply to list of all relevant measurements, including IMAGE_POOR & IMAGE_FEW */
	image_bad = IMAGE_BAD;
	IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL;
	m = catalog[i].averageT[j].measureOffset;
	N = 0;
	for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
	  // skip measurements that do not match the current photcode
	  int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
	  if (ecode != thisCode) { continue; }

	  // SKIP gpc1 stack data
	  if (isGPC1stack(catalog[i].measureT[m].photcode)) continue;
	  
	  // SKIP gpc1 forced-warp data
	  if (isGPC1warp(catalog[i].measureT[m].photcode)) continue;
	  
	  // NOTE: we do not skip MEAS_BAD because this measurement is just an internal assessment of the outliers
	  Mcal  = getMcal  (m, i, flatcorr, catalog);
	  if (isnan(Mcal)) continue;
	  Mmos  = getMmos  (m, i);
	  if (isnan(Mmos)) continue;
	  Mgrid = getMgrid (m, i);
	  if (isnan(Mgrid)) continue;

	  Msys = PhotSysTiny (&catalog[i].measureT[m], &catalog[i].averageT[j], &catalog[i].secfilt[j*Nsecfilt], MAG_CLASS_PSF);
	  list[N] = Msys - Mcal - Mmos - Mgrid;
	  dlist[N] = MAX (catalog[i].measureT[m].dM, MIN_ERROR);
	  ilist[N] = m;
	  N++;
	  Nave ++;
	}
	if (N < TOOFEW) continue;

	/* mark bad measures (> 3 sigma deviant) */
	for (k = 0; k < N; k++) {
	  // treat the scatter of the star as a systematic term.  this is a bit of an
	  // over-estimage (a perfect Gauss distribution with perfect errors would have
	  // mySigma = sqrt(2) too large)
	  float mySigma = hypot (stats.sigma, dlist[k]);
	  if (fabs (list[k] - stats.median) > NSIGMA_REJECT*mySigma) {
	    catalog[i].measureT[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM;
	    if (final) {
	      // for the final pass, we have a duplicate set of values in measure and measureT
	      catalog[i].measure[ilist[k]].dbFlags |= ID_MEAS_POOR_PHOTOM;
	    }
	    Ndel ++;
	  }
	}
	IMAGE_BAD = image_bad;
      }
    }
  }
  if (VERBOSE) fprintf (stderr, OFF_T_FMT" measures marked poor, "OFF_T_FMT" total\n", Ndel, Nave);

  free (list);
  free (dlist);
  free (ilist);
  free (tlist);
}

StatType statsStarN (Catalog *catalog, int Ncatalog, int Nsec, int seccode, FlatCorrectionTable *flatcorr) {

  off_t j, k, m, Ntot;
  int i, n, N;
  double *list, *dlist;
  float Mcal, Mmos, Mgrid;
  StatType stats;
  // int N1, N2, N3, N4, N0;
  // N1 = N2 = N3 = N4 = N0 = 0;

  int Nsecfilt = GetPhotcodeNsecfilt ();

  Ntot = 0;
  for (i = 0; i < Ncatalog; i++) {
    Ntot += catalog[i].Naverage; 
  }

  ALLOCATE (list, double, Ntot);
  ALLOCATE (dlist, double, Ntot);

  n = 0;
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {

      /* calculate the average value for a single star */
      if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) { continue;  }
      m = catalog[i].averageT[j].measureOffset;

      N = 0;
      for (k = 0; k < catalog[i].averageT[j].Nmeasure; k++, m++) {
	int ecode = GetPhotcodeEquivCodebyCode (catalog[i].measureT[m].photcode);
	if (ecode != seccode) { continue;}
	Mcal = getMcal  (m, i, flatcorr, catalog);
	if (isnan(Mcal)) { continue;}
	Mmos = getMmos  (m, i);
	if (isnan(Mmos)) { continue; }
	Mgrid = getMgrid (m, i);
	if (isnan(Mgrid)) { continue;}
	N++;
      }
      
      list[n] = N;
      dlist[n] = 1;
      n++;
    }
  }

  liststats_setmode (&stats, "MEAN");

  liststats (list, dlist, NULL, n, &stats);
  free (list);
  free (dlist);
  return (stats);
}

// stats for a single secfilt at a time
StatType statsStarX (Catalog *catalog, int Ncatalog, int Nsec) {

  off_t j, Ntot;
  int i, n;
  double *list, *dlist;
  StatType stats;

  int Nsecfilt = GetPhotcodeNsecfilt ();

  Ntot = 0;
  for (i = 0; i < Ncatalog; i++) {
    Ntot += catalog[i].Naverage; 
  }

  ALLOCATE (list, double, Ntot);
  ALLOCATE (dlist, double, Ntot);

  n = 0;
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {

      /* calculate the average value for a single star */
      if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  
	
      float Mchisq = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq;
      if (isnan(Mchisq)) continue;
      list[n] = Mchisq;
      dlist[n] = 1;
      n++;
    }
  }

  liststats_setmode (&stats, "MEAN");

  liststats (list, dlist, NULL, n, &stats);
  free (list);
  free (dlist);
  return (stats);
}

StatType statsStarS (Catalog *catalog, int Ncatalog, int Nsec) {

  int i, n;
  off_t j, Ntot;
  double *list, *dlist;
  float dM;
  StatType stats;

  int Nsecfilt = GetPhotcodeNsecfilt ();

  Ntot = 0;
  for (i = 0; i < Ncatalog; i++) {
    Ntot += catalog[i].Naverage; 
  }

  ALLOCATE (list, double, Ntot);
  ALLOCATE (dlist, double, Ntot);

  n = 0;
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {

      /* calculate the average value for a single star */
      if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  

      dM = catalog[i].secfilt[Nsecfilt*j+Nsec].dM;
      if (isnan(dM)) continue;
      list[n] = dM;
      dlist[n] = 1;
      n++;
    }
  }

  liststats_setmode (&stats, "MEAN");

  liststats (list, dlist, NULL, n, &stats);
  free (list);
  free (dlist);
  return (stats);
}

void plot_stars (Catalog *catalog, int Ncatalog) {

  int i, bin;
  off_t j;
  float dMrel;
  double *xlist, *Mlist;
  Graphdata graphdata;

# define NBIN 200
  ALLOCATE (xlist, double, NBIN);
  ALLOCATE (Mlist, double, NBIN);

  int Nsecfilt = GetPhotcodeNsecfilt ();

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

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

    for (i = 0; i < NBIN; i++) xlist[i] = 0.00025*i;
    bzero (Mlist, NBIN*sizeof(double));
    for (i = 0; i < Ncatalog; i++) {
      for (j = 0; j < catalog[i].Naverage; j++) {
	if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;  
	dMrel = catalog[i].secfilt[Nsecfilt*j+Nsec].dM;
	bin = dMrel / 0.00025;
	bin = MAX (0, MIN (NBIN-1, bin));
	Mlist[bin] += 1.0;
      }
    }

    plot_defaults (&graphdata);
    graphdata.style = 1;
    plot_list (&graphdata, xlist, Mlist, NBIN, "dMrel hist", "%s.dMhist.png", OUTROOT);
  }
  free (xlist);
  free (Mlist);
}

void plot_chisq (Catalog *catalog, int Ncatalog) {

  off_t j, Ntotal;
  int i, N, value;
  double *xlist, *ylist;
  Graphdata graphdata;

  int Nsecfilt = GetPhotcodeNsecfilt ();

  Ntotal = 0;
  for (i = 0; i < Ncatalog; i++) Ntotal += catalog[i].Naverage;

  ALLOCATE (xlist, double, Ntotal);
  ALLOCATE (ylist, double, Ntotal);

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

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

    N = 0;
    for (i = 0; i < Ncatalog; i++) {
      for (j = 0; j < catalog[i].Naverage; j++) {
	if (catalog[i].secfilt[Nsecfilt*j+Nsec].flags & STAR_BAD) continue;
	xlist[N] = catalog[i].secfilt[Nsecfilt*j+Nsec].M;
	value    = catalog[i].secfilt[Nsecfilt*j+Nsec].Mchisq;
	if (isnan(value)) continue;
	ylist[N] = value;
	N++;
      }
    }

    plot_defaults (&graphdata);
    graphdata.ymin = -3.0;
    plot_list (&graphdata, xlist, ylist, N, "chisq", "%s.chisq.png", OUTROOT);
  }

  free (xlist);
  free (ylist);
}

void plot_star_coords (Catalog *catalog, int Ncatalog) {

  int i;
  off_t j, N;
  double *xlist, *ylist;
  // double Xmin, Ymin, Xmax, Ymax;
  Graphdata graphdata;

  N = 0; 
  for (i = 0; i < Ncatalog; i++) {
    N += catalog[i].Naverage;
  }
  ALLOCATE (xlist, double, N);
  ALLOCATE (ylist, double, N);

  N = 0;
  // Xmin = Ymin = +360.0;
  // Xmax = Ymax = -360.0;
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {
      xlist[N] = catalog[i].averageT[j].R;
      ylist[N] = catalog[i].averageT[j].D;
      N++;
    }
  }
  plot_defaults (&graphdata);
  plot_list (&graphdata, xlist, ylist, N, "coords", "%s.coords.png", OUTROOT);

  free (xlist);
  free (ylist);
}

