# include "relphot.h"

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

  /* 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*PhotNsec);
  ALLOCATE (subcatalog[0].measure, Measure, NMEASURE);
  Nmeasure = Naverage = 0;

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

    /* XXX this limitation is absurd
    if (catalog[0].average[i].R < fullregion[0].RA[0]) continue;
    if (catalog[0].average[i].R > fullregion[0].RA[1]) continue;
    if (catalog[0].average[i].D < fullregion[0].DEC[0]) continue;
    if (catalog[0].average[i].D > fullregion[0].DEC[1]) continue;
    */

    /* start with all stars good */
    subcatalog[0].average[Naverage] = catalog[0].average[i];
    subcatalog[0].average[Naverage].offset = Nmeasure;
    for (j = 0; j < PhotNsec; j++) {
      subcatalog[0].secfilt[PhotNsec*Naverage+j] = catalog[0].secfilt[PhotNsec*i+j];
    }

    if (RESET) {
      float *p;
      p = (PhotPrimary) ? &subcatalog[0].average[Naverage].M : &subcatalog[0].secfilt[PhotNsec*Naverage+PhotSec].M_PS;
      *p = NO_MAG;
      p = (PhotPrimary) ? &subcatalog[0].average[Naverage].dM : &subcatalog[0].secfilt[PhotNsec*Naverage+PhotSec].dM_PS;
      *p = NO_MAG;
      subcatalog[0].average[Naverage].code &= ~ID_STAR_FEW;
      subcatalog[0].average[Naverage].code &= ~ID_STAR_POOR;
    }

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

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

      /* select measurements by photcode */
      ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[offset].source);
      if (ecode != photcode[0].code) continue;

      /* select measurements by time */
      if (TimeSelect) {
	if (catalog[0].measure[offset].t < TSTART) continue;
	if (catalog[0].measure[offset].t > TSTOP) continue;
      }

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

      /* select measurements by mag limit */
      mag = PhotCat (&catalog[0].measure[offset]);
      if (mag > MAG_LIM) continue;

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

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

      subcatalog[0].measure[Nmeasure].flags &= ~ID_MEAS_SKIP;
      subcatalog[0].measure[Nmeasure]        = catalog[0].measure[offset];
      subcatalog[0].measure[Nmeasure].averef = Naverage;
      if (RESET) { 
	subcatalog[0].measure[Nmeasure].Mcal_PS = 0;
	subcatalog[0].measure[Nmeasure].flags &= 0xff00;
	subcatalog[0].measure[Nmeasure].flags &= ~ID_MEAS_POOR;
	subcatalog[0].measure[Nmeasure].flags &= ~ID_MEAS_AREA;
	subcatalog[0].measure[Nmeasure].flags &= ~ID_MEAS_NOCAL;
      }
      Nmeasure ++;
      Nm ++;
      if (Nmeasure == NMEASURE) {
	NMEASURE += 1000;
	REALLOCATE (subcatalog[0].measure, Measure, NMEASURE);
      }
    }
    if (Nm < 2) { /* enough measurements in band? */
      Nmeasure -= Nm;
      continue; 
    }
    subcatalog[0].average[Naverage].Nm = Nm;
    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 (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 (TRUE);
}
