Index: /trunk/Ohana/src/relastro/Makefile
===================================================================
--- /trunk/Ohana/src/relastro/Makefile	(revision 39238)
+++ /trunk/Ohana/src/relastro/Makefile	(revision 39239)
@@ -10,4 +10,6 @@
 MAN	=	$(HOME)/doc
 INC	= 	$(HOME)/include
+TESTDIR	=	$(HOME)/test
+TESTBIN	=	$(HOME)/test
 include ../../Makefile.Common
 
@@ -21,6 +23,4 @@
 
 install: $(DESTBIN)/relastro $(DESTBIN)/relastro_client 
-
-# $(DESTBIN)/testparallax
 
 RELASTRO = \
@@ -161,2 +161,16 @@
 $(TESTPAR): $(INC)/relastro.h
 $(BIN)/testparallax.$(ARCH): $(TESTPAR)
+
+fitpm: $(BIN)/fitpm.$(ARCH)
+
+FITPM = \
+$(SRC)/fitpm.$(ARCH).o 	      \
+$(SRC)/ParFactor.$(ARCH).o           \
+$(SRC)/FitAstromOps.$(ARCH).o        \
+$(SRC)/FitPM.$(ARCH).o               \
+$(SRC)/FitPM_IRLS.$(ARCH).o               \
+$(SRC)/mkpolyterm.$(ARCH).o            \
+$(SRC)/fitpoly.$(ARCH).o
+
+$(FITPM): $(INC)/relastro.h
+$(BIN)/fitpm.$(ARCH): $(FITPM)
Index: /trunk/Ohana/src/relastro/src/FitAstromOps.c
===================================================================
--- /trunk/Ohana/src/relastro/src/FitAstromOps.c	(revision 39238)
+++ /trunk/Ohana/src/relastro/src/FitAstromOps.c	(revision 39239)
@@ -2,4 +2,64 @@
 
 // I am modifying FitPM with an eye to (a) threaded operations and (b) bootstrap resampling tests.
+
+int FitAstromPoints_Project (FitStats *fitStats, double *Tmean, double *Trange, double *parRange) {
+
+  int k;
+
+  int Npoints = fitStats->Npoints;
+  FitAstromPoint *points = fitStats->points;
+
+  // find Tmin & Tmax from the list of accepted measurements
+  double Tmin  = points[0].T;
+  double Tmax  = points[0].T;
+  double pRmin = +2.0;
+  double pRmax = -2.0;
+  double pDmin = +2.0;
+  double pDmax = -2.0;
+
+  *Tmean = 0.0;
+
+  double Tsum = 0.0;
+  double Wsum = 0.0;
+  for (k = 0; k < Npoints; k++) {
+    Tmin = MIN(Tmin, points[k].T);
+    Tmax = MAX(Tmax, points[k].T);
+
+    float wx = 1.0 / SQ(points[k].dX);
+
+    Tsum += points[k].T * wx;
+    Wsum += wx;
+
+    // at this point, T is in years since J2000
+    ParFactor (&points[k].pR, &points[k].pD, points[k].R, points[k].D, points[k].T);
+    pRmin = MIN (pRmin, points[k].pR);
+    pRmax = MAX (pRmax, points[k].pR);
+    pDmin = MIN (pDmin, points[k].pD);
+    pDmax = MAX (pDmax, points[k].pD);
+  }
+  *Trange = Tmax - Tmin;
+
+  // mean epoch
+  *Tmean = Tsum / Wsum;
+
+  // for HIGH_SPEED, just use the center of the range
+  if (RELASTRO_OP == OP_HIGH_SPEED) {
+    *Tmean = 0.5*(Tmax - Tmin);
+  }
+
+  *parRange = hypot (pRmax - pRmin, pDmax - pDmin);
+
+  /* we need to do the fit in a locally linear space; choose a ref coordinate */
+  fitStats->coords.crval1 = points[0].R;
+  fitStats->coords.crval2 = points[0].D;
+
+  // project all of the R,D coordinates to a plane centered on this coordinate. set
+  // the times to be relative to Tmean
+  for (k = 0; k < Npoints; k++) {
+    RD_to_XY (&points[k].X, &points[k].Y, points[k].R, points[k].D, &fitStats->coords);
+    points[k].T -= *Tmean;
+  }	  
+  return TRUE;
+}
 
 FitStats *FitStatsInit (int Nmax, int Nboot) {
Index: /trunk/Ohana/src/relastro/src/UpdateObjects.c
===================================================================
--- /trunk/Ohana/src/relastro/src/UpdateObjects.c	(revision 39238)
+++ /trunk/Ohana/src/relastro/src/UpdateObjects.c	(revision 39239)
@@ -625,64 +625,4 @@
 }
 
