# include "addstar.h"

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].dM > SNLIMIT) continue;
    if (XMAX && (instars[j].X > XMAX)) continue;
    if (XMIN && (instars[j].X < XMIN)) continue;
    if (YMAX && (instars[j].Y > YMAX)) continue;
    if (YMIN && (instars[j].Y < YMIN)) continue;
    stars[N] = instars[j];

    XY_to_RD (&stars[N].R, &stars[N].D, stars[N].X, stars[N].Y, &image[0].coords);
    stars[N].R = ohana_normalize_angle (stars[N].R);
    stars[N].found = -1;
    stars[N].code = image[0].photcode;

    RMIN = MIN (RMIN, stars[N].R);
    RMAX = MAX (RMAX, stars[N].R);
    DMIN = MIN (DMIN, stars[N].D);
    DMAX = MAX (DMAX, stars[N].D);
    /** additional quantities to supply to Stars based on the image data **/

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

    if ((stars[N].M > 25.0) && (stars[N].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].X, stars[N].Y);
      dMx = scat_subpix (stars[N].X, stars[N].Y);
      if (!isnan(stars[N].dM)) {
	stars[N].dM = hypot (stars[N].dM, dMx);
      }
    }

    if (!isnan(stars[N].M)) {
      stars[N].M += MTIME - dMs;
    }
    if (!isnan(stars[N].Map)) {
      stars[N].Map += MTIME - dMs;
    }
    
    // XXX currently, this ID is internal only; 
    // should we use the psphot / other external ID, if available?
    stars[N].detID = N; // sequence number within image
    stars[N].imageID = imageID; // does this need to be updated?

    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, int *Nstars, Stars *instars, 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);
}
