# include "relastro.h"
# define PAR_TOOFEW 5

int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m, int applyGalaxyOffset);
int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats);

static off_t Nmax;
static double *X, *dX;
static double *Y, *dY;
static double *R, *dR;
static double *D, *dD;
static double *pX;
static double *pY;
static double *T;
static double *dT;
static double *C_blue;
static double *C_red;

static Coords coords;

static time_t T2000;

void initFitStats (FitStats *fitStats) {
  fitStats->Nave = 0;  
  fitStats->Npm = 0;   
  fitStats->Npar = 0;  
  fitStats->Nskip = 0; 
  fitStats->Noffset = 0;
  return;
}

void sumFitStats (FitStats *srcFitStats, FitStats *tgtFitStats) {
  tgtFitStats->Nave    += srcFitStats->Nave    ;  
  tgtFitStats->Npm     += srcFitStats->Npm     ;   
  tgtFitStats->Npar    += srcFitStats->Npar    ;  
  tgtFitStats->Nskip   += srcFitStats->Nskip   ; 
  tgtFitStats->Noffset += srcFitStats->Noffset ;
  return;
}

void initObjectData (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].average[j].Nmeasure);
    }
  }

  ALLOCATE (R, double, MAX (1, Nmax));
  ALLOCATE (D, double, MAX (1, Nmax));
  ALLOCATE (T, double, MAX (1, Nmax));
  ALLOCATE (X, double, MAX (1, Nmax));
  ALLOCATE (Y, double, MAX (1, Nmax));

  ALLOCATE (dR, double, MAX (1, Nmax));
  ALLOCATE (dD, double, MAX (1, Nmax));
  ALLOCATE (dT, double, MAX (1, Nmax));
  ALLOCATE (dX, double, MAX (1, Nmax));
  ALLOCATE (dY, double, MAX (1, Nmax));

  ALLOCATE (pX, double, MAX (1, Nmax));
  ALLOCATE (pY, double, MAX (1, Nmax));

  ALLOCATE (C_blue, double, MAX (1, Nmax));
  ALLOCATE (C_red,  double, MAX (1, Nmax));

  /* project coordinates to a plane centered on the object with units of arcsec */
  InitCoords (&coords, "DEC--SIN");
  coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0;

  // use J2000 as a reference time
  T2000 = ohana_date_to_sec ("2000/01/01,12:00:00");
}  

void freeObjectData () {

  free (R);
  free (D);
  free (T);
  free (X);
  free (Y);

  free (dR);
  free (dD);
  free (dT);
  free (dX);
  free (dY);

  free (pX);
  free (pY);

  free (C_blue);
  free (C_red);
}  

