# include <ohana.h>

int main (int argc, char **argv) {

  float *ra1, *ra2, *dec1, *dec2;
  float x, y, x2, y2, xy, R, D, Rx, Ry, Dx, Dy, N;
  float x3, x2y, x4, Rx2;
  float DEC, RA, DX, DY, RX, RY, d2RA, d2DEC, dRA, dDEC;
  float Sx2, Sy2, Sxy, Srx, Sry, Sdx, Sdy;
  int i, Npairs, NPAIR;

  NPAIR = 100;
  ALLOCATE (ra1, float, NPAIR);
  ALLOCATE (ra2, float, NPAIR);
  ALLOCATE (dec1, float, NPAIR);
  ALLOCATE (dec2, float, NPAIR);

/*
  ALLOCATE (a, double *, 3);
  ALLOCATE (a[0], double, 3);
  ALLOCATE (a[1], double, 3);
  ALLOCATE (a[2], double, 3);
  ALLOCATE (a[3], double, 4);
  ALLOCATE (b, double *, 3);
  ALLOCATE (b[0], double, 1);
  ALLOCATE (b[1], double, 1);
  ALLOCATE (b[2], double, 1);
  ALLOCATE (b[3], double, 1);
*/

  if (argc > 1) {
    fprintf (stderr, "USAGE: list_astro x\nTakes list from the stdin in the");
    fprintf (stderr, " format:\n  X Y RA DEC\nReturns:\n  RA, RX, RY, dRA\n  DEC, DX, DY, dDEC\n");
    exit(0);
  }
    
  for (i = 0; fscanf (stdin, "%f %f %f %f", &ra1[i], &dec1[i], &ra2[i], &dec2[i]) != EOF; i++) {
    if (i == NPAIR - 1) {
      NPAIR += 50;
      REALLOCATE (ra1, float, NPAIR);
      REALLOCATE (ra2, float, NPAIR);
      REALLOCATE (dec1, float, NPAIR);
      REALLOCATE (dec2, float, NPAIR);
    }
  }

  
  
  N = x = y = x2 = y2 = xy = R = D = Rx = Ry = Dx = Dy = 0;
  Npairs = i;
  for (i = 0; i < Npairs; i++) {
    x  += ra1[i];
    y  += dec1[i];
    x2 += ra1[i]*ra1[i];
    y2 += dec1[i]*dec1[i];
    xy += ra1[i]*dec1[i];
    x3 += ra1[i]*ra1[i]*ra1[i];
    x2y += ra1[i]*ra1[i]*dec1[i];
    x4 += ra1[i]*ra1[i]*ra1[i]*ra1[i];
    R  += ra2[i];
    D  += dec2[i];
    Rx += ra2[i]*ra1[i];
    Ry += ra2[i]*dec1[i];
    Rx2 += ra2[i]*ra1[i]*ra1[i];
    Dx += dec2[i]*ra1[i];
    Dy += dec2[i]*dec1[i];
    N  += 1.0;
  }

/*
  a[0][0] = N;
  a[0][1] = a[1][0] = x;
  a[0][2] = a[2][0] = y;
  a[0][3] = a[3][0] = x2;
  a[1][1] = x2;
  a[1][2] = a[2][1] = xy;
  a[1][3] = a[3][1] = x3;
  a[2][2] = y2;
  a[2][3] = a[3][2] = x2y;
  a[3][3] = x4;
  b[0][0] = R;
  b[1][0] = Rx;
  b[2][0] = Ry;
  b[1][0] = Rx2;
*/

  Sx2 = x2 - x*x/N;
  Sy2 = y2 - y*y/N;
  Sxy = xy - x*y/N;
  Srx = Rx - R*x/N;
  Sry = Ry - R*y/N;
  Sdx = Dx - D*x/N;
  Sdy = Dy - D*y/N;
  
  RX = (Srx*Sy2 - Sry*Sxy) / (Sx2*Sy2 - Sxy*Sxy);
  RY = (Sry*Sx2 - Srx*Sxy) / (Sx2*Sy2 - Sxy*Sxy);
  RA = R/N - RX*x/N - RY*y/N;

/*
  gaussj (a, 4, b, 1);

  RA = b[0][0];
  RX = b[1][0];
  RY = b[2][0];
  RX2 = b[3][0];
*/
  
  DX  = (Sdx*Sy2 - Sdy*Sxy) / (Sx2*Sy2 - Sxy*Sxy);
  DY  = (Sdy*Sx2 - Sdx*Sxy) / (Sx2*Sy2 - Sxy*Sxy);
  DEC = D/N - DX*x/N - DY*y/N;

  d2RA = d2DEC = N = 0;
  N = x = y = x2 = y2 = xy = R = D = Rx = Ry = Dx = Dy = 0;
  for (i = 0; i < Npairs; i++) {
    d2RA += pow((ra2[i] - RA - RX*ra1[i] - RY*dec1[i]), 2.0);
    d2DEC += pow((dec2[i] - DEC - DX*ra1[i] - DY*dec1[i]), 2.0);
    N += 1.0;
  }

  dRA  = 3600.0*sqrt(d2RA / (N - 3.0));
  dDEC = 3600.0*sqrt(d2DEC / (N - 3.0));

  fprintf (stdout, " %12.9f %12.9e %12.9e %12.9f\n", RA, RX, RY, dRA);
  fprintf (stdout, " %12.9f %12.9e %12.9e %12.9f\n", DEC, DX, DY, dDEC);
  exit (0);
}


