# include "photdbc.h"

Catalog *subcatalog (Catalog *catalog) {
  
  int i, j, offset, ecode;
  int NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm, Nsecfilt;
  int mag;
  Catalog *subcatalog;

  ALLOCATE (subcatalog, Catalog, 1);

  Nsecfilt = GetPhotcodeNsecfilt ();

  subcatalog[0] = catalog[0];

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

  for (i = 0; i < catalog[0].Naverage; i++) {
    /* exclude stars with too few measurements */
    if (catalog[0].average[i].Nm < NMEAS_MIN) continue; 

    /* exclude stars not in range */
    if (catalog[0].average[i].R < RMIN) continue;
    if (catalog[0].average[i].R > RMAX) continue;
    if (catalog[0].average[i].D < DMIN) continue;
    if (catalog[0].average[i].D > DMAX) 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];
    }

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

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

      /* exclude measurements by photcode 
      if (DropPhotcode) {
	ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[offset].source);
	if (ecode == photcode[0].code) continue;
      }*/

      /* select measurements by quality 
      if (DropNonStellar && catalog[0].measure[offset].dophot != 1) continue;
      */

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

      /* select measurements by mag limit */
      mag = iPhotInst (&catalog[0].measure[offset]);
      if (mag < MMIN) continue;
      if (mag > MMAX) continue;

      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);
      }
    }
    if (Nm < NMEAS_MIN) { /* enough measurements meet criteria? */
      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;

  ALLOCATE (subcatalog[0].missing, Missing, 1);
  subcatalog[0].Nmissing = 0;

  if (catalog[0].Naverage > 0) free (catalog[0].average);
  if (catalog[0].Nmeasure > 0) free (catalog[0].measure);
  if (catalog[0].Naverage > 0) free (catalog[0].secfilt);

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

}