// This function operates on both Measure and MeasureTiny.  In the big stages, this should
// be called with just MeasureTiny set and Measure == NULL
int UpdateObjects (Catalog *catalog, int Ncatalog, int Nloop) {

  initObjectData (catalog, Ncatalog);

  // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
  int Nsecfilt = GetPhotcodeNsecfilt ();
  if (Ncatalog) {
    assert (catalog[0].Nsecfilt == Nsecfilt);
  }

  FitStats sumStatsChips; initFitStats (&sumStatsChips);
  FitStats sumStatsStack; initFitStats (&sumStatsStack);

  int i;
  for (i = 0; i < Ncatalog; i++) {

    if (VERBOSE2) fprintf (stderr, "astrometrize catalog %d : "OFF_T_FMT" ave, "OFF_T_FMT" meas\n", i,  catalog[i].Naverage,  catalog[i].Nmeasure);

    FitStats fitStatsChips; initFitStats (&fitStatsChips);
    FitStats fitStatsStack; initFitStats (&fitStatsStack);

    off_t j;
    for (j = 0; j < catalog[i].Naverage; j++) {
      /* calculate the average value of R,D for a single star */
      off_t m = catalog[i].average[j].measureOffset;
      MeasureTiny *measure = &catalog[i].measureT[m];
      Measure *measureBig = catalog[i].measure ? &catalog[i].measure[m] : NULL;
      Average *average = &catalog[i].average[j];
      SecFilt *secfilt = &catalog[i].secfilt[j*Nsecfilt];

      UpdateObjects_Stack(average, secfilt, measure, measureBig, Nsecfilt, &fitStatsStack);
      UpdateObjects_Chips(average, secfilt, measure, measureBig, Nsecfilt, &fitStatsChips, i, m, Nloop);
    }
    if (VERBOSE2) fprintf (stderr, "catalog %d : chips "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n",  i,  fitStatsChips.Nave,  fitStatsChips.Npm,  fitStatsChips.Npar,  fitStatsChips.Nskip, fitStatsChips.Noffset);
    if (VERBOSE2) fprintf (stderr, "catalog %d : stack "OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par : Nskip "OFF_T_FMT", Noffset "OFF_T_FMT"\n",  i,  fitStatsStack.Nave,  fitStatsStack.Npm,  fitStatsStack.Npar,  fitStatsStack.Nskip, fitStatsStack.Noffset);
    sumFitStats (&fitStatsChips, &sumStatsChips);
    sumFitStats (&fitStatsStack, &sumStatsStack);
  }
  freeObjectData ();

  if (VERBOSE && (Ncatalog > 1)) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n",  (sumStatsChips.Nave + sumStatsChips.Npm + sumStatsChips.Npar),  sumStatsChips.Nave,  sumStatsChips.Npm,  sumStatsChips.Npar,  sumStatsChips.Nskip, sumStatsChips.Noffset);
  if (VERBOSE && (Ncatalog > 1)) fprintf (stderr, "fitted "OFF_T_FMT" objects ("OFF_T_FMT" ave, "OFF_T_FMT" pm, "OFF_T_FMT" par), skipped "OFF_T_FMT", "OFF_T_FMT" have too large an offset\n",  (sumStatsStack.Nave + sumStatsStack.Npm + sumStatsStack.Npar),  sumStatsStack.Nave,  sumStatsStack.Npm,  sumStatsStack.Npar,  sumStatsStack.Nskip, sumStatsStack.Noffset);
  return (TRUE);
}

// This function operates on both Measure and MeasureTiny.  In the big stages, this should
// be called with just MeasureTiny set and Measure == NULL
int UpdateObjects_Chips (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats, int i, off_t m, int applyGalaxyOffset) {

  int setRefColor = areImagesMatched();

  /* calculate the average value of R,D for a single star */

  PMFit fit;    memset (&fit,    0, sizeof(fit));
  PMFit fitAve; memset (&fitAve, 0, sizeof(fitAve)); fitAve.chisq = NAN;
  PMFit fitPM;  memset (&fitPM,  0, sizeof(fitPM));  fitPM.chisq = NAN;
  PMFit fitPAR; memset (&fitPAR, 0, sizeof(fitPAR)); fitPAR.chisq = NAN;

  // if we fail to fit the astrometry for some reason, we need to set/reset these
  average[0].flags |= ID_STAR_NO_ASTROM;
  average[0].ChiSqAve  = NAN;
  average[0].ChiSqPM   = NAN;
  average[0].ChiSqPar  = NAN;
  average[0].Npos = 0;

  // an object with no measurements is externally supplied
  if (average[0].Nmeasure == 0) return TRUE;

  int NcBlue = 0;
  int NcRed = 0;
  int N = 0;

  int mode = FIT_MODE; // start with the globally-defined fit mode

  int XVERB = FALSE;
  XVERB |= (average[0].objID == OBJ_ID_SRC) && (average[0].catID == CAT_ID_SRC);
  XVERB |= (average[0].objID == OBJ_ID_DST) && (average[0].catID == CAT_ID_DST);

  // find the basic properties of the detections for this object (Tmin, Tmax, Tmean)
  off_t k;
  for (k = 0; k < average[0].Nmeasure; k++) {

    if (XVERB) {
      char *date = ohana_sec_to_date (measure[k].t);
      int dbFlagsBig = measureBig ? measureBig[k].dbFlags : 0;
      fprintf (stderr, OFF_T_FMT" %f %f %s : 0x%08x : 0x%08x\n",  k, measure[k].R, measure[k].D, date, measure[k].dbFlags, dbFlagsBig);
      free (date);
    }

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

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

    // reset the bit to note that a detection was used (or not)
    measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
    if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; }

    // does the measurement pass the supplied filtering constraints?
    // MeasFilterTestTiny does not test psfQF
    // exclude bad detections based on: photcodes, psfQF, time range, photflags & astromBadMask, mag_inst
    int keepMeasure = measureBig ? MeasFilterTest(&measureBig[k], FALSE) : MeasFilterTestTiny(&measure[k], FALSE);
    if (!keepMeasure) {
      continue;
    }

    double Ri = getMeanR (&measure[k], average, secfilt);
    double Di = getMeanD (&measure[k], average, secfilt);

    // if we are correcting for the Galaxy Motion Model, only should apply it here 
    // (a) when we are working to correct the images (mean R,D assumed to be at J2000) and
    // (b) if we think the measure R,D is already at the image epoch position
    if (USE_GALAXY_MODEL && applyGalaxyOffset) {
      Ri -= measure[k].RoffGAL / 3600.0;
      Di -= measure[k].DoffGAL / 3600.0;
    }

    // XXX add in dR,dD GAL here

    // mark (as POOR) any measurements which are deviant from the mean by > ExcludeBogusRadius
    if (ExcludeBogus) {
      coords.crval1 = average[0].R;
      coords.crval2 = average[0].D;
      double Xi, Yi;
      RD_to_XY (&Xi, &Yi, Ri, Di, &coords);
      double radius = hypot(Xi, Yi);
      if (radius > ExcludeBogusRadius) {
	measure[k].dbFlags |= ID_MEAS_POOR_ASTROM;
	if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_POOR_ASTROM; }
	continue;
      }
      measure[k].dbFlags &= ~ID_MEAS_POOR_ASTROM;
      if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_POOR_ASTROM; }
    }

    // outlier rejection
    if (FALSE && FlagOutlier && (measure[k].dbFlags & ID_MEAS_POOR_ASTROM)) {
      continue;
    }

    R[N] = Ri;
    D[N] = Di;

    // measure[k].t is UNIX seconds, T2000 is UNIX seconds for J2000.
    // T[] is time in years since J2000 (jd = 2451545)
    T[N] = (measure[k].t - T2000) / (86400*365.25) ; // time relative to J2000 in years

    // dX, dY : error in arcsec -- 
    dX[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_RA);
    dY[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_DEC);

    // allow a given photcode or measurement to be
    // ignored if the error is NAN (for photcode, set astromErrSys to NaN)
    if (isnan(dX[N])) continue;
    if (isnan(dY[N])) continue;

    // add systematic error in quadrature, if desired
    // only do this after the fit has converged (or you will never improve the poor images)
    // if (INCLUDE_SYS_ERR) {
    // float dRsys = FromShortPixels(measure[k].dRsys);
    // dX[N] = hypot(dX[N], dRsys);
    // dY[N] = hypot(dY[N], dRsys);
    // }

    // dX[N] = 0.1;
    // dY[N] = 0.1;

    dT[N] = measure[k].dt;

    // XXX this is (slightly) inconsistent: dX,dY are the X and Y direction errors in
    // arcseconds.  dR, dD are the errors in those directions in degrees.  IF we have
    // non-circular errors (different values for X and Y), then dR and dD will be
    // incorrect: they would need to be rotated to take out the position angle
    dR[N] = dX[N] / 3600.0;
    dD[N] = dY[N] / 3600.0;

    if (setRefColor) {
      float colorBlue = getColorBlue (m+k, i);
      if (!isnan(colorBlue)) {
	C_blue[NcBlue] = colorBlue;
	NcBlue++;
      }
      float colorRed = getColorRed (m+k, i);
      if (!isnan(colorRed)) {
	C_red[NcRed] = colorRed;
	NcRed++;
      }
    }

    measure[k].dbFlags |= ID_MEAS_USED_OBJ;
    if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_USED_OBJ; }

    N++;
  } // loop over measurements : average[0].Nmeasure 

  if (N < 1) { 
    if (isfinite(average[0].Rstk) && isfinite(average[0].Dstk)) {
      average[0].R  = average[0].Rstk;
      average[0].D  = average[0].Dstk;
      average[0].dR = average[0].dRstk;
      average[0].dD = average[0].dDstk;
    }
    return FALSE;
  }

  // if we have too few good detections for the desired fit, or too limited a
  // baseline, use a fit with fewer parameters.  XXX if we have too few measurements
  // for even the average position, consider including the lower-quality detections?

  // find Tmin & Tmax from the list of accepted measurements
  double Tmean = 0.0;
  double Tmin = T[0];
  double Tmax = T[0];
  for (k = 0; k < N; k++) {
    Tmin = MIN(Tmin, T[k]);
    Tmax = MAX(Tmax, T[k]);
    Tmean += T[k];
  }
  double Trange = Tmax - Tmin;

  if (RELASTRO_OP == OP_HIGH_SPEED) {
    Tmean = 0.5*(Tmax - Tmin);
  } else {
    Tmean /= (float) N;
  }

  /* we need to do the fit in a locally linear space; choose a ref coordinate */
  coords.crval1 = R[0];
  coords.crval2 = D[0];

  // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq

  // project all of the R,D coordinates to a plane centered on this coordinate. set
  // the times to be relative to Tmean (this is required for parallax as well)
  for (k = 0; k < N; k++) {
    RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords);
    T[k] -= Tmean;
    if (XVERB) {
      fprintf (stderr, OFF_T_FMT" %f %f %f  %f %f +/- %f %f\n",  k, T[k], R[k], D[k], X[k], Y[k], dX[k], dY[k]);
    }
  }	  

  // *** first fit for the proper motion (skip fit if Trange or Npts is too small) ***
  if ((mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) {
    if (Trange < PM_DT_MIN) {
      mode = FIT_AVERAGE;
      goto skipPM;
    }
    if (N <= PM_TOOFEW) {
      mode = FIT_AVERAGE;
      goto skipPM;
    }

    FitPM (&fitPM, X, dX, Y, dY, T, N, XVERB);

    if (XVERB) fprintf (stderr, "fitted PM:  %f - %f : %f %f : %f %f : %f vs %f\n", Tmin, Tmax, fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.chisq, fitAve.chisq);

    // project Ro, Do back to RA,DEC
    XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &coords);
    if (XVERB) fprintf (stderr, "project: %f %f : %f %f : %f\n", fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.p);
    if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");
    // XXX : does this make sense at 0,360 boundary?

    fitPM.p  = fitPM.dp  = 0.0;
    average[0].flags |= ID_STAR_FIT_PM;
    fitStats->Npm ++;

    // XXX a hard-wired hack...
    if ((fabs(fitPM.uR) > 2.0) || (fabs(fitPM.uD) > 2.0)) {
      mode = FIT_AVERAGE;
      average[0].flags |= ID_STAR_BAD_PM;
    }
  }
  
