# include "relphot.h"

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 StatDataSet arrays only once (or once per thread).  this
// function finds the largest value of Nmeasure so we can allocate that max size when
// needed
static int Nmax; 
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 SetMrelInfoReset (SetMrelInfo *results) {
  results->Nfew  = 0;
  results->Ncode  = 0;
  results->Nsys  = 0;
  results->Nbad  = 0;
  results->Ncal  = 0;
  results->Nmos  = 0;
  results->Ngrid = 0;
}

static int NsynthNames = 6;
static float synthMaxMags[] = {13.5, 13.5, 13.5, 13.0, 12.0, 13.5};
static char *synthNames[]   = { "g",  "r",  "i",  "z",  "y", "w"};

void SetMrelInfoResetObject (SetMrelInfo *results) {

  int i;
  for (i = 0; i < results->Nsecfilt; i++) {
    results->havePS1[i] = FALSE;
    results->haveSYN[i] = FALSE;
    results->needSYN[i] = FALSE;
    results->measSYN[i] = -1;

    results->psfQfMax[i]  = NAN;
    results->psfQfPerfMax[i]  = NAN;

    results->Nmeas[i] = 0;
    results->NmeasGood[i] = 0;
    results->Next[i] = 0;
    results->NexpPS1[i] = 0;
    results->havePS1[i] = 0;
    results->haveUbercal[i] = 0;

    results->Nstargal   = 0;
    results->Nphotflags = 0;

    results->tessID[i]    = 0;
    results->projID[i]    = 0;
    results->skycellID[i] = 0;

    results->minUbercalDist[i] = 1000;

    results->psfData[i].Nlist = 0;
    results->aperData[i].Nlist = 0;
    results->kronData[i].Nlist = 0;
  }
}

void SetMrelInfoInit (SetMrelInfo *results, int allocLists) {

  int i;
  SetMrelInfoReset (results);
  if (allocLists) {
    results->Nsecfilt = GetPhotcodeNsecfilt ();

    results->psfData  = StatDataSetAlloc (results->Nsecfilt, Nmax);
    results->aperData = StatDataSetAlloc (results->Nsecfilt, Nmax);
    results->kronData = StatDataSetAlloc (results->Nsecfilt, Nmax);

    ALLOCATE (results->psfqf_list,     double, Nmax);
    ALLOCATE (results->psfqfperf_list, double, Nmax);
    ALLOCATE (results->stargal_list,   double, Nmax);
    ALLOCATE (results->photflag_list,  uint32_t, Nmax);

    ALLOCATE (results->havePS1,   int, results->Nsecfilt);
    ALLOCATE (results->haveSYN,   int, results->Nsecfilt);
    ALLOCATE (results->needSYN,   int, results->Nsecfilt);
    ALLOCATE (results->measSYN,   int, results->Nsecfilt);
    ALLOCATE (results->minSYN,  float, results->Nsecfilt);

    ALLOCATE (results->psfQfMax,     float, results->Nsecfilt);
    ALLOCATE (results->psfQfPerfMax, float, results->Nsecfilt);

    ALLOCATE (results->Nmeas,          int, results->Nsecfilt);
    ALLOCATE (results->NmeasGood,      int, results->Nsecfilt);
    ALLOCATE (results->Next,           int, results->Nsecfilt);
    ALLOCATE (results->NexpPS1,        int, results->Nsecfilt);
    ALLOCATE (results->haveUbercal,    int, results->Nsecfilt);

    ALLOCATE (results->tessID,         int, results->Nsecfilt);
    ALLOCATE (results->projID,         int, results->Nsecfilt);
    ALLOCATE (results->skycellID,      int, results->Nsecfilt);

    ALLOCATE (results->minUbercalDist, float, results->Nsecfilt);

    SetMrelInfoResetObject (results);

    for (i = 0; i < NsynthNames; i++) {
      int photcode = GetPhotcodeCodebyName (synthNames[i]);
      if (!photcode) continue;
      int Nsec = GetPhotcodeNsec (photcode);
      if (Nsec < 0) continue;
      results->minSYN[Nsec] = synthMaxMags[i];
    }
  } else {
    results->Nsecfilt = 0;
    results->psfData  = NULL;
    results->aperData = NULL;
    results->kronData = NULL;

    results->psfqf_list     = NULL;
    results->psfqfperf_list = NULL;
    results->stargal_list   = NULL;
    results->photflag_list  = NULL;

    results->havePS1 = NULL;
    results->haveSYN = NULL;
    results->needSYN = NULL;
    results->measSYN = NULL;
    results->minSYN  = NULL;

    results->psfQfMax = NULL;
    results->psfQfPerfMax = NULL;

    results->Nmeas = NULL;
    results->NmeasGood = NULL;
    results->Next = NULL;
    results->NexpPS1 = NULL;
    results->haveUbercal = NULL;

    results->tessID = NULL;
    results->projID = NULL;
    results->skycellID = NULL;

    results->minUbercalDist = NULL;
  }
}

void SetMrelInfoFree (SetMrelInfo *results) {
  StatDataSetFree (results->psfData,  results->Nsecfilt);
  StatDataSetFree (results->aperData, results->Nsecfilt);
  StatDataSetFree (results->kronData, results->Nsecfilt);
  FREE (results->psfqf_list);
  FREE (results->psfqfperf_list);
  FREE (results->stargal_list);
  FREE (results->photflag_list);

  FREE (results->havePS1);
  FREE (results->haveSYN);
  FREE (results->needSYN);
  FREE (results->measSYN);
  FREE (results->minSYN );

  FREE (results->psfQfMax);
  FREE (results->psfQfPerfMax);

  FREE (results->Nmeas);
  FREE (results->NmeasGood);
  FREE (results->Next);
  FREE (results->NexpPS1);
  FREE (results->haveUbercal);

  FREE (results->tessID);
  FREE (results->projID);
  FREE (results->skycellID);

  FREE (results->minUbercalDist);
}

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++) {
    setMrel_catalog_alt (catalog, i, FALSE, flatcorr, &results, Nsecfilt);
    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);
  SetMrelInfoFree (&summary);

  return (TRUE);
}

// setMrelOutput is NOT threaded because we read/write catalogs one at a time (for now)
int setMrelOutput (Catalog *catalog, int Ncatalog, FlatCorrectionTable *flatcorr) {

  int i;

  int Nsecfilt = GetPhotcodeNsecfilt ();

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

  for (i = 0; i < Ncatalog; i++) {
    setMrel_catalog_alt (catalog, i, TRUE, flatcorr, &results, Nsecfilt);
    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);
  SetMrelInfoFree (&summary);
  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);
    SetMrelInfoFree (&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);

  SetMrelInfoFree (&summary);
  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;
      SetMrelInfoFree (&results);
      return NULL;
    }

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

    setMrel_catalog_alt (catalog, i, FALSE, flatcorr, &results, Nsecfilt);
    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; }


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_OBJ_POOR;
	  Ndel ++;
	  if (dM > MaxScatter)   { Nscat ++; }
	  if (isnan(Mchisq))     { Nnan ++; }
	  if (Mchisq > MaxChisq) { Nchi ++; }
	} else {
	  catalog[i].secfilt[Nsecfilt*j+Nsec].flags &= ~ID_OBJ_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 (> 5 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 ++;
	  } else {
	    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;
	    }
	  }
	}
	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((double)(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);
}

