# include "absphot.h"

int subset_catalog (Catalog *subcatalog, Catalog *catalog) {
  
  // XXX define this : obs - ref vs color(obs)
  // target is the ref photcode (derived from obs photcode table entry)
  PhotCode *target; 

  int i, j, offset, averef, ecode;
  int NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm;
  float mag;

  /* we are moving only the subset of measurements from catalog[0] to subcatalog[0] */
  if (subcatalog[0].Naverage == 0) {
    Naverage = 0;
    Nmeasure = 0;
    NAVERAGE = 50;
    NMEASURE = 1000;
    ALLOCATE (subcatalog[0].average, Average, NAVERAGE);
    ALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*PhotNsec);
    ALLOCATE (subcatalog[0].measure, Measure, NMEASURE);
  } else {
    Naverage = subcatalog[0].Naverage;
    Nmeasure = subcatalog[0].Nmeasure;
    NAVERAGE = Naverage + 50;
    NMEASURE = Nmeasure + 1000;
    ALLOCATE (subcatalog[0].average, Average, NAVERAGE);
    ALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*PhotNsec);
    ALLOCATE (subcatalog[0].measure, Measure, NMEASURE);
  }

  /* exclude stars not in range or with too few measurements */
  for (i = 0; i < catalog[0].Nmeasure; i++) {

    // find desired target photcode measurements
    if (catalog[0].measure[i].photcode != target[0].code) continue;

    // some possible exclusions based on the target measurements
    // mag = PhotCat (&catalog[0].measure[i]);
    // if (mag > MAG_LIM) continue;

    /* select measurements by measurement error */
    // if ((SIGMA_LIM > 0) && (catalog[0].measure[i].dM > SIGMA_LIM)) continue;

    averef = catalog[0].measure[i].averef;

    // copy the average and secfilt data for this average object
    // XXX potentially skip this source on some basis here
    subcatalog[0].average[Naverage] = catalog[0].average[averef];
    subcatalog[0].average[Naverage].offset = Nmeasure;
    for (j = 0; j < PhotNsec; j++) {
      subcatalog[0].secfilt[PhotNsec*Naverage+j] = catalog[0].secfilt[PhotNsec*averef+j];
    }

    // copy the target photcode measurement
    // XXX how do we manage two target photcode measurements for the same star?
    // looks like we create two output average stars, keeping the average 
    // quantities for both.

    subcatalog[0].average[Naverage].Nm = 1;

    subcatalog[0].measure[Nmeasure] = catalog[0].measure[i];
    subcatalog[0].measure[Nmeasure].averef = Naverage;

    Nmeasure ++;
    if (Nmeasure == NMEASURE) {
      NMEASURE += 1000;
      REALLOCATE (subcatalog[0].measure, Measure, NMEASURE);
    }

    Naverage ++;
    if (Naverage == NAVERAGE) {
      NAVERAGE += 50;
      REALLOCATE (subcatalog[0].average, Average, NAVERAGE);
      REALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*PhotNsec);
    }
  }
  REALLOCATE (subcatalog[0].average, Average, MAX (Naverage, 1));
  REALLOCATE (subcatalog[0].measure, Measure, MAX (Nmeasure, 1));
  REALLOCATE (subcatalog[0].secfilt, SecFilt, PhotNsec*MAX (Naverage, 1));
  subcatalog[0].Naverage = Naverage;
  subcatalog[0].Nmeasure = Nmeasure;
  subcatalog[0].Nsecfilt = catalog[0].Nsecfilt;

  if (VERBOSE) {
    fprintf (stderr, "%d: using %d stars (%d measures) for catalog\n", i, 
	     subcatalog[0].Naverage, subcatalog[0].Nmeasure);
  }
  return (TRUE);
}
