# include "relastro.h"

// NOTE: we only measure the systematic floor of the astrometric scatter per stack, no change to the calibration
int UpdateStacks (Catalog *catalog, int Ncatalog) {

  if (SKIP_PS1_STACK) return TRUE;

  off_t Nimage;
  Image *image = getimages (&Nimage, NULL);

  off_t Nskip = 0;
  off_t Nmeas = 0;

  // measure the scatter for each stack
  for (off_t i = 0; i < Nimage; i++) {

    /* skip all except stack images */
    if (!isGPC1stack(image[i].photcode)) continue;

    /* convert measure coordinates to raw entries */
    off_t Nraw;
    StarData *raw = getImageRaw (catalog, Ncatalog, i, &Nraw, MODE_SIMPLE);
    if (!raw) {
      Nskip ++;
      continue;
    }
    if (Nraw <= IMFIT_TOO_FEW) {
      Nskip ++;
      free (raw);
      continue;
    }

    /* convert average coordinates to ref entries */
    off_t Nref;
    StarData *ref = getImageRef (catalog, Ncatalog, i, &Nref, MODE_SIMPLE);
    if (!ref) {
      Nskip ++;
      free (raw);
      continue;
    }

    // note that Nraw & Nref must be equal: if not, we made a programming error in one of these two functions.
    assert (Nraw == Nref);

    int Nstat;
    float dLsig, dMsig, dRsig;
    GetScatterRawRef(&dLsig, &dMsig, &dRsig, &Nstat, raw, ref, Nraw, IMFIT_SYS_SIGMA_LIM);

    // XXX: I need to convert dLsig, dMsig from degrees to pixels
    dLsig *= 3600.0;
    dMsig *= 3600.0;

    image[i].dXpixSys = dLsig;
    image[i].dYpixSys = dMsig;
    image[i].nFitAstrom = Nstat;
  }

  fprintf (stderr, "UpdateStacks: %d measured, %d skipped\n", (int) Nmeas, (int) Nskip);
  return (TRUE);
}

// XXX 2 hardwired hacks in this file: 1) photcode hardwired for GPC1 stacks, 2) platescale
