# include "addstar.h"

Stars *FilterStars (Stars *instars, Image *image) {

  int j, N;
  float MTIME, dMs;
  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);
    while (stars[N].R <    0.0) stars[N].R += 360.0;
    while (stars[N].R >= 360.0) stars[N].R -= 360.0;
    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;

    stars[N].M       = MIN (stars[N].M + MTIME, NO_MAG);
    stars[N].dM      = MIN (stars[N].dM, NO_ERR);

    stars[N].Mgal    = MIN (stars[N].Mgal + MTIME, NO_MAG);
    
    if (SUBPIX) {
      dMs = get_subpix (stars[N].X, stars[N].Y);
      stars[N].M    -= dMs;
      stars[N].Mgal -= dMs;
      stars[N].Map  -= dMs;
      dMs = scat_subpix (stars[N].X, stars[N].Y);
      stars[N].dM = hypot (stars[N].dM, dMs);
    }
    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);
}
