# include "relphot.h"

// we've just reloaded the data from disk; we now need to apply the Image/Mosaic/Grid
// calibrations determined by the rest of the program.  We also need to set the final
// output dbFlags values

void setMrelFinal (Catalog *catalog, FlatCorrectionTable *flatcorr, int simpleAverage) {

  off_t i;
  int ecode;

  int Nsecfilt = GetPhotcodeNsecfilt ();

  /* if we reset the catalog, reset all the current measurements */
  if (RESET) {
    // flags used by the photometry analysis (excluding UBERCAL)
    unsigned int PHOTOM_FLAGS = 
      ID_MEAS_NOCAL |          // detection ignored for this analysis (photcode, time range)
      ID_MEAS_POOR_PHOTOM |    // detection is photometry outlier 
      ID_MEAS_SKIP_PHOTOM |    // detection was ignored for photometry measurement 
      ID_MEAS_AREA |	       // detetion was outside acceptable area of device
      ID_MEAS_SYNTH_MAG |      // magnitude is synthetic
      ID_MEAS_STACK_PRIMARY |  // this stack measurement is in the primary skycell
      ID_MEAS_STACK_PHOT_SRC;  // this measurement supplied the stack photometry

    // ID_MEAS_PHOTOM_UBERCAL -- externally-supplied zero point from ubercal analysis
    // this is set by 'setphot', do not reset here

    // flags used by the photometry analysis (excluding UBERCAL)
    // unsigned int secfiltFlags = 
    //   ID_PHOTOM_PASS_0 | // average measured at pass 0
    //   ID_PHOTOM_PASS_1 | // average measured at pass 1
    //   ID_PHOTOM_PASS_2 | // average measured at pass 2
    //   ID_PHOTOM_PASS_3 | // average measured at pass 3
    //   ID_PHOTOM_PASS_4 | // average measured at pass 3
    //   ID_SECF_USE_SYNTH | // average measured at pass 3
    //   ID_SECF_USE_UBERCAL | // average measured at pass 3
    //   ID_SECF_OBJ_EXT; // average measured at pass 3
    
    int Ns;
    for (Ns = 0; Ns < Nphotcodes; Ns++) {
	
      int thisCode = photcodes[Ns][0].code;
      int Nsec = GetPhotcodeNsec(thisCode);

      for (i = 0; i < catalog[0].Naverage; i++) {
	off_t N = Nsecfilt*i+Nsec;
	dvo_secfilt_init (&catalog[0].secfilt[N]);

	off_t m = catalog[0].average[i].measureOffset;
	off_t j;
	for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) {
	
	  // skip measurements that do not match the current photcode
	  ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].photcode);
	  if (ecode != thisCode) { continue; }

	  /* select measurements by time */
	  if (TimeSelect) {
	    if (catalog[0].measure[m].t < TSTART) continue;
	    if (catalog[0].measure[m].t > TSTOP) continue;
	  }
	
	  // only reset Mcal for measures with a matching image
	  // do not reset Mcal for ubercal images unless explicitly requested

	  if (catalog[0].measure[m].dbFlags & ID_MEAS_PHOTOM_UBERCAL) {
	    if (!KEEP_UBERCAL) {
	      catalog[0].measure[m].Mcal = 0.0;
	      catalog[0].measure[m].dbFlags &= ~ID_MEAS_PHOTOM_UBERCAL;
	    } 
	  } else {
	    if (RESET_ZEROPTS && (getImageEntry (m, 0) >= 0)) {
	      catalog[0].measure[m].Mcal = 0.0;
	    }
	  }
	  catalog[0].measure[m].dbFlags &= ~PHOTOM_FLAGS;
	}
      }
    }
  }

  // this sets flags in the measureT element, not the measure element
  setExclusions (catalog, 1, VERBOSE);  /* mark by area */

  /* set catalog[0].found[i] = FALSE */
  ALLOCATE (catalog[0].found_t, off_t, MAX (1, Nsecfilt*catalog[0].Naverage));
  ALLOCATE (catalog[0].foundWarp_t, off_t, MAX (1, Nsecfilt*catalog[0].Naverage));
  for (i = 0; i < Nsecfilt*catalog[0].Naverage; i++) {
    catalog[0].found_t[i] = FALSE;
    catalog[0].foundWarp_t[i] = FALSE;
  }

  for (i = 0; i < catalog[0].Naverage; i++) {
    catalog[0].average[i].psfQF     = NAN;	// force recalculation below
    catalog[0].average[i].psfQFperf = NAN;	// force recalculation below
    catalog[0].average[i].stargal   = NAN;	// force recalculation below
  }

  // XXX make this optional? (do not clean for -averages?)
  if (!simpleAverage) clean_measures (catalog, 1, TRUE, flatcorr);    /* mark outliers ID_MEAS_POOR_PHOTOM */
  for (i = 0; i < 5; i++) {
    skip_measurements (catalog, i, flatcorr);       /* set ID_MEAS_SKIP for measures to be skipped */
    setMrelOutput  (catalog, 1, i, flatcorr);    /* set Mrel using remaining measures */
  }
  if (!simpleAverage) setMcalOutput (catalog, 1, flatcorr);

  /* clear ID_STAR_POOR, ID_STAR_FEW values before writing ??? */
  /* ID_MEAS_NOCAL is an internal bit, so it should be cleared */
  // XXX eventually: move these internal flags to an internal-only catalog->FOO element (like found)
  for (i = 0; i < catalog[0].Naverage; i++) {
    catalog[0].average[i].flags &= ~ID_STAR_FEW;
    catalog[0].average[i].flags &= ~ID_STAR_POOR;
    off_t j;
    for (j = 0; j < Nsecfilt; j++) {
	catalog[0].secfilt[i*Nsecfilt+j].flags &= ~ID_STAR_FEW; 
	catalog[0].secfilt[i*Nsecfilt+j].flags &= ~ID_STAR_POOR; 
    }
    off_t m = catalog[0].average[i].measureOffset;
    for (j = 0; j < catalog[0].average[i].Nmeasure; j++, m++) {
      catalog[0].measure[m].dbFlags &= ~ID_MEAS_NOCAL;
    }
  }
}

