Index: /trunk/Ohana/src/tools/Makefile
===================================================================
--- /trunk/Ohana/src/tools/Makefile	(revision 41683)
+++ /trunk/Ohana/src/tools/Makefile	(revision 41684)
@@ -61,4 +61,25 @@
 $(PROGRAMS) : % : $(BIN)/%.$(ARCH)
 
+$(BIN)/mpcorb_predict.$(ARCH) : $(SRC)/mpcorb_predict.$(ARCH).o
+	@if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi
+	$(CC) $(FULL_CFLAGS) -o $@ $^ $(FULL_LDFLAGS) -lsla
+	@echo "compiled $*"
+	@echo ""
+
+$(DESTBIN)/mpcorb_predict: $(BIN)/mpcorb_predict.$(ARCH)
+	@if [ ! -d $(DESTBIN) ]; then mkdir -p $(DESTBIN); fi
+	rm -f $(DESTBIN)/mpcorb_predict
+	cp $(BIN)/mpcorb_predict.$(ARCH) $(DESTBIN)/mpcorb_predict
+	@echo "installed mpcorb_predict"
+	@echo ""
+
+$(DESTLIB)/libsla.a:
+	echo "missing SLALIB libsla.a: install in $(DESTLIB)"
+
+$(DESTINC)/slalib.h:
+	echo "missing SLALIB slalib.h: install in $(DESTINC)"
+
+mpcorb_predict: $(DESTBIN)/mpcorb_predict $(DESTLIB)/libsla.a $(DESTINC)/slalib.h
+
 %.clean:
 	rm -f $(SRC)/$*.$(ARCH).o
