# include "relastro.h"
int plotChipFits (double *Ro, double *Do, char *mode, int Nimage);
int saveCenter (Image *image, double *Ro, double *Do, int im);

int UpdateChips (Catalog *catalog, int Ncatalog) {

  int Nskip, Nmosaic, NnewFit, NoldFit;

  /* we can measure new image parameters for each non-mosaic chip independently */
  off_t i, Nimage, Nraw, Nref, nFitAstr;
  Image *image;
  StarData *raw, *ref;
  Coords *oldCoords;
  float dXpixSys, dYpixSys;
  double *Ro, *Do;
  char *mode;

  Nskip = Nmosaic = NnewFit = NoldFit = 0;

  image = getimages (&Nimage, NULL);

  // save fit results for summary plot
  ALLOCATE (Ro, double, Nimage);
  ALLOCATE (Do, double, Nimage);
  ALLOCATE (mode, char, Nimage);

  for (i = 0; i < Nimage; i++) {

    /* skip all except WRP images */
    if (strcmp(&image[i].coords.ctype[4], "-WRP")) {
      Nmosaic ++;
      mode[i] = 0;
      continue;
    }

    /* convert measure coordinates to raw entries */
    raw = getImageRaw (catalog, Ncatalog, i, &Nraw, MODE_MOSAIC);
    if (!raw) {
      Nskip ++;
      mode[i] = 0;
      continue;
    }

    /* convert average coordinates to ref entries */
    ref = getImageRef (catalog, Ncatalog, i, &Nref, MODE_MOSAIC);
    if (!ref) {
      Nskip ++;
      mode[i] = 0;
      continue;
    }

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

    // save these in case of failure
    saveCoords (&image[i].coords, i);
    dXpixSys = image[i].dXpixSys;
    dYpixSys = image[i].dYpixSys;
    nFitAstr = image[i].nFitAstrom;

    // FitChip does iterative, clipped fitting
    // fprintf (stderr, "image "OFF_T_FMT" : Nstars: "OFF_T_FMT"\n",  i,  Nraw);
    if (!FitChip (raw, ref, Nraw, &image[i])) {
      if (VERBOSE) fprintf (stderr, "reject fit for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name,  i,  Nraw);

      // restore status quo ante
      oldCoords = getCoords (i);
      memcpy (&image[i].coords, oldCoords, sizeof(Coords));
      image[i].dXpixSys = dXpixSys;
      image[i].dYpixSys = dYpixSys; 
      image[i].nFitAstrom = nFitAstr;

      saveCenter (image, &Ro[i], &Do[i], i);
      mode[i] = 1;
      NoldFit ++;
      free (raw);
      free (ref);
      continue;
    }

    if (!checkStarMap (i)) {
      if (VERBOSE) fprintf (stderr, "fit diverges too much for image %s ("OFF_T_FMT") : Nstars: "OFF_T_FMT"\n", image[i].name,  i,  Nraw);
      // restore status quo ante
      oldCoords = getCoords (i);
      memcpy (&image[i].coords, oldCoords, sizeof(Coords));
      image[i].dXpixSys = dXpixSys;
      image[i].dYpixSys = dYpixSys; 
      image[i].nFitAstrom = nFitAstr;

      saveCenter (image, &Ro[i], &Do[i], i);
      mode[i] = 2;
      image[i].flags |= ID_IMAGE_ASTROM_POOR;
      NoldFit ++;
      free (raw);
      free (ref);
      continue;
    } 

    saveCenter (image, &Ro[i], &Do[i], i);
    mode[i] = 3;
    NnewFit ++;
    free (raw);
    free (ref);
  }

  plotChipFits (Ro, Do, mode, Nimage);

  fprintf (stderr, "UpdateChips: %d fitted, %d keep old, %d skipped, %d mosaic (skipped)\n", NnewFit, NoldFit, Nskip, Nmosaic);
  return (TRUE);
}

int saveCenter (Image *image, double *Ro, double *Do, int im) {

  Mosaic *mosaic;
  Coords *moscoords, *imcoords;
  double X, Y, L, M, P, Q, R, D;

  moscoords = NULL;
  if (!strcmp(&image[im].coords.ctype[4], "-WRP")) {
    mosaic = getMosaicForImage (im);
    if (mosaic == NULL) return FALSE;  // if we cannot find the associated image, skip it
    moscoords = &mosaic[0].coords;
  }
  imcoords = &image[im].coords;
  
  if (!strcmp(&image[im].coords.ctype[4], "-WRP")) {
    X = 0.5*image[im].NX; 
    Y = 0.5*image[im].NY;
  } else {
    X = 0.0; 
    Y = 0.0;
  }

  if (moscoords == NULL) {
    // this is a Simple image (not a mosaic)
    // note that for a Simple image, L,M = P,Q
    XY_to_LM (&L, &M, X, Y, imcoords);
    LM_to_RD (&R, &D, L, M, imcoords);
  } else {
    XY_to_LM (&L, &M, X, Y, imcoords);
    XY_to_LM (&P, &Q, L, M, moscoords);
    LM_to_RD (&R, &D, P, Q, moscoords);
  }

  double Rmid;
  if (UserCatalog) {
      Rmid = UserCatalogRA;
  } else {
      Rmid = 0.5*(UserPatch.Rmin + UserPatch.Rmax);
  }

  R = ohana_normalize_angle_to_midpoint (R, Rmid);

  *Ro = R;
  *Do = D;

  return (TRUE);
}
 
int plotChipFits (double *Ro, double *Do, char *mode, int Nimage) {

  static int kapa = -1; 

  int i, N;
  double Rmin, Rmax, Dmin, Dmax;
  float *xvec, *yvec;
  Graphdata graphdata;

  if (!relastroGetVisual()) return (TRUE);

  if (kapa == -1) {
    kapa = KapaOpenNamedSocket("kapa", "relastro");
    if (kapa == -1) {
      fprintf (stderr, "can't open kapa window\n");
      return FALSE;
    }
  }

  Rmin = +720;
  Rmax = -720;
  Dmin = +90;
  Dmax = -90;

  // find the R, D range
  for (i = 0; i < Nimage; i++) {
    if (!mode[i]) continue;

    Rmin = MIN(Rmin, Ro[i]);
    Rmax = MAX(Rmax, Ro[i]);
    Dmin = MIN(Dmin, Do[i]);
    Dmax = MAX(Dmax, Do[i]);
  }

  ALLOCATE (xvec, float, Nimage);
  ALLOCATE (yvec, float, Nimage);

  bzero (&graphdata, sizeof(Graphdata));
  plot_defaults (&graphdata);
  graphdata.xmin = Rmin;
  graphdata.xmax = Rmax;
  graphdata.ymin = Dmin;
  graphdata.ymax = Dmax;
  graphdata.style = 2;
  graphdata.size = 1;

  KapaSetFont (kapa, "helvetica", 14);
  KapaSetLimits (kapa, &graphdata);
  KapaBox (kapa, &graphdata);

  // *** good images ***
  N = 0;
  for (i = 0; i < Nimage; i++) {
    if (mode[i] != 3) continue;
    xvec[N] = Ro[i];
    yvec[N] = Do[i];
    N++;
  }
  graphdata.ptype = 7;
  graphdata.color = KapaColorByName("black");
  KapaPrepPlot (kapa, N, &graphdata);
  KapaPlotVector (kapa, N, xvec, "x");
  KapaPlotVector (kapa, N, yvec, "y");
  
  // *** reject fit ***
  N = 0;
  for (i = 0; i < Nimage; i++) {
    if (mode[i] != 1) continue;
    xvec[N] = Ro[i];
    yvec[N] = Do[i];
    N++;
  }
  graphdata.ptype = 3;
  graphdata.color = KapaColorByName("red");
  KapaPrepPlot (kapa, N, &graphdata);
  KapaPlotVector (kapa, N, xvec, "x");
  KapaPlotVector (kapa, N, yvec, "y");
  
  // *** divergent fit ***
  N = 0;
  for (i = 0; i < Nimage; i++) {
    if (mode[i] != 2) continue;
    xvec[N] = Ro[i];
    yvec[N] = Do[i];
    N++;
  }
  graphdata.ptype = 2;
  graphdata.color = KapaColorByName("blue");
  KapaPrepPlot (kapa, N, &graphdata);
  KapaPlotVector (kapa, N, xvec, "x");
  KapaPlotVector (kapa, N, yvec, "y");
  
  free (xvec);
  free (yvec);

  return (TRUE);
}

// XXX if (!FindMosaicForImage (image, Nimage, i)) { }