/* ID_MEAS_SKIP marks measurements which were not used to calculate Mrel */
// For each star & photcode, we are selecting the measurements to use in the average
// photometry value.  This function is called several times, lowering the bar to
// acceptance on each pass. The rejections at the different passes are:

// 0 : only use measurements thought to be GOOD (photflags not POOR or BAD)
// 1 : accept measurements thought to be POOR based on photflags
// 2 : accept the measurements marked as outliers; accept images marked as outliers
// 3 : accept measurements thought to be BAD based on photflags
// 4 : accept the measurements outside of the instrumental magnitude limits (eg, SAT)

// XXX if I'm reading this correctly, passes 0, 1, 2 are identical & accept all
// non-outliers with Imag range

// int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure) {
// 
//   off_t k;
// 
//   int Nsecfilt = GetPhotcodeNsecfilt ();
// 
//   off_t m = average[0].measureOffset;
// 
//   for (k = 0; k < average[0].Nmeasure; k++, m++) {
//     fprintf (stderr, "meas: %08x\n", measure[m].dbFlags);
//   }
// 
//   int Ns;
//   for (Ns = 0; Ns < Nsecfilt; Ns++) {
//     fprintf (stderr, "secf: %08x\n", secfilt[Ns].flags);
//   }
//   return 1;
// }

int print_measure_set (Average *average, SecFilt *secfilt, Measure *measure);

