# include "photdbc.h"

// copy a catalog to a new subcatalog, applying some filters
// the supplied subcatalog must already be locked, opened, and created
int make_subcatalog (Catalog *subcatalog, Catalog *catalog) {
  
  int i, j, offset;
  int NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm, Nsecfilt;
  double mag, minMag;

  Nsecfilt = GetPhotcodeNsecfilt ();

  /* we are moving only the subset of measurements from catalog[0] to subcatalog[0] */
  NAVERAGE = 50;
  NMEASURE = 1000;
  Nmeasure = Naverage = 0;
  REALLOCATE (subcatalog[0].average, Average, NAVERAGE);
  REALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt);
  REALLOCATE (subcatalog[0].measure, Measure, NMEASURE);

  for (i = 0; i < catalog[0].Naverage; i++) {
    // perform exclusions based on average properties

    // exclude stars with too few measurements
    if (NMEAS_MIN && (catalog[0].average[i].Nm < NMEAS_MIN)) continue; 

    /* assign average and secfilt values */
    subcatalog[0].average[Naverage] = catalog[0].average[i];
    subcatalog[0].average[Naverage].offset = Nmeasure;
    for (j = 0; j < Nsecfilt; j++) {
      subcatalog[0].secfilt[Nsecfilt*Naverage+j] = catalog[0].secfilt[Nsecfilt*i+j];
    }

    minMag = 32;
    Nm = 0;
    for (j = 0; j < catalog[0].average[i].Nm; j++) {

      offset = catalog[0].average[i].offset + j;

      # if 0
      if (DropPhotcode) {
	ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[offset].source);
	if (ecode == photcode[0].code) continue;
      }
      if (DropNonStellar && catalog[0].measure[offset].dophot != 1) continue;
      # endif

      /* exclude measurements by measurement error */
      if (SIGMA_MAX && (catalog[0].measure[offset].dM_PS > SIGMA_MAX)) continue;

      /* select measurements by mag limit */
      if (ExcludeByInstMag) {
	mag = PhotInst (&catalog[0].measure[offset]);
	if (mag < INST_MAG_MIN) continue;
	if (mag > INST_MAG_MAX) continue;
      }

      if (ExcludeByMaxMinMag) {
	mag = PhotSys (&catalog[0].measure[offset], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]);
	minMag = MIN (minMag, mag);
      }

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

    // exclude faint objects
    if (ExcludeByMaxMinMag && (minMag > MAX_MIN_MAG)) {
      Nmeasure -= Nm;
      continue; 
    }

    // after measurement exclusion, exclude stars with too few measurements
    if (Nm < NMEAS_MIN) {
      Nmeasure -= Nm;
      continue; 
    }
    subcatalog[0].average[Naverage].Nn = 0;
    subcatalog[0].average[Naverage].Nm = Nm;
    Naverage ++;
    if (Naverage == NAVERAGE) {
      NAVERAGE += 50;
      REALLOCATE (subcatalog[0].average, Average, NAVERAGE);
      REALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt);
    }
  }
  REALLOCATE (subcatalog[0].average, Average, MAX (Naverage, 1));
  REALLOCATE (subcatalog[0].measure, Measure, MAX (Nmeasure, 1));
  REALLOCATE (subcatalog[0].secfilt, SecFilt, Nsecfilt*MAX (Naverage, 1));
  subcatalog[0].Naverage = Naverage;
  subcatalog[0].Nmeasure = Nmeasure;
  subcatalog[0].Nsecfilt = catalog[0].Nsecfilt;

  // XXX for now, don't copy the missing entries (these should be re-computed)
  ALLOCATE (subcatalog[0].missing, Missing, 1);
  subcatalog[0].Nmissing = 0;

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