# include "addstar.h"

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

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 (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);
}
  
SaveCalibration (float M, float dM, float Mr, float Mc, float A, int N, float ra, float dec, float x, float y) {

  Mobs[Ncal] = M;
  dMobs[Ncal] = dM;
  Mref[Ncal] = Mr;
  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 (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);
  }
}

FindCalibration (Image *image) {

  int i, MaxN, *Nlist, Nkeep;
  float N, M1, M2, Klam, Clam, Xlam, Mabs, dM, dm, *Dmag;
  float dAs, dMs, dYs;
  
  /* 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);
  Nkeep = 0;
  Klam = thiscode[0].K;
  Clam = thiscode[0].C;
  Xlam = thiscode[0].X;
  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;

    /* subpix correction : make parameters global! */
    dAs = 0;
    if (Y[i] < 400) dAs = 0.06 - 0.00015*Y[i];
    dYs = Y[i] - (int) (Y[i]) + 0.55;
    dMs = dAs * sin (6.28*dYs);
    /* Mabs -= dMs; */

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

    /* skip stars brighter than 8.0 */
    /* if (Mabs < 8.0) continue; */
    
    Dmag[Nkeep] = (Mabs - 0.001*Mref[i]);
    Nkeep ++;
  }

  if (Nkeep < 3) {
    fprintf (stderr, "too few stars\n");
    image[0].Mcal = 10000;
    image[0].dMcal = 10000;
    return;
  }
  fsort (Dmag, Nkeep);

  /* accumulate delta mags (25% - 75% of list) */
  M1 = 0.0;
  M2 = 0.0;
  N  = 0.0;
  for (i = 0.25*Nkeep; i <= 0.75*Nkeep; i++) {

    /* a straight mean works better, given the bad photometry of skyprobe */
    M1 += Dmag[i];
    M2 += SQ (Dmag[i]);
    N  += 1.0; 

  }

  if (N > 1) {
    M1 = M1 / N;
    M2 = sqrt (fabs(M2/N - M1*M1));
    fprintf (stderr, "N: %.0f, mean: %f, stdev: %f, precision: %f\n", N, M1, M2, M2 / sqrt (N));
    image[0].Mcal = 1000 * M1;
    image[0].dMcal = 1000 * M2 / sqrt (N);
  } else {
    fprintf (stderr, "too few stars\n");
    image[0].Mcal = 10000;
    image[0].dMcal = 10000;
  }
}