skipPM:
  // fit the parallax + proper-motion model
  // NOTE : we only fit PAR if we have already fitted for proper motion. if we do not fit PM or we fail
  // to fit PM, we do not attempt PAR.  thus failure to fit PAR falls back to PM-only
  if (mode == FIT_PM_AND_PAR) {
    if (Trange < PM_DT_MIN) {
      mode = FIT_PM_ONLY;
      goto skipPAR;
    }
    if (N <= PAR_TOOFEW) {
      mode = FIT_PM_ONLY;
      goto skipPAR;
    }
    float pXmin = +2.0;
    float pXmax = -2.0;
    float pYmin = +2.0;
    float pYmax = -2.0;
    for (k = 0; k < N; k++) {
      ParFactor (&pX[k], &pY[k], R[k], D[k], T[k], Tmean);
      pXmin = MIN (pXmin, pX[k]);
      pXmax = MAX (pXmax, pX[k]);
      pYmin = MIN (pYmin, pY[k]);
      pYmax = MAX (pYmax, pY[k]);
    }
    float dXRange = pXmax - pXmin;
    float dYRange = pYmax - pYmin;
    float parRange = hypot (dXRange, dYRange);
	
    if (parRange < PAR_FACTOR_MIN) {
      mode = FIT_PM_ONLY;
      goto skipPAR;
    }

    FitPMandPar (&fitPAR, X, dX, Y, dY, T, pX, pY, N, XVERB);
    if (XVERB) fprintf (stderr, "fitted PM+PAR:  %f - %f : %f %f : %f %f : %f %f : %f vs %f vs %f\n", Tmin, Tmax, fitPAR.Ro, fitPAR.Do, fitPAR.uR, fitPAR.uD, fitPAR.p, fitPAR.dp, fitPAR.chisq, fitPM.chisq, fitAve.chisq);

    XY_to_RD (&fitPAR.Ro, &fitPAR.Do, fitPAR.Ro, fitPAR.Do, &coords);
    average[0].flags |= ID_STAR_FIT_PAR;
    fitStats->Npar ++;

    if (fabs(fitPM.Ro) < 0.01) fprintf (stderr, "watch out for 0,360 boundary\n");

    // XXX a hard-wired hack...
    if ((fabs(fitPAR.uR) > 2.0) || (fabs(fitPAR.uD) > 2.0)) {
      mode = FIT_PM_ONLY;
    }
  }	  

