# include "relastro.h"
void dump_measures (Average *average, Measure *measure); // in ImageOps.c

// this function is called by UpdateObjectOffsets or relastro_images.c for reset
int UpdateMeasures (Catalog *catalog, int Ncatalog) {

  off_t i, j, Nimage;
  Image *image;
  double DPOS_MAX_ASEC;

  int badImage = 
    ID_IMAGE_ASTROM_NOCAL | 
    ID_IMAGE_ASTROM_POOR | 
    ID_IMAGE_ASTROM_FAIL | 
    ID_IMAGE_ASTROM_SKIP | 
    ID_IMAGE_ASTROM_FEW;

  if (RESET_BAD_IMAGES) badImage = 0;

  image = getimages (&Nimage, NULL);

  // track measurements which are far from their original position (poor at pole)
  int NoffRAave = 0;  int NoffRAori = 0; 
  int NoffDECave = 0; int NoffDECori = 0;

  for (i = 0; i < Ncatalog; i++) {
    myAssert (!catalog[i].Nmeasure || catalog[i].measureT, "programming error");
    for (j = 0; j < catalog[i].Nmeasure; j++) {
      MeasureTiny *measureT = &catalog[i].measureT[j];
      Measure     *measureB = catalog[i].measure ? &catalog[i].measure[j] : NULL;

      off_t im = getImageByID (measureT->imageID);
      if (im < 0) continue; // detections without imageIDs are not associated with a known image (eg, refs)

      // check that we have the right image:
      myAssert (measureT->photcode == image[im].photcode, "image photcode mismatch");
      myAssert (measureT->t >= image[im].tzero, "image time mismatch (1)");
      myAssert (measureT->t <= image[im].tzero + image[im].NX*image[im].trate/5000.0, "image time mismatch (2)");
      // I'm allowing the trate*NX to be a factor of 2x too small

      // skip measurements on images that have failed solutions (divergent or otherwise)
      // XXX apply this or not??
      if (FALSE && (image[im].flags & badImage)) continue;

      Coords *moscoords = image[im].coords.mosaic;
      Coords *imcoords = &image[im].coords;

      if (moscoords) {
	DPOS_MAX_ASEC = 3600.0*DPOS_MAX*hypot(moscoords->cdelt1, moscoords->cdelt2);
      } else {
	DPOS_MAX_ASEC = 3600.0*DPOS_MAX*hypot(imcoords->cdelt1, imcoords->cdelt2);
      }

      double X = measureT->Xccd;
      double Y = measureT->Yccd;
      if (USE_FIXED_PIXCOORDS) {
	if (isfinite(measureT->Xfix) && isfinite(measureT->Yfix)) {
	  float dX = measureT->Xfix - measureT->Xccd;
	  float dY = measureT->Yfix - measureT->Yccd;
	  if (hypot(dX,dY) < 2.0) {
	    X = measureT->Xfix;
	    Y = measureT->Yfix;
	  } 
	}
      }

      double R, D;
      XY_to_RD (&R, &D, X, Y, imcoords);

      int n = measureT[0].averef;
      Average *average = &catalog[i].average[n];

      double dR = 3600.0*(average[0].R - R);
      double dD = 3600.0*(average[0].D - D);

      // make sure detection is on the same side of the 0,360 boundary as the average
      // this will give some funny results withing ~1 arcsec of the pol
      if (dR > +180.0*3600.0) {
	// average on high end of boundary, move star up
	R += 360.0;
	dR = 3600.0*(average[0].R - R);
      }
      if (dR < -180.0*3600.0) {
	// average on low end of boundary, move star down
	R -= 360.0;
	dR = 3600.0*(average[0].R - R);
      }

      float csdec = cos(average[0].D * RAD_DEG);

      // complain if the new location is far from the average location
      // NOTE: This should never happen, or our StarMap tests are not working
      if (fabs(dR*csdec) > 3.0*ADDSTAR_RADIUS) {
	NoffRAave ++;
	if (VERBOSE2 && catalog[i].measure) {
	  fprintf (stderr, "measurement is far from average location (R): %f %f (%f %f %f)\n", average[0].R, average[0].D, dR, dR*csdec, dD);
	  dump_measures (&average[0], catalog[i].measure);
	}
      }
      if (fabs(dD) > 3.0*ADDSTAR_RADIUS) {
	NoffDECave ++;
	if (VERBOSE2 && catalog[i].measure) {
	  fprintf (stderr, "measurement is far from average location (D): %f %f (%f %f)\n", average[0].R, average[0].D, dR, dD);
	  dump_measures (&average[0], catalog[i].measure);
	}
      }

      // complain if the new location is far from the old location
      if (fabs(csdec*(measureT[0].R - R)) > DPOS_MAX_ASEC) {
	NoffRAori ++;
	if (VERBOSE2 && catalog[i].measure) {
	  fprintf (stderr, "measurement moves far from original location (R): %f %f (%f %f %f %f)\n", average[0].R, average[0].D, measureT[0].R, dR, csdec*(measureT[0].R - R), dD);
	  dump_measures (&average[0], catalog[i].measure);
	}
      }
      if (fabs(measureT[0].D - D) > DPOS_MAX_ASEC) {
	NoffDECori ++;
	if (VERBOSE2 && catalog[i].measure) {
	  fprintf (stderr, "measurement moves far from original location (D): %f %f (%f %f %f)\n", average[0].R, average[0].D, dR, measureT[0].D, dD);
	  dump_measures (&average[0], catalog[i].measure);
	}
      }
      measureT->R = R;
      measureT->D = D;
      if (measureB) {
	measureB->R = R;
	measureB->D = D;
      }
    }
  }
  // printNcatTotal ();

  int Noff = NoffRAave + NoffDECave + NoffRAori + NoffDECori;
  if (VERBOSE && (Noff > 0)) fprintf (stderr, "Noff ave RA %d, Noff ave DEC %d, Noff ori RA %d, Noff ori DEC %d\n", NoffRAave, NoffDECave, NoffRAori, NoffDECori);

  return (TRUE);
}

// this function operates on Measure, not MeasureTiny
int UpdateMeasuresOld (Catalog *catalog, int Ncatalog) {

  off_t i, Nimage;
  Image *image;

  int badImage = 
    ID_IMAGE_ASTROM_NOCAL | 
    ID_IMAGE_ASTROM_POOR | 
    ID_IMAGE_ASTROM_FAIL | 
    ID_IMAGE_ASTROM_SKIP | 
    ID_IMAGE_ASTROM_FEW;

  if (RESET_BAD_IMAGES) badImage = 0;

  image = getimages (&Nimage, NULL);

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

    /* skip DIS images (since they have no associated detections) */
    if (!strcmp(&image[i].coords.ctype[4], "-DIS")) continue;

    // skip images that have failed solutions (divergent or otherwise)
    if (image[i].flags & badImage) continue;

    /* convert measure coordinates to raw entries */
    fixImageRaw (catalog, Ncatalog, i);
  }

  // printNcatTotal ();

  return (TRUE);
}

