# include "astro.h"
# define TIME {gettimeofday (&now, (struct timeval *) NULL); \
fprintf (stderr, "start time: %f\n", now.tv_sec+now.tv_usec*1e-6);}

void find_shift (stars1, stars2, N1, N2, Rx, Ry, n, S, dx, dy)
Stars stars1[], stars2[];
int N1, N2;
double Rx, Ry;
double *dx, *dy;
double *S;
int *n;
{

  double *DX, *DX2, *DY, *DY2, *N;
  double Xrange, Yrange, Xzero, Yzero;
  double dX, dY, sx, sy, s, dN, mean, sigma;
  double *buf;
  int i, j, dxi, dyi;
  struct timeval now;  
  double fx1, fx2, fy1, fy2;
  int k;

  *S = 1000000;
  *n = *dx = *dy = 0;
  NPIX = 1 + 30 * sqrt(N1*N2) / MIN (N1, N2);
  mean = 2*N1*N2 / (NPIX*NPIX);
  sigma = sqrt (mean);
  ALLOCATE (N,    double, NPIX*NPIX);
  ALLOCATE (DX,   double, NPIX*NPIX);
  ALLOCATE (DY,   double, NPIX*NPIX);
  ALLOCATE (DX2,  double, NPIX*NPIX);
  ALLOCATE (DY2,  double, NPIX*NPIX);
  bzero (N,   NPIX*NPIX*sizeof(double));
  bzero (DX,  NPIX*NPIX*sizeof(double));
  bzero (DY,  NPIX*NPIX*sizeof(double));
  bzero (DX2, NPIX*NPIX*sizeof(double));
  bzero (DY2, NPIX*NPIX*sizeof(double));
  ranges (stars1, stars2, N1, N2, &Xrange, &Yrange, &Xzero, &Yzero, Rx, Ry);
  fx1 = (NPIX - 1) / Xrange;
  fx2 = -fx1 * Xzero;
  fy1 = (NPIX - 1) / Yrange;
  fy2 = -fy1 * Yzero;
  for (i = 0; i < N1; i++) {
    for (j = 0; j < N2; j++) {
      dX = stars1[i].X - stars2[j].X;
      dY = stars1[i].Y - stars2[j].Y;
      dxi = fx1*dX + fx2;
      dyi = fy1*dY + fy2;
      k = dxi*NPIX + dyi;
      N[k]   += 1.0;
      DX[k]  += dX;
      DY[k]  += dY;
      DX2[k] += dX*dX;
      DY2[k] += dY*dY;
    }
  }

  for (i = 0; i < NPIX; i++) {
    for (j = 0; j < NPIX; j++) {
      if (N[i*NPIX + j] < 4) 
	continue;
      if (N[i*NPIX + j] - mean < 5*sigma) 
	continue;
		 
      sx = DX2[i*NPIX + j] / (N[i*NPIX + j]) - SQ(DX[i*NPIX + j]) / SQ(N[i*NPIX + j]);
      sy = DY2[i*NPIX + j] / (N[i*NPIX + j]) - SQ(DY[i*NPIX + j]) / SQ(N[i*NPIX + j]);
      /* sx, sy are the positional scatter in the bin */
      s = sqrt (sx + sy + 1e-11) / SQ(N[i*NPIX + j] - 3);
      if (N[i*NPIX + j] > mean + 10*sigma) {
	if (s < *S) { 
	  *S = s;
	  *n = N[i*NPIX + j];
	  *dx = DX[i*NPIX + j] / (double) N[i*NPIX + j];
	  *dy = DY[i*NPIX + j] / (double) N[i*NPIX + j];
	}
	/* fprintf (stderr, "point: %f %f  %f %f\n", s, N[i*NPIX + j],  DX[i*NPIX + j] / (double) N[i*NPIX + j], DY[i*NPIX + j] / (double) N[i*NPIX + j]);  */
      }
    }
  }

  free (N);
  free (DX);
  free (DY);
  free (DX2);
  free (DY2);

}