skipPAR:
  {
    // ALWAYS fit the average model
    StatType statsR, statsD;
    liststats_pos (X, dX, N, &statsR, XVERB); // WARNING: this function modifies R (do not use after here)
    liststats_pos (Y, dY, N, &statsD, XVERB); // WARNING: this function modifies D (do not use after here)

    // project Ro, Do back to RA,DEC
    XY_to_RD (&fitAve.Ro, &fitAve.Do, statsR.mean, statsD.mean, &coords);
    if (XVERB) fprintf (stderr, "average: %f %f\n", fitAve.Ro, fitAve.Do);

    fitAve.dRo = statsR.sigma;
    fitAve.dDo = statsD.sigma;

    fitAve.chisq = (N > 1) ? 0.5 * (statsR.chisq + statsD.chisq) : NAN;
    fitAve.Nfit = N;

    fitAve.uR = fitAve.duR = 0.0;
    fitAve.uD = fitAve.duD = 0.0;
    fitAve.p  = fitAve.dp  = 0.0;
    average[0].flags |= ID_STAR_FIT_AVE;
    fitStats->Nave ++;
  }

  if (setRefColor) {
    float colorMedian;
    dsort (C_blue, NcBlue);
    colorMedian = (NcBlue > 0) ? C_blue[(int)(0.5*NcBlue)] : NAN;
    average[0].refColorBlue = colorMedian;
    dsort (C_red, NcRed);
    colorMedian = (NcRed > 0) ? C_red[(int)(0.5*NcRed)] : NAN;
    average[0].refColorRed = colorMedian;
  }

  /* choose the result based on the chisq values */
  // XXXX for now, just use the mode as the result:
  int result = mode;

  switch (result) {
    case FIT_AVERAGE:
      average[0].flags |= ID_STAR_USE_AVE;
      fit = fitAve;
      break;
    case FIT_PM_ONLY:
      average[0].flags |= ID_STAR_USE_PM;
      fit = fitPM;
      break;
    case FIT_PM_AND_PAR:
      average[0].flags |= ID_STAR_USE_PAR;
      fit = fitPAR;
      break;
  }
  if (XVERB) {
    fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f) plx=(%f +/- %f)\n",
	     average[0].R, 
	     average[0].D, 
	     fit.Ro, fit.Do, 
	     3600*(average[0].R - fit.Ro), 
	     3600*(average[0].D - fit.Do),
	     fit.uR, fit.uD, fit.p, fit.dp);
    fprintf (stderr, "-----\n");
  }

  // make sure that the fit succeeded
  int status = TRUE;
  status &= finite(fit.Ro);
  status &= finite(fit.Do);
  status &= finite(fit.dRo);
  status &= finite(fit.dDo);
  status &= finite(fit.uR);
  status &= finite(fit.uD);
  status &= finite(fit.duR);
  status &= finite(fit.duD);
  status &= finite(fit.p);
  status &= finite(fit.dp);
  if (!status) {
    fitStats->Nskip ++;
    return FALSE; // XXX ??
  }

  // what is the offset relative to the mean fit position?
  coords.crval1 = average[0].R;
  coords.crval2 = average[0].D;
  if (isnan(coords.crval1)) {
    return (FALSE);
  }
  if (isnan(coords.crval2)) {
    return (FALSE);
  }

  double dXoff, dYoff;
  RD_to_XY (&dXoff, &dYoff, fit.Ro, fit.Do, &coords);
  float dPos = hypot (dXoff, dYoff);
  if (dPos > MaxMeanOffset) {
    if (fitStats->Noffset < 100) {
      fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", coords.crval1, coords.crval2, fit.Ro, fit.Do, dXoff, dYoff);
    }
    fitStats->Noffset ++;
    return FALSE; // XXX ??
  }

  if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f) chisq=(%f, %f, %f)\n",
		      average[0].R,
		      average[0].D,
		      fit.Ro, fit.Do,
		      3600*(average[0].R - fit.Ro),
		      3600*(average[0].D - fit.Do),
		      average[0].uR,
		      average[0].uD,
		      fitAve.chisq, fitPM.chisq, fitPAR.chisq);

  average[0].R  	= fit.Ro; // RA in degrees
  average[0].D  	= fit.Do; // DEC in degrees
  average[0].dR 	= fit.dRo; // RA scatter in arcsec
  average[0].dD 	= fit.dDo; // DEC scatter in arcsec

  average[0].uR         = fit.uR; // RA proper motion in arcsec/year
  average[0].uD         = fit.uD; // DEC proper motion in arcsec/year
  average[0].duR        = fit.duR; // RA proper motion error in arcsec/year
  average[0].duD        = fit.duD; // DEC proper motion error in arcsec/year

  average[0].P          = fit.p; // parallax in arcsec
  average[0].dP         = fit.dp; // parallax error in arcsec

  average[0].ChiSqAve   = fitAve.chisq;
  average[0].ChiSqPM    = fitPM.chisq;
  average[0].ChiSqPar   = fitPAR.chisq;

  average[0].Tmean      = (Tmean * 86400 * 365.25) + T2000;
  average[0].Trange     = (Trange * 86400 * 365.25);
  average[0].Npos       = fit.Nfit;

  // unset the NO_ASTROM bit (not(NO_ASTROM) == HAVE_ASTROM)
  average[0].flags &= ~ID_STAR_NO_ASTROM;

  return (TRUE);
}

