# include "astro.h"

int line_fit (stars1, stars2, N1, N2, dx, dy, Rx, Ry,
		X_O, X_X, X_Y, dX, Y_O, Y_X, Y_Y, dY)
Stars stars1[], stars2[];
int N1, N2;
double dx, dy, Rx, Ry;
double *X_O, *X_X, *X_Y, *dX;
double *Y_O, *Y_X, *Y_Y, *dY;
{

  int i, j, iteration;
  double Xmin, Ymin, Xmax, Ymax;
  double x, y, x2, y2, xy, N;
  double X, Y, X2, Y2, Xy, Yx;
  double Sx2, Sy2, Sxy, SXx, SXy, SYx, SYy;
  double Dx, Dy, DD, d2X, d2Y, factor, dX_Y, dY_X;
  FILE *f;

  factor = 1.0;
  for (iteration = 0; iteration < 10; iteration ++) {
    x = y = x2 = y2 = xy = X = Y = X2 = Xy = Yx = Y2 = N = 0;
    for (i = 0; i < N1; i++) {
      for (j = 0; j < N2; j++) {
	if (iteration == 0) {
	  /* first iteration, more lax, use 2*RADIUS */
	  Dx = 0.5*(stars1[i].X - Rx*stars2[j].X - dx) / Rx;
	  Dy = 0.5*(stars1[i].Y - Ry*stars2[j].Y - dy) / Ry;
	}
	else {
	  Dx = (stars1[i].X - (*X_O) - (*X_X)*stars2[j].X - (*X_Y)*stars2[j].Y) / Rx;
	  Dy = (stars1[i].Y - (*Y_O) - (*Y_X)*stars2[j].X - (*Y_Y)*stars2[j].Y) / Ry;
	}
	DD = factor * hypot (Dx, Dy);
	if (DD < RADIUS) {
	  x  += stars2[j].X;
	  y  += stars2[j].Y;
	  x2 += stars2[j].X*stars2[j].X;
	  y2 += stars2[j].Y*stars2[j].Y;
	  xy += stars2[j].X*stars2[j].Y;
	  X  += stars1[i].X;
	  Y  += stars1[i].Y;
	  X2 += stars1[i].X*stars2[j].X;
	  Xy += stars1[i].X*stars2[j].Y;
	  Yx += stars1[i].Y*stars2[j].X;
	  Y2 += stars1[i].Y*stars2[j].Y;
	  N  += 1.0;
	}
      }
    }

    fprintf (stderr, "%d matched stars, iteration %d\n", (int) N, iteration + 1);
    if (N < 3) {
      fprintf (stderr, "too few stars, skipping\n");
      *X_O = *X_X = *X_Y = *Y_O = *Y_X = *Y_Y = 0;
      return (N);
    }
    Sx2 = x2 - x*x/N;
    Sy2 = y2 - y*y/N;
    Sxy = xy - x*y/N;
    SXx = X2 - X*x/N;
    SXy = Xy - X*y/N;
    SYx = Yx - Y*x/N;
    SYy = Y2 - Y*y/N;
    
    *X_X = (SXx*Sy2 - SXy*Sxy) / (Sx2*Sy2 - Sxy*Sxy);
    *X_Y = (SXy*Sx2 - SXx*Sxy) / (Sx2*Sy2 - Sxy*Sxy);
    *X_O = X/N - (*X_X)*x/N - (*X_Y)*y/N;
    
    *Y_X = (SYx*Sy2 - SYy*Sxy) / (Sx2*Sy2 - Sxy*Sxy);
    *Y_Y = (SYy*Sx2 - SYx*Sxy) / (Sx2*Sy2 - Sxy*Sxy);
    *Y_O = Y/N - (*Y_X)*x/N - (*Y_Y)*y/N;
    /* make small course correction, keeping CCD rectilinear: */
/*    dX_Y = (*X_Y + *Y_X*(*Y_Y / *X_X)); */
/*    dY_X = 0.5*(*Y_X + *X_Y*(*X_Y / *Y_Y)); */
/*    *X_Y -= dX_Y;  */
/*    *Y_X -= dY_X; */
    factor = factor * 1.1;
  }
  factor = factor / 1.1;
  
  d2X = d2Y = N = 0;
  Ymin = Xmin = 10000;
  Xmax = Ymax = 0;
  for (i = 0; i < N1; i++) {
    for (j = 0; j < N2; j++)  {
      Dx = (stars1[i].X - (*X_O) - (*X_X)*stars2[j].X - (*X_Y)*stars2[j].Y) / Rx;
      Dy = (stars1[i].Y - (*Y_O) - (*Y_X)*stars2[j].X - (*Y_Y)*stars2[j].Y) / Ry;
      DD = factor * hypot (Dx, Dy);
      if (DD < RADIUS) {
	d2X += Dx*Dx;
	d2Y += Dy*Dy;
	N += 1.0;
	Xmin = MIN (Xmin, stars2[j].X);
	Xmax = MAX (Xmax, stars2[j].X);
	Ymin = MIN (Ymin, stars2[j].Y);
	Ymax = MAX (Ymax, stars2[j].Y);
      }
    }
  }
  fprintf (stderr, "NPOINTS: %f\n", N);
  fprintf (stderr, "Xrange: %f -- %f\n", Xmin, Xmax);
  fprintf (stderr, "Yrange: %f -- %f\n", Ymin, Ymax);
  (*dX) = sqrt(d2X / (N - 3.0));
  (*dY) = sqrt(d2Y / (N - 3.0));

  fprintf (stderr, "%f %f %f %f\n", *X_O, *X_X, *X_Y, *dX);
  fprintf (stderr, "%f %f %f %f\n", *Y_O, *Y_X, *Y_Y, *dY);

  return ((int)N);

}
