# include "relastro.h"

static int Nmax;
static double *X, *dX;
static double *Y, *dY;
static double *R, *dR;
static double *D, *dD;
static double *pX;
static double *pY;
static double *T;
static double *dT;

void initObjectData (Catalog *catalog, int Ncatalog) {

  int i, j;
  
  Nmax = 0;
  for (i = 0; i < Ncatalog; i++) {
    for (j = 0; j < catalog[i].Naverage; j++) {
      Nmax = MAX (Nmax, catalog[i].average[j].Nm);
    }
  }

  ALLOCATE (R, double, MAX (1, Nmax));
  ALLOCATE (D, double, MAX (1, Nmax));
  ALLOCATE (T, double, MAX (1, Nmax));
  ALLOCATE (X, double, MAX (1, Nmax));
  ALLOCATE (Y, double, MAX (1, Nmax));

  ALLOCATE (dR, double, MAX (1, Nmax));
  ALLOCATE (dD, double, MAX (1, Nmax));
  ALLOCATE (dT, double, MAX (1, Nmax));
  ALLOCATE (dX, double, MAX (1, Nmax));
  ALLOCATE (dY, double, MAX (1, Nmax));

  ALLOCATE (pX, double, MAX (1, Nmax));
  ALLOCATE (pY, double, MAX (1, Nmax));
}  