Index: /trunk/Ohana/src/tools/src/mpcorb_predict.c
===================================================================
--- /trunk/Ohana/src/tools/src/mpcorb_predict.c	(revision 41684)
+++ /trunk/Ohana/src/tools/src/mpcorb_predict.c	(revision 41684)
@@ -0,0 +1,915 @@
+# include <ohana.h>
+# include <gfitsio.h>
+# include <convert.h>
+# include <slalib.h>
+
+/* convert an MPCORB.DAT file to a useful format that can be loaded by mknight */
+
+/* The MPCORB file is NOT fixed bytes/row.  We need to handle fractional lines at the end
+ * of each read block.
+ */
+
+/* read in chunks of ~64MB */
+# define NBYTE 0x400000
+# define DEBUG 0
+# define VERBOSE 1
+# define TESTING 0
+
+typedef struct {
+  char   ID[8];
+  double epoch;
+  double mean_anomaly;
+  double perihelion;
+  double ascend_node;
+  double inclination;
+  double eccentricity;
+  double semimajor_a;
+  double mjdMin;
+  double mjdMax;
+  double RatMin;
+  double RatMax;
+  double DatMin;
+  double DatMax;
+} Planets;
+
+int Shutdown (char *format, ...);
+int mpcorb_parseline (char *line, Planets *planet, int Nmax);
+Planets *mpcorb_read_text (char *filename, int *nplanets);
+Planets *mpcorb_read_fits (char *filename, int *nplanets);
+void     mpcorb_save_fits (char *filename, Planets *planets, int Nplanets);
+
+void mpcorb_predict (Planets *planet, double mjdObs, double *Robs, double *Dobs, int FullCalc);
+
+void mpcorb_trange (int argc, char **argv);
+void mpcorb_moment (int argc, char **argv);
+
+void mpcorb_testing (int argc, char **argv);
+
+void mpcorb_predict_test ();
+void mpcorb_predict_test_calc (Planets *planet, double mjdObs, double *Robs, double *Dobs);
+
+# define PS1_LONGITUDE -156.2559041668965 /* degrees East */
+# define PS1_LATITUDE    20.7070994446013 /* degrees North */
+# define PS1_ALTITUDE  3067.7 /* meters */
+
+int main (int argc, char **argv) {
+
+  if (TESTING) mpcorb_testing (argc, argv);
+
+  if ((argc != 9) && (argc != 10)) goto usage;
+
+  if (!strcasecmp (argv[1], "trange")) mpcorb_trange (argc, argv);
+  if (!strcasecmp (argv[1], "moment")) mpcorb_moment (argc, argv);
+
+usage:
+  fprintf (stderr, "USAGE: %s trange (MPCORB.DAT) (Tmin) (Tmax) (Rmin) (Rmax) (Dmin) (Dmax) (range.fits)\n", argv[0]);
+  fprintf (stderr, "USAGE: %s moment (range.fits) (Tobs) (Rmin) (Rmax) (Dmin) (Dmax) (output)\n", argv[0]);
+
+  fprintf (stderr, " trange : generate a table of asteroid orbital elements which fall within the region at Tmin or Tmax (in MJD)\n");
+  fprintf (stderr, " moment : generate a table of asteroid positions which fall within the region at Tobs (in MJD)\n");
+
+  exit (2);
+}
+
+// generate a table of asteroid orbital elements which fall within the region at Tmin or Tmax (in MJD)
+void mpcorb_trange (int argc, char **argv) {
+
+  if (argc != 10) Shutdown ("USAGE: %s trange (MPCORB.DAT) (Tmin) (Tmax) (Rmin) (Rmax) (Dmin) (Dmax) (range.fits)\n", argv[0]);
+
+  char  *filename = argv[2];
+  double mjdMin   = atof (argv[3]);
+  double mjdMax   = atof (argv[4]);
+  double Rmin     = atof (argv[5]);
+  double Rmax     = atof (argv[6]);
+  double Dmin     = atof (argv[7]);
+  double Dmax     = atof (argv[8]);
+  char  *output   = argv[9];
+
+  int Nplanets = 0;
+  Planets *planets = mpcorb_read_text (filename, &Nplanets);
+  if (VERBOSE) fprintf (stderr, "loaded %d planets\n", Nplanets);
+
+  int NplanetsSave = 0;
+  int NPLANETSSAVE = 1000;
+  ALLOCATE_PTR (planetsSave, Planets, NPLANETSSAVE);
+
+  // transform all objects and identify those in the target region:
+  for (int i = 0; i < Nplanets; i++) {
+    if (i % 1000 == 0) fprintf (stderr, ".");
+
+    // is the object in the region for the start of the period?
+    double RatMin, DatMin;
+    mpcorb_predict (&planets[i], mjdMin, &RatMin, &DatMin, FALSE);
+    int inRangeMin = TRUE;
+    inRangeMin = inRangeMin && (RatMin > Rmin);
+    inRangeMin = inRangeMin && (RatMin < Rmax);
+    inRangeMin = inRangeMin && (DatMin > Dmin);
+    inRangeMin = inRangeMin && (DatMin < Dmax);
+    
+    // is the object in the region for the end of the period?
+    double RatMax, DatMax;
+    mpcorb_predict (&planets[i], mjdMax, &RatMax, &DatMax, FALSE);
+    int inRangeMax = TRUE;
+    inRangeMax = inRangeMax && (RatMax > Rmin);
+    inRangeMax = inRangeMax && (RatMax < Rmax);
+    inRangeMax = inRangeMax && (DatMax > Dmin);
+    inRangeMax = inRangeMax && (DatMax < Dmax);
+
+    // NOTE: if an object is moving so fast that it traverses the entire region, we will
+    // miss it.  I'm going to ignore this case for now.
+    if (!inRangeMin && !inRangeMax) continue;
+
+    planetsSave[NplanetsSave] = planets[i];
+    planetsSave[NplanetsSave].mjdMin = mjdMin;
+    planetsSave[NplanetsSave].RatMin = RatMin;
+    planetsSave[NplanetsSave].DatMin = DatMin;
+    planetsSave[NplanetsSave].mjdMax = mjdMax;
+    planetsSave[NplanetsSave].RatMax = RatMax;
+    planetsSave[NplanetsSave].DatMax = DatMax;
+    
+    NplanetsSave++;
+    CHECK_REALLOCATE (planetsSave, Planets, NPLANETSSAVE, NplanetsSave, 1000);
+  }
+  fprintf (stderr, "\n");
+
+  mpcorb_save_fits (output, planetsSave, NplanetsSave);
+
+  exit (0);
+}
+
+// generate a table of asteroid positions which fall within the region at Tobs (in MJD)
+void mpcorb_moment (int argc, char **argv) {
+
+  if (argc != 9) Shutdown ("USAGE: %s moment (range.fits) (Tobs) (Rmin) (Rmax) (Dmin) (Dmax) (output)\n", argv[0]);
+
+  char  *filename = argv[2];
+  double mjdObs   = atof (argv[3]);
+  double Rmin     = atof (argv[4]);
+  double Rmax     = atof (argv[5]);
+  double Dmin     = atof (argv[6]);
+  double Dmax     = atof (argv[7]);
+  char  *output   = argv[8];
+
+  FILE *fout = fopen (output, "w");
+  if (!fout) Shutdown ("ERROR: unable to open file for output %s\n", output);
+
+  int Nplanets = 0;
+  Planets *planets = mpcorb_read_fits (filename, &Nplanets);
+  if (VERBOSE) fprintf (stderr, "loaded %d planets\n", Nplanets);
+
+  // transform all objects and identify those in the target region:
+  for (int i = 0; i < Nplanets; i++) {
+    if (i % 1000 == 0) fprintf (stderr, ".");
+
+    // is the object in the region for the start of the period?
+    // first, use the fast, but inaccurate, calculation
+    double Robs, Dobs;
+    mpcorb_predict (&planets[i], mjdObs, &Robs, &Dobs, FALSE);
+    if (Robs < Rmin) continue;
+    if (Robs > Rmax) continue;
+    if (Dobs < Dmin) continue;
+    if (Dobs > Dmax) continue;
+    
+    // if it is in the region, use the more accurate calculation
+    mpcorb_predict (&planets[i], mjdObs, &Robs, &Dobs, TRUE);
+    fprintf (fout, "%8s %12.6f %12.6f\n", planets[i].ID, Robs, Dobs);
+  }
+  fclose (fout);
+  fprintf (stderr, "\n");
+
+  exit (0);
+}
+
+// STATUS is value expected for success
+# define CHECK_STATUS(STATUS,MSG,...) if (!(STATUS)) { Shutdown (MSG, __VA_ARGS__); }
+
+// write minor planet orbital elements (plus (R,D,T)[min,max]) to a FITS table
+void mpcorb_save_fits (char *filename, Planets *planets, int Nplanets) {
+
+  Header header;
+  Header theader;
+  Matrix matrix;
+  FTable ftable;
+
+  gfits_init_header (&header);
+  header.extend = TRUE;
+  gfits_create_header (&header);
+  gfits_create_matrix (&header, &matrix);
+
+  gfits_create_table_header (&theader, "BINTABLE", "MPCORB");
+
+  gfits_define_bintable_column (&theader, "8A", "ID",           "name",         "none",     1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "EPOCH",        "epoch",        "MJD",      1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "MEAN_ANOMALY", "mean_anomaly", "degrees",  1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "PERIHELION",   "perihelion",   "degrees",  1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "ASCEND_NODE",  "ascend_node",  "degrees",  1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "INCLINATION",  "inclination",  "degrees",  1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "ECCENTRICITY", "eccentricity", "unitless", 1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "SEMIMAJOR_A",  "semimajor_a",  "AU",       1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "MJD_MIN",      "mjdMin",       "MJD",      1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "MJD_MAX",      "mjdMax",       "MJD",      1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "RAT_MIN",      "RatMin",       "degrees",  1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "RAT_MAX",      "RatMax",       "degrees",  1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "DAT_MIN",      "DatMin",       "degrees",  1.0, 0.0);
+  gfits_define_bintable_column (&theader,  "D", "DAT_MAX",      "DatMax",       "degrees",  1.0, 0.0);
+
+  // generate the output array that carries the data
+  gfits_create_table (&theader, &ftable);
+
+  // create intermediate storage arrays
+  ALLOCATE_PTR (ID,                     char,   Nplanets*8);
+  ALLOCATE_PTR (epoch,                  double, Nplanets);
+  ALLOCATE_PTR (mean_anomaly,	        double, Nplanets);
+  ALLOCATE_PTR (perihelion,	        double, Nplanets);
+  ALLOCATE_PTR (ascend_node,	        double, Nplanets);
+  ALLOCATE_PTR (inclination,	        double, Nplanets);
+  ALLOCATE_PTR (eccentricity,	        double, Nplanets);
+  ALLOCATE_PTR (semimajor_a,	        double, Nplanets);
+  ALLOCATE_PTR (mjdMin, 	        double, Nplanets);
+  ALLOCATE_PTR (mjdMax, 	        double, Nplanets);
+  ALLOCATE_PTR (RatMin, 	        double, Nplanets);
+  ALLOCATE_PTR (RatMax, 	        double, Nplanets);
+  ALLOCATE_PTR (DatMin, 	        double, Nplanets);
+  ALLOCATE_PTR (DatMax, 	        double, Nplanets);
+
+  // set intermediate storage arrays
+  for (int i = 0; i < Nplanets; i++) {
+    epoch[i]           = planets[i].epoch;
+    mean_anomaly[i]    = planets[i].mean_anomaly;
+    perihelion[i]      = planets[i].perihelion;
+    ascend_node[i]     = planets[i].ascend_node;
+    inclination[i]     = planets[i].inclination;
+    eccentricity[i]    = planets[i].eccentricity;
+    semimajor_a[i]     = planets[i].semimajor_a;
+    mjdMin[i]          = planets[i].mjdMin;
+    mjdMax[i]          = planets[i].mjdMax;
+    RatMin[i]          = planets[i].RatMin;
+    RatMax[i]          = planets[i].RatMax;
+    DatMin[i]          = planets[i].DatMin;
+    DatMax[i]          = planets[i].DatMax;
+
+    memcpy(&ID[i*8], planets[i].ID, 8);
+  }
+
+  // add the columns to the output array
+  gfits_set_bintable_column (&theader, &ftable, "ID",           ID,           Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "EPOCH",        epoch,        Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "MEAN_ANOMALY", mean_anomaly, Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "PERIHELION",   perihelion,   Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "ASCEND_NODE",  ascend_node,  Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "INCLINATION",  inclination,  Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "ECCENTRICITY", eccentricity, Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "SEMIMAJOR_A",  semimajor_a,  Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "MJD_MIN",      mjdMin,       Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "MJD_MAX",      mjdMax,       Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "RAT_MIN",      RatMin,       Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "RAT_MAX",      RatMax,       Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "DAT_MIN",      DatMin,       Nplanets);
+  gfits_set_bintable_column (&theader, &ftable, "DAT_MAX",      DatMax,       Nplanets);
+
+  // free intermediate storage arrays
+  free (ID);
+  free (epoch);
+  free (mean_anomaly);
+  free (perihelion);
+  free (ascend_node);
+  free (inclination);
+  free (eccentricity);
+  free (semimajor_a);
+  free (mjdMin);
+  free (mjdMax);
+  free (RatMin);
+  free (RatMax);
+  free (DatMin);
+  free (DatMax);
+
+  // open the output file
+  FILE *f = fopen (filename, "w");
+  if (!f) Shutdown ("ERROR: cannot open mpcorb file for output %s\n", filename);
+
+  int status;
+  status = gfits_fwrite_header  (f, &header);
+  CHECK_STATUS (status, "ERROR: cannot write header for mpcorbs %s\n", filename);
+
+  status = gfits_fwrite_matrix  (f, &matrix);
+  CHECK_STATUS (status, "ERROR: cannot write matrix for mpcorbs %s\n", filename);
+
+  status = gfits_fwrite_Theader (f, &theader);
+  CHECK_STATUS (status, "ERROR: cannot write table header for mpcorbs %s\n", filename);
+
+  status = gfits_fwrite_table  (f, &ftable);
+  CHECK_STATUS (status, "ERROR: cannot write table data for mpcorbs %s\n", filename);
+
+  gfits_free_header (&header);
+  gfits_free_matrix (&matrix);
+  gfits_free_header (&theader);
+  gfits_free_table (&ftable);
+
+  int fd = fileno (f);
+
+  status = fflush (f);
+  CHECK_STATUS (!status, "ERROR: cannot flush file mpcorbs %s\n", filename);
+
+  status = fsync (fd);
+  CHECK_STATUS (!status, "ERROR: cannot flush file mpcorbs %s\n", filename);
+
+  status = fclose (f);
+  CHECK_STATUS (!status, "ERROR: problem closing mpcorbs file %s\n", filename);
+
+  return;
+}
+
+# define GET_COLUMN(OUT,NAME,TYPE) \
+  TYPE *OUT = gfits_get_bintable_column_data (&theader, &ftable, NAME, type, &Nrow, &Ncol); \
+  myAssert (!strcmp(type, #TYPE), "wrong column type");
+
+// load planet data from a fits table
+Planets *mpcorb_read_fits (char *filename, int *nplanets) {
+
+  int Ncol;
+  off_t Nrow;
+  Header header;
+  Header theader;
+  Matrix matrix;
+  FTable ftable;
+
+  FILE *f = fopen (filename, "r");
+  if (!f) Shutdown ("ERROR: cannot open mpcorb fits file %s\n", filename);
+
+  // load in PHU segment (ignore)
+  if (!gfits_fread_header (f, &header)) Shutdown ("can't read mpcorb fits file header %s\n", filename);
+  if (!gfits_fread_matrix (f, &matrix, &header)) Shutdown ("can't read mpcorb fits file matrix %s\n", filename);
+  ftable.header = &theader;
+
+  // load data for this header 
+  if (!gfits_load_header (f, &theader)) Shutdown ("can't read mpcorb fits file table header %s\n", filename);
+  if (!gfits_fread_ftable_data (f, &ftable, FALSE)) Shutdown ("can't read mpcorb fits file table data %s\n", filename);
+
+  // we read the entire block of data, then extract the columns, then set the Planets structure values.
+  // I free the FITS table data after extracting the colums to avoid having 3 copies in memory.
+
+  char type[16]; // used by the functions called by GET_COLUMN
+
+  GET_COLUMN (ID,           "ID",            char);
+  GET_COLUMN (epoch,        "EPOCH",         double);
+  GET_COLUMN (mean_anomaly, "MEAN_ANOMALY",  double);
+  GET_COLUMN (perihelion,   "PERIHELION",    double);
+  GET_COLUMN (ascend_node,  "ASCEND_NODE",   double);
+  GET_COLUMN (inclination,  "INCLINATION",   double);
+  GET_COLUMN (eccentricity, "ECCENTRICITY",  double);
+  GET_COLUMN (semimajor_a,  "SEMIMAJOR_A",   double);
+  GET_COLUMN (mjdMin,       "MJD_MIN",       double);
+  GET_COLUMN (mjdMax,       "MJD_MAX",       double);
+  GET_COLUMN (RatMin,       "RAT_MIN",       double);
+  GET_COLUMN (RatMax,       "RAT_MAX",       double);
+  GET_COLUMN (DatMin,       "DAT_MIN",       double);
+  GET_COLUMN (DatMax,       "DAT_MAX",       double);
+
+  // free the memory associated with the FITS files
+  gfits_free_header (&header);
+  gfits_free_matrix (&matrix);
+  gfits_free_header (&theader);
+  gfits_free_table (&ftable);
+
+  ALLOCATE_PTR (planets, Planets, Nrow);
+  for (int i = 0; i < Nrow; i++) {
+    planets[i].epoch        = epoch[i];
+    planets[i].mean_anomaly = mean_anomaly[i];
+    planets[i].perihelion   = perihelion[i];
+    planets[i].ascend_node  = ascend_node[i];
+    planets[i].inclination  = inclination[i];
+    planets[i].eccentricity = eccentricity[i];
+    planets[i].semimajor_a  = semimajor_a[i];
+    planets[i].mjdMin       = mjdMin[i];
+    planets[i].mjdMax       = mjdMax[i];
+    planets[i].RatMin       = RatMin[i];
+    planets[i].RatMax       = RatMax[i];
+    planets[i].DatMin       = DatMin[i];
+    planets[i].DatMax       = DatMax[i];
+
+    memcpy(planets[i].ID, &ID[i*8], 8);
+    planets[i].ID[7] = 0;
+  }
+  fprintf (stderr, "loaded data for %lld planets\n", (long long) Nrow);
+
+  free (ID          );
+  free (epoch       );
+  free (mean_anomaly);
+  free (perihelion  );
+  free (ascend_node );
+  free (inclination );
+  free (eccentricity);
+  free (semimajor_a );
+  free (mjdMin      );
+  free (mjdMax      );
+  free (RatMin      );
+  free (RatMax      );
+  free (DatMin      );
+  free (DatMax      );
+
+  *nplanets = Nrow;
+  return planets;
+}
+
+void mpcorb_testing (int argc, char **argv) {
+
+  if (argc != 3) Shutdown ("USAGE: %s (MPCORB.DAT) (epoch)\n", argv[0]);
+
+  char  *filename = argv[1];
+  double mjdObs   = atof (argv[2]);
+
+  int Nplanets = 0;
+  Planets *planets = mpcorb_read_text (filename, &Nplanets);
+  if (VERBOSE) fprintf (stderr, "loaded %d planets\n", Nplanets);
+
+  if (TESTING == 1) { mpcorb_predict_test(planets, Nplanets); exit (0); }
+
+  if (TESTING == 2) {
+    double Robs, Dobs;
+    for (int i = 0; i < 10; i++) {
+      mpcorb_predict (&planets[i], mjdObs, &Robs, &Dobs, TRUE);
+      char ra_string[64], de_string[64];
+      hms_format (ra_string, 64, Robs/15.0);
+      hms_format (de_string, 64, Dobs);
+      fprintf (stderr, "result: %12.6f %12.6f : %s %s\n", Robs, Dobs, ra_string, de_string);
+    }
+    exit (0);
+  }
+  fprintf (stderr, "ERROR: invalid value of TESTING : %d\n", TESTING);
+  exit (2);
+}
+    
+// Supply this test with the full list of minor planets from MPCORB.DAT.
+// It will compare the positions for the first 10 minor planets with the
+// reference data downloaded from the MPC for:
+// MJD = 59396.0 @ (Long,Lat,Alt) = (0.0, 0.0, 0.0)
+
+void mpcorb_predict_test (Planets *planets, int Nplanets) {
+  
+  int    Nref = 10;
+  double Rref[10] = { 49.14708333, 352.90125000, 251.95416667, 175.74666667, 266.25541667, 298.38000000,  52.84458333, 121.94916667, 192.29500000, 144.20583333 };
+  double Dref[10] = { 11.48944444,   8.68722222,  -4.19666667,   9.37111111, -16.97583333,  -8.30500000,  22.13472222,  21.50888889,  -1.52000000,  11.57138889 };
+
+  for (int i = 0; (i < Nref) && (i < Nplanets); i++) {
+
+    if (VERBOSE) fprintf (stderr, "%8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.7f %12.7f\n",
+	     planets[i].ID,
+	     planets[i].epoch,
+	     planets[i].mean_anomaly,
+	     planets[i].perihelion,
+	     planets[i].ascend_node,
+	     planets[i].inclination,
+	     planets[i].eccentricity,
+	     planets[i].semimajor_a);
+    
+    // generate predictions for (Long,Lat,Alt) = (0.0, 0.0, 0.0)
+    double Rpre, Dpre;
+    mpcorb_predict_test_calc (&planets[i], 59396.0, &Rpre, &Dpre);
+
+    fprintf (stderr, "result: %8s : %12.6f %12.6f : %6.2f %6.2f\n", planets[i].ID, Rpre, Dpre, 3600*(Rref[i]-Rpre)*cos(Dref[i]*RAD_DEG), 3600*(Dref[i]-Dpre));
+  }
+}
+
+// attempt to predict as a test for Long,Lat = 0,0, with all other corrections
+// this attempt uses slaPertel to calculate perturbed versions of the orbital elements
+// NOTE this test seems to have worked: the calculated positions of the first 10 asteroids reported
+// by MPC match to within 2 arcsec.  This is much larger than the Planetary abberation,
+// so slaPlante must already be accounting for that effect.
+
+void mpcorb_predict_test_calc (Planets *planet, double mjdObs, double *Robs, double *Dobs) {
+
+  int jstat;
+  
+  double dTT = slaDtt (mjdObs);
+  double mjdObsTT = mjdObs + dTT / 86400.0;
+
+  Planets tmpPlanet;
+
+  // correct the orbital elements to the perturbed version
+  slaPertel (2,
+	     planet->epoch,
+	     mjdObsTT,
+	     planet->epoch,
+	     planet->inclination*RAD_DEG,
+	     planet->ascend_node*RAD_DEG,
+	     planet->perihelion*RAD_DEG,
+	     planet->semimajor_a,
+	     planet->eccentricity,
+	     planet->mean_anomaly*RAD_DEG,
+	     &tmpPlanet.epoch,
+	     &tmpPlanet.inclination,  // returned as radians and used as radians below
+	     &tmpPlanet.ascend_node,  // returned as radians and used as radians below
+	     &tmpPlanet.perihelion,   // returned as radians and used as radians below
+	     &tmpPlanet.semimajor_a,
+	     &tmpPlanet.eccentricity,
+	     &tmpPlanet.mean_anomaly, // returned as radians and used as radians below
+	     &jstat);
+
+  // this function returns the topocentric apparent R,D.  I assume it
+  // has no atm corrections.
+
+  double Rtopo, Dtopo, dist;
+  slaPlante (mjdObsTT,
+	     0.0, // PS1_LONGITUDE*RAD_DEG,
+	     0.0, // PS1_LATITUDE*RAD_DEG,
+	     2,
+	     tmpPlanet.epoch,
+	     tmpPlanet.inclination,
+	     tmpPlanet.ascend_node,
+	     tmpPlanet.perihelion,
+	     tmpPlanet.semimajor_a,
+	     tmpPlanet.eccentricity,
+	     tmpPlanet.mean_anomaly,
+	     0.0, &Rtopo, &Dtopo, &dist, &jstat);
+
+  if (VERBOSE) fprintf (stderr, "Plante   : %12.6f %12.6f\n", Rtopo*DEG_RAD, Dtopo*DEG_RAD);
+
+  // NOTE: the atm effects seem to be large and incompatible with MPC outputs.
+  // Setting the pressure to 0.0 seems to avoid those corrections.
+
+  double Rap, Dap;
+  slaOap ("R",
+	  Rtopo,
+	  Dtopo,
+	  mjdObs,
+	  -0.4,
+	  0.0, // LONGITUDE
+	  0.0, // LATITUDE
+	  0.0, // ALTITUDE,
+	  0.0, // xp
+	  0.0, // yp
+	  273.15, // T(K)
+	  0.0, // P(mB) p = 1013.25 * exp ( -hm / ( 29.3 * tsl ) )
+	  0.0, // RH(%)
+	  0.7, // Wavelength (um)
+	  0.0, // lapse rate
+	  &Rap, &Dap);
+
+  if (VERBOSE) fprintf (stderr, "Oap   : %12.6f %12.6f : %6.2f %6.2f\n", Rap*DEG_RAD, Dap*DEG_RAD, 3600*(Rap-Rtopo)*DEG_RAD*cos(Dtopo), 3600*(Dap-Dtopo)*DEG_RAD);
+
+  double Rmean, Dmean;
+  slaAmp (Rap, Dap, mjdObsTT, 2000.0, &Rmean, &Dmean);
+
+  if (VERBOSE) fprintf (stderr, "Amp   : %12.6f %12.6f : %6.2f %6.2f\n", Rmean*DEG_RAD, Dmean*DEG_RAD, 3600*(Rmean-Rap)*DEG_RAD*cos(Dap), 3600*(Dmean-Dap)*DEG_RAD);
+
+  if (0) {
+    // slaAmp applies both precession/nutation as well as abberation
+    // but, I think abberation is not appropriate to solar system bodies.
+
+    // here, I instead calculate just the precession and nutation with slalib
+
+    double Rpre, Dpre;
+    Rpre = Rap;
+    Dpre = Dap;
+    slaPreces ("FK5", 2021.5, 2000.000000, &Rpre, &Dpre);
+    if (VERBOSE) fprintf (stderr, "Prec     : %12.6f %12.6f : %6.2f %6.2f\n", Rpre*DEG_RAD, Dpre*DEG_RAD, 3600*(Rpre-Rap)*DEG_RAD*cos(Dap), 3600*(Dpre-Dap)*DEG_RAD);
+  }
+
+  *Robs = Rmean * DEG_RAD;
+  *Dobs = Dmean * DEG_RAD;
+}
+
+// mjdObs is NOT in TT
+void mpcorb_predict (Planets *planet, double mjdObs, double *Robs, double *Dobs, int FullCalc) {
+
+  int jstat;
+  
+  double dTT = slaDtt (mjdObs);
+  double mjdObsTT = mjdObs + dTT / 86400.0;
+
+  Planets tmpPlanet;
+
+  // correct the orbital elements to the perturbed version
+  if (FullCalc) { slaPertel (2,
+	     planet->epoch,
+	     mjdObsTT,
+	     planet->epoch,
+	     planet->inclination*RAD_DEG,
+	     planet->ascend_node*RAD_DEG,
+	     planet->perihelion*RAD_DEG,
+	     planet->semimajor_a,
+	     planet->eccentricity,
+	     planet->mean_anomaly*RAD_DEG,
+	     &tmpPlanet.epoch,
+	     &tmpPlanet.inclination,  // returned as radians and used as radians below
+	     &tmpPlanet.ascend_node,  // returned as radians and used as radians below
+	     &tmpPlanet.perihelion,   // returned as radians and used as radians below
+	     &tmpPlanet.semimajor_a,
+	     &tmpPlanet.eccentricity,
+	     &tmpPlanet.mean_anomaly, // returned as radians and used as radians below
+	     &jstat);
+  } else {
+    tmpPlanet = *planet;
+  }
+
+  // this function returns the topocentric apparent R,D.  I assume it
+  // has no atm corrections.
+
+  double Rtopo, Dtopo, dist;
+  slaPlante (mjdObsTT,
+	     PS1_LONGITUDE*RAD_DEG,
+	     PS1_LATITUDE*RAD_DEG,
+	     2,
+	     planet->epoch, 
+	     planet->inclination*RAD_DEG,
+	     planet->ascend_node*RAD_DEG,
+	     planet->perihelion*RAD_DEG,
+	     planet->semimajor_a,
+	     planet->eccentricity,
+	     planet->mean_anomaly*RAD_DEG,
+	     0.0, &Rtopo, &Dtopo, &dist, &jstat);
+
+  // NOTE: the atm effects seem to be large and incompatible with MPC outputs.
+  // Setting the pressure to 0.0 seems to avoid those corrections.
+
+  double Rap, Dap;
+  if (FullCalc) {
+    slaOap ("R",
+	  Rtopo,
+	  Dtopo,
+	  mjdObs,
+	  -0.4,
+	  PS1_LONGITUDE*RAD_DEG,
+	  PS1_LATITUDE*RAD_DEG,
+	  PS1_ALTITUDE,
+	  0.0, // xp
+	  0.0, // yp
+	  273.15, // T(K)
+	  0.0, // P(mB) p = 1013.25 * exp ( -hm / ( 29.3 * tsl ) )
+	  0.0, // RH(%)
+	  0.7, // Wavelength (um)
+	  0.0065, // lapse rate
+	  &Rap, &Dap);
+  } else { Rap = Rtopo; Dap = Dtopo; }
+
+  double Rmean, Dmean;
+  if (FullCalc) {
+    slaAmp (Rap, Dap, mjdObsTT, 2000.0, &Rmean, &Dmean);
+  } else {Rmean = Rap; Dmean = Dap; }
+
+  *Robs = Rmean * DEG_RAD;
+  *Dobs = Dmean * DEG_RAD;
+}
+
+# define MPC_DOUBLE(FIELD,START,END) {FIELD = strtod(&line[START-1], NULL);}
+
+// part the data in a single line of the MPCORB.DAT file
+int mpcorb_parseline (char *line, Planets *planet, int Nmax) {
+
+  // fields I need:
+
+  // epoch       : 21 - 25
+  // inclination : 60 - 68
+  // longitude of ascending node (J2000) : 49 - 57
+  // argument of perihelion (J2000) : 38 - 46
+  // semimajor axis (AU) : 93 - 103
+  // eccentricity : 71 - 79
+  // Mean anomaly : 27 - 35
+  
+  // these fields are saved in the data file with whitespace between, so I can
+  // set the byte after the range to NULL and scanf the values.  or just use strtod as needed:
+
+  strncpy (planet->ID, line, 7); planet->ID[7] = 0;
+
+  MPC_DOUBLE (planet->mean_anomaly, 27, 35);
+  MPC_DOUBLE (planet->perihelion,   38, 46);
+  MPC_DOUBLE (planet->ascend_node,  49, 57);
+  MPC_DOUBLE (planet->inclination,  60, 68);
+  MPC_DOUBLE (planet->eccentricity, 71, 79);
+  MPC_DOUBLE (planet->semimajor_a,  93, 103);
+
+  // the epoch is stored in a compressed format:
+  int year = 1800 + 100*(line[20] - 'I') + 10*(line[21] - '0') + line[22] - '0';
+  char mchar = line[23];
+  int month = (mchar > 64) ? mchar - 'A' + 10 : mchar - '0';
+  char dchar = line[24];
+  int day = (dchar > 64) ? dchar - 'A' + 10 : dchar - '0';
+
+  double jd = day - 32075 + (int)(1461*(year + 4800 + (int)((month-14)/12))/4)
+    + (int)(367*(month - 2 - (int)((month - 14)/12)*12)/12)
+    - (int)(3*(int)((year + 4900 + (int)((month - 14)/12))/100)/4) - 0.5;
+
+  planet->epoch = jd - 2400000.5;;
+  return TRUE;
+}
+
+// read the full MPCORB file, storing desired information in Planets structure
+Planets *mpcorb_read_text (char *filename, int *nplanets) {
+
+  ALLOCATE_PTR (buffer, char, NBYTE);
+
+  // scan through the entire MPCORB file
+  FILE *f = fopen (filename, "r");
+  if (f == NULL) Shutdown ("can't open data file: %s", filename);
+
+  int Nbyte  = 0;  // number of bytes read on each pass
+  int Nextra = 0;  // number excess bytes from last partial row
+  int Nline  = 0;  // number of lines read so far in file
+
+  int Nplanets = 0;
+  int NPLANETS = 10000;
+  ALLOCATE_PTR (planets, Planets, NPLANETS);
+
+  while ((Nbyte = fread (&buffer[Nextra], 1, NBYTE-Nextra, f)) != 0) {
+    if (Nbyte == -1) Shutdown ("error reading from data file %s", filename);
+    if (DEBUG) fprintf (stderr, "read %d bytes", Nbyte);
+
+    Nbyte += Nextra;
+
+    if (VERBOSE) fprintf (stderr, "read .. ");
+
+    /* find bounds on first complete line */
+    char *p = buffer;
+    char *q = memchr (p, '\n', Nbyte);
+    if (q == NULL) Shutdown ("incomplete line at end of file\n");
+    int offset = p - buffer; // offset within this scan
+
+    // scan through entire buffer for star coords
+    while (1) {
+
+      Nline ++;
+      if ((Nline > 43) && (q - p > 5)) {
+	mpcorb_parseline (p, &planets[Nplanets], Nbyte - offset);
+	Nplanets ++;
+	CHECK_REALLOCATE (planets, Planets, NPLANETS, Nplanets, 10000);
+      }
+
+      /* start of the next line */
+      p = q + 1;
+      offset = p - buffer; // offset within this scan
+      if (offset == Nbyte) {
+	// last buffer read is a complete line
+	Nextra = 0;
+	break;
+      }
+
+      /* end of the next line */
+      q = memchr (p, '\n', Nbyte - offset);
+      if (q == NULL) {
+	// last, incomplete line in buffer
+	Nextra = Nbyte - offset;
+	break;
+      } 
+    }
+
+    if (VERBOSE) fprintf (stderr, "scan %d planets (%s).. ", Nplanets, planets[Nplanets-1].ID);
+
+    // at end, p points at the start of last, partial line
+    if (Nextra) memmove (buffer, p, Nextra);
+
+    // TESTING:
+    if (TESTING && (Nplanets > 50)) goto skip;
+
+    // XX temporary:
+    // if (Nplanets > 100000) goto skip;
+
+  }
+  if (VERBOSE) fprintf (stderr, "\n");
+  
+skip:
+
+  fclose (f);
+  free (buffer);
+
+  *nplanets = Nplanets;
+  return planets;
+}
+
+// print an error message and exit
+int Shutdown (char *format, ...) {  
+  va_list argp;
+  char formatplus[1024];
+  
+  strcpy (formatplus, format);
+  strcat (formatplus, "\n");
+
+  va_start (argp, format);
+  vfprintf (stderr, formatplus, argp);
+  va_end (argp);
+
+  fprintf (stderr, "ERROR: mpcorb_convert halted\n");
+  exit (1);
+}
+
+// The following test case provided by Serge can be perfectly reproduced, but only if I
+// use Dtt for the epoch, and this is clearly wrong.
+void mpcorb_predict_test_from_serge () {
+
+  double epochUT = 57435.5;
+
+  Planets testPlanet;
+  // testPlanet.epoch        = epochUT + slaDtt(epochUT) / 86400.0; // XXX this is what the epoch should be
+  testPlanet.epoch        = slaDtt(epochUT); // XXX this is the value that works with the test values
+  testPlanet.inclination  = 0.043306854729735354 * DEG_RAD;
+  testPlanet.ascend_node  = 4.728915154005972    * DEG_RAD;
+  testPlanet.perihelion   = 2.2756608097151414   * DEG_RAD;
+  testPlanet.mean_anomaly = 2.888163341284418    * DEG_RAD;
+  testPlanet.eccentricity = 0.09528339999999975;
+  testPlanet.semimajor_a  = 2.782049599999998;
+		
+  double Robs, Dobs;
+  mpcorb_predict (&testPlanet, epochUT, &Robs, &Dobs, TRUE);
+
+  fprintf (stderr, "%f vs %f : %f\n",  3.92105439135726600*DEG_RAD, Robs, 3600*(3.92105439135726600*DEG_RAD - Robs));
+  fprintf (stderr, "%f vs %f : %f\n", -0.33499850519808244*DEG_RAD, Dobs, 3600*(-0.33499850519808244*DEG_RAD - Dobs));
+  exit (0);
+}
+
+// The following test code uses unperturbed orbital elements and gets the
+// wrong answer by 20 - 60 arcsec for the top 10 minor planets.
+void mpcorb_predict_t0 (Planets *planet, double mjdObs, double *Robs, double *Dobs) {
+
+  int jstat;
+  
+  double dTT = slaDtt (mjdObs);
+  double mjdObsTT = mjdObs + dTT / 86400.0;
+
+  // this function returns the topocentric apparent R,D.  I assume it
+  // has no atm corrections.
+
+  double Rrad, Drad, dist;
+  slaPlante (mjdObsTT,
+	     0.0, // PS1_LONGITUDE*RAD_DEG,
+	     0.0, // PS1_LATITUDE*RAD_DEG,
+	     2,
+	     planet->epoch,
+	     planet->inclination*RAD_DEG,
+	     planet->ascend_node*RAD_DEG,
+	     planet->perihelion*RAD_DEG,
+	     planet->semimajor_a,
+	     planet->eccentricity,
+	     planet->mean_anomaly*RAD_DEG,
+	     0.0, &Rrad, &Drad, &dist, &jstat);
+
+  fprintf (stderr, "Plante   : %12.6f %12.6f\n", Rrad*DEG_RAD, Drad*DEG_RAD);
+
+  // I'm not certain if I should use realistic (STD) atm parameters or not.
+  // If I set them to zero, is the correction not applied?
+ 
+  double Rap, Dap;
+  slaOap ("R",
+	  Rrad,
+	  Drad,
+	  mjdObs,
+	  -0.4,
+	  0.0, // LONGITUDE
+	  0.0, // LATITUDE
+	  0.0, // ALTITUDE,
+	  0.0, // xp
+	  0.0, // yp
+	  273.15, // T(K)
+	  1013.25, // P(mB) p = 1013.25 * exp ( -hm / ( 29.3 * tsl ) )
+	  0.0, // RH(%)
+	  0.7, // Wavelength (um)
+	  0.0065, // lapse rate
+	  &Rap, &Dap);
+
+  // XXX the atm effects seem to be large and incompatible with MPC outputs
+  // fprintf (stderr, "Oap (STD): %12.6f %12.6f\n", Rap*DEG_RAD, Dap*DEG_RAD);
+
+  slaOap ("R",
+	  Rrad,
+	  Drad,
+	  mjdObs,
+	  -0.4,
+	  0.0, // LONGITUDE
+	  0.0, // LATITUDE
+	  0.0, // ALTITUDE,
+	  0.0, // xp
+	  0.0, // yp
+	  273.15, // T(K)
+	  0.0, // P(mB) p = 1013.25 * exp ( -hm / ( 29.3 * tsl ) )
+	  0.0, // RH(%)
+	  0.7, // Wavelength (um)
+	  0.0, // lapse rate
+	  &Rap, &Dap);
+
+  fprintf (stderr, "Oap (0.0): %12.6f %12.6f : %6.2f %6.2f\n", Rap*DEG_RAD, Dap*DEG_RAD, 3600*(Rap-Rrad)*DEG_RAD*cos(Drad), 3600*(Dap-Drad)*DEG_RAD);
+
+  double Rmean, Dmean;
+  slaAmp (Rap, Dap, mjdObsTT, 2000.0, &Rmean, &Dmean);
+
+  fprintf (stderr, "Amp1     : %12.6f %12.6f : %6.2f %6.2f\n", Rmean*DEG_RAD, Dmean*DEG_RAD, 3600*(Rmean-Rap)*DEG_RAD*cos(Dap), 3600*(Dmean-Dap)*DEG_RAD);
+
+  // apply Amp to output from Plante
+  double Ralt, Dalt;
+  slaAmp (Rrad, Drad, mjdObsTT, 2000.0, &Ralt, &Dalt);
+
+  fprintf (stderr, "Amp2     : %12.6f %12.6f : %6.2f %6.2f\n", Ralt*DEG_RAD, Dalt*DEG_RAD, 3600*(Ralt-Rrad)*DEG_RAD*cos(Drad), 3600*(Dalt-Drad)*DEG_RAD);
+
+
+  // slaAmp applies both precession/nutation as well as abberation
+  // but, I think abberation is not appropriate to solar system bodies.
+
+  // here, I instead calculate just the precession and nutation with slalib
+
+  double Rpre, Dpre;
+  Rpre = Rap;
+  Dpre = Dap;
+  slaPreces ("FK5", 2021.5, 2000.000000, &Rpre, &Dpre);
+  fprintf (stderr, "Prec     : %12.6f %12.6f : %6.2f %6.2f\n", Rpre*DEG_RAD, Dpre*DEG_RAD, 3600*(Rpre-Rap)*DEG_RAD*cos(Dap), 3600*(Dpre-Dap)*DEG_RAD);
+
+  *Robs = Rmean * DEG_RAD;
+  *Dobs = Dmean * DEG_RAD;
+}
+