# include "getstar.h"

/* add selected catalog objects to the output catalog */
int select_by_region (Catalog *output, Catalog *catalog, SkyRegion *region, int start, int end) {

  int i, j, n, Nm, offset, m, Nsecfilt, code, Nsec;
  int Nave, NAVE, Nmeas, NMEAS;
  double R, D;
  float mag;

  Nsecfilt = output[0].Nsecfilt;

  if (output == NULL) Shutdown ("output not defined");

  /* identify selection criteria */
  if (VERBOSE) fprintf (stderr, "extracting for time range %d to %d\n", start, end);
  if (VERBOSE) fprintf (stderr, "extracting for region %f,%f to %f,%f\n", 
			region[0].Rmin, region[0].Dmin, region[0].Rmax, region[0].Dmax);

  if (output[0].header.buffer != NULL) {
    gfits_modify (&output[0].header, "RA0",  "%lf", 1, region[0].Rmin);
    gfits_modify (&output[0].header, "DEC0", "%lf", 1, region[0].Dmin);
    gfits_modify (&output[0].header, "RA1",  "%lf", 1, region[0].Rmax);
    gfits_modify (&output[0].header, "DEC1", "%lf", 1, region[0].Dmax);
  }

  code = photcode[0].code;
  Nsec = GetPhotcodeNsec (code);

  /* select the average objects in this region */
  Nave = output[0].Naverage;
  NAVE = output[0].Naverage + 1000;
  REALLOCATE (output[0].average, Average, NAVE);
  REALLOCATE (output[0].secfilt, SecFilt, NAVE*Nsecfilt);

  Nmeas = output[0].Nmeasure;
  NMEAS = output[0].Nmeasure + 1000;
  REALLOCATE (output[0].measure, Measure, NMEAS);

  for (i = 0; i < catalog[0].Naverage; i++) {
    n = catalog[0].measure[i].averef;

    R = catalog[0].average[i].R;
    D = catalog[0].average[i].D;
    
    if (R < region[0].Rmin) continue;
    if (R > region[0].Rmax) continue;
    if (D < region[0].Dmin) continue;
    if (D > region[0].Dmax) continue;

    if (MagLimitUse) {
      mag = NAN;
      if (Nsec != -1) {
	mag = catalog[0].secfilt[i*Nsecfilt + Nsec].M;
      } else {
	offset = catalog[0].average[i].measureOffset;
	for (m = 0; m < catalog[0].average[i].Nmeasure; m++) {
	  if (catalog[0].measure[offset + m].photcode == code) {
	    mag = PhotRel (&catalog[0].measure[offset + m], &catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt]);
	    break;
	  }
	}
      }
      if (isnan(mag) || (mag > MagLimitValue)) continue;
    }

    output[0].average[Nave] = catalog[0].average[i];
    output[0].average[Nave].measureOffset = Nmeas;
    for (j = 0; j < Nsecfilt; j++) {
      output[0].secfilt[Nsecfilt*Nave + j] = catalog[0].secfilt[Nsecfilt*i + j];
    }

    Nm = 0;
    offset = catalog[0].average[i].measureOffset;

    for (j = 0; j < catalog[0].average[i].Nmeasure; j++) {
      output[0].measure[Nmeas] = catalog[0].measure[offset + j];
      output[0].measure[Nmeas].averef = Nave;
      Nmeas ++;
      
      CHECK_REALLOCATE (output[0].measure, Measure, NMEAS, Nmeas, 1000);
    }      

    Nave ++;
    if (Nave == NAVE) {
      NAVE += 1000;
      REALLOCATE (output[0].average, Average, NAVE);
      REALLOCATE (output[0].secfilt, SecFilt, NAVE*Nsecfilt);
    }
  }
  output[0].Naverage = Nave;
  output[0].Nmeasure = Nmeas;
  output[0].Nsecf_mem = Nave*Nsecfilt;

  fprintf (stderr, "output catalog has %d stars (%d measures, %d secfilt)\n",
	   output[0].Naverage, output[0].Nmeasure, output[0].Nsecfilt);
  return (TRUE);
}
