# include "relastro.h"

void setMrelFinal (Catalog *catalog) {

  int i, j, m, ecode;
  float *p;
  short *q;

  /* if we reset the catalog, reset all the current measurements */
  if (RESET) {

    for (i = 0; i < catalog[0].Naverage; i++) {
      p = (PhotPrimary) ? &catalog[0].average[i].M  : &catalog[0].secfilt[PhotNsec*i+PhotSec].M;
      *p = NO_MAG;
      p = (PhotPrimary) ? &catalog[0].average[i].dM : &catalog[0].secfilt[PhotNsec*i+PhotSec].dM;
      *p = NO_MAG;
      q = (PhotPrimary) ? &catalog[0].average[i].Xm : &catalog[0].secfilt[PhotNsec*i+PhotSec].Xm;
      *q = NO_MAG;

      m = catalog[0].average[i].offset;
      for (j = 0; j < catalog[0].average[i].Nm; j++, m++) {
	
	/* select measurements by photcode */
	ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].photcode);
	if (ecode != photcode[0].code) continue;
	
	/* select measurements by time */
	if (TimeSelect) {
	  if (catalog[0].measure[m].t < TSTART) continue;
	  if (catalog[0].measure[m].t > TSTOP) continue;
	}
	
	catalog[0].measure[m].Mcal = 0;
	catalog[0].measure[m].flags &= 0xff00;
	catalog[0].measure[m].flags &= ~ID_MEAS_POOR;
	catalog[0].measure[m].flags &= ~ID_MEAS_SKIP;
	catalog[0].measure[m].flags &= ~ID_MEAS_AREA;
	catalog[0].measure[m].flags &= ~ID_MEAS_NOCAL;
      }
    }
  }

  setExclusions (catalog, 1);  /* mark by area */

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

  for (i = 0; i < 5; i++) {
    skip_measurements (catalog, i);       /* set ID_MEAS_SKIP for measures to be skipped */
    setMrelOutput  (catalog, 1, FALSE);   /* set Mrel using all measures not skipped */
    clean_measures (catalog, 1, TRUE);    /* mark outliers ID_MEAS_POOR */
    setMrelOutput  (catalog, 1, TRUE);    /* set Mrel using remaining measures */
  }
  setMcalOutput (catalog, 1);

  /* clear ID_STAR_POOR, ID_STAR_FEW, ID_MEAS_NOCAL values before writing ??? */
  for (i = 0; i < catalog[0].Naverage; i++) {
    catalog[0].average[i].code &= ~ID_STAR_FEW;
    catalog[0].average[i].code &= ~ID_STAR_POOR;
    m = catalog[0].average[i].offset;
    for (j = 0; j < catalog[0].average[i].Nm; j++, m++) {
      catalog[0].measure[m].flags &= ~ID_MEAS_NOCAL;
    }
  }
}

/* ID_MEAS_SKIP marks measurements which were not used to calculate Mrel */
void skip_measurements (Catalog *catalog, int pass) {

  int i, k, m, ecode, d1, d2;
  int Ntot, Ntry, Nkeep, Nskip;
  float mag;

  Ntot = Ntry = Nskip = Nkeep = 0;

  /* allow measures from images marked POOR and FEW */
  if (pass >= 3) IMAGE_BAD = ID_IMAGE_NOCAL;
  
  /* allow measures marked as outliers (POOR) and off image region (AREA) */
  if (pass >= 3) MEAS_BAD = ID_MEAS_NOCAL | ID_MEAS_SKIP;

  /* mark measures which should be ignored on second pass */
  for (i = 0; i < catalog[0].Naverage; i++) {
    Ntot += catalog[0].average[i].Nm;
    if (catalog[0].found[i]) continue;

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

      /* clear SKIP for all measures at first */
      catalog[0].measure[m].flags &= ~ID_MEAS_SKIP;

      /** never use these measurements (wrong photcode, bad time range) */
      /* skipped via NOCAL, don't mark as skipped */
      ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].photcode);
      if (ecode != photcode[0].code) continue;

      /* skip measurements by time range */
      if (TimeSelect) {
	if (catalog[0].measure[m].t < TSTART) goto skip;
	if (catalog[0].measure[m].t > TSTOP) goto skip;
      }

      /* skip measurements with sat. dophot values */
      if ((pass < 4) && (catalog[0].measure[m].dophot == 10)) goto skip;

      /* skip measurements from BAD images and mosaics */
      if ((getMcal  (m, 0)) == NO_MAG) goto skip;
      if ((getMmos  (m, 0)) == NO_MAG) goto skip;

      /* skip measurements by inst mag limit */
      if ((pass < 4) && ImagSelect) {
	mag = PhotInst (&catalog[0].measure[m]);
	if (mag < ImagMin) goto skip;
	if (mag > ImagMax) goto skip;
      }

      d1 = (catalog[0].measure[m].dophot == 1);
      d2 = (catalog[0].measure[m].dophot == 2);

      /* skip measurements by measurement error */
      if ((pass < 2) && !(d1 || d2)) goto skip;

      /* if ((pass < 1) && (catalog[0].measure[m].dM > SIGMA_LIM)) goto skip; */
      if ((pass < 1) && !d1) goto skip;
      Nkeep ++;

      continue;

    skip:
      catalog[0].measure[m].flags |= ID_MEAS_SKIP;
      Nskip ++;
    }
  }
  fprintf (stderr, "pass %d, Ntot: %d, Ntry: %d, Nskip: %d, Nkeep: %d\n",
	   pass, Ntot, Ntry, Nskip, Nkeep);
}

/* dophot types:

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

