# include "addstar.h"

static int Ncal, NCAL;
static float *Mobs, *dMobs, *Mref, *dMref, *Cref, *Airm;
static float *RA, *DEC, *X, *Y;
int *Nstar;

void InitCalibration () {

    fprintf (stderr, "calibrating the image...  %d\n", CalReference);
    Ncal = 0;
    NCAL = 1000;
    ALLOCATE (Mobs,  float,  NCAL);
    ALLOCATE (dMobs, float, NCAL);
    ALLOCATE (Mref,  float,  NCAL);
    ALLOCATE (dMref,  float,  NCAL);
    ALLOCATE (Cref,  float,  NCAL);
    ALLOCATE (Airm,  float,  NCAL);
    ALLOCATE (Nstar, int, NCAL);
    ALLOCATE (RA, float, NCAL);
    ALLOCATE (DEC, float, NCAL);
    ALLOCATE (X, float, NCAL);
    ALLOCATE (Y, float, NCAL);
}
  
void SaveCalibration (float M, float dM, float Mr, float dMr, float Mc, float A, int N, float ra, float dec, float x, float y) {

  Mobs[Ncal] = M;
  dMobs[Ncal] = dM;
  Mref[Ncal] = Mr;
  dMref[Ncal] = dMr;
  Cref[Ncal] = Mc;
  Airm[Ncal] = A;
  Nstar[Ncal] = N;

  RA[Ncal] = ra;
  DEC[Ncal] = dec;
  X[Ncal] = x;
  Y[Ncal] = y;
  Ncal ++;

  if (Ncal == NCAL) {
    NCAL += 1000;
    REALLOCATE (dMobs, float, NCAL);
    REALLOCATE (Mobs,  float, NCAL);
    REALLOCATE (Mref,  float, NCAL);
    REALLOCATE (dMref,  float, NCAL);
    REALLOCATE (Cref,  float, NCAL);
    REALLOCATE (Airm,  float, NCAL);
    REALLOCATE (Nstar, int,   NCAL);
    REALLOCATE (RA, float, NCAL);
    REALLOCATE (DEC, float, NCAL);
    REALLOCATE (X, float, NCAL);
    REALLOCATE (Y, float, NCAL);
  }
}

void FindCalibration (Image *image) {

  int i, MaxN, *Nlist, Nkeep;
  float N, M1, M2, Klam, Clam, Xlam, Mabs, *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 = 10000;
    image[0].dMcal = 10000;
    return;
  }
  /* create a hash array from Nstar[] entries */
  ALLOCATE (Nlist, int, MAX (1, 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;
  Klam = thiscode[0].K;
  Clam = thiscode[0].C;
  Xlam = thiscode[0].X[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;
    Mabs = 0.001*(Mobs[i] + Klam*(Airm[i] - 1000) + Clam + Xlam*(Mref[i] - Cref[i])) - ZeroPt;

    /* note: subpix correction is applied in gstars */

    if (DUMP_MATCHES) 
      fprintf (stdout, "%d  %6.3f %6.3f %6.3f  %10.6f %10.6f  %7.2f %7.2f  %7.2f  %7.2f\n", 
	       i, Mabs, 0.001*Mref[i], 0.001*Cref[i], RA[i], DEC[i], X[i], Y[i], Airm[i], 0.001*dMobs[i]);

    /* skip stars brighter than 8.0 */
    if (Mabs > 9.0) continue;
    if (Mabs < 5.0) continue;
    
    dMr = MAX (0.005, 0.001*dMref[i]);
    dMo = MAX (0.005, 0.001*dMobs[i]);
      
    Dmag[Nkeep]  = (Mabs - 0.001*Mref[i]);
    dDmag[Nkeep] = (dMr*dMr + dMo*dMo);
    Nkeep ++;
  }

  if (Nkeep < 5) {
    fprintf (stderr, "too few stars\n");
    image[0].Mcal = 10000;
    image[0].dMcal = 10000;
    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 */

  /* find lo (Dmed - 0.1) + hi (Dmed + 0.1) 
  lo = -1; hi = -1;
  for (i = 0; (lo == -1) && (i < Nkeep); i++) { if (Dmed - Dmag[i] < 0.1) lo = i; }
  for (i = Nkeep - 1; (hi == -1) && (i >= 0); i--) { if (Dmag[i] - Dmed < 0.1) hi = i; }
  if (lo == -1) lo = 0;
  if (hi == -1) hi = Nkeep - 1;
  */

  /* 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 = 1000 * M1;
    image[0].dMcal = 1000 * M2 / sqrt (N);
    image[0].Mxxxx = N;
  } else {
    fprintf (stderr, "too few stars\n");
    image[0].Mcal = 10000;
    image[0].dMcal = 10000;
    image[0].Mxxxx = 0;
  }
}