int gaussj (double **a, int n, double **b, int m) {

  int *indxc,*indxr,*ipiv;
  int i, icol, irow, j, k, l, ll;
  double big,dum,pivinv;
  
  ALLOCATE (indxc, int, n);
  ALLOCATE (indxr, int, n);
  ALLOCATE (ipiv, int, n);
  for (j = 0; j < n; j++) 
    ipiv[j] = 0;

  for (i = 0; i < n; i++) {
    big = 0.0;
    for (j = 0; j < n; j++) {
      if (ipiv[j] != 1) {
	for (k = 0; k < n; k++) {
	  if (ipiv[k] == 0) {
	    if (fabs (a[j][k]) >= big) {
	      big  = fabs (a[j][k]);
	      irow = j;
	      icol = k;
	    }
	  } 
	  else 
	    if (ipiv[k] > 1) {
	      fprintf (stderr, "GAUSSJ: Singular Matrix! (1)\n");
	      return (0);
	    }
	}
      }
    }
    ipiv[icol]++;
    if (irow != icol) {
      for (l = 0; l < n; l++) 
	SWAP (a[irow][l], a[icol][l]);
      for (l = 0; l < m; l++) 
	SWAP (b[irow][l], b[icol][l]);
    }
    indxr[i] = irow;
    indxc[i] = icol;
    if (a[icol][icol] == 0.0) {
      fprintf (stderr, "GAUSSJ: Singular Matrix! (2)\n");
      return (0);
    }
    pivinv = 1.0 / a[icol][icol];
    a[icol][icol] = 1.0;
    for (l = 0; l < n; l++) 
      a[icol][l] *= pivinv;
    for (l = 0; l < m; l++) 
      b[icol][l] *= pivinv;
    for (ll = 0; ll < n; ll++) {
      if (ll != icol) {
	dum = a[ll][icol];
	a[ll][icol] = 0.0;
	for (l = 0; l < n; l++) 
	  a[ll][l] -= a[icol][l]*dum;
	for (l = 0; l < m; l++) 
	  b[ll][l] -= b[icol][l]*dum;
      }
    }
  }

  for (l = n - 1; l >= 0; l--) {
    if (indxr[l] != indxc[l])
      for (k = 0; k < n; k++)
	SWAP (a[k][indxr[l]], a[k][indxc[l]]);
  }
  free (ipiv);
  free (indxr);
  free (indxc);
  return (1);
}

