Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitPM.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitPM.c	(revision 12047)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/FitPM.c	(revision 12047)
@@ -0,0 +1,101 @@
+# include "relastro.h"
+
+/* do we want an init function which does the alloc and a clear function to free? */
+
+FitPM (double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts) {
+
+  double **A, **B;
+
+  ALLOCATE (A, double *, 5);
+  ALLOCATE (B, double *, 5);
+  for (i = 0; i < 5; i++) {
+    ALLOCATE (A[i], double, 5);
+    ALLOCATE (B[i], double, 1);
+    bzero (A[i], 5*sizeof(double));
+    bzero (B[i], 1*sizeof(double));
+  }
+
+  for (i = 0; i < Npts; i++) {
+    /* handle case where dX or dY = 0.0 */
+    wx = 1.0 / SQ(dX[i]);
+    wy = 1.0 / SQ(dY[i]);
+
+    Wx += wx;
+    Wy += wy;
+
+    Tx += T[i]*wx;
+    Ty += T[i]*wy;
+    
+    Tx2 += SQ(T[i])*wx;
+    Ty2 += SQ(T[i])*wy;
+    
+    PR += pR[i]*wx;
+    PD += pD[i]*wy;
+    
+    PRT += pR[i]*T[i]*wx;
+    PDT += pD[i]*T[i]*wy;
+    
+    PRX += pR[i]*X[i]*wx;
+    PDY += pD[i]*Y[i]*wy;
+    
+    PR2 += SQ(pR[i])*wx;
+    PD2 += SQ(pD[i])*wy;
+
+    Xs += X[i]*wx;
+    Ys += Y[i]*wy;
+
+    XT += X[i]*T[i]*wx;
+    YT += Y[i]*T[i]*wy;
+  }
+
+  A[0][0] = Wx;
+  A[0][1] = Tx;
+  A[0][4] = PR;
+
+  A[1][0] = Tx;
+  A[1][1] = Tx2;
+  A[1][4] = PRT;
+
+  A[2][2] = Wy;
+  A[2][3] = Ty;
+  A[2][4] = PD;
+
+  A[3][2] = Ty;
+  A[3][3] = Ty2;
+  A[3][4] = PDT;
+
+  A[4][0] = PR;
+  A[4][1] = PRT;
+  A[4][2] = PD;
+  A[4][3] = PDT;
+  A[4][4] = PR2 + PD2;
+
+  B[0][0] = Xs;
+  B[1][0] = XT;
+  B[2][0] = Ys;
+  B[3][0] = YT;
+  B[4][0] = PRX + PDY;
+
+  dgaussj (A, 5, B, 1);
+
+  Ro = B[0];
+  uR = B[1];
+  Do = B[2];
+  uD = B[3];
+  p  = B[4];
+  
+  dRo = sqrt(A[0][0]);
+  duR = sqrt(A[1][1]);
+  dDo = sqrt(A[2][2]);
+  duD = sqrt(A[3][3]);
+  dp  = sqrt(A[4][4]);
+  
+  for (i = 0; i < 5; i++) {
+    free (A[i]);
+    free (B[i]);
+  }
+  free (A);
+  free (B);
+
+  return (fit);
+}
Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateObjects.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateObjects.c	(revision 12046)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/UpdateObjects.c	(revision 12047)
@@ -2,6 +2,9 @@
 
 static int Nmax;
-static double *listR, *listD;
-static double *dlistR, *dlistD;
+static double *listX, *dlistX;
+static double *listY, *dlistY;
+static double *listR, *dlistR;
+static double *listD, *dlistD;
+static double *listT, *dlistT;
 
 void initObjectData (Catalog *catalog, int Ncatalog) {
@@ -25,4 +28,15 @@
   float chisq;
   StatType statsR, statsD;
+  Coords coords;
+
+  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");
 
   for (i = 0; i < Ncatalog; i++) {
@@ -39,4 +53,5 @@
 	listR[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
 	listD[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
+	listT[N] = catalog[i].measure[m].t;
 
 	/* the astrometric errors are not being carried yet (but should be!) */
@@ -44,4 +59,6 @@
 	dlistR[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
 	dlistD[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
+	dlistT[N] = 1.0; /* XXX hard-wire this for now */
+
 	N++;
       }
@@ -52,18 +69,35 @@
       }	
 
+      /* we need to do the fit in a locally linear space; choose a ref coordinate */
+      coords.crval1 = listR[0];
+      coords.crval2 = listD[0];
+      
+      /* project all of the R,D coordinates to a plane centered on this coordinate */
+      /* calculate pR[i], pD[i] for each point */
+      for (k = 0; k < N; k++) {
+	  RD_to_XY (&listX[k], &listY[k], listR[k], listD[k], &coords);
+	  ParFactor (&pX[k], &pY[k], listR[k], listD[k]);
+      }	  
+
       /* in here, we should be fitting the parallax and proper-motion components */
+      /* XXX if I fit R,D,T as above, I will be getting the wrong answers near the pole */
+      if (FitProperMotion) {
+	  fit = FitPM (listX, dlistX, listY, dlistY, listT, pR, pD, N);
+      }
 
-      liststats (listR, dlistR, N, &statsR);
-      liststats (listR, dlistR, N, &statsD);
+      if (AverageCoords) {
+	  liststats (listR, dlistR, N, &statsR);
+	  liststats (listD, dlistD, N, &statsD);
 
-      catalog[i].average[j].R = statsR.mean;
-      catalog[i].average[j].dR = statsR.sigma;
+	  catalog[i].average[j].R = statsR.mean;
+	  catalog[i].average[j].dR = statsR.sigma;
 
-      catalog[i].average[j].D = statsD.mean;
-      catalog[i].average[j].dD = statsD.sigma;
+	  catalog[i].average[j].D = statsD.mean;
+	  catalog[i].average[j].dD = statsD.sigma;
 
-      chisq = 0.5*(statsR.chisq + statsD.chisq);
+	  chisq = 0.5*(statsR.chisq + statsD.chisq);
+      }
+
       catalog[i].average[j].Xp = (statsR.Nmeas > 1) ? 100.0*log10(chisq) : NO_MAG;
-
     }
   }
@@ -71,2 +105,27 @@
 }
 
+/* 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?
+
+*/
Index: /branches/dvo-mods-2007-02/Ohana/src/relastro/src/dvo_astrom_ops.c
===================================================================
--- /branches/dvo-mods-2007-02/Ohana/src/relastro/src/dvo_astrom_ops.c	(revision 12047)
+++ /branches/dvo-mods-2007-02/Ohana/src/relastro/src/dvo_astrom_ops.c	(revision 12047)
@@ -0,0 +1,41 @@
+# include "relastro.h"
+
+double getMeanR (Measure *measure, Average *average, SecFilt *secfilt) {
+
+    double ra;
+
+    /* the measure carries the instantaneous mean position at the epoch t */ 
+    ra = average[0].R - measure[0].dR / 3600.0;
+
+    /* possible corrections to mean ra:
+
+    - proper-motion and parallax
+    - abberation
+    - precession and nutation, etc
+    - refraction
+    - DCR
+
+    */
+
+    return (ra);
+}
+
+double getMeanD (Measure *measure, Average *average, SecFilt *secfilt) {
+
+    double dec;
+
+    /* the measure carries the instantaneous mean position at the epoch t */ 
+    dec = average[0].D - measure[0].dD / 3600.0;
+
+    /* possible corrections to mean ra:
+
+    - proper-motion and parallax
+    - abberation
+    - precession and nutation, etc
+    - refraction
+    - DCR
+
+    */
+
+    return (dec);
+}