# include "addstar.h"

static int InternalCal;
static int Ncal, NCAL, *Nstar;
static float *Mobs, *dMobs, *Mref, *dMref, *Minst;

void InitCalibration (int mode) {

  InternalCal = mode;

  fprintf (stderr, "calibrating the image...\n");
  Ncal = 0;
  NCAL = 1000;
  ALLOCATE (Mobs,  float,  NCAL);
  ALLOCATE (dMobs, float, NCAL);
  ALLOCATE (Mref,  float,  NCAL);
  ALLOCATE (dMref, float,  NCAL);
  ALLOCATE (Minst, float,  NCAL);
  ALLOCATE (Nstar, int, NCAL);
}
  
void SaveCalibration (float Mo, float dMo, float Mr, float dMr, float Mi, int N) {

  Mobs[Ncal]  = Mo;
  dMobs[Ncal] = dMo;
  Mref[Ncal]  = Mr;
  dMref[Ncal] = dMr;
  Minst[Ncal] = Mi;
  Nstar[Ncal] = N;
  Ncal ++;

  if (Ncal == NCAL) {
    NCAL += 1000;
    REALLOCATE (Mobs,  float, NCAL);
    REALLOCATE (dMobs, float, NCAL);
    REALLOCATE (Mref,  float, NCAL);
    REALLOCATE (dMref, float, NCAL);
    REALLOCATE (Minst, float, NCAL);
    REALLOCATE (Nstar, int,   NCAL);
  }
}

/* use the linked list to navigate the measures; safe for unsorted measures */
void AddToCalibration (Average *average, Measure *measure, Measure *new, int *next, int Nstar) {

  int i, j, m, found0, found1, found2;
  float CalM0, CalM1, CalM2, dCalM;
  float Mcal, color, factor;
  short CalC0, CalC1, CalC2;

  PhotCode *mycode;  // photcode of this measurement
  PhotCode *incode;  // mycode.equiv (internal reference)
  PhotCode *excode;  // incode.equiv (external reference)

  found0 = found1 = found2 = FALSE;
  CalM0 = CalM1 = CalM2 = dCalM = NO_MAG;

  // we have two options here: 
  //  - calibrate to internal system (Mcal)
  //  - calibrate to external system (Mref)

  mycode = GetPhotcodebyCode (new[0].source);
  incode = GetPhotcodebyCode (mycode[0].equiv);
  excode = GetPhotcodebyCode (incode[0].equiv);

  if (InternalCal) {
    CalC0 = incode[0].code;
  } else {
    CalC0 = excode[0].code;
  }
  CalC1 = mycode[0].c1;
  CalC2 = mycode[0].c2;

  m = average[0].offset;
  for (i = 0; i < average[0].Nm; i++) {
    if (measure[m].source == CalC0) { 
      found0 = TRUE; 
      CalM0  = measure[m].M_PS; 
      dCalM  = measure[m].dM_PS; 
    }
    if (measure[m].source == CalC1) { 
      found1 = TRUE; 
      CalM1  = measure[m].M_PS; 
    }
    if (measure[m].source == CalC2) { 
      found2 = TRUE; 
      CalM2  = measure[m].M_PS; 
    }
    if (found0 && found1 && found2) {
      Mcal   = new[0].M_PS + 0.001*mycode[0].C + mycode[0].K*(new[0].airmass_PS - 1.0) - ZeroPt;
      color  = CalM1 - CalM2 - 0.001*mycode[0].dX;
      factor = color;
      for (j = 0; j < mycode[0].Nc; j++) {
	Mcal += mycode[0].X[j]*factor;
	factor *= color;
      }
      if (!InternalCal) {
	Mcal += 0.001*incode[0].C;
      }
      // if we want to apply a Mcal -> Mref color correction, we need the additional color term
      SaveCalibration (Mcal, new[0].dM_PS, CalM0, dCalM, new[0].M_PS - ZeroPt - new[0].dt_PS, Nstar);
      return;
    }
    m = next[m];
  }
  return;
}

void FindCalibration (Image *image) {

  int i, MaxN, *Nlist, Nkeep;
  float N, M1, M2, *Dmag, *dDmag;
  float dMo, dMr, Mw, Dmed, W1, W2, NSigma;

  /* reject multiple matched-stars */
  /* find maximum value of Nstar[] */
  MaxN = -1;
  for (i = 0; i < Ncal; i++) {
    MaxN = MAX (Nstar[i], MaxN);
  }
  if (MaxN == -1) {
    fprintf (stderr, "no clean stars\n");
    image[0].Mcal_PS = 10.000;
    image[0].dMcal_PS = 10.000;
    return;
  }
  /* create a hash array from Nstar[] entries */
  ALLOCATE (Nlist, int, MaxN + 1);
  bzero (Nlist, MAX (0, MaxN*sizeof(int)));
  for (i = 0; i < Ncal; i++) {
    Nlist[Nstar[i]] ++;
  }
  
  /* accumulate delta mags */
  ALLOCATE (Dmag, float, Ncal);
  ALLOCATE (dDmag, float, Ncal);
  Nkeep = 0;
  for (i = 0; i < Ncal; i++) {
    /* if this entry has too many (or two few?) matches, skip it */
    if (Nlist[Nstar[i]] != 1) continue;

    /* clip by instrumental magnitude */
    if (Minst[i] > CAL_INSTMAG_MAX) continue;
    if (Minst[i] < CAL_INSTMAG_MIN) continue;
    
    /* XXX EAM: note the artificial 0.005 dmag here */
    dMr = MAX (0.005, dMref[i]);
    dMo = MAX (0.005, dMobs[i]);

    Dmag[Nkeep] = (Mobs[i] - Mref[i]);
    dDmag[Nkeep] = (dMr*dMr + dMo*dMo);
    Nkeep ++;
  }

  if (Nkeep < 5) {
    fprintf (stderr, "too few stars\n");
    image[0].Mcal_PS = 10.000;
    image[0].dMcal_PS = 10.000;
    return;
  }
  fsort2 (Dmag, dDmag, Nkeep);

  /* take sort list of Dmag, find median */
  Dmed = Dmag[(int)(0.5*Nkeep)];

  /* exclude points with abs(Dmag - Dmed) / dDmag > 2.5 */

  /* accumulate delta mags (25% - 75% of clipped range) */
  W1 = 0.0;
  W2 = 0.0;
  M1 = 0.0;
  M2 = 0.0;
  N  = 0.0;
  for (i = 0; i < Nkeep; i++) {
    NSigma = fabs (Dmag[i] - Dmed) / sqrt (dDmag[i]);
    if (NSigma > 2.5) continue;
    W1 += Dmag[i] / dDmag[i];
    W2 += 1 / dDmag[i];
    M1 += Dmag[i];
    M2 += SQ (Dmag[i]);
    N  += 1.0; 
  }

  if (N > 1) {
    M1 = M1 / N;
    M2 = sqrt (fabs(M2/N - M1*M1));
    Mw = W1 / W2;
    fprintf (stderr, "N: %.0f, mean: %f, wt mean: %f, stdev: %f, precision: %f\n", N, M1, Mw, M2, M2 / sqrt (N));
    image[0].Mcal_PS = M1;
    image[0].dMcal_PS = M2 / sqrt (N);
    image[0].Mxxxx = N;
  } else {
    fprintf (stderr, "too few stars\n");
    image[0].Mcal_PS = 10.000;
    image[0].dMcal_PS = 10.000;
    image[0].Mxxxx = 0;
  }
}
