# include "addstar.h"
# define EXTERNAL_ID TRUE

/* exclude some detections based on various filter options:
 * *** S/N limit
 * *** Acceptable window on detector
 * *** bad photFlags
 */

// the imageID supplied here is the sequence **within this set**
// this value is updated based on the image table later
Stars *FilterStars (Stars *instars, Image *image, unsigned int imageID) {

  int j, N;
  float MTIME, dMs, dMx;
  Stars *stars;
  float RMIN, RMAX, DMIN, DMAX;

  /* correct instrumental mags for exposure time */
  MTIME = (image[0].exptime > 0) ? 2.500*log10(image[0].exptime) : 0.0;

  RMIN = 360.0;
  RMAX =   0.0;
  DMIN = +90.0;
  DMAX = -90.0;

  /* modify resulting star list */
  ALLOCATE (stars, Stars, image[0].nstar);
  for (N = j = 0; j < image[0].nstar; j++) {

    /* allow for some dynamic filtering of star list */
    if (SNLIMIT && instars[j].measure.dM > SNLIMIT) continue;
    if (XMAX && (instars[j].measure.Xccd > XMAX)) continue;
    if (XMIN && (instars[j].measure.Xccd < XMIN)) continue;
    if (YMAX && (instars[j].measure.Yccd > YMAX)) continue;
    if (YMIN && (instars[j].measure.Yccd < YMIN)) continue;
    if (PHOTFLAG_EXCLUDE && (instars[j].measure.photFlags & PHOTFLAG_EXCLUDE)) continue;

    stars[N] = instars[j];

    XY_to_RD (&stars[N].average.R, &stars[N].average.D, stars[N].measure.Xccd, stars[N].measure.Yccd, &image[0].coords);
    stars[N].average.R = ohana_normalize_angle (stars[N].average.R);

    stars[N].measure.photcode = image[0].photcode;

    // determine the full coverage of this set of measurements
    RMIN = MIN (RMIN, stars[N].average.R);
    RMAX = MAX (RMAX, stars[N].average.R);
    DMIN = MIN (DMIN, stars[N].average.D);
    DMAX = MAX (DMAX, stars[N].average.D);

    /** additional quantities to supply to Stars based on the image data **/

    /* calculate accurate per-star airmass and azimuth */
    stars[N].measure.airmass = airmass (image[0].secz, stars[N].average.R, stars[N].average.D, image[0].sidtime, image[0].latitude);
    stars[N].measure.az      = azimuth (image[0].sidtime - stars[N].average.R, stars[N].average.D, image[0].latitude);
    stars[N].measure.Mcal    = image[0].Mcal;
    stars[N].measure.t       = image[0].tzero + 1e-4*stars[N].measure.Yccd*image[0].trate;  /* trate is in 0.1 msec / row */
    stars[N].measure.dt      = MTIME;

    // watch out for any strange values:
    if ((stars[N].measure.M > 25.0) && (stars[N].measure.M < 32.0)) {
      fprintf (stderr, "*");
    }
    // stars->M is either NAN or a valid inst magnitude
    // stars->dM is either NAN or a valid error

    dMs  = 0.0;
    dMx = 0.0;
    if (SUBPIX) {
      dMs =  get_subpix (stars[N].measure.Xccd, stars[N].measure.Yccd);
      dMx = scat_subpix (stars[N].measure.Xccd, stars[N].measure.Yccd);
      if (!isnan(stars[N].measure.dM)) {
	stars[N].measure.dM = hypot (stars[N].measure.dM, dMx);
      }
    }

    if (!isnan(stars[N].measure.M)) {
      stars[N].measure.M   += MTIME - dMs;
    }
    if (!isnan(stars[N].measure.Map)) {
      stars[N].measure.Map += MTIME - dMs;
    }
    
    // the external ID is supplied, but do we trust it?
    if (!EXTERNAL_ID) {
      stars[N].measure.detID = N; // sequence number within image
    }

    if (PSPS_ID) {
      double mjd;
      mjd = ohana_sec_to_mjd (image[0].tzero);
      stars[N].measure.extID = CreatePSPSDetectionID(mjd, image[0].ccdnum, stars[N].measure.detID);
    } else {
      stars[N].measure.extID = 0;
    }

    stars[N].measure.imageID = imageID; // this value is updated in UpdateImageIDs

    N ++;
  }
  image[0].nstar = N;
  REALLOCATE (stars, Stars, image[0].nstar);
  free (instars);

  if (VERBOSE) fprintf (stderr, "read %d stars from target file\n", image[0].nstar);
  if (VERBOSE) fprintf (stderr, "stars cover region %f,%f - %f,%f\n", RMIN, DMIN, RMAX, DMAX);
  return (stars);
}

Stars *MergeStars (Stars *stars, unsigned int *Nstars, Stars *instars, unsigned int Ninstars) {

  int i, j;

  if (stars == NULL) {
    ALLOCATE (stars, Stars, Ninstars);
  } else {
    REALLOCATE (stars, Stars, *Nstars + Ninstars);
  }

  for (j = 0, i = *Nstars; i < *Nstars + Ninstars; i++, j++) {
    stars[i] = instars[j];
  }
  
  *Nstars += Ninstars;
  return (stars);
}