int UpdateObjects (Catalog *catalog, int Ncatalog) {

  int i, j, k, m, N;
  StatType statsR, statsD;
  Coords coords;
  PMFit fit;
  time_t To;
  int mode, Nave, Npm, Npar, Nskip;
  double Tmin, Tmax;
  float errorScale;
  PhotCode *code;

  initObjectData (catalog, Ncatalog);

  /* project coordinates to a plane centered on the object with units of arcsec */
  coords.crval1 = 0;
  coords.crval2 = 0;
  coords.crpix1 = 0;
  coords.crpix2 = 0;
  coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0;
  coords.pc1_1  = coords.pc2_2 = 1.0;
  coords.pc1_2  = coords.pc2_1 = 0.0;
  coords.Npolyterms = 1;
  strcpy (coords.ctype, "RA---SIN");

  // use J2000 as a reference time
  To = ohana_date_to_sec ("2000/01/01");
  Nave = Npar = Npm = 0;

  for (i = 0; i < Ncatalog; i++) {

    if (VERBOSE) fprintf (stderr, "astrometrize catalog %d : %d ave, %d meas\n", i, catalog[i].Naverage, catalog[i].Nmeasure);

    Nskip = 0;
    for (j = 0; j < catalog[i].Naverage; j++) {

      /* calculate the average value of R,D for a single star */
      if (catalog[i].average[j].code & STAR_BAD) {
	Nskip ++;
	continue;  
      }

      N = 0;
      m = catalog[i].average[j].offset;

      Tmin = Tmax = (catalog[i].measure[m].t - To) / (86400*365.25);
      mode = FIT_MODE;

      for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
	if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
	
	R[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
	D[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
	T[N] = (catalog[i].measure[m].t - To) / (86400*365.25) ; // time relative to J2000 in years

	Tmin = MIN(Tmin, T[N]);
	Tmax = MAX(Tmax, T[N]);

	/* the astrometric errors are not being carried yet (but should be!) */
	/* we use the photometric mag error as a weighting term */

	code = GetPhotcodebyCode (catalog[0].measure[m].photcode);
	errorScale = code[0].astromScale;
	dR[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR) * errorScale;
	dD[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR) * errorScale;
	dT[N] = catalog[i].measure[m].dt;

	N++;
      }

      if ((Tmax - Tmin) < PM_DT_MIN) mode = FIT_AVERAGE;
      if ((mode == FIT_PM_ONLY) && (N < PM_TOOFEW)) mode = FIT_AVERAGE;

      // XXX This criterion needs to be better considered: adjust to match Ndof
      if (N < POS_TOOFEW) { /* too few measurements */
	catalog[i].average[j].code |= ID_STAR_FEW;
	continue;
      } else {
	catalog[i].average[j].code &= ~ID_STAR_FEW;
      }	

      /* we need to do the fit in a locally linear space; choose a ref coordinate */
      coords.crval1 = R[0];
      coords.crval2 = D[0];
      
      /* project all of the R,D coordinates to a plane centered on this coordinate */
      for (k = 0; k < N; k++) {
	RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords);
	dX[k] =  dR[k];
	dY[k] =  dD[k];
	// fprintf (stderr, "%d %f %f %f  %f %f\n", k, T[k], R[k], D[k], X[k], Y[k]);
      }	  

      /* fit the model components as needed */
      switch (mode) {
	case FIT_AVERAGE:
	  liststats (R, dR, N, &statsR);
	  liststats (D, dD, N, &statsD);

	  fit.Ro = statsR.mean;
	  fit.dRo = 3600.0*statsR.sigma;

	  fit.Do = statsD.mean;
	  fit.dDo = 3600.0*statsD.sigma;

	  fit.chisq = 0.5*(statsR.chisq + statsD.chisq);
	  fit.Nfit = N;
	  Nave ++;
	  break;

	case FIT_PM_ONLY:
	  FitPM (&fit, X, dX, Y, dY, T, N);
	  // fprintf (stderr, "fitted:  %f - %f : %f %f : %f %f : %f\n", Tmin, Tmax, fit.Ro, fit.Do, fit.uR, fit.uD, fit.p);
	  // project Ro, Do back to RA,DEC
	  XY_to_RD (&fit.Ro, &fit.Do, fit.Ro, fit.Do, &coords);
	  // fprintf (stderr, "project: %f %f : %f %f : %f\n", fit.Ro, fit.Do, fit.uR, fit.uD, fit.p);
	  // continue;
	  Npm ++;
	  break;

	case FIT_PAR_ONLY:
	  for (k = 0; k < N; k++) {
	    ParFactor (&pX[k], &pY[k], R[k], D[k], T[k]);
	  }
	  FitPar (&fit, X, dX, Y, dY, pX, pY, N);

	  // project Ro, Do back to RA,DEC
	  XY_to_RD (&fit.Ro, &fit.Do, fit.Ro, fit.Do, &coords);
	  Npar ++;
	  break;

	case FIT_PM_AND_PAR:
	  for (k = 0; k < N; k++) {
	    ParFactor (&pX[k], &pY[k], R[k], D[k], T[k]);
	  }
	  FitPMandPar (&fit, X, dX, Y, dY, T, pX, pY, N);
	  XY_to_RD (&fit.Ro, &fit.Do, fit.Ro, fit.Do, &coords);
	  Npar ++;
	  break;

	default:
	  fprintf (stderr, "programming error at %s, %d", __FILE__, __LINE__);
	  exit (2);
      }	  

      if (0 && (j < 100)) {
	fprintf (stderr, "%f %f -> %f %f (%f,%f)\n",
		 catalog[i].average[j].R, 
		 catalog[i].average[j].D, 
		 fit.Ro, fit.Do, 
		 3600*(catalog[i].average[j].R - fit.Ro), 
		 3600*(catalog[i].average[j].D - fit.Do));
      }

      // the measure fields must be updated before the average fields
      m = catalog[i].average[j].offset;
      for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
	if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
	setMeanR (fit.Ro, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
	setMeanD (fit.Do, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
      }      

      catalog[i].average[j].R  	= fit.Ro; // RA in degrees
      catalog[i].average[j].D  	= fit.Do; // DEC in degrees
      catalog[i].average[j].dR 	= fit.dRo; // RA scatter in arcsec
      catalog[i].average[j].dD 	= fit.dDo; // DEC scatter in arcsec

      catalog[i].average[j].uR  = fit.uR; // RA proper motion in arcsec/year
      catalog[i].average[j].uD  = fit.uD; // DEC proper motion in arcsec/year
      catalog[i].average[j].duR = fit.duR; // RA proper motion error in arcsec/year
      catalog[i].average[j].duD = fit.duD; // DEC proper motion error in arcsec/year

      catalog[i].average[j].P   = fit.p; // parallax in arcsec
      catalog[i].average[j].dP  = fit.dp; // parallax error in arcsec

      catalog[i].average[j].Xp  = (fit.Nfit > 1) ? 100.0*log10(fit.chisq) : NAN_S_SHORT;
    }

    if (VERBOSE) fprintf (stderr, "catalog %d : %d ave, %d pm, %d par : Nskip % d\n", i, Nave, Npm, Npar, Nskip);
  }

  if (VERBOSE) fprintf (stderr, "fitted %d objects (%d ave, %d pm, %d par)\n", Nave + Npm + Npar, Nave, Npm, Npar);
  return (TRUE);
}

/* fitting proper-motion and parallax:

given a source at position r,d, at a time t, we need to calculate a vector (pr,pd)

let x,y be the coordinate in the linearized frame with y parallel to DEC lines

L,B are the ecliptic longitude and latitude of the object, 
dL and dB are the offsets in the L and B directions

dL = sin(t - topp)
dB = cos(t - topp)*sin(B)

these need to be rotated to the R,D frame to yield pR,pD.  Then, the equation of motion
for the source in the x,y frame is:

x = Ro + uR * (t - to) + p * pR 
y = Do + uD * (t - to) + p * pD

the unknowns in these equations are Ro, uR, Do, uD, and p

XXX think through the concepts for the pole a bit better.  all objects near the pole 
move the same way with the same phase.  choose a projection center and define dL,dB relative 
to that center point coordinate system?

*/