// This function operates on both Measure and MeasureTiny.  In the big stages, this should
// be called with just MeasureTiny set and Measure == NULL
int UpdateObjects_Stack (Average *average, SecFilt *secfilt, MeasureTiny *measure, Measure *measureBig, int Nsecfilt, FitStats *fitStats) {

  off_t k;

  // set the default values
  average[0].Rstk  = NAN; // RA in degrees
  average[0].Dstk  = NAN; // DEC in degrees
  average[0].dRstk = NAN; // RA scatter in arcsec
  average[0].dDstk = NAN; // DEC scatter in arcsec

  /* calculate the average value of R,D for a single star */
  PMFit fitAve;
  memset (&fitAve, 0, sizeof(fitAve));
  fitAve.chisq = NAN;

  if (average[0].Nmeasure == 0) return TRUE;

  int N = 0;

  int XVERB = FALSE;
  XVERB |= (average[0].objID == OBJ_ID_SRC) && (average[0].catID == CAT_ID_SRC);
  XVERB |= (average[0].objID == OBJ_ID_DST) && (average[0].catID == CAT_ID_DST);

  // find the basic properties of the detections for this object (Tmin, Tmax, Tmean)
  for (k = 0; k < average[0].Nmeasure; k++) {

    if (XVERB) {
      char *date = ohana_sec_to_date (measure[k].t);
      int dbFlagsBig = measureBig ? measureBig[k].dbFlags : 0;
      fprintf (stderr, "stack: "OFF_T_FMT" %f %f %s : 0x%08x : 0x%08x\n",  k, measure[k].R, measure[k].D, date, measure[k].dbFlags, dbFlagsBig);
      free (date);
    }

    // SKIP everything except gpc1 stack data
    if (!isGPC1stack(measure[k].photcode)) continue;

    // exclude bad detections based on: photcodes, psfQF, time range, photflags & astromBadMask, mag_inst
    int keepMeasure = measureBig ? MeasFilterTest(&measureBig[k], FALSE) : MeasFilterTestTiny(&measure[k], FALSE);
    if (!keepMeasure) {
      continue;
    }

    R[N] = getMeanR (&measure[k], average, secfilt);
    D[N] = getMeanD (&measure[k], average, secfilt);

    // dX, dY : error in arcsec -- 
    dX[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_RA);
    dY[N] = GetAstromErrorTiny (&measure[k], ERROR_MODE_DEC);

    // allow a given photcode or measurement to be
    // ignored if the error is NAN (for photcode, set astromErrSys to NaN)
    if (isnan(dX[N])) continue;
    if (isnan(dY[N])) continue;

    // XXX this is (slightly) inconsistent: dX,dY are the X and Y direction errors in
    // arcseconds.  dR, dD are the errors in those directions in degrees.  IF we have
    // non-circular errors (different values for X and Y), then dR and dD will be
    // incorrect: they would need to be rotated to take out the position angle
    dR[N] = dX[N] / 3600.0;
    dD[N] = dY[N] / 3600.0;

    // XXX use a different flag for stack measurements?
    // measure[k].dbFlags |= ID_MEAS_USED_OBJ;
    // if (measureBig) { measureBig[k].dbFlags |= ID_MEAS_USED_OBJ; }

    N++;
  } // loop over measurements : average[0].Nmeasure 

  // if we have too few good detections for the desired fit, or too limited a
  // baseline, use a fit with fewer parameters.  XXX if we have too few measurements
  // for even the average position, consider including the lower-quality detections?

  // too few measurements for average position (require 2 values)
  if (N < 1) return FALSE; // XXX ?? 

  // find the mean position
  StatType statsR, statsD;
  liststats_pos (R, dR, N, &statsR, XVERB); // WARNING: this function modifies R (do not use after here)
  liststats_pos (D, dD, N, &statsD, XVERB); // WARNING: this function modifies D (do not use after here)

  fitAve.Ro = statsR.mean;
  fitAve.dRo = 3600.0*statsR.sigma;

  fitAve.Do = statsD.mean;
  fitAve.dDo = 3600.0*statsD.sigma;

  fitAve.chisq = 0.5 * (statsR.chisq + statsD.chisq);
  fitAve.Nfit = N;

  // XXX choose stack flag? average[0].flags |= ID_STAR_FIT_AVE;
  fitStats->Nave ++;

  if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f)\n",
		      average[0].R, 
		      average[0].D, 
		      fitAve.Ro, fitAve.Do, 
		      3600*(average[0].R - fitAve.Ro), 
		      3600*(average[0].D - fitAve.Do));

  // make sure that the fit succeeded
  int status = TRUE;
  status &= finite(fitAve.Ro);
  status &= finite(fitAve.Do);
  status &= finite(fitAve.dRo);
  status &= finite(fitAve.dDo);
  if (!status) {
    fitStats->Nskip ++;
    return FALSE;
  }

  // what is the offset relative to the mean fit position?
  coords.crval1 = average[0].R;
  coords.crval2 = average[0].D;

  double dXoff, dYoff;
  RD_to_XY (&dXoff, &dYoff, fitAve.Ro, fitAve.Do, &coords);
  float dPos = hypot (dXoff, dYoff);
  if (dPos > MaxMeanOffset) {
    if (fitStats->Noffset < 100) {
      fprintf (stderr, "(%f,%f) -> (%f,%f) (%f,%f)\n", coords.crval1, coords.crval2, fitAve.Ro, fitAve.Do, dXoff, dYoff);
    }
    fitStats->Noffset ++;
    return FALSE;
  }

  // set the stack position values
  average[0].Rstk  = fitAve.Ro; // RA in degrees
  average[0].Dstk  = fitAve.Do; // DEC in degrees
  average[0].dRstk = fitAve.dRo; // RA scatter in arcsec
  average[0].dDstk = fitAve.dDo; // DEC scatter in arcsec

  return (TRUE);
}



/* fitting proper-motion and parallax:

   given a source at position r,d, at a time t, we need to calculate a vector (pr,pd)

   let x,y be the coordinate in the linearized frame with y parallel to DEC lines

   L,B are the ecliptic longitude and latitude of the object, 
   dL and dB are the offsets in the L and B directions

   dL = sin(t - topp)
   dB = cos(t - topp)*sin(B)

   these need to be rotated to the R,D frame to yield pR,pD.  Then, the equation of motion
   for the source in the x,y frame is:

   x = Ro + uR * (t - to) + p * pR 
   y = Do + uD * (t - to) + p * pD

   the unknowns in these equations are Ro, uR, Do, uD, and p

   XXX think through the concepts for the pole a bit better.  all objects near the pole 
   move the same way with the same phase.  choose a projection center and define dL,dB relative 
   to that center point coordinate system?

*/
