# include "relastro.h"

// XXX make these user parameters
# define FIT_CHIP_MAX_ERROR 0.05
# define FIT_CHIP_NITER     3
# define FIT_CHIP_NSIGMA    3.0

// XXX we should test if the fit is sufficiently constrained across the chip, or if the
// new positions deviate too much from the old positions

// XXX add visualization tools: per chip residual plots

// XXX save the fit[0].Npts value in the image table?

// XXX save measurements of the fit quality (scatter, chisq) in the image table

void FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords) {

  int i, Nscatter, Niter, skip;
  CoordFit *fit;
  double dL, dM, dR, dRmax, *values;

  ALLOCATE (values, double, Nmatch);

  for (Niter = 0; Niter < FIT_CHIP_NITER; Niter ++) {

    // measure the scatter distribution (use only the bright end detections)
    for (i = Nscatter = 0; i < Nmatch; i++) {
      if (raw[i].mask) continue;
      if (raw[i].dMag > FIT_CHIP_MAX_ERROR) continue;

      dL = raw[i].L - ref[i].L;
      dM = raw[i].M - ref[i].M;
      dR = hypot (dL, dM);
    
      values[Nscatter] = dR;
      Nscatter++;
    }

    // for a 2D Gaussian, 40% of the points are within 1 sigma; dRmax is ~ 3 sigma
    dsort (values, Nscatter);
    dRmax = FIT_CHIP_NSIGMA*values[(int)(0.40*Nscatter)];

    // fit the requested order polynomial
    if (CHIPORDER > 0) {
      coords[0].Npolyterms = CHIPORDER;
    }
    fit = fit_init (coords[0].Npolyterms);

    // generate the fit matches
    for (i = 0; i < Nmatch; i++) {
      if (raw[i].mask) continue;

      // only keep objects within dRmax
      dL = raw[i].L - ref[i].L;
      dM = raw[i].M - ref[i].M;
      dR = hypot (dL, dM);
      if (dR > dRmax) continue;
    
      fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, raw[i].dPos);
    }

    // check if the fit has enough data points for the polynomial order
    skip = FALSE;
    switch (coords[0].Npolyterms) {
      case 0:
      case 1:
	skip = (fit[0].Npts < 8);
	break;
      case 2:
	skip = (fit[0].Npts < 11);
	break;
      case 3:
	skip = (fit[0].Npts < 15);
	break;
      default:
	fprintf (stderr, "invalid chip order %d\n", coords[0].Npolyterms);
	abort ();
    }
    if (skip) {
      fprintf (stderr, "insufficient measurements (%d) for requested order (%d)\n", fit[0].Npts, coords[0].Npolyterms);
      fit_free (fit);
      free (values);
      return;
    }

    fprintf (stderr, "scatter limit: %f based on %d detections; using %d of %d for fit\n", dRmax, Nscatter, fit[0].Npts, Nmatch);

    // measure the fit, update the coords & object coordinates
    fit_eval (fit);
    fit_apply_coords (fit, coords);
    fit_free (fit);

    for (i = 0; i < Nmatch; i++) {
      XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, coords);
    }
  }

  free (values);
  return;
}

/* in the mosaic case, we have four coord systems of interest:
   R,D : the sky
   P,Q : the tangent plane
   L,M : the focal plane
   X,Y : the chip

   R,D -> P,Q (projection)
   P,Q -> L,M (polynomial transformation : DIS)
   L,M -> X,Y (polynomial transformation : WRP)

Some details about this function:

- the initial value of the clipping radius is set based on the distribution of the dR
  values for the bright sources.

- NITER clipping passes are performed

- the per-star astrometric errors are included in the fit.  The photcode table controls
  whether only the reported positional errors are used, or if the magnitudes errors are
  scaled to determine the error, and if there is a systematic floor for all sources.

- input detections are pre-filtered on the basis of the photFlags, etc, in bcatalog

- outlying detections are clipped in the iterative passes, but the clipped detections are
  not recorded

*/

/* example using fit_apply() :

  f = fopen ("test3.dat", "w");

  // apply new coords to raw (X,Y -> L,M)
  for (i = 0; i < Nmatch; i++) {
    fprintf (f, "%f %f  %f %f  ", raw[i].X, raw[i].Y, raw[i].L, raw[i].M);
    fit_apply (newfit, &L1, &M1, raw[i].X - coords[0].crpix1, raw[i].Y - coords[0].crpix2);
    fprintf (f, "%f %f\n", L1, M1);
  }
  fclose (f);
*/