void skip_measurements (Catalog *catalog, int pass, FlatCorrectionTable *flatcorr) {

  off_t i, k, m;
  off_t Ntot, Ntry, Nkeep, Nskip;
  float mag;

  int Nsecfilt = GetPhotcodeNsecfilt ();

  Ntot = Ntry = Nskip = Nkeep = 0;

  /* allow measures from images marked POOR and FEW */
  if (pass >= 3) IMAGE_BAD = ID_IMAGE_PHOTOM_NOCAL;
  
  /* allow measures marked as outliers (POOR) and off image region (AREA) */
  if (pass >= 3) {
    MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_SKIP_PHOTOM;
  } else {
    MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_PHOTOM | ID_MEAS_SKIP_PHOTOM | ID_MEAS_AREA;
  }

  /* mark measures which should be ignored on this pass */
  for (i = 0; i < catalog[0].Naverage; i++) {
    Ntot += catalog[0].average[i].Nmeasure;

    if (FALSE && (catalog[0].average[i].objID == 0x46a4) && (catalog[0].average[i].catID == 0xf40e)) {
      fprintf (stderr, "test obj\n");
      print_measure_set (&catalog[0].average[i], &catalog[0].secfilt[i*Nsecfilt], catalog[0].measure);
    }

    // mark measurements for each secfilt separately
    int Ns;
    for (Ns = 0; Ns < Nphotcodes; Ns++) {
      
      int thisCode = photcodes[Ns][0].code;
      int Nsec = GetPhotcodeNsec(thisCode);
      
      /* star/photcodes already calibrated */
      if (catalog[0].found_t[Nsecfilt*i+Nsec]) continue;  

      m = catalog[0].average[i].measureOffset;
      for (k = 0; k < catalog[0].average[i].Nmeasure; k++, m++) {
	Ntry++;

	// skip measurements not related to this photcode
	PhotCode *code = GetPhotcodebyCode (catalog[0].measure[m].photcode);
	if (!code) continue;
	if (code->equiv != thisCode) continue;

	// skip measurements by time range (mark as skipped for this secfilt)
	// XXX note that this is a bit dangerous : some objects may never get calibrated 
	if (TimeSelect) {
	  if (catalog[0].measure[m].t < TSTART) goto skip;
	  if (catalog[0].measure[m].t > TSTOP) goto skip;
	}
	
	// we now have a measurement relavant to this photcode (and time range)

	// clear SKIP for all measures at first 
	catalog[0].measureT[m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;
	catalog[0].measure [m].dbFlags &= ~ID_MEAS_SKIP_PHOTOM;

	// skip measurements from BAD images and mosaics (not REF mags)
	// do NOT skip measurements without a matching image (REF mags)
	off_t Nim = getImageEntry (m, 0);
	if (Nim > -1) {
	  if (isnan(getMcal (m, 0, flatcorr, catalog))) goto skip;
	  if (isnan(getMmos (m, 0))) goto skip;
	}

	// PASS 4 : skip measurements by inst mag limit (not REF mags)
	if ((pass < 4) && ImagSelect) {
	  if (Nim > -1) {
	    mag = PhotInst (&catalog[0].measure[m], MAG_CLASS_PSF);
	    if (mag < ImagMin) goto skip;
	    if (mag > ImagMax) goto skip;
	  }
	}

	// PASS 3 : accept bad measurements (eg, SAT, CR), internal outliers
	if (pass < 3) {
	  if (catalog[0].measure[m].photFlags & code->photomBadMask) goto skip;
	  if ((catalog[0].measure[m].photcode > 10000) && (catalog[0].measure[m].photcode < 10500)) {
	    if (catalog[0].measure[m].psfQF < 0.85) goto skip;
	  }
	}
	
	// PASS 2 : (no additional cut)

	// PASS 1 : accept poor measurements as well (eg, POOR FIT, etc)
	if ((pass < 1) && (catalog[0].measure[m].photFlags & code->photomPoorMask)) goto skip;

	// PASS 0 : skip poor, outlier, bad measurements as well (eg, POOR FIT, etc)

	Nkeep ++;
	
	continue;

      skip:
	catalog[0].measure [m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
	catalog[0].measureT[m].dbFlags |= ID_MEAS_SKIP_PHOTOM;
	Nskip ++;
      }
    }
  }
  if (VERBOSE) fprintf (stderr, "pass %d, Ntot: "OFF_T_FMT", Ntry: "OFF_T_FMT", Nskip: "OFF_T_FMT", Nkeep: "OFF_T_FMT"\n",
			pass, Ntot, Ntry, Nskip, Nkeep);
}

/* dophot types:

1 - star
2 - galaxy
3 - star
4 - weak
5 - edge
7 - weak
9 - weak
10 - sat
*/