-int FitAstromPoints_Project (FitStats *fitStats, double *Tmean, double *Trange, double *parRange) {
-
-  int k;
-
-  int Npoints = fitStats->Npoints;
-  FitAstromPoint *points = fitStats->points;
-
-  // find Tmin & Tmax from the list of accepted measurements
-  double Tmin  = points[0].T;
-  double Tmax  = points[0].T;
-  double pRmin = +2.0;
-  double pRmax = -2.0;
-  double pDmin = +2.0;
-  double pDmax = -2.0;
-
-  *Tmean = 0.0;
-
-  double Tsum = 0.0;
-  double Wsum = 0.0;
-  for (k = 0; k < Npoints; k++) {
-    Tmin = MIN(Tmin, points[k].T);
-    Tmax = MAX(Tmax, points[k].T);
-
-    float wx = 1.0 / SQ(points[k].dX);
-
-    Tsum += points[k].T * wx;
-    Wsum += wx;
-
-    // at this point, T is in years since J2000
-    ParFactor (&points[k].pR, &points[k].pD, points[k].R, points[k].D, points[k].T);
-    pRmin = MIN (pRmin, points[k].pR);
-    pRmax = MAX (pRmax, points[k].pR);
-    pDmin = MIN (pDmin, points[k].pD);
-    pDmax = MAX (pDmax, points[k].pD);
-  }
-  *Trange = Tmax - Tmin;
-
-  // mean epoch
-  *Tmean = Tsum / Wsum;
-
-  // for HIGH_SPEED, just use the center of the range
-  if (RELASTRO_OP == OP_HIGH_SPEED) {
-    *Tmean = 0.5*(Tmax - Tmin);
-  }
-
-  *parRange = hypot (pRmax - pRmin, pDmax - pDmin);
-
-  /* we need to do the fit in a locally linear space; choose a ref coordinate */
-  fitStats->coords.crval1 = points[0].R;
-  fitStats->coords.crval2 = points[0].D;
-
-  // project all of the R,D coordinates to a plane centered on this coordinate. set
-  // the times to be relative to Tmean
-  for (k = 0; k < Npoints; k++) {
-    RD_to_XY (&points[k].X, &points[k].Y, points[k].R, points[k].D, &fitStats->coords);
-    points[k].T -= *Tmean;
-  }	  
-  return TRUE;
-}
-
 int CatalogMaxNmeasure (Catalog *catalog, int Ncatalog) {
 
Index: /trunk/Ohana/src/relastro/src/fitpm.c
===================================================================
--- /trunk/Ohana/src/relastro/src/fitpm.c	(revision 39239)
+++ /trunk/Ohana/src/relastro/src/fitpm.c	(revision 39239)
@@ -0,0 +1,205 @@
+# include "relastro.h"
+
+double rnd_gauss (double mean, double sigma);
+void gauss_init (int Nbin);
+double gaussian (double x, double mean, double sigma);
+int mkstar (FitStats *fitStats, double Ro, double Do, double uR, double uD, double plx, int Npoints, int Nbad);
+
+# define MJD_MIN_MJD 55197   /* 2010/01/01,00:00:00 */
+# define MJD_MAX_MJD 57023   /* 2015/01/01,00:00:00 */
+# define MJD_MIN_YRS 10.00   /* 2010/01/01,00:00:00 */
+# define MJD_MAX_YRS 14.9993 /* 2015/01/01,00:00:00 */
+# define MJD_REF_YRS 12.4148 /* 2012/06/01,00:00:00 */
+
+# define POS_ERROR 0.010 /* arcsec */
+# define OUT_ERROR 0.500 /* arcsec */
+
+# define N_POINTS 100
+# define N_OUTLIERS 10
+
+# define dcos(THETA) cos(RAD_DEG*THETA)
+# define dsin(THETA) sin(RAD_DEG*THETA)
+
+int main (int argc, char **argv) {
+  
+  // plan_tests (14);
+
+  // diag ("relastro fitpm tests");
+  
+  // init the random seed
+  { 
+    struct timeval now;
+    gettimeofday (&now, NULL);
+    long A = now.tv_sec + now.tv_usec * 1000000;
+    srand48(A);
+  }
+  gauss_init (2048);
+
+  FitStats *fitStats = FitStatsInit (N_POINTS + N_OUTLIERS, 0);
+
+  int i;
+
+  // XXX try a single star to see if it looks good:
+  if (0) {
+    mkstar (fitStats, 3.0, 2.0, 0.2, 0.4, 0.1, N_POINTS, N_OUTLIERS);
+    for (i = 0; i < fitStats->Npoints; i++) {
+      fprintf (stdout, "%f %f\n", 3600*(fitStats->points[i].R - 3.0), 3600*(fitStats->points[i].D - 2.0));
+    }
+    exit (0);
+  }
+
+  int Nstars = 1000;
+  for (i = 0; i < Nstars; i++) {
+
+    double Ro  = 360.0*(drand48() - 0.5);
+    double Phi = 2.0*(drand48() - 0.5);
+    double Do  = DEG_RAD * asin(Phi);
+
+    double uR  = 0.5*(drand48() - 0.5);
+    double uD  = 0.5*(drand48() - 0.5);
+    double plx = 0.0*(drand48() + 0.0);
+
+    mkstar (fitStats, Ro, Do, uR, uD, plx, N_POINTS, N_OUTLIERS);
+
+    double Tmean, Trange, parRange;
+    FitAstromPoints_Project (fitStats, &Tmean, &Trange, &parRange);
+
+    FitAstromResult fitPM;    
+    FitAstromResultInit (&fitPM);
+
+    if (0) {
+      FitPM (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints);
+    } else {
+      FitPM_IRLS (&fitPM, fitStats->fitdataPM, fitStats->points, fitStats->Npoints, 0);
+    }
+
+    XY_to_RD (&fitPM.Ro, &fitPM.Do, fitPM.Ro, fitPM.Do, &fitStats->coords);
+
+    fprintf (stdout, "%12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f %12.8f %8.6f %8.6f %8.6f | %12.8f\n",
+	     Ro, Do, uR, uD, plx, 
+	     fitPM.Ro, fitPM.Do, fitPM.uR, fitPM.uD, fitPM.p, Tmean); 
+
+    // ok (TRUE, "fitted star");
+  }
+  // return exit_status();
+  exit (0);
+}
+
+int mkstar (FitStats *fitStats, double Ro, double Do, double uR, double uD, double plx, int Npoints, int Nbad) {
+  
+  int i;
+
+  int Ntotal = Npoints + Nbad;
+
+  FitAstromPoint *points = fitStats->points;
+
+  for (i = 0; i < Npoints; i++) {
+    FitAstromPointInit (&points[i]);
+
+    // ParFactor expects a time which is in years since J2000 (MJD_..._YRS is so defined)
+    points[i].T = MJD_MIN_YRS + drand48()*(MJD_MAX_YRS - MJD_MIN_YRS);
+    
+    double pR, pD;
+    ParFactor (&pR, &pD, Ro, Do, points[i].T);
+
+    double dR = rnd_gauss(0.0, POS_ERROR);
+    double dD = rnd_gauss(0.0, POS_ERROR);
+
+    points[i].R = Ro + (dR + plx*pR + uR*(points[i].T - MJD_REF_YRS))/3600/dcos(Do);
+    points[i].D = Do + (dD + plx*pD + uD*(points[i].T - MJD_REF_YRS))/3600;
+
+    points[i].dX = POS_ERROR;
+    points[i].dY = POS_ERROR;
+  }
+
+  for (i = Npoints; i < Ntotal; i++) {
+    FitAstromPointInit (&points[i]);
+
+    // ParFactor expects a time which is in years since J2000 (MJD_..._YRS is so defined)
+    points[i].T = MJD_MIN_YRS + drand48()*(MJD_MAX_YRS - MJD_MIN_YRS);
+    
+    double pR, pD;
+    ParFactor (&pR, &pD, Ro, Do, points[i].T);
+
+    double dR = OUT_ERROR*(drand48() - 0.5);
+    double dD = OUT_ERROR*(drand48() - 0.5);
+
+    points[i].R = Ro + (dR + plx*pR + uR*(points[i].T - MJD_REF_YRS))/3600/dcos(Do);
+    points[i].D = Do + (dD + plx*pD + uD*(points[i].T - MJD_REF_YRS))/3600;
+
+    points[i].dX = POS_ERROR;
+    points[i].dY = POS_ERROR;
+  }
+  fitStats->Npoints = Ntotal;
+
+  return TRUE;
+}
+
+static int Ngaussint = 0;
+static double *gaussint;
+
+extern double drand48();
+
+double gaussian (double x, double mean, double sigma) {
+
+  double f;
+
+  f = exp (-0.5 * SQ(x - mean) / SQ(sigma)) / sqrt(2 * M_PI * SQ(sigma));
+
+  return (f);
+
+}
+
+/* integrate a gaussian from -5 sigma to +5 sigma */
+void gauss_init (int Nbin) {
+ 
+  int i;
+  double val, x, dx, dx1, dx2, dx3, df;
+  double mean, sigma;
+ 
+  /* no need to generate this if it already exists */
+  if (Ngaussint == Nbin) return;
+
+  // A = time(NULL);
+  // // XXX this is expensive if called a lot (1 sec min)
+  // // for (B = 0; A == time(NULL); B++);
+  // B = A + 10000;
+  // srand48(B);
+ 
+  Ngaussint = Nbin;
+  ALLOCATE (gaussint, double, Ngaussint + 1);
+
+  val = 0;
+  dx = 1.0 / Ngaussint;
+  dx1 = dx / 3.0;
+  dx2 = 2.0*dx/3.0;
+  dx3 = dx;
+  mean = 0.0;
+  sigma = 1.0;
+ 
+  for (i = 0, x = -7.0; (i < Ngaussint) && (x < 7.0); x += dx)  {
+    df = (3.0*gaussian(x    , mean, sigma) + 
+          9.0*gaussian(x+dx1, mean, sigma) +
+          9.0*gaussian(x+dx2, mean, sigma) + 
+          3.0*gaussian(x+dx3, mean, sigma)) * (dx1/8.0);
+    val += df;
+    if (val > (i + 0.5) / (double) Ngaussint) {
+      gaussint[i] = x + dx / 2.0;
+      i++;
+    }
+  }
+}
+
+double rnd_gauss (double mean, double sigma) {
+ 
+  int i;
+  double y;
+ 
+  y = drand48();
+  i = Ngaussint*y;
+  y = gaussint[i]*sigma + mean;
+ 
+  return (y);
+ 
+}
+