Index: trunk/Ohana/src/relastro/Makefile
===================================================================
--- trunk/Ohana/src/relastro/Makefile	(revision 12067)
+++ trunk/Ohana/src/relastro/Makefile	(revision 12332)
@@ -15,5 +15,5 @@
 FULL_CFLAGS   = $(BASE_CFLAGS)
 FULL_CPPFLAGS = $(BASE_CPPFLAGS)
-FULL_LDFLAGS  = -ldvo -lFITS -lohana $(BASE_LDFLAGS)
+FULL_LDFLAGS  = -lkapa -ldvo -lFITS -lohana $(BASE_LDFLAGS)
 
 relastro: $(BIN)/relastro.$(ARCH)
@@ -21,35 +21,42 @@
 
 RELASTRO = \
-$(SRC)/relastro.$(ARCH).o	 \
-$(SRC)/initialize.$(ARCH).o	 \
 $(SRC)/ConfigInit.$(ARCH).o	 \
-$(SRC)/args.$(ARCH).o		 \
-$(SRC)/liststats.$(ARCH).o	 \
-$(SRC)/load_images.$(ARCH).o	 \
-$(SRC)/name_region.$(ARCH).o	 \
-$(SRC)/find_images.$(ARCH).o	 \
-$(SRC)/find_regions.$(ARCH).o	 \
-$(SRC)/get_regions.$(ARCH).o	 \
-$(SRC)/load_catalogs.$(ARCH).o	 \
-$(SRC)/gcatalog.$(ARCH).o 	 \
-$(SRC)/bcatalog.$(ARCH).o	 \
-$(SRC)/GridOps.v2.$(ARCH).o	 \
-$(SRC)/StarOps.$(ARCH).o	 \
+$(SRC)/FitChip.$(ARCH).o         \
+$(SRC)/FitMosaic.$(ARCH).o       \
+$(SRC)/FitPM.$(ARCH).o           \
+$(SRC)/FitPMandPar.$(ARCH).o     \
+$(SRC)/FitSimple.$(ARCH).o       \
 $(SRC)/ImageOps.$(ARCH).o	 \
 $(SRC)/MosaicOps.$(ARCH).o	 \
-$(SRC)/global_stats.$(ARCH).o	 \
+$(SRC)/ParFactor.$(ARCH).o       \
+$(SRC)/SetSignals.$(ARCH).o 	 \
+$(SRC)/Shutdown.$(ARCH).o 	 \
+$(SRC)/UpdateChips.$(ARCH).o     \
+$(SRC)/UpdateMosaic.$(ARCH).o    \
+$(SRC)/UpdateObjects.$(ARCH).o   \
+$(SRC)/UpdateSimple.$(ARCH).o    \
+$(SRC)/args.$(ARCH).o		 \
+$(SRC)/bcatalog.$(ARCH).o	 \
+$(SRC)/dvo_astrom_ops.$(ARCH).o  \
+$(SRC)/fitpoly.$(ARCH).o         \
+$(SRC)/initialize.$(ARCH).o	 \
+$(SRC)/liststats.$(ARCH).o	 \
+$(SRC)/load_catalogs.$(ARCH).o	 \
+$(SRC)/load_images.$(ARCH).o	 \
+$(SRC)/mkpolyterm.$(ARCH).o      \
 $(SRC)/plot_scatter.$(ARCH).o	 \
 $(SRC)/plotstuff.$(ARCH).o	 \
-$(SRC)/misc.$(ARCH).o		 \
+$(SRC)/select_images.$(ARCH).o	 \
+$(SRC)/sort.$(ARCH).o	         \
+$(SRC)/relastro.$(ARCH).o	 \
 $(SRC)/reload_catalogs.$(ARCH).o \
-$(SRC)/image-db.$(ARCH).o \
-$(SRC)/free_catalogs.$(ARCH).o   \
-$(SRC)/check_permissions.$(ARCH).o   \
+$(SRC)/save_catalogs.$(ARCH).o   \
+$(SRC)/setExclusions.$(ARCH).o 	 \
+$(SRC)/write_coords.$(ARCH).o
+
+DROP = \
 $(SRC)/setMrelFinal.$(ARCH).o 	 \
-$(SRC)/write_coords.$(ARCH).o 	 \
-$(SRC)/setExclusions.$(ARCH).o 	 \
-$(SRC)/Shutdown.$(ARCH).o 	 \
-$(SRC)/SetSignals.$(ARCH).o 	 \
-$(SRC)/wcatalog.$(ARCH).o
+$(SRC)/StarOps.$(ARCH).o	 
+
 
 $(RELASTRO): $(INC)/relastro.h
Index: trunk/Ohana/src/relastro/doc/parallax.pro
===================================================================
--- trunk/Ohana/src/relastro/doc/parallax.pro	(revision 12067)
+++ trunk/Ohana/src/relastro/doc/parallax.pro	(revision 12332)
@@ -332,2 +332,17 @@
 end
 
+
+macro parfactor
+
+  $RA = $1
+  $Dec = $2
+  $e = 23 + 27/60
+  
+  create s 0 360
+  set pR = dcos($e)*d
+
+  set pR =  +(dcos($e)*dsin(s)*dcos($RA) - dcos(s)*dsin($RA))
+  set pD =  -(dcos($e)*dsin(s)*dsin($RA) + dcos(s)*dcos($RA))*dsin($Dec) + dsin($e)*dsin(s)*dcos($Dec)
+
+  lim -1.1 1.1 -1.1 1.1; clear; box; plot pR pD
+end
Index: trunk/Ohana/src/relastro/include/relastro.h
===================================================================
--- trunk/Ohana/src/relastro/include/relastro.h	(revision 12332)
+++ trunk/Ohana/src/relastro/include/relastro.h	(revision 12332)
@@ -0,0 +1,285 @@
+# include <ohana.h>
+# include <dvo.h>
+# include <kapa.h>
+# include <signal.h>
+
+typedef struct {
+  double R, D;  /* Sky Coords    - degrees */
+  double P, Q;  /* Tangent Plane - pixels  */
+  double L, M;  /* Focal Plane   - pixels  */
+  double X, Y;  /* Chip Coords   - pixels  */
+  double Mag, dMag;
+  int mask;
+} StarData;
+
+// structure to hold coordinate fitting terms
+typedef struct {
+    int Npts;
+    int Nterms;
+    int Norder;
+    int Nsums;
+    int Nelems;
+    double **sum;
+    double **xsum;
+    double **ysum;
+    double **xfit;
+    double **yfit;
+} CoordFit;
+
+typedef struct {
+  double Ro, dRo;
+  double Do, dDo;
+  
+  double uR, duR;
+  double uD, duD;
+
+  double p, dp;
+
+  double chisq;
+  int Nfit;
+} PMFit;
+
+typedef struct {
+  unsigned int start;
+  unsigned int stop;
+  float Mcal;
+  float dMcal;
+  short Xm;
+  float secz;
+  char code;
+  Coords coords;
+} Mosaic; 
+
+typedef struct {
+  double median;
+  double mean;
+  double sigma;
+  double error;
+  double chisq;
+  double min;
+  double max;
+  double total;
+  int    Nmeas;
+} StatType;
+
+/* global variables set in parameter file */
+char   ImageCat[256];
+char   ImageTemplate[256];
+char   CatTemplate[256];
+char   GSCFILE[256];
+char   CATDIR[256];
+char   CATMODE[16];    /* raw, mef, split, mysql */
+char   CATFORMAT[16];  /* internal, elixir, loneos, panstarrs */
+char   CameraConfig[256];
+char   SKY_TABLE[256];
+int    SKY_DEPTH;  /** XXX EAM : depth of catalog tables, fix usage */
+
+double MAG_LIM;
+double SIGMA_LIM;
+double IMAGE_SCATTER;
+double IMAGE_OFFSET;
+double STAR_SCATTER;
+double STAR_CHISQ;
+double MIN_ERROR;
+
+int    VERBOSE;
+
+int    NLOOP;
+int    RESET;
+int    UPDATE;
+int    PLOTSTUFF;
+int    SAVEPLOT;
+int    SHOW_PARAMS;
+char   MOSAICNAME[256];
+char   STATMODE[32];
+int    STAR_BAD;
+int    MEAS_BAD;
+int    STAR_TOOFEW;
+int    IMAGE_TOOFEW;
+double IMAGE_GOOD_FRACTION;
+int    IMAGE_BAD;
+int    FREEZE_IMAGES;
+int    USE_GRID;
+int    PLOTDELAY;
+
+int    RELPHOT_GRID_X;
+int    RELPHOT_GRID_Y;
+int    RELPHOT_GRID_BINNING;
+
+char          *PHOTCODE_LIST;
+int           Nphotcodes;
+PhotCode     **photcodes;
+int            PhotPrimary;
+int            PhotNsec;
+
+int AreaSelect;
+double AreaXmin, AreaXmax, AreaYmin, AreaYmax;
+
+int ImagSelect, ImagMin, ImagMax;
+
+double  PlotMmin, PlotMmax, PlotdMmin, PlotdMmax;
+enum {black, white, red, orange, yellow, green, blue, indigo, violet};
+
+int DophotSelect, DophotValue;
+
+int TimeSelect;
+time_t TSTART, TSTOP;
+
+int FIT_MODE;
+enum {FIT_NONE, FIT_AVERAGE, FIT_PM_ONLY, FIT_PM_AND_PAR};
+
+int FIT_TARGET;
+enum {TARGET_NONE, TARGET_OBJECTS, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS};
+
+SkyRegion UserPatch;
+
+int DoUpdateObjects;
+int DoUpdateSimple;
+int DoUpdateChips;
+int DoUpdateMosaics;
+
+/*** relphot prototypes ***/
+void          ConfigInit          PROTO((int *argc, char **argv));
+void          GetConfig           PROTO((char *config, char *field, char *format, int N, void *ptr));
+char         *GetPhotnamebyCode   PROTO((PhotCodeData *photcodes, int code));
+void          InterpolateGrid     PROTO((float *buffer, int Nx, int Ny, Coords *ccd, Coords *gcoords));
+int          *SelectRefMosaic     PROTO((Mosaic **refmosaic, int *Nimage));
+int           args                PROTO((int argc, char **argv));
+int           bcatalog            PROTO((Catalog *subcatalog, Catalog *catalog));
+void          clean_images        PROTO(());
+void          clean_measures      PROTO((Catalog *catalog, int Ncatalog, int final));
+void          clean_mosaics       PROTO(());
+void          clean_stars         PROTO((Catalog *catalog, int Ncatalog));
+int           corner_check        PROTO((double *x1, double *y1, double *x2, double *y2));
+void          dumpGrid            PROTO(());
+void          dump_grid           PROTO(()); 
+int           edge_check          PROTO((double *x1, double *y1, double *x2, double *y2));
+void          findImages          PROTO((Catalog *catalog, int Ncatalog));
+int           findMosaics         PROTO((Catalog *catalog, int Ncatalog));
+Image        *find_images         PROTO((FITS_DB *db, GSCRegion *region, int Nregion, int *Nimage, int **LineNum));
+void set_db (FITS_DB *in);
+int Shutdown (char *format, ...);
+void TrapSignal (int sig);
+void SetProtect (int mode);
+int SetSignals ();
+
+GSCRegion    *find_regions        PROTO((Image *image, int Nimage, int *Nregions, GSCRegion *fullregion));
+void          freeGridBins        PROTO((int Ncatalog));
+void          freeImageBins       PROTO((int Ncatalog));
+void          freeMosaicBins      PROTO((int Ncatalog));
+void          free_catalogs       PROTO((Catalog *catalog, int Ncatalog));
+int           gcatalog            PROTO((Catalog *catalog, int FINAL));
+Coords       *getCoords           PROTO((int meas, int cat));
+float         getMcal             PROTO((int meas, int cat));
+float         getMgrid            PROTO((int meas, int cat));
+float         getMmos             PROTO((int meas, int cat));
+float         getMrel             PROTO((Catalog *catalog, int meas, int cat));
+GSCRegion    *get_regions         PROTO((double minRa, double maxRa, double minDec, double maxDec, int *Nregions));
+void          getfullregion       PROTO((Image *image, int Nimage, GSCRegion *fullregion));
+Image        *getimage            PROTO((int N));
+Image        *getimages           PROTO((int *N));
+void          global_stats        PROTO((Catalog *catalog, int Ncatalog));
+void          initGrid            PROTO((int dX, int dY));
+void          initGridBins        PROTO((Catalog *catalog, int Ncatalog));
+void          initImageBins       PROTO((Catalog *catalog, int Ncatalog));
+void          initImages          PROTO((Image *input, int N));
+void          initMosaicBins      PROTO((Catalog *catalog, int Ncatalog));
+void          initMosaicGrid      PROTO((Image *image, int Nimage));
+void          initMosaics         PROTO((Image *image, int Nimage));
+void          initMrel            PROTO((Catalog *catalog, int Ncatalog));
+void          initialize          PROTO((int argc, char **argv));
+void          initstats           PROTO((char *mode));
+int           liststats           PROTO((double *value, double *dvalue, int N, StatType *stats));
+Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int subselect));
+SkyList      *load_images         PROTO((FITS_DB *db, SkyRegion *region));
+Image        *select_images       PROTO((SkyList *skylist, Image *timage, int Ntimage, int **LineNumber, int *Nimage));
+
+void check_permissions (char *basefile);
+void lock_image_db (FITS_DB *db, char *filename);
+void unlock_image_db (FITS_DB *db);
+void create_image_db (FITS_DB *db);
+
+int           main                PROTO((int argc, char **argv));
+void          mark_images         PROTO((Image *image, int Nimage, Image *timage, int Ntimage));
+void          matchImage          PROTO((Catalog *catalog, int meas, int cat));
+void          matchMosaics        PROTO((Catalog *catalog, int meas, int cat));
+GSCRegion    *name_region         PROTO((char *name, int *Nregions));
+double        opening_angle       PROTO((double x1, double y1, double x2, double y2, double x3, double y3));
+void          plot_chisq          PROTO((Catalog *catalog, int Ncatalog));
+void          plot_defaults       PROTO((Graphdata *graphdata));
+void          plot_grid           PROTO((Catalog *catalog));
+void          plot_images         PROTO(());
+void          plot_list           PROTO((Graphdata *graphdata, double *xlist, double *ylist, int N, char *label, char *file));
+void          plot_mosaic_fields  PROTO((Catalog *catalog));
+void          plot_mosaics        PROTO(());
+void          plot_scatter        PROTO((Catalog *catalog, int Ncatalog));
+void          plot_star_coords    PROTO((Catalog *catalog, int Ncatalog));
+void          plot_stars          PROTO((Catalog *catalog, int Ncatalog));
+void          reload_catalogs     PROTO((SkyList *skylist));
+int           setExclusions       PROTO((Catalog *catalog, int Ncatalog));
+void          setMcal             PROTO((Catalog *catalog, int Poor));
+void          setMcalFinal        PROTO(());
+int           setMcalOutput       PROTO((Catalog *catalog, int Ncatalog));
+void          setMgrid            PROTO((Catalog *catalog));
+int           setMmos             PROTO((Catalog *catalog, int Poor));
+int           setMrel             PROTO((Catalog *catalog, int Ncatalog));
+void          setMrelFinal        PROTO((Catalog *catalog));
+int           setMrelOutput       PROTO((Catalog *catalog, int Ncatalog, int mark));
+void          set_ZP              PROTO((double ZERO));
+int           setrefcode          PROTO((Image *image, int Nimage)); 
+void          skip_measurements   PROTO((Catalog *catalog, int pass));
+void          sortA               PROTO((double *X, int N));
+void          sortB               PROTO((double *X, double *Y, int N));
+void          sortC               PROTO((double *X, double *Y, double *F1, double *F2, int N));
+void          sortD               PROTO((double *X, double *Y, double *Z, int N));
+StatType      statsImageM         PROTO((Catalog *catalog));
+StatType      statsImageN         PROTO((Catalog *catalog));
+StatType      statsImageX         PROTO((Catalog *catalog));
+StatType      statsImagedM        PROTO((Catalog *catalog));
+StatType      statsMosaicM        PROTO((Catalog *catalog));
+StatType      statsMosaicN        PROTO((Catalog *catalog));
+StatType      statsMosaicX        PROTO((Catalog *catalog));
+StatType      statsMosaicdM       PROTO((Catalog *catalog));
+StatType      statsStarN          PROTO((Catalog *catalog, int Ncatalog));
+StatType      statsStarS          PROTO((Catalog *catalog, int Ncatalog));
+StatType      statsStarX          PROTO((Catalog *catalog, int Ncatalog));
+void          wcatalog            PROTO((Catalog *catalog));
+void          wimages             PROTO(());
+void          write_coords        PROTO((Header *header, Coords *coords));
+
+
+double **array_init (int Nx, int Ny);
+void array_free (double **array, int Nx);
+CoordFit *fit_init (int order);
+void fit_free (CoordFit *fit);
+void fit_add (CoordFit *fit, double x1, double y1, double x2, double y2, double wt);
+void fit_eval (CoordFit *fit);
+double **poly2d_dx (double **poly, int Nx, int Ny);
+double **poly2d_dy (double **poly, int Nx, int Ny);
+double **poly2d_copy (double **poly, int Nx, int Ny);
+double poly2d_eval (double **poly, int Nx, int Ny, double x, double y);
+void fit_apply_coords (CoordFit *fit, Coords *coords);
+int CoordsGetCenter (CoordFit *fit, double tol, double *xo, double *yo);
+CoordFit *CoordsSetCenter (CoordFit *input, double Xo, double Yo);
+void FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
+void FitMosaic (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
+void FitSimple (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
+void initObjectData (Catalog *catalog, int Ncatalog);
+int UpdateObjects (Catalog *catalog, int Ncatalog);
+int UpdateSimple (Catalog *catalog, int Ncatalog);
+int UpdateChips (Catalog *catalog, int Ncatalog);
+int UpdateMosaic (Catalog *catalog, int Ncatalog);
+int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon);
+int ParFactor (double *pR, double *pD, double R, double D, time_t T);
+int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, time_t *T, int Npts);
+int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, time_t *T, double *pR, double *pD, int Npts);
+
+Image *getMosaicForImage (int N);
+
+StarData *getImageRef (Catalog *catalog, int Ncatalog, int im, int *Nstars, int isMosaic);
+StarData *getImageRaw (Catalog *catalog, int Ncatalog, int im, int *Nstars, int isMosaic);
+
+double getMeanR (Measure *measure, Average *average, SecFilt *secfilt);
+double getMeanD (Measure *measure, Average *average, SecFilt *secfilt);
+int setMeanR (double ra_fit, Measure *measure, Average *average, SecFilt *secfilt);
+int setMeanD (double dec_fit, Measure *measure, Average *average, SecFilt *secfilt);
Index: trunk/Ohana/src/relastro/include/relphot.h
===================================================================
--- trunk/Ohana/src/relastro/include/relphot.h	(revision 12067)
+++ 	(revision )
@@ -1,212 +1,0 @@
-# include <ohana.h>
-# include <dvo.h>
-# include <kapa.h>
-# include <signal.h>
-
-/* # define GRID_V1 */
-# define GRID_V2
-
-# if (0)
-typedef struct {
-  double xmin, xmax, ymin, ymax;
-  int style, ptype, ltype, etype, color;
-  double lweight, size;
-} Graphdata;
-# endif
-
-typedef struct {
-  unsigned int start;
-  unsigned int stop;
-  float Mcal;
-  float dMcal;
-  short Xm;
-  float secz;
-  char code;
-  Coords coords;
-} Mosaic; 
-
-typedef struct {
-  double median;
-  double mean;
-  double sigma;
-  double error;
-  double chisq;
-  double min;
-  double max;
-  double total;
-  int    Nmeas;
-} StatType;
-
-/* global variables set in parameter file */
-char   ImageCat[256];
-char   ImageTemplate[256];
-char   CatTemplate[256];
-char   GSCFILE[256];
-char   CATDIR[256];
-char   CATMODE[16];    /* raw, mef, split, mysql */
-char   CATFORMAT[16];  /* internal, elixir, loneos, panstarrs */
-char   CameraConfig[256];
-
-double MAG_LIM;
-double SIGMA_LIM;
-double IMAGE_SCATTER;
-double IMAGE_OFFSET;
-double STAR_SCATTER;
-double STAR_CHISQ;
-double MIN_ERROR;
-
-int    VERBOSE;
-
-int    NLOOP;
-int    RESET;
-int    UPDATE;
-int    PLOTSTUFF;
-int    SAVEPLOT;
-int    SHOW_PARAMS;
-char   MOSAICNAME[256];
-char   STATMODE[32];
-int    STAR_BAD;
-int    MEAS_BAD;
-int    STAR_TOOFEW;
-int    IMAGE_TOOFEW;
-double IMAGE_GOOD_FRACTION;
-int    IMAGE_BAD;
-int    FREEZE_IMAGES;
-int    USE_GRID;
-int    PLOTDELAY;
-
-int    RELPHOT_GRID_X;
-int    RELPHOT_GRID_Y;
-int    RELPHOT_GRID_BINNING;
-
-PhotCode      *photcode;
-int            PhotPrimary;
-int            PhotNsec;
-int            PhotSec;
-
-int AreaSelect;
-double AreaXmin, AreaXmax, AreaYmin, AreaYmax;
-
-int ImagSelect, ImagMin, ImagMax;
-
-double  PlotMmin, PlotMmax, PlotdMmin, PlotdMmax;
-enum {black, white, red, orange, yellow, green, blue, indigo, violet};
-
-int TimeSelect;
-time_t TSTART, TSTOP;
-
-# ifdef GRID_V1
-int setGridMeasure (int meas, int cat, double X, double Y);
-# endif
-
-# ifdef GRID_V2
-int setGridMeasure (int meas, int cat, double X, double Y, int ccdnum);
-# endif
-
-/*** relphot prototypes ***/
-void          ConfigInit          PROTO((int *argc, char **argv));
-void          GetConfig           PROTO((char *config, char *field, char *format, int N, void *ptr));
-char         *GetPhotnamebyCode   PROTO((PhotCodeData *photcodes, int code));
-void          InterpolateGrid     PROTO((float *buffer, int Nx, int Ny, Coords *ccd, Coords *gcoords));
-int          *SelectRefMosaic     PROTO((Mosaic **refmosaic, int *Nimage));
-int           args                PROTO((int argc, char **argv));
-int           bcatalog            PROTO((Catalog *subcatalog, Catalog *catalog, GSCRegion *fullregion));
-void          clean_images        PROTO(());
-void          clean_measures      PROTO((Catalog *catalog, int Ncatalog, int final));
-void          clean_mosaics       PROTO(());
-void          clean_stars         PROTO((Catalog *catalog, int Ncatalog));
-int           corner_check        PROTO((double *x1, double *y1, double *x2, double *y2));
-void          dumpGrid            PROTO(());
-void          dump_grid           PROTO(()); 
-int           edge_check          PROTO((double *x1, double *y1, double *x2, double *y2));
-void          findImages          PROTO((Catalog *catalog, int Ncatalog));
-int           findMosaics         PROTO((Catalog *catalog, int Ncatalog));
-Image        *find_images         PROTO((FITS_DB *db, GSCRegion *region, int Nregion, int *Nimage, int **LineNum));
-void set_db (FITS_DB *in);
-int Shutdown (char *format, ...);
-void TrapSignal (int sig);
-void SetProtect (int mode);
-int SetSignals ();
-
-GSCRegion    *find_regions        PROTO((Image *image, int Nimage, int *Nregions, GSCRegion *fullregion));
-void          freeGridBins        PROTO((int Ncatalog));
-void          freeImageBins       PROTO((int Ncatalog));
-void          freeMosaicBins      PROTO((int Ncatalog));
-void          free_catalogs       PROTO((Catalog *catalog, int Ncatalog));
-int           gcatalog            PROTO((Catalog *catalog, int FINAL));
-Coords       *getCoords           PROTO((int meas, int cat));
-float         getMcal             PROTO((int meas, int cat));
-float         getMgrid            PROTO((int meas, int cat));
-float         getMmos             PROTO((int meas, int cat));
-float         getMrel             PROTO((Catalog *catalog, int meas, int cat));
-GSCRegion    *get_regions         PROTO((double minRa, double maxRa, double minDec, double maxDec, int *Nregions));
-void          getfullregion       PROTO((Image *image, int Nimage, GSCRegion *fullregion));
-Image        *getimage            PROTO((int N));
-Image        *getimages           PROTO((int *N));
-void          global_stats        PROTO((Catalog *catalog, int Ncatalog));
-void          initGrid            PROTO((int dX, int dY));
-void          initGridBins        PROTO((Catalog *catalog, int Ncatalog));
-void          initImageBins       PROTO((Catalog *catalog, int Ncatalog));
-void          initImages          PROTO((Image *input, int N));
-void          initMosaicBins      PROTO((Catalog *catalog, int Ncatalog));
-void          initMosaicGrid      PROTO((Image *image, int Nimage));
-void          initMosaics         PROTO((Image *image, int Nimage));
-void          initMrel            PROTO((Catalog *catalog, int Ncatalog));
-void          initialize          PROTO((int argc, char **argv));
-void          initstats           PROTO((char *mode));
-int           liststats           PROTO((double *value, double *dvalue, int N, StatType *stats));
-Catalog      *load_catalogs       PROTO((GSCRegion *region, int Nregion, GSCRegion *fullregion));
-GSCRegion    *load_images         PROTO((FITS_DB *db, char *seed, int *nregion, GSCRegion *fullregion));
-
-void check_permissions (char *basefile);
-void lock_image_db (FITS_DB *db, char *filename);
-void unlock_image_db (FITS_DB *db);
-void create_image_db (FITS_DB *db);
-
-int           main                PROTO((int argc, char **argv));
-void          mark_images         PROTO((Image *image, int Nimage, Image *timage, int Ntimage));
-void          matchImage          PROTO((Catalog *catalog, int meas, int cat));
-void          matchMosaics        PROTO((Catalog *catalog, int meas, int cat));
-GSCRegion    *name_region         PROTO((char *name, int *Nregions));
-double        opening_angle       PROTO((double x1, double y1, double x2, double y2, double x3, double y3));
-void          plot_chisq          PROTO((Catalog *catalog, int Ncatalog));
-void          plot_defaults       PROTO((Graphdata *graphdata));
-void          plot_grid           PROTO((Catalog *catalog));
-void          plot_images         PROTO(());
-void          plot_list           PROTO((Graphdata *graphdata, double *xlist, double *ylist, int N, char *label, char *file));
-void          plot_mosaic_fields  PROTO((Catalog *catalog));
-void          plot_mosaics        PROTO(());
-void          plot_scatter        PROTO((Catalog *catalog, int Ncatalog));
-void          plot_star_coords    PROTO((Catalog *catalog, int Ncatalog));
-void          plot_stars          PROTO((Catalog *catalog, int Ncatalog));
-void          reload_catalogs     PROTO((GSCRegion *region, int Nregion, GSCRegion *fullregion));
-int           setExclusions       PROTO((Catalog *catalog, int Ncatalog));
-void          setMcal             PROTO((Catalog *catalog, int Poor));
-void          setMcalFinal        PROTO(());
-int           setMcalOutput       PROTO((Catalog *catalog, int Ncatalog));
-void          setMgrid            PROTO((Catalog *catalog));
-int           setMmos             PROTO((Catalog *catalog, int Poor));
-int           setMrel             PROTO((Catalog *catalog, int Ncatalog));
-void          setMrelFinal        PROTO((Catalog *catalog));
-int           setMrelOutput       PROTO((Catalog *catalog, int Ncatalog, int mark));
-void          set_ZP              PROTO((double ZERO));
-int           setrefcode          PROTO((Image *image, int Nimage)); 
-void          skip_measurements   PROTO((Catalog *catalog, int pass));
-void          sortA               PROTO((double *X, int N));
-void          sortB               PROTO((double *X, double *Y, int N));
-void          sortC               PROTO((double *X, double *Y, double *F1, double *F2, int N));
-void          sortD               PROTO((double *X, double *Y, double *Z, int N));
-StatType      statsImageM         PROTO((Catalog *catalog));
-StatType      statsImageN         PROTO((Catalog *catalog));
-StatType      statsImageX         PROTO((Catalog *catalog));
-StatType      statsImagedM        PROTO((Catalog *catalog));
-StatType      statsMosaicM        PROTO((Catalog *catalog));
-StatType      statsMosaicN        PROTO((Catalog *catalog));
-StatType      statsMosaicX        PROTO((Catalog *catalog));
-StatType      statsMosaicdM       PROTO((Catalog *catalog));
-StatType      statsStarN          PROTO((Catalog *catalog, int Ncatalog));
-StatType      statsStarS          PROTO((Catalog *catalog, int Ncatalog));
-StatType      statsStarX          PROTO((Catalog *catalog, int Ncatalog));
-void          wcatalog            PROTO((Catalog *catalog));
-void          wimages             PROTO(());
-void          write_coords        PROTO((Header *header, Coords *coords));
Index: trunk/Ohana/src/relastro/src/ConfigInit.c
===================================================================
--- trunk/Ohana/src/relastro/src/ConfigInit.c	(revision 12067)
+++ trunk/Ohana/src/relastro/src/ConfigInit.c	(revision 12332)
@@ -1,3 +1,3 @@
-# include "relphot.h"
+# include "relastro.h"
 
 void ConfigInit (int *argc, char **argv) {
@@ -5,5 +5,6 @@
   double ZERO_POINT;
   char  *config, *file;
-  char   PhotCodeFile[256];
+  char CatdirPhotcodeFile[256];
+  char MasterPhotcodeFile[256];
 
   /*** load configuration info ***/
@@ -32,4 +33,6 @@
   ScanConfig (config, "CATMODE",                "%s",  0, CATMODE);
   ScanConfig (config, "CATFORMAT",              "%s",  0, CATFORMAT);
+  ScanConfig (config, "PHOTCODE_FILE",         	"%s",  0, MasterPhotcodeFile);
+
   sprintf (ImageCat, "%s/Images.dat", CATDIR);
 
@@ -41,5 +44,4 @@
   }
 
-  GetConfig (config, "PHOTCODE_FILE",          "%s",  0, PhotCodeFile);
   GetConfig (config, "ZERO_PT",                "%lf", 0, &ZERO_POINT);
 
@@ -52,6 +54,8 @@
   if (*CATFORMAT == 0) strcpy (CATFORMAT, "ELIXIR");
 
-  if (!LoadPhotcodes (PhotCodeFile)) {
-    fprintf (stderr, "error loading photcodes\n");
+  /* XXX this does not yet write out the master photcode table */
+  sprintf (CatdirPhotcodeFile, "%s/Photcodes.dat", CATDIR);
+  if (!LoadPhotcodes (CatdirPhotcodeFile, MasterPhotcodeFile)) {
+    fprintf (stderr, "error loading photcode table %s or master file %s\n", CatdirPhotcodeFile, MasterPhotcodeFile);
     exit (1);
   }
Index: trunk/Ohana/src/relastro/src/FitChip.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitChip.c	(revision 12332)
+++ trunk/Ohana/src/relastro/src/FitChip.c	(revision 12332)
@@ -0,0 +1,36 @@
+# include "relastro.h"
+
+void FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords) {
+
+  int i;
+  CoordFit *fit;
+
+  fit = fit_init (coords[0].Npolyterms);
+  for (i = 0; i < Nmatch; i++) {
+    if (raw[i].mask) continue;
+    fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, 1.0);
+  }
+  fit_eval (fit);
+  fit_apply_coords (fit, coords);
+  fit_free (fit);
+  /* FitChip and FitSimple update the coords in different ways? maybe not... */
+}
+
+/* XXX the supporting fitting functions are not quite right
+   I probably need to apply the projection and then do the fit 
+   to L,M.  I also need to incorporate the transformations 
+   determined for psastro */
+
+
+/* in the simple case, we only have three coord systems of interest:
+   R,D : the sky
+   P,Q : the tangent plane
+   X,Y : the chip
+
+   R,D -> P,Q (projection)
+   P,Q -> X,Y (polynomial transformation)
+*/
+
+/* XXX I'm not using the errors at all : this could at least be done with the dMag values */
+
+   
Index: trunk/Ohana/src/relastro/src/FitMosaic.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitMosaic.c	(revision 12332)
+++ trunk/Ohana/src/relastro/src/FitMosaic.c	(revision 12332)
@@ -0,0 +1,36 @@
+# include "relastro.h"
+
+void FitMosaic (StarData *raw, StarData *ref, int Nmatch, Coords *coords) {
+
+  int i;
+  CoordFit *fit;
+
+  fit = fit_init (coords[0].Npolyterms);
+  for (i = 0; i < Nmatch; i++) {
+    if (raw[i].mask) continue;
+    fit_add (fit, raw[i].L, raw[i].M, ref[i].P, ref[i].Q, 1.0);
+  }
+  fit_eval (fit);
+  fit_apply_coords (fit, coords);
+  fit_free (fit);
+  /* FitChip and FitSimple update the coords in different ways? maybe not... */
+}
+
+/* XXX the supporting fitting functions are not quite right
+   I probably need to apply the projection and then do the fit 
+   to L,M.  I also need to incorporate the transformations 
+   determined for psastro */
+
+
+/* in the simple case, we only have three coord systems of interest:
+   R,D : the sky
+   P,Q : the tangent plane
+   X,Y : the chip
+
+   R,D -> P,Q (projection)
+   P,Q -> X,Y (polynomial transformation)
+*/
+
+/* XXX I'm not using the errors at all : this could at least be done with the dMag values */
+
+   
Index: trunk/Ohana/src/relastro/src/FitPM.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitPM.c	(revision 12332)
+++ trunk/Ohana/src/relastro/src/FitPM.c	(revision 12332)
@@ -0,0 +1,72 @@
+# include "relastro.h"
+
+/* do we want an init function which does the alloc and a clear function to free? */
+int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, time_t *T, int Npts) {
+
+  int i;
+
+  double **A, **B;
+  double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT;
+
+  /* do I need to do this as 2 2x2 matrix equations? */
+  A = array_init (4, 4);
+  B = array_init (4, 1);
+
+  Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0;
+  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;
+    
+    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[1][0] = Tx;
+  A[1][1] = Tx2;
+
+  A[2][2] = Wy;
+  A[2][3] = Ty;
+
+  A[3][2] = Ty;
+  A[3][3] = Ty2;
+
+  B[0][0] = Xs;
+  B[1][0] = XT;
+  B[2][0] = Ys;
+  B[3][0] = YT;
+
+  dgaussj (A, 4, B, 1);
+
+  fit[0].Ro = B[0][0];
+  fit[0].uR = B[1][0];
+  fit[0].Do = B[2][0];
+  fit[0].uD = B[3][0];
+  fit[0].p  = 0.0;
+  
+  fit[0].dRo = sqrt(A[0][0]);
+  fit[0].duR = sqrt(A[1][1]);
+  fit[0].dDo = sqrt(A[2][2]);
+  fit[0].duD = sqrt(A[3][3]);
+  fit[0].dp  = 0.0;
+  
+  array_free (A, 4);
+  array_free (B, 4);
+
+  return (TRUE);
+}
Index: trunk/Ohana/src/relastro/src/FitPMandPar.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitPMandPar.c	(revision 12332)
+++ trunk/Ohana/src/relastro/src/FitPMandPar.c	(revision 12332)
@@ -0,0 +1,98 @@
+# include "relastro.h"
+
+/* do we want an init function which does the alloc and a clear function to free? */
+int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, time_t *T, double *pR, double *pD, int Npts) {
+
+  int i;
+
+  double **A, **B;
+  double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT;
+  double PR, PD, PRT, PDT, PRX, PDY, PR2, PD2;
+
+  A = array_init (5, 5);
+  B = array_init (5, 1);
+
+  PR = PD = PRT = PDT = PRX = PDY = PR2 = PD2 = 0.0;
+  Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0;
+  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);
+
+  fit[0].Ro = B[0][0];
+  fit[0].uR = B[1][0];
+  fit[0].Do = B[2][0];
+  fit[0].uD = B[3][0];
+  fit[0].p  = B[4][0];
+  
+  fit[0].dRo = sqrt(A[0][0]);
+  fit[0].duR = sqrt(A[1][1]);
+  fit[0].dDo = sqrt(A[2][2]);
+  fit[0].duD = sqrt(A[3][3]);
+  fit[0].dp  = sqrt(A[4][4]);
+  
+  array_free (A, 5);
+  array_free (B, 5);
+
+  /* get the chisq from the matrix values */
+
+  return (TRUE);
+}
Index: trunk/Ohana/src/relastro/src/FitSimple.c
===================================================================
--- trunk/Ohana/src/relastro/src/FitSimple.c	(revision 12332)
+++ trunk/Ohana/src/relastro/src/FitSimple.c	(revision 12332)
@@ -0,0 +1,35 @@
+# include "relastro.h"
+
+void FitSimple (StarData *raw, StarData *ref, int Nmatch, Coords *coords) {
+
+  int i;
+  CoordFit *fit;
+
+  fit = fit_init (coords[0].Npolyterms);
+  for (i = 0; i < Nmatch; i++) {
+    if (raw[i].mask) continue;
+    fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, 1.0);
+  }
+  fit_eval (fit);
+  fit_apply_coords (fit, coords);
+  fit_free (fit);
+}
+
+/* XXX the supporting fitting functions are not quite right
+   I probably need to apply the projection and then do the fit 
+   to L,M.  I also need to incorporate the transformations 
+   determined for psastro */
+
+
+/* in the simple case, we only have three coord systems of interest:
+   R,D : the sky
+   P,Q : the tangent plane
+   X,Y : the chip
+
+   R,D -> P,Q (projection)
+   P,Q -> X,Y (polynomial transformation)
+*/
+
+/* XXX I'm not using the errors at all : this could at least be done with the dMag values */
+
+   
Index: trunk/Ohana/src/relastro/src/ImageOps.c
===================================================================
--- trunk/Ohana/src/relastro/src/ImageOps.c	(revision 12067)
+++ trunk/Ohana/src/relastro/src/ImageOps.c	(revision 12332)
@@ -1,3 +1,3 @@
-# include "relphot.h"
+# include "relastro.h"
 
 static unsigned int *start;
@@ -13,4 +13,13 @@
 static int          Nimage;
 
+Image *getimages (int *N) {
+  *N = Nimage;
+  return (image);
+}
+
+Image *getimage (int N) {
+  return (&image[N]);
+}
+
 void initImages (Image *input, int N) {
 
@@ -68,13 +77,11 @@
 }
 
-/* select all image equivalent to the current photcode */
+/* match measurements to images */
 void findImages (Catalog *catalog, int Ncatalog) {
 
-  int i, j, ecode;
+  int i, j;
 
   for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Nmeasure; j++) {
-      ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[j].source);
-      if (photcode[0].code != ecode) continue;
       matchImage (catalog, j, i);
     }
@@ -82,42 +89,17 @@
 }
 
+/* modify this function to use the measure->imageID field */
 void matchImage (Catalog *catalog, int meas, int cat) {
 
-  int i, ave, ccdnum;
-  double ra, dec, X, Y;
-  char *pname, *filter, *p, base[256];
+  int i;
   Measure *measure;
   
   measure = &catalog[cat].measure[meas];
   for (i = 0; i < Nimage; i++) {
-    if (image[0].source == -1) continue;
-    if (measure[0].source != image[i].source) continue;
+    if (image[0].photcode == -1) continue;
+    if (measure[0].photcode != image[i].photcode) continue;
     if (measure[0].t < start[i]) continue;
     if (measure[0].t > stop[i]) continue;
     
-# ifdef GRID_V2 /* this section is added to support GridOps.v2.c */
-    if (USE_GRID) {
-
-      /* identify the ccd on the basis of the photcode name */
-      pname = GetPhotcodeNamebyCode (image[i].source);
-      filter = photcode[0].name;
-      sprintf (base, "%s.%s.", MOSAICNAME, filter);
-      if (strncmp (pname, base, strlen (base))) continue;
-      p = pname + strlen(base);
-      if (*p == 0) continue;
-      ccdnum = atoi (p);
-      /* some ambiguity here between seq number and id number */
-
-      /* add this measurement to the grid cell for this chip */
-      ave = measure[0].averef;
-      ra  = catalog[cat].average[ave].R - measure[0].dR_PS / 3600.0;
-      dec = catalog[cat].average[ave].D - measure[0].dD_PS / 3600.0;
-       
-      /* X,Y always positive-definite in range 0,0 - dX, dY */
-      RD_to_XY (&X, &Y, ra, dec, &image[i].coords);
-      setGridMeasure (meas, cat, X, Y, ccdnum);
-    }
-# endif
-
     bin[cat][meas] = i;
 
@@ -133,18 +115,4 @@
     return;
   }
-  /*  fprintf (stderr, "can't find source image for this measurement: %d (%d)\n", measure[0].t, measure[0].source); */
-}
-
-float getMcal (int meas, int cat) {
-
-  int i;
-  float value;
-
-  i = bin[cat][meas];
-  if (i == -1) return (NO_MAG);
-
-  if (image[i].code & IMAGE_BAD)  return (NO_MAG);  
-  value = image[i].Mcal_PS;
-  return (value);
 }
 
@@ -158,128 +126,4 @@
 }
 
-/* determine Mcal values for all images */
-void setMcal (Catalog *catalog, int PoorImages) {
-
-  int i, j, m, c, n, N, Nmax, mark, bad;
-  float Msys, Mrel, Mmos, Mgrid;
-  double *list, *dlist;
-  StatType stats;
-
-  if (FREEZE_IMAGES) return;
-
-  if (PoorImages) {
-    IMAGE_BAD = STAR_BAD = MEAS_BAD = 0;
-  }
-
-  Nmax = 0;
-  for (i = 0; i < Nimage; i++) {
-    Nmax = MAX (Nmax, Nlist[i]);
-  }
-  ALLOCATE (list, double, Nmax);
-  ALLOCATE (dlist, double, Nmax);
-
-  for (i = 0; i < Nimage; i++) {
-    
-    /* on PoorImages run, skip good images */
-    if (PoorImages) {
-      bad = image[i].code & (ID_IMAGE_FEW | ID_IMAGE_POOR | ID_IMAGE_SKIP);
-      if (!bad) continue;
-    }      
-
-    N = 0;
-    for (j = 0; j < Nlist[i]; j++) {
-      
-      m = mlist[i][j];
-      c = clist[i][j];
-      
-      if (catalog[c].measure[m].flags & MEAS_BAD) continue;
-      if ((Mmos  = getMmos  (m, c)) == NO_MAG) continue;
-      if ((Mgrid = getMgrid (m, c)) == NO_MAG) continue;
-      if ((Mrel  = getMrel  (catalog, m, c)) == NO_MAG) continue;
-      
-      n = catalog[c].measure[m].averef;
-      Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
-      list[N] = Msys - Mrel - Mmos - Mgrid;
-      dlist[N] = MAX (catalog[c].measure[m].dM_PS, MIN_ERROR);
-      N++;
-    }
-    /* Nlist[i] is all measurements, N is good measurements */
-
-    /* too few good measurements or too many bad measurements */
-    if (!PoorImages) {
-      mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*Nlist[i]);
-      if (mark) {
-	image[i].code |= ID_IMAGE_FEW;
-      } else {
-	image[i].code &= ~ID_IMAGE_FEW;
-      }      
-    }
-    
-    liststats (list, dlist, N, &stats);
-    image[i].Mcal_PS  = stats.mean;
-    image[i].dMcal_PS = stats.sigma;
-    image[i].Xm    = 100.0*log10(stats.chisq);
-  }
-  free (list);
-  free (dlist);
-  if (PoorImages) {
-    IMAGE_BAD = ID_IMAGE_POOR | ID_IMAGE_FEW | ID_IMAGE_SKIP;
-    STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
-    MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR | ID_MEAS_SKIP | ID_MEAS_AREA;
-  }
-  return;
-}
-
-/* mark image if: abs(Mcal) too large, dMcal too large */
-void clean_images () {
-
-  int i, N, mark, Nmark;
-  double *mlist, *slist, *dlist;
-  double MaxOffset, MaxScatter, MedOffset;
-  StatType stats;
-
-  if (FREEZE_IMAGES) return;
-
-  if (VERBOSE) fprintf (stderr, "marking poor images\n");
-
-  ALLOCATE (mlist, double, Nimage);
-  ALLOCATE (slist, double, Nimage);
-  ALLOCATE (dlist, double, Nimage);
-
-  for (i = N = 0; i < Nimage; i++) {
-    if (image[i].code & IMAGE_BAD) continue;
-    mlist[N] = fabs (image[i].Mcal_PS);
-    slist[N] = image[i].dMcal_PS;
-    dlist[N] = 1;
-    N++;
-  }
-  initstats ("MEAN");
-  liststats (mlist, dlist, N, &stats);
-  MaxOffset = MAX (IMAGE_OFFSET, 3*stats.sigma);
-  MedOffset = stats.median;
-  liststats (slist, dlist, N, &stats);
-  MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median);
-  fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset);
-  
-  Nmark = 0;
-  for (i = 0; i < Nimage; i++) {
-    mark = FALSE;
-    image[i].code &= ~ID_IMAGE_POOR;
-    mark = (image[i].dMcal_PS > MaxScatter) || (fabs(image[i].Mcal_PS - MedOffset) > MaxOffset);
-    if (mark) { 
-      Nmark ++;
-      image[i].code |= ID_IMAGE_POOR;
-    } else {
-      image[i].code &= ~ID_IMAGE_POOR;
-    }
-  }
-
-  fprintf (stderr, "%d images marked poor\n", Nmark);
-  initstats (STATMODE);
-  free (mlist);
-  free (slist);
-  free (dlist);
-}
-
 void plot_images () {
 
@@ -296,7 +140,7 @@
   /**** dMcal vs airmass ****/
   for (i = 0; i < Nimage; i++) {
-    Mlist[i] = image[i].Mcal_PS;
-    dlist[i] = image[i].dMcal_PS;
-    xlist[i] = image[i].secz_PS;
+    Mlist[i] = image[i].Mcal;
+    dlist[i] = image[i].dMcal;
+    xlist[i] = image[i].secz;
   }
 
@@ -316,5 +160,5 @@
   bzero (Mlist, NBIN*sizeof(double));
   for (i = 0; i < Nimage; i++) {
-    bin = image[i].dMcal_PS / 0.00025;
+    bin = image[i].dMcal / 0.00025;
     bin = MAX (0, MIN (NBIN - 1, bin));
     Mlist[bin] += 1.0;
@@ -330,134 +174,88 @@
 }
 
-StatType statsImageN (Catalog *catalog) {
-
-  int i, j, m, c, n, N;
-  double *list, *dlist;
-  StatType stats;
-
-  bzero (&stats, sizeof (StatType));
-  if (FREEZE_IMAGES) return (stats);
-
-  ALLOCATE (list, double, Nimage);
-  ALLOCATE (dlist, double, Nimage);
-
-  n = 0;
-  for (i = 0; i < Nimage; i++) {
-    if (image[i].code & IMAGE_BAD)  continue;
-
-    N = 0;
-    for (j = 0; j < Nlist[i]; j++) {
-
-      m = mlist[i][j];
-      c = clist[i][j];
-
-      if (getMcal  (m, c) == NO_MAG) continue;
-      if (getMmos  (m, c) == NO_MAG) continue;
-      if (getMgrid (m, c) == NO_MAG) continue;
-      N++;
+/* XXX re-write to use image->nstars if that is possible */
+StarData *getImageRaw (Catalog *catalog, int Ncatalog, int im, int *Nstars, int isMosaic) {
+
+  int i, m, c;
+
+  Image *mosaic;
+  Coords *moscoords;
+  StarData *raw;
+  
+  ALLOCATE (raw, StarData, Nlist[im]);
+
+  mosaic = NULL;
+  moscoords = NULL;
+  if (isMosaic) {
+    mosaic = getMosaicForImage (im);
+    moscoords = &mosaic[0].coords;
+  }
+
+  for (i = 0; i < Nlist[im]; i++) {
+    m = mlist[im][i];
+    c = clist[im][i];
+
+    /* apply the current image transformation or use the current value of R+dR, D+dD? */
+    raw[i].X = catalog[c].measure[m].Xccd;
+    raw[i].Y = catalog[c].measure[m].Yccd;
+    
+    raw[i].Mag  = catalog[c].measure[m].M;
+    raw[i].dMag = catalog[c].measure[m].dM;
+
+    /* note that for a Simple image, L,M = P,Q */
+    XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords);
+    if (isMosaic) {
+      XY_to_LM (&raw[i].P, &raw[i].Q, raw[i].L, raw[i].M, moscoords);
+      LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, moscoords);
+    } else {
+      raw[i].P = raw[i].L;
+      raw[i].Q = raw[i].M;
+      LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, &image[im].coords);
     }
-    list[n] = N;
-    dlist[n] = 1;
-    n++;
-  }
-
-  liststats (list, dlist, n, &stats);
-  free (list);
-  free (dlist);
-  return (stats);
-}
-
-StatType statsImageX (Catalog *catalog) {
-
-  int i, n;
-  double *list, *dlist;
-  StatType stats;
-
-  bzero (&stats, sizeof (StatType));
-  if (FREEZE_IMAGES) return (stats);
-
-  ALLOCATE (list, double, Nimage);
-  ALLOCATE (dlist, double, Nimage);
-
-  n = 0;
-  for (i = 0; i < Nimage; i++) {
-
-    if (image[i].code & IMAGE_BAD)  continue;
-
-    list[n] = pow (10.0, 0.01*image[i].Xm);
-    dlist[n] = 1;
-    n++;
-  }
-
-  liststats (list, dlist, n, &stats);
-  free (list);
-  free (dlist);
-  return (stats);
-}
-
-StatType statsImageM (Catalog *catalog) {
-
-  int i, n;
-  double *list, *dlist;
-  StatType stats;
-
-  bzero (&stats, sizeof (StatType));
-  if (FREEZE_IMAGES) return (stats);
-
-  ALLOCATE (list, double, Nimage);
-  ALLOCATE (dlist, double, Nimage);
-
-  n = 0;
-  for (i = 0; i < Nimage; i++) {
-
-    if (image[i].code & IMAGE_BAD)  continue;
-
-    list[n] = image[i].Mcal_PS;
-    dlist[n] = 1;
-    n++;
-  }
-
-  liststats (list, dlist, n, &stats);
-  free (list);
-  free (dlist);
-  return (stats);
-}
-
-StatType statsImagedM (Catalog *catalog) {
-
-  int i, n;
-  double *list, *dlist;
-  StatType stats;
-
-  bzero (&stats, sizeof (StatType));
-  if (FREEZE_IMAGES) return (stats);
-
-  ALLOCATE (list, double, Nimage);
-  ALLOCATE (dlist, double, Nimage);
-
-  n = 0;
-  for (i = 0; i < Nimage; i++) {
-
-    if (image[i].code & IMAGE_BAD)  continue;
-
-    list[n] = image[i].dMcal_PS;
-    dlist[n] = 1;
-    n++;
-  }
-
-  liststats (list, dlist, n, &stats);
-  free (list);
-  free (dlist);
-  return (stats);
-}
-
-Image *getimages (int *N) {
-
-  *N = Nimage;
-  return (image);
-}
-
-Image *getimage (int N) {
-  return (&image[N]);
-}
-
+  }
+  
+  *Nstars = Nlist[im];
+  return (raw);
+}
+
+StarData *getImageRef (Catalog *catalog, int Ncatalog, int im, int *Nstars, int isMosaic) {
+
+  int i, m, c, n;
+
+  Image *mosaic;
+  Coords *moscoords;
+  StarData *ref;
+  
+  ALLOCATE (ref, StarData, Nlist[im]);
+
+  mosaic = NULL;
+  moscoords = NULL;
+  if (isMosaic) {
+    mosaic = getMosaicForImage (im);
+    moscoords = &mosaic[0].coords;
+  }
+
+  for (i = 0; i < Nlist[im]; i++) {
+    m = mlist[im][i];
+    c = clist[im][i];
+    n = catalog[c].measure[m].averef;
+
+    /* apply the current image transformation or use the current value of R+dR, D+dD? */
+    ref[i].R = catalog[c].average[n].R;
+    ref[i].D = catalog[c].average[n].D;
+
+    /* note that for a Simple image, L,M = P,Q */
+    RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, &image[im].coords);
+    if (isMosaic) {
+      LM_to_XY (&ref[i].M, &ref[i].L, ref[i].P, ref[i].Q, moscoords);
+      LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, moscoords);
+    } else {
+      ref[i].L = ref[i].P;
+      ref[i].M = ref[i].Q;
+      LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords);
+    }
+  }
+  
+  *Nstars = Nlist[im];
+  return (ref);
+}
Index: trunk/Ohana/src/relastro/src/MosaicOps.c
===================================================================
--- trunk/Ohana/src/relastro/src/MosaicOps.c	(revision 12067)
+++ trunk/Ohana/src/relastro/src/MosaicOps.c	(revision 12332)
@@ -1,3 +1,3 @@
-# include "relphot.h"
+# include "relastro.h"
 
 Image *getimages (int *N);
@@ -16,4 +16,5 @@
 
 /* find mosaic frames (unique time periods & photcode name matches mosaic) */
+/* XXX what is a mosaic?  do we need to use the 'DIS/WRP' info to track? */
 void initMosaics (Image *image, int Nimage) {
 
@@ -36,5 +37,5 @@
 
     /* select valid mosaic images by photcode */
-    pname = GetPhotcodeNamebyCode (image[i].source);
+    pname = GetPhotcodeNamebyCode (image[i].photcode);
     status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME));
     if (status) continue;
@@ -69,5 +70,5 @@
     mosaic[Nmosaic].Xm    = 0.0;
     mosaic[Nmosaic].code  = image[i].code;
-    mosaic[Nmosaic].secz  = image[i].secz_PS;
+    mosaic[Nmosaic].secz  = image[i].secz;
 
     /* add image to mosaic image list */
@@ -86,104 +87,5 @@
     }
   }
-
-  initMosaicGrid (image, Nimage);
   return;
-}
-
-void initMosaicGrid (Image *image, int Nimage) {
-
-  /* find max dR, dD range for all mosaics */
-  /* define mosaic.coords to cover dR, dD */
-  /* send results to initGridBins */
-
-  int i, j, m, NX, NY;
-  int dXmax, dYmax;
-  double dS, dX, dY;
-  double R, D, Rmin, Rmax, Dmin, Dmax;
-  double Mcal, dMcal, Xm;
-
-  dXmax = dYmax = 0.0;
-  for (i = 0; i < Nmosaic; i++) {
-    Dmin = Rmin =  360.0;
-    Dmax = Rmax = -360.0;
-    dS = 0.0;
-    Mcal = dMcal = Xm = 0;
-    for (j = 0; j < Nimlist[i]; j++) {
-      m = imlist[i][j];
-      NX = image[m].NX;
-      NY = image[m].NY;
-      dS += hypot (image[m].coords.cdelt1*image[m].coords.pc1_1, image[m].coords.cdelt1*image[m].coords.pc2_1);
-      XY_to_RD (&R, &D, 0.0, 0.0, &image[m].coords);
-      Rmin = MIN (Rmin, R);
-      Rmax = MAX (Rmax, R);
-      Dmin = MIN (Dmin, D);
-      Dmax = MAX (Dmax, D);
-      XY_to_RD (&R, &D, (double) NX, 0.0, &image[m].coords);
-      Rmin = MIN (Rmin, R);
-      Rmax = MAX (Rmax, R);
-      Dmin = MIN (Dmin, D);
-      Dmax = MAX (Dmax, D);
-      XY_to_RD (&R, &D, (double) NX, (double) NY, &image[m].coords);
-      Rmin = MIN (Rmin, R);
-      Rmax = MAX (Rmax, R);
-      Dmin = MIN (Dmin, D);
-      Dmax = MAX (Dmax, D);
-      XY_to_RD (&R, &D, 0.0, (double) NY, &image[m].coords);
-      Rmin = MIN (Rmin, R);
-      Rmax = MAX (Rmax, R);
-      Dmin = MIN (Dmin, D);
-      Dmax = MAX (Dmax, D);
-      Mcal += image[m].Mcal_PS;
-      dMcal += image[m].dMcal_PS;
-      Xm += image[m].Xm;
-      /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */
-      image[m].Mcal_PS = 0.0;
-      image[m].dMcal_PS = image[m].Xm = NO_MAG;
-    }
-    dS /= Nimlist[i];
-    strcpy (mosaic[i].coords.ctype, "RA---TAN");
-    mosaic[i].coords.crval1 = Rmin;
-    mosaic[i].coords.crval2 = Dmin;
-    mosaic[i].coords.crpix1 = 0.0;
-    mosaic[i].coords.crpix2 = 0.0;
-    mosaic[i].coords.cdelt1 = dS;
-    mosaic[i].coords.cdelt2 = dS;
-    mosaic[i].coords.pc1_1  = 1.0;
-    mosaic[i].coords.pc2_2  = 1.0;
-    mosaic[i].coords.pc1_2  = 0.0;
-    mosaic[i].coords.pc2_1  = 0.0;
-    RD_to_XY (&dX, &dY, Rmax, Dmax, &mosaic[i].coords);
-
-    mosaic[i].Mcal  = Mcal / Nimlist[i];
-    mosaic[i].dMcal = dMcal / Nimlist[i];
-    mosaic[i].Xm    = Xm / Nimlist[i];
-  }
-  if (!USE_GRID) return;
-
-  dXmax = MAX (dXmax, dX);
-  dYmax = MAX (dYmax, dY);
-  initGrid (dXmax, dYmax);
-  return;
-}
-
-void setMcalFinal () {
-
-  int i, j, im, Nimage;
-  Image *image;
-
-  if (!MOSAICNAME[0]) return;
-
-  image = getimages (&Nimage);
-
-  for (i = 0; i < Nmosaic; i++) {
-    for (j = 0; j < Nimlist[i]; j++) {
-      im = imlist[i][j];
-      image[im].Mcal_PS = mosaic[i].Mcal;
-      image[im].dMcal_PS = mosaic[i].dMcal;
-      image[im].Xm = mosaic[i].Xm;
-      image[im].code |= (mosaic[i].code & ID_IMAGE_FEW);
-      image[im].code |= (mosaic[i].code & ID_IMAGE_POOR);
-    }
-  }      
 }
 
@@ -238,5 +140,5 @@
 int findMosaics (Catalog *catalog, int Ncatalog) {
   
-  int i, j, ecode;
+  int i, j;
 
   if (!MOSAICNAME[0]) return (FALSE);
@@ -247,6 +149,4 @@
 	if (catalog[i].measure[j].t > TSTOP) continue;
       }
-      ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[j].source);
-      if (photcode[0].code != ecode) continue;
       matchMosaics (catalog, j, i);
     }
@@ -255,4 +155,5 @@
 }
 
+/* modify this function to use the measure->imageID field ? */
 void matchMosaics (Catalog *catalog, int meas, int cat) {
 
@@ -263,16 +164,4 @@
     if (catalog[cat].measure[meas].t > mosaic[i].stop) continue;
     
-# ifdef GRID_V1
-    if (USE_GRID) {
-      ave = catalog[cat].measure[meas].averef;
-      ra  = catalog[cat].average[ave].R_PS - catalog[cat].measure[meas].dR_PS / 3600.0;
-      dec = catalog[cat].average[ave].D_PS - catalog[cat].measure[meas].dD_PS / 3600.0;
-
-      /* X,Y always positive-definite in range 0,0 - dX, dY */
-      RD_to_XY (&X, &Y, ra, dec, &mosaic[i].coords);
-      setGridMeasure (meas, cat, X, Y);
-    }
-# endif
-
     bin[cat][meas] = i;
 
@@ -288,273 +177,5 @@
     return;
   }
-  fprintf (stderr, "missed measurement\n");
   return;
-}
-
-float getMmos (int meas, int cat) {
-
-  int i;
-  float value;
-
-  if (!MOSAICNAME[0]) return (0);
-  i = bin[cat][meas];
-  if (i == -1) return (NO_MAG);
-
-  if (mosaic[i].code & IMAGE_BAD) return (NO_MAG);  
-  value = mosaic[i].Mcal;
-  return (value);
-}
-
-int setMmos (Catalog *catalog, int PoorImages) {
-
-  int i, j, m, c, n, N, Nmax, mark, bad;
-  float Msys, Mrel, Mcal, Mgrid;
-  double *list, *dlist, *Mlist, *dMlist;
-  StatType stats;
-  Image *image;
-
-  if (!MOSAICNAME[0]) return (FALSE);
-
-  image = getimages (&N);
-
-  if (PoorImages) {
-    IMAGE_BAD = 0;
-  }
-
-  Nmax = 0;
-  for (i = 0; i < Nmosaic; i++) {
-    Nmax = MAX (Nmax, Nlist[i]);
-  }
-  ALLOCATE (list, double, Nmax);
-  ALLOCATE (dlist, double, Nmax);
-  ALLOCATE (Mlist, double, Nmax);
-  ALLOCATE (dMlist, double, Nmax);
-
-  for (i = 0; i < Nmosaic; i++) {
-    
-    /* on PoorImages run, skip good images */
-    if (PoorImages) {
-      bad = mosaic[i].code & (ID_IMAGE_FEW | ID_IMAGE_POOR | ID_IMAGE_SKIP);
-      if (!bad) continue;
-    }      
-
-    N = 0;
-    for (j = 0; j < Nlist[i]; j++) {
-      
-      m = mlist[i][j];
-      c = clist[i][j];
-      
-      if (catalog[c].measure[m].flags & MEAS_BAD) continue;
-      if ((Mcal  = getMcal  (m, c)) == NO_MAG) continue;
-      if ((Mgrid = getMgrid (m, c)) == NO_MAG) continue;
-      if ((Mrel  = getMrel  (catalog, m, c)) == NO_MAG) continue;
-      
-      n = catalog[c].measure[m].averef;
-      Msys = PhotSys (&catalog[c].measure[m], &catalog[c].average[n], &catalog[c].secfilt[n*PhotNsec]);
-      list[N]  = Msys - Mrel - Mcal - Mgrid;
-      dlist[N] = MAX (catalog[c].measure[m].dM_PS, MIN_ERROR);
-      Mlist[N] = Msys;
-      dMlist[N] = list[N];
-      N++;
-    }
-    /* Nlist[i] is all measurements, N is good measurements */
-
-    /* too few good measurements or too many bad measurements (skip in PoorImages run) */
-
-    if (!PoorImages) {
-      mark = (N < IMAGE_TOOFEW) || (N < IMAGE_GOOD_FRACTION*Nlist[i]);
-      if (mark) {
-	fprintf (stderr, "marked image %s (%d), (%d < %d) || (%d < %f*%d)\n", image[imlist[i][0]].name, i, N, IMAGE_TOOFEW, N, IMAGE_GOOD_FRACTION, Nlist[i]);
-	mosaic[i].code |= ID_IMAGE_FEW;
-      } else {
-	mosaic[i].code &= ~ID_IMAGE_FEW;
-      }
-    }
-    liststats (list, dlist, N, &stats);
-    if (PoorImages) fprintf (stderr, "Mmos: %f %f %d %d\n", stats.mean, stats.sigma, stats.Nmeas, N);
-    mosaic[i].Mcal  = stats.mean;
-    mosaic[i].dMcal = stats.sigma;
-    mosaic[i].Xm    = 100.0*log10(stats.chisq);
-  }
-  free (list);
-  free (dlist);
-  free (Mlist);
-  free (dMlist);
-
-  if (PoorImages) {
-    IMAGE_BAD = ID_IMAGE_POOR | ID_IMAGE_FEW | ID_IMAGE_SKIP;
-    STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
-    MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR | ID_MEAS_SKIP | ID_MEAS_AREA;
-  }
-  return (TRUE);
-}
-  
-StatType statsMosaicM (Catalog *catalog) {
-
-  int i, n;
-  double *list, *dlist;
-  StatType stats;
-
-  bzero (&stats, sizeof (StatType));
-  if (!MOSAICNAME[0]) return (stats);
-
-  ALLOCATE (list, double, Nmosaic);
-  ALLOCATE (dlist, double, Nmosaic);
-
-  n = 0;
-  for (i = 0; i < Nmosaic; i++) {
-    if (mosaic[i].code & IMAGE_BAD) continue;
-    list[n] = mosaic[i].Mcal;
-    dlist[n] = 1;
-    n++;
-  }
-
-  liststats (list, dlist, n, &stats);
-  free (list);
-  free (dlist);
-  return (stats);
-}
-
-StatType statsMosaicdM (Catalog *catalog) {
-
-  int i, n;
-  double *list, *dlist;
-  StatType stats;
-
-  bzero (&stats, sizeof (StatType));
-  if (!MOSAICNAME[0]) return (stats);
-
-  ALLOCATE (list, double, Nmosaic);
-  ALLOCATE (dlist, double, Nmosaic);
-
-  n = 0;
-  for (i = 0; i < Nmosaic; i++) {
-
-    if (mosaic[i].code & IMAGE_BAD) continue;
-    list[n] = mosaic[i].dMcal;
-    dlist[n] = 1;
-    n++;
-  }
-
-  liststats (list, dlist, n, &stats);
-  free (list);
-  free (dlist);
-  return (stats);
-}
-
-StatType statsMosaicN (Catalog *catalog) {
-
-  int i, j, m, c, n, N;
-  double *list, *dlist;
-  StatType stats;
-
-  bzero (&stats, sizeof (StatType));
-  if (!MOSAICNAME[0]) return (stats);
-
-  ALLOCATE (list, double, Nmosaic);
-  ALLOCATE (dlist, double, Nmosaic);
-
-  n = 0;
-  for (i = 0; i < Nmosaic; i++) {
-    if (mosaic[i].code & IMAGE_BAD)  continue;
-
-    N = 0;
-    for (j = 0; j < Nlist[i]; j++) {
-
-      m = mlist[i][j];
-      c = clist[i][j];
-
-      if (getMcal  (m, c) == NO_MAG) continue;
-      if (getMgrid (m, c) == NO_MAG) continue;
-      if (getMrel  (catalog, m, c) == NO_MAG) continue;
-      N++;
-    }
-    list[n] = N;
-    dlist[n] = 1;
-    n++;
-  }
-  fprintf (stderr, "Nmosaic: %d, n: %d\n", Nmosaic, n);
-
-  liststats (list, dlist, n, &stats);
-  free (list);
-  free (dlist);
-  return (stats);
-}
-
-StatType statsMosaicX (Catalog *catalog) {
-
-  int i, n;
-  double *list, *dlist;
-  StatType stats;
-
-  bzero (&stats, sizeof (StatType));
-  if (!MOSAICNAME[0]) return (stats);
-
-  ALLOCATE (list, double, Nmosaic);
-  ALLOCATE (dlist, double, Nmosaic);
-
-  n = 0;
-  for (i = 0; i < Nmosaic; i++) {
-
-    if (mosaic[i].code & IMAGE_BAD) continue;
-    list[n] = pow(10.0, 0.01*mosaic[i].Xm);
-    dlist[n] = 1;
-    n++;
-  }
-
-  liststats (list, dlist, n, &stats);
-  free (list);
-  free (dlist);
-  return (stats);
-}
-
-/* mark mosaic if: abs(Mcal - <Mcal>) too large, dMcal too large */
-void clean_mosaics () {
-
-  int i, N, mark, Nmark;
-  double *mlist, *slist, *dlist;
-  double MaxOffset, MedOffset, MaxScatter;
-  StatType stats;
-
-  if (!MOSAICNAME[0]) return;
-
-  if (VERBOSE) fprintf (stderr, "marking poor mosaics\n");
-
-  ALLOCATE (mlist, double, Nmosaic);
-  ALLOCATE (slist, double, Nmosaic);
-  ALLOCATE (dlist, double, Nmosaic);
-
-  for (i = N = 0; i < Nmosaic; i++) {
-    if (mosaic[i].code & IMAGE_BAD) continue;
-    mlist[N] = mosaic[i].Mcal;
-    slist[N] = mosaic[i].dMcal;
-    dlist[N] = 1;
-    N++;
-  }
-  initstats ("MEAN");
-  liststats (mlist, dlist, N, &stats);
-  MaxOffset = MAX (IMAGE_OFFSET, 2*stats.sigma);
-  MedOffset = stats.median;
-  liststats (slist, dlist, N, &stats);
-  MaxScatter = MAX (IMAGE_SCATTER, 2*stats.median);
-  fprintf (stderr, "Mrel: %f, dMrel: %f, Max Scatter: %f, Max Offset: %f\n", MedOffset, stats.median, MaxScatter, MaxOffset);
-  
-  Nmark = 0;
-  for (i = 0; i < Nmosaic; i++) {
-    mark = FALSE;
-    mark = (mosaic[i].dMcal > MaxScatter) || (fabs(mosaic[i].Mcal - MedOffset) > MaxOffset);
-    if (mark) { 
-      Nmark ++;
-      mosaic[i].code |= ID_IMAGE_POOR;
-    } else {
-      mosaic[i].code &= ~ID_IMAGE_POOR;
-    }
-  }
-
-  fprintf (stderr, "%d mosaics marked poor\n", Nmark);
-  initstats (STATMODE);
-  free (mlist);
-  free (slist);
-  free (dlist);
 }
 
@@ -591,6 +212,6 @@
 
       ave = catalog[c].measure[m].averef;
-      xlist[N] = catalog[c].average[ave].R - catalog[c].measure[m].dR_PS / 3600.0;
-      ylist[N] = catalog[c].average[ave].D - catalog[c].measure[m].dD_PS / 3600.0;
+      xlist[N] = catalog[c].average[ave].R - catalog[c].measure[m].dR / 3600.0;
+      ylist[N] = catalog[c].average[ave].D - catalog[c].measure[m].dD / 3600.0;
       N++;
     }
@@ -655,4 +276,12 @@
 }
 
+Image *getMosaicForImage (int Nim) {
+
+    Image *mosaic;
+
+    mosaic = NULL;
+    return mosaic;
+}
+
 int *SelectRefMosaic (Mosaic **refmosaic, int *Nimage) {
 
Index: trunk/Ohana/src/relastro/src/ParFactor.c
===================================================================
--- trunk/Ohana/src/relastro/src/ParFactor.c	(revision 12332)
+++ trunk/Ohana/src/relastro/src/ParFactor.c	(revision 12332)
@@ -0,0 +1,62 @@
+# include "relastro.h"
+
+# if (0)
+/* Low precision formulae for the sun, from Almanac p. C24 (1990) */
+/* ra and dec are returned as decimal hours and decimal degrees. */
+void lpsun (double jd, double *ra, double *dec) {
+
+  double n, L, g, lambda,epsilon,alpha,delta,x,y,z;
+
+  n = jd - J2000;
+  L = 280.460 + 0.9856474 * n;
+  g = (357.528 + 0.9856003 * n)/DEG_IN_RADIAN;
+  lambda = (L + 1.915 * sin(g) + 0.020 * sin(2. * g))/DEG_IN_RADIAN;
+  epsilon = (23.439 - 0.0000004 * n)/DEG_IN_RADIAN;
+
+  // this is the conversion from ecliptic to celestial coords
+  x = cos(lambda);
+  y = cos(epsilon)*sin(lambda);
+  z = sin(epsilon)*sin(lambda);
+
+  *ra = (atan_circ(x,y))*HRS_IN_RADIAN;
+  *dec = (asin(z))*DEG_IN_RADIAN;
+}
+# endif
+
+/* code borrowed from Skycalc : fix this stuff XXX */
+/* Low precision formulae for the sun, from Almanac p. C24 (1990) */
+int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon) {
+
+  double n, L, g;
+
+# define J2000 2451545.       /* Julian date at standard epoch */
+
+  n = jd - J2000;
+  L = 280.460 + 0.9856474 * n;
+  g = (357.528 + 0.9856003 * n)*RAD_DEG;
+  *lambda = L + 1.915 * sin(g) + 0.020 * sin(2. * g); // longitude in degrees
+  *beta = 0.0;					  // approx latitude
+  *epsilon = (23.439 - 0.0000004 * n);		  // obliquity of ecliptic in degrees
+  return TRUE;
+}
+
+/* given RA, DEC, Time, calculate the parallax factor */
+int ParFactor (double *pR, double *pD, double R, double D, time_t T) {
+
+  double jd;
+  double L, B, E, e, s, r, d;
+
+  /* given a time T in UNIX seconds, determine the solar longitude S */
+
+  jd = sec_to_jd (T);
+  sun_ecliptic (jd, &L, &B, &E);
+
+  e = E * RAD_DEG;
+  s = L * RAD_DEG;
+  r = R * RAD_DEG;
+  d = D * RAD_DEG;
+  
+  *pR =  +(cos(e)*sin(s)*cos(r) - cos(s)*sin(r));
+  *pD =  -(cos(e)*sin(s)*sin(r) + cos(s)*cos(r))*sin(d) + sin(e)*sin(s)*cos(d);
+  return TRUE;
+}
Index: trunk/Ohana/src/relastro/src/SetSignals.c
===================================================================
--- trunk/Ohana/src/relastro/src/SetSignals.c	(revision 12067)
+++ trunk/Ohana/src/relastro/src/SetSignals.c	(revision 12332)
@@ -1,3 +1,3 @@
-# include "relphot.h"
+# include "relastro.h"
 
 static int Protect = FALSE;
Index: trunk/Ohana/src/relastro/src/Shutdown.c
===================================================================
--- trunk/Ohana/src/relastro/src/Shutdown.c	(revision 12067)
+++ trunk/Ohana/src/relastro/src/Shutdown.c	(revision 12332)
@@ -1,3 +1,3 @@
-# include "relphot.h"
+# include "relastro.h"
 
 static FITS_DB *db;
@@ -27,2 +27,4 @@
 }
 
+
+/* XXX this is probably not needed anymore : just protect the write statements */
Index: trunk/Ohana/src/relastro/src/StarOps.c
===================================================================
--- trunk/Ohana/src/relastro/src/StarOps.c	(revision 12067)
+++ trunk/Ohana/src/relastro/src/StarOps.c	(revision 12332)
@@ -1,3 +1,3 @@
-# include "relphot.h"
+# include "relastro.h"
 
 static int Nmax;
@@ -28,5 +28,5 @@
   if (catalog[cat].average[ave].code & STAR_BAD) return (NO_MAG);  
  
-  value = PhotPrimary ? catalog[cat].average[ave].M : catalog[cat].secfilt[PhotNsec*ave+PhotSec].M_PS;
+  value = catalog[cat].secfilt[PhotNsec*ave+PhotSec].M;
   return (value);
 }
@@ -54,5 +54,5 @@
 	Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
 	list[N] = Msys - Mcal - Mmos - Mgrid;
-	dlist[N] = MAX (catalog[i].measure[m].dM_PS, MIN_ERROR);
+	dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
 	N++;
       }
@@ -65,13 +65,7 @@
       liststats (list, dlist, N, &stats);
 
-      if (PhotPrimary) {
-	catalog[i].average[j].M = stats.mean;
-	catalog[i].average[j].dM = stats.sigma;
-	catalog[i].average[j].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
-      } else {
-	catalog[i].secfilt[PhotNsec*j+PhotSec].M_PS  = stats.mean;
-	catalog[i].secfilt[PhotNsec*j+PhotSec].dM_PS = stats.sigma;
-	catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
-      }      
+      catalog[i].secfilt[PhotNsec*j+PhotSec].M  = stats.mean;
+      catalog[i].secfilt[PhotNsec*j+PhotSec].dM = stats.sigma;
+      catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
     }
   }
@@ -112,5 +106,5 @@
 	Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
 	list[N] = Msys - Mcal - Mmos - Mgrid;
-	dlist[N] = MAX (catalog[i].measure[m].dM_PS, MIN_ERROR);
+	dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
 	N++;
       }
@@ -121,13 +115,7 @@
 
       /* use sigma or error in dM for output? */
-      if (PhotPrimary) {
-	catalog[i].average[j].M = stats.mean;
-	catalog[i].average[j].dM = MAX (stats.sigma, stats.error);
-	catalog[i].average[j].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
-      } else {
-	catalog[i].secfilt[PhotNsec*j+PhotSec].M_PS  = stats.mean;
-	catalog[i].secfilt[PhotNsec*j+PhotSec].dM_PS = MAX (stats.error, stats.sigma);
-	catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
-      }      
+      catalog[i].secfilt[PhotNsec*j+PhotSec].M  = stats.mean;
+      catalog[i].secfilt[PhotNsec*j+PhotSec].dM = MAX (stats.error, stats.sigma);
+      catalog[i].secfilt[PhotNsec*j+PhotSec].Xm = (stats.Nmeas > 1) ? 100.0*log10(stats.chisq) : NO_MAG;
     }
   }
@@ -156,5 +144,5 @@
 	if ((Mmos  = getMmos  (m, i)) == NO_MAG) continue;
 	if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue;
-	catalog[i].measure[m].Mcal_PS = Mcal + Mmos + Mgrid;
+	catalog[i].measure[m].Mcal = Mcal + Mmos + Mgrid;
       }
     }
@@ -183,9 +171,9 @@
     for (j = 0; j < catalog[i].Naverage; j++) {
       if (catalog[i].average[j].code & STAR_BAD) continue;
-      Xm = PhotPrimary ? catalog[i].average[j].Xm : catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
+      Xm = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
       if (Xm == -1) continue;
       Chisq = pow (10.0, 0.01*Xm);
       xlist[Ntot] = Chisq;
-      slist[Ntot] = PhotPrimary ? catalog[i].average[j].dM : catalog[i].secfilt[PhotNsec*j+PhotSec].dM_PS;
+      slist[Ntot] = catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
       dlist[Ntot] = 1;
       Ntot ++;
@@ -203,6 +191,6 @@
   for (i = 0; i < Ncatalog; i++) {
     for (j = 0; j < catalog[i].Naverage; j++) {
-      dM = PhotPrimary ? catalog[i].average[j].dM : catalog[i].secfilt[PhotNsec*j+PhotSec].dM_PS;
-      Xm = PhotPrimary ? catalog[i].average[j].Xm    : catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
+      dM = catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
+      Xm = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
       Chisq = pow (10.0, 0.01*Xm);
       mark = (dM > MaxScatter) || (Xm == NO_MAG) || (Chisq > MaxChisq);
@@ -269,5 +257,5 @@
 	Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
 	list[N] = Msys - Mcal - Mmos - Mgrid;
-	dlist[N] = MAX (catalog[i].measure[m].dM_PS, MIN_ERROR);
+	dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
 	N++;
       }
@@ -301,5 +289,5 @@
 	Msys = PhotSys (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
 	list[N] = Msys - Mcal - Mmos - Mgrid;
-	dlist[N] = MAX (catalog[i].measure[m].dM_PS, MIN_ERROR);
+	dlist[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
 	ilist[N] = m;
 	N++;
@@ -387,5 +375,5 @@
       if (catalog[i].average[j].code & STAR_BAD) continue;  
 
-      Xm = PhotPrimary ? catalog[i].average[j].Xm : catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
+      Xm = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
       if (Xm == NO_MAG) continue;
       list[n] = pow (10.0, 0.01*Xm);
@@ -423,5 +411,5 @@
       if (catalog[i].average[j].code & STAR_BAD) continue;  
 
-      dM = PhotPrimary ? catalog[i].average[j].dM : catalog[i].secfilt[PhotNsec*j+PhotSec].dM_PS;
+      dM = catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
       list[n] = dM;
       dlist[n] = 1;
@@ -452,5 +440,5 @@
     for (j = 0; j < catalog[i].Naverage; j++) {
       if (catalog[i].average[j].code & STAR_BAD) continue;  
-      dMrel = PhotPrimary ? catalog[i].average[j].dM : catalog[i].secfilt[PhotNsec*j+PhotSec].dM_PS;
+      dMrel = catalog[i].secfilt[PhotNsec*j+PhotSec].dM;
       bin = dMrel / 0.00025;
       bin = MAX (0, MIN (NBIN-1, bin));
@@ -483,6 +471,6 @@
     for (j = 0; j < catalog[i].Naverage; j++) {
       if (catalog[i].average[j].code & STAR_BAD) continue;
-      xlist[N] = PhotPrimary ? catalog[i].average[j].M : catalog[i].secfilt[PhotNsec*j+PhotSec].M_PS;
-      value    = PhotPrimary ? catalog[i].average[j].Xm   : catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
+      xlist[N] = catalog[i].secfilt[PhotNsec*j+PhotSec].M;
+      value    = catalog[i].secfilt[PhotNsec*j+PhotSec].Xm;
       if (value == NO_MAG) continue;
       ylist[N] = 0.01*value;
Index: trunk/Ohana/src/relastro/src/UpdateChips.c
===================================================================
--- trunk/Ohana/src/relastro/src/UpdateChips.c	(revision 12332)
+++ trunk/Ohana/src/relastro/src/UpdateChips.c	(revision 12332)
@@ -0,0 +1,36 @@
+# include "relastro.h"
+
+int UpdateChips (Catalog *catalog, int Ncatalog) {
+
+  /* we can measure new image parameters for each non-mosaic chip independently */
+  int i, Nimage, Nstars;
+  Image *image;
+  StarData *raw, *ref;
+
+  image = getimages (&Nimage);
+
+  for (i = 0; i < Nimage; i++) {
+
+    /* skip all except WRP images */
+    if (strcmp(&image[i].coords.ctype[4], "-WRP")) continue;
+
+    /* XXX I probably need modes for getImageRaw/Ref to distinguish SIMPLE, CHIP, MOSAIC */
+
+    /* convert measure coordinates to raw entries */
+    raw = getImageRaw (catalog, Ncatalog, i, &Nstars, TRUE);
+
+    /* convert average coordinates to ref entries */
+    ref = getImageRef (catalog, Ncatalog, i, &Nstars, TRUE);
+
+    /* XXX grab this code from mosastro */
+    /* XXX check that it matches psastro */
+    /* XXX this may need a 'mosaic' and 'non-mosaic' version */
+    FitChip (raw, ref, Nstars, &image[i].coords);
+
+    free (raw);
+    free (ref);
+  }
+
+  return (TRUE);
+}
+
Index: trunk/Ohana/src/relastro/src/UpdateImages.c
===================================================================
--- trunk/Ohana/src/relastro/src/UpdateImages.c	(revision 12067)
+++ trunk/Ohana/src/relastro/src/UpdateImages.c	(revision 12332)
@@ -42,6 +42,6 @@
 	/* the astrometric errors are not being carried yet (but should be!) */
 	/* we use the photometric mag error as a weighting term */
-	dlistR[N] = MAX (catalog[i].measure[m].dM_PS, MIN_ERROR);
-	dlistD[N] = MAX (catalog[i].measure[m].dM_PS, MIN_ERROR);
+	dlistR[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
+	dlistD[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
 	N++;
       }
Index: trunk/Ohana/src/relastro/src/UpdateMosaic.c
===================================================================
--- trunk/Ohana/src/relastro/src/UpdateMosaic.c	(revision 12332)
+++ trunk/Ohana/src/relastro/src/UpdateMosaic.c	(revision 12332)
@@ -0,0 +1,36 @@
+# include "relastro.h"
+
+int UpdateMosaic (Catalog *catalog, int Ncatalog) {
+
+  /* we can measure new image parameters for each mosaic independently */
+  int i, Nimage, Nstars;
+  Image *image;
+  StarData *raw, *ref;
+
+  image = getimages (&Nimage);
+
+  for (i = 0; i < Nimage; i++) {
+
+    /* skip all except DIS images */
+    if (strcmp(&image[i].coords.ctype[4], "-DIS")) continue;
+
+    /* XXX I probably need modes for getImageRaw/Ref to distinguish SIMPLE, CHIP, MOSAIC */
+
+    /* convert measure coordinates to raw entries */
+    raw = getImageRaw (catalog, Ncatalog, i, &Nstars, TRUE);
+
+    /* convert average coordinates to ref entries */
+    ref = getImageRef (catalog, Ncatalog, i, &Nstars, TRUE);
+
+    /* XXX grab this code from mosastro */
+    /* XXX check that it matches psastro */
+    /* XXX this may need a 'mosaic' and 'non-mosaic' version */
+    FitMosaic (raw, ref, Nstars, &image[i].coords);
+
+    free (raw);
+    free (ref);
+  }
+
+  return (TRUE);
+}
+
Index: trunk/Ohana/src/relastro/src/UpdateObjects.c
===================================================================
--- trunk/Ohana/src/relastro/src/UpdateObjects.c	(revision 12067)
+++ trunk/Ohana/src/relastro/src/UpdateObjects.c	(revision 12332)
@@ -2,6 +2,12 @@
 
 static int Nmax;
-static double *listR, *listD;
-static double *dlistR, *dlistD;
+static double *X, *dX;
+static double *Y, *dY;
+static double *R, *dR;
+static double *D, *dD;
+static double *pX;
+static double *pY;
+static time_t *T;
+static double *dT;
 
 void initObjectData (Catalog *catalog, int Ncatalog) {
@@ -16,6 +22,18 @@
   }
 
-  ALLOCATE (list, double, MAX (1, Nmax));
-  ALLOCATE (dlist, double, MAX (1, Nmax));
+  ALLOCATE (R, double, MAX (1, Nmax));
+  ALLOCATE (D, double, MAX (1, Nmax));
+  ALLOCATE (T, time_t, MAX (1, Nmax));
+  ALLOCATE (X, double, MAX (1, Nmax));
+  ALLOCATE (Y, double, MAX (1, Nmax));
+
+  ALLOCATE (dR, double, MAX (1, Nmax));
+  ALLOCATE (dD, double, MAX (1, Nmax));
+  ALLOCATE (dT, double, MAX (1, Nmax));
+  ALLOCATE (dX, double, MAX (1, Nmax));
+  ALLOCATE (dY, double, MAX (1, Nmax));
+
+  ALLOCATE (pX, double, MAX (1, Nmax));
+  ALLOCATE (pY, double, MAX (1, Nmax));
 }  
 
@@ -23,6 +41,19 @@
 
   int i, j, k, m, N;
-  float chisq;
   StatType statsR, statsD;
+  Coords coords;
+  PMFit fit;
+
+  initObjectData (catalog, Ncatalog);
+
+  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++) {
@@ -31,38 +62,102 @@
       /* calculate the average value of R,D for a single star */
       if (catalog[i].average[j].code & STAR_BAD) continue;  
-      m = catalog[i].average[j].offset;
 
       N = 0;
+      m = catalog[i].average[j].offset;
       for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
 	if (catalog[i].measure[m].flags & MEAS_BAD) continue;
 	
-	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]);
+	R[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
+	D[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
+	T[N] = catalog[i].measure[m].t;
 
 	/* the astrometric errors are not being carried yet (but should be!) */
 	/* we use the photometric mag error as a weighting term */
-	dlistR[N] = MAX (catalog[i].measure[m].dM_PS, MIN_ERROR);
-	dlistD[N] = MAX (catalog[i].measure[m].dM_PS, MIN_ERROR);
+	dR[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
+	dD[N] = MAX (catalog[i].measure[m].dM, MIN_ERROR);
+	dT[N] = catalog[i].measure[m].dt;
+
 	N++;
       }
+
+      // XXX This criterion needs to be better considered: adjust to match Ndof
       if (N < STAR_TOOFEW) { /* too few measurements */
 	catalog[i].average[j].code |= ID_STAR_FEW;
+	continue;
       } else {
 	catalog[i].average[j].code &= ~ID_STAR_FEW;
       }	
 
-      /* in here, we should be fitting the parallax and proper-motion components */
+      /* we need to do the fit in a locally linear space; choose a ref coordinate */
+      coords.crval1 = R[0];
+      coords.crval2 = D[0];
+      
+      /* project all of the R,D coordinates to a plane centered on this coordinate */
+      for (k = 0; k < N; k++) {
+	  RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords);
+      }	  
 
-      liststats (listR, dlistR, N, &statsR);
-      liststats (listR, dlistR, N, &statsD);
+      /* fit the model components as needed */
+      switch (FIT_MODE) {
+	case FIT_AVERAGE:
+	  liststats (R, dR, N, &statsR);
+	  liststats (D, dD, N, &statsD);
 
-      catalog[i].average[j].R = statsR.mean;
-      catalog[i].average[j].dR = statsR.sigma;
+	  fit.Ro = statsR.mean;
+	  fit.dRo = statsR.sigma;
 
-      catalog[i].average[j].D = statsD.mean;
-      catalog[i].average[j].dD = statsD.sigma;
+	  fit.Do = statsD.mean;
+	  fit.dDo = statsD.sigma;
 
-      chisq = 0.5*(statsR.chisq + statsD.chisq);
-      catalog[i].average[j].Xp = (statsR.Nmeas > 1) ? 100.0*log10(chisq) : NO_MAG;
+	  fit.chisq = 0.5*(statsR.chisq + statsD.chisq);
+	  fit.Nfit = N;
+	  break;
+
+	case FIT_PM_ONLY:
+	  FitPM (&fit, X, dX, Y, dY, T, N);
+	  break;
+
+	case FIT_PM_AND_PAR:
+	  for (k = 0; k < N; k++) {
+	    ParFactor (&pX[k], &pY[k], R[k], D[k], T[k]);
+	  }
+	  FitPMandPar (&fit, X, dX, Y, dY, T, pX, pY, N);
+	  break;
+	default:
+	  fprintf (stderr, "programming error at %s, %s", __FILE__, __LINE__);
+	  exit (2);
+      }	  
+
+      if (0 && (j < 100)) {
+	  fprintf (stderr, "%f %f -> %f %f (%f,%f)\n",
+		   catalog[i].average[j].R, 
+		   catalog[i].average[j].D, 
+		   fit.Ro, fit.Do, 
+		   3600*(catalog[i].average[j].R - fit.Ro), 
+		   3600*(catalog[i].average[j].D - fit.Do));
+      }
+
+      // the measure fields must be updated before the average fields
+      m = catalog[i].average[j].offset;
+      for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
+	if (catalog[i].measure[m].flags & MEAS_BAD) continue;
+	setMeanR (fit.Ro, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
+	setMeanD (fit.Do, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]);
+      }      
+
+      catalog[i].average[j].R  	= fit.Ro;
+      catalog[i].average[j].D  	= fit.Do;
+      catalog[i].average[j].dR 	= fit.dRo;
+      catalog[i].average[j].dD 	= fit.dDo;
+
+      catalog[i].average[j].uR  = fit.uR;
+      catalog[i].average[j].uD  = fit.uD;
+      catalog[i].average[j].duR = fit.duR;
+      catalog[i].average[j].duD = fit.duD;
+
+      catalog[i].average[j].P   = fit.p;
+      catalog[i].average[j].dP  = fit.dp;
+
+      catalog[i].average[j].Xp  = (fit.Nfit > 1) ? 100.0*log10(fit.chisq) : NO_MAG;
 
     }
@@ -71,2 +166,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: trunk/Ohana/src/relastro/src/UpdateSimple.c
===================================================================
--- trunk/Ohana/src/relastro/src/UpdateSimple.c	(revision 12332)
+++ trunk/Ohana/src/relastro/src/UpdateSimple.c	(revision 12332)
@@ -0,0 +1,37 @@
+# include "relastro.h"
+
+int UpdateSimple (Catalog *catalog, int Ncatalog) {
+
+  /* we can measure new image parameters for each non-mosaic chip independently */
+  int i, Nimage, Nstars;
+  Image *image;
+  StarData *raw, *ref;
+
+  image = getimages (&Nimage);
+
+  for (i = 0; i < Nimage; i++) {
+
+    /* skip WRP and DIS images */
+    if (!strcmp(&image[i].coords.ctype[4], "-WRP")) continue;
+    if (!strcmp(&image[i].coords.ctype[4], "-DIS")) continue;
+
+    /* XXX I probably need modes for getImageRaw/Ref to distinguish SIMPLE, CHIP, MOSAIC */
+
+    /* convert measure coordinates to raw entries */
+    raw = getImageRaw (catalog, Ncatalog, i, &Nstars, FALSE);
+
+    /* convert average coordinates to ref entries */
+    ref = getImageRef (catalog, Ncatalog, i, &Nstars, FALSE);
+
+    /* XXX grab this code from mosastro */
+    /* XXX check that it matches psastro */
+    /* XXX this may need a 'mosaic' and 'non-mosaic' version */
+    FitSimple (raw, ref, Nstars, &image[i].coords);
+
+    free (raw);
+    free (ref);
+  }
+
+  return (TRUE);
+}
+
Index: trunk/Ohana/src/relastro/src/args.c
===================================================================
--- trunk/Ohana/src/relastro/src/args.c	(revision 12067)
+++ trunk/Ohana/src/relastro/src/args.c	(revision 12332)
@@ -1,3 +1,3 @@
-# include "relphot.h"
+# include "relastro.h"
 void usage (void);
 
@@ -6,4 +6,54 @@
   int N;
   double trange;
+
+  /* possible operations */
+  FIT_TARGET = TARGET_NONE;
+  FIT_MODE = FIT_AVERAGE;
+  if ((N = get_argument (argc, argv, "-update-objects"))) {
+    remove_argument (N, &argc, argv);
+    FIT_TARGET = TARGET_OBJECTS;
+
+    // check for object fitting modes (not valid for images)
+    if ((N = get_argument (argc, argv, "-pm"))) {
+	remove_argument (N, &argc, argv);
+	FIT_MODE = FIT_PM_ONLY;
+    }
+    if ((N = get_argument (argc, argv, "-pmpar"))) {
+	remove_argument (N, &argc, argv);
+	FIT_MODE = FIT_PM_AND_PAR;
+    }
+  }
+  if ((N = get_argument (argc, argv, "-update-simple"))) {
+    remove_argument (N, &argc, argv);
+    FIT_TARGET = TARGET_SIMPLE;
+  }
+  if ((N = get_argument (argc, argv, "-update-chips"))) {
+    remove_argument (N, &argc, argv);
+    FIT_TARGET = TARGET_CHIPS;
+  }
+  if ((N = get_argument (argc, argv, "-update-mosaics"))) {
+    remove_argument (N, &argc, argv);
+    FIT_TARGET = TARGET_MOSAICS;
+  }
+  if (FIT_TARGET == TARGET_NONE) usage();
+
+  /* specify portion of the sky : allow default of all sky? */
+  UserPatch.Rmin = 0;
+  UserPatch.Rmax = 360;
+  UserPatch.Dmin = -90;
+  UserPatch.Dmax = +90;
+  if ((N = get_argument (argc, argv, "-region"))) {
+    remove_argument (N, &argc, argv);
+    UserPatch.Rmin = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    UserPatch.Rmax = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    UserPatch.Dmin = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    UserPatch.Dmax = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  } else {
+    usage ();
+  }
 
   /* define time */
@@ -34,21 +84,9 @@
   }
 
-  /* specify portion of the sky */
-  UserPatch.Rmin = 0;
-  UserPatch.Rmax = 360;
-  UserPatch.Dmin = -90;
-  UserPatch.Dmax = +90;
-  UserPatchSelect = FALSE;
-  if ((N = get_argument (argc, argv, "-region"))) {
-    remove_argument (N, &argc, argv);
-    UserPatch.Rmin = atof (argv[N]);
-    remove_argument (N, &argc, argv);
-    UserPatch.Rmax = atof (argv[N]);
-    remove_argument (N, &argc, argv);
-    UserPatch.Dmin = atof (argv[N]);
-    remove_argument (N, &argc, argv);
-    UserPatch.Dmax = atof (argv[N]);
-    remove_argument (N, &argc, argv);
-    UserPatchSelect = TRUE;
+  PHOTCODE_LIST = NULL;
+  if ((N = get_argument (argc, argv, "-photcode"))) {
+    remove_argument (N, &argc, argv);
+    PHOTCODE_LIST = strcreate(argv[N]);
+    remove_argument (N, &argc, argv);
   }
 
@@ -79,11 +117,4 @@
   }
 
-  NLOOP = 8;
-  if ((N = get_argument (argc, argv, "-n"))) {
-    remove_argument (N, &argc, argv);
-    NLOOP = atof (argv[N]);
-    remove_argument (N, &argc, argv);
-  }
-
   RESET = FALSE;
   if ((N = get_argument (argc, argv, "-reset"))) {
@@ -104,4 +135,5 @@
   }
 
+  /* XXX update these for relevant plots */
   PlotMmin = 10.0; PlotMmax = 20.0; PlotdMmin = -1.0; PlotdMmax = 1.0;
   if ((N = get_argument (argc, argv, "-plrange"))) {
@@ -117,28 +149,5 @@
   }
 
-  /* group images by mosaic, find Mmos */
-  MOSAICNAME[0] = 0;
-  if ((N = get_argument (argc, argv, "-mosaic"))) {
-    remove_argument (N, &argc, argv);
-    strcpy (MOSAICNAME, argv[N]);
-    remove_argument (N, &argc, argv);
-  }
-
-  FREEZE_IMAGES = FALSE;
-  if ((N = get_argument (argc, argv, "-imfreeze"))) {
-    remove_argument (N, &argc, argv);
-    FREEZE_IMAGES = TRUE;
-  }
-
-  USE_GRID = FALSE;
-  if ((N = get_argument (argc, argv, "-grid"))) {
-    remove_argument (N, &argc, argv);
-    USE_GRID = TRUE;
-    if (!MOSAICNAME[0]) {
-      fprintf (stderr, "-grid is only valid with -mosaic\n");
-      exit (2);
-    }
-  }
-
+  /* XXX update this */
   MIN_ERROR = 0.001;
   if ((N = get_argument (argc, argv, "-minerror"))) {
@@ -172,41 +181,32 @@
     ImagSelect = TRUE;
   }
-
-  /* possible operations */
-  UpdateObjects = FALSE;
-  if ((N = get_argument (argc, argv, "-update-objects"))) {
-    remove_argument (N, &argc, argv);
-    UpdateObjects = TRUE;
-  }
-  UpdateImages = FALSE;
-  if ((N = get_argument (argc, argv, "-update-images"))) {
-    remove_argument (N, &argc, argv);
-    UpdateImages = TRUE;
-  }
-
-  if ( UserPatchSelect && (argc != 2)) usage ();
-  if (!UserPatchSelect && (argc != 3)) usage ();
-
+  
+  /* XXX drop this? */
+  DophotSelect = FALSE;
+  if ((N = get_argument (argc, argv, "-dophot"))) {
+    remove_argument (N, &argc, argv);
+    DophotValue = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    DophotSelect = TRUE;
+  }
+
+  if (argc != 1) usage ();
   return TRUE;
 }
 
 void usage () {
-  fprintf (stderr, "ERROR: USAGE: relphot (region) (photcode)\n");
-  fprintf (stderr, "       or:    relphot (photcode) -region RA RA DEC DEC\n");
-  fprintf (stderr, "  options: \n");
+  fprintf (stderr, "ERROR: USAGE: relastro -region RA RA DEC DEC\n");
+  fprintf (stderr, "  working options: \n");
   fprintf (stderr, "  -time (start) (stop)\n");
   fprintf (stderr, "  -v\n");
   fprintf (stderr, "  -plot\n");
   fprintf (stderr, "  -plotdelay (seconds)\n");
-  fprintf (stderr, "  -statmode (mode)\n");
-  fprintf (stderr, "  -n (nloop)\n");
-  fprintf (stderr, "  -reset\n");
   fprintf (stderr, "  -update\n");
   fprintf (stderr, "  -params\n");
-  fprintf (stderr, "  -mosaic (mosaic)\n");
-  fprintf (stderr, "  -imfreeze\n");
-  fprintf (stderr, "  -grid\n");
+  fprintf (stderr, "  -reset\n");
   fprintf (stderr, "  -area Xmin Xmax Ymin Ymax\n");
-  fprintf (stderr, "  -instmag min max\n");
+  fprintf (stderr, "  -instmag min max\n\n");
+  fprintf (stderr, "  planned options: \n");
+  fprintf (stderr, "  -photcode code[,code,code,..]\n");
   fprintf (stderr, "  \n");
   exit (2);
Index: trunk/Ohana/src/relastro/src/bcatalog.c
===================================================================
--- trunk/Ohana/src/relastro/src/bcatalog.c	(revision 12067)
+++ trunk/Ohana/src/relastro/src/bcatalog.c	(revision 12332)
@@ -1,7 +1,7 @@
-# include "relphot.h"
+# include "relastro.h"
 
 int bcatalog (Catalog *subcatalog, Catalog *catalog) {
   
-  int i, j, offset, ecode;
+  int i, j, k, offset, found;
   int NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm;
   float mag;
@@ -19,11 +19,4 @@
     if (catalog[0].average[i].Nm < 2) continue; 
 
-    /* XXX this limitation is absurd
-    if (catalog[0].average[i].R < fullregion[0].RA[0]) continue;
-    if (catalog[0].average[i].R > fullregion[0].RA[1]) continue;
-    if (catalog[0].average[i].D < fullregion[0].DEC[0]) continue;
-    if (catalog[0].average[i].D > fullregion[0].DEC[1]) continue;
-    */
-
     /* start with all stars good */
     subcatalog[0].average[Naverage] = catalog[0].average[i];
@@ -34,9 +27,5 @@
 
     if (RESET) {
-      float *p;
-      p = (PhotPrimary) ? &subcatalog[0].average[Naverage].M : &subcatalog[0].secfilt[PhotNsec*Naverage+PhotSec].M_PS;
-      *p = NO_MAG;
-      p = (PhotPrimary) ? &subcatalog[0].average[Naverage].dM : &subcatalog[0].secfilt[PhotNsec*Naverage+PhotSec].dM_PS;
-      *p = NO_MAG;
+      // XXX reset the ra,dec coords?
       subcatalog[0].average[Naverage].code &= ~ID_STAR_FEW;
       subcatalog[0].average[Naverage].code &= ~ID_STAR_POOR;
@@ -48,7 +37,12 @@
       offset = catalog[0].average[i].offset + j;
 
-      /* select measurements by photcode */
-      ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[offset].source);
-      if (ecode != photcode[0].code) continue;
+      /* select measurements by photcode, if specified */
+      if (Nphotcodes > 0) {
+	found = FALSE;
+	for (k = 0; (k < Nphotcodes) && !found; k++) {
+	  if (photcodes[k][0].code == catalog[0].measure[offset].photcode) found = TRUE;
+	}
+	if (!found) continue;
+      }
 
       /* select measurements by time */
@@ -59,5 +53,5 @@
 
       /* select measurements by quality */
-      if (catalog[0].measure[offset].dophot != 1) continue;
+      if (DophotSelect && (catalog[0].measure[offset].dophot != DophotValue)) continue;
 
       /* select measurements by mag limit */
@@ -66,5 +60,5 @@
 
       /* select measurements by measurement error */
-      if (catalog[0].measure[offset].dM_PS > SIGMA_LIM) continue;
+      if ((SIGMA_LIM > 0) && (catalog[0].measure[offset].dM > SIGMA_LIM)) continue;
 
       /* select measurements by mag limit */
@@ -79,5 +73,5 @@
       subcatalog[0].measure[Nmeasure].averef = Naverage;
       if (RESET) { 
-	subcatalog[0].measure[Nmeasure].Mcal_PS = 0;
+	subcatalog[0].measure[Nmeasure].Mcal = 0;
 	subcatalog[0].measure[Nmeasure].flags &= 0xff00;
 	subcatalog[0].measure[Nmeasure].flags &= ~ID_MEAS_POOR;
@@ -111,8 +105,4 @@
   subcatalog[0].Nsecfilt = catalog[0].Nsecfilt;
 
-  if (catalog[0].Naverage > 0) free (catalog[0].average);
-  if (catalog[0].Nmeasure > 0) free (catalog[0].measure);
-  if (catalog[0].Naverage > 0) free (catalog[0].secfilt);
-
   if (VERBOSE) {
     fprintf (stderr, "%d: using %d stars (%d measures) for catalog\n", i, 
Index: trunk/Ohana/src/relastro/src/dvo_astrom_ops.c
===================================================================
--- trunk/Ohana/src/relastro/src/dvo_astrom_ops.c	(revision 12332)
+++ trunk/Ohana/src/relastro/src/dvo_astrom_ops.c	(revision 12332)
@@ -0,0 +1,91 @@
+# 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);
+}
+
+int setMeanR (double ra_fit, Measure *measure, Average *average, SecFilt *secfilt) {
+
+  /* math to get from new fitted position to new measure offset
+     ra_obs = average[0].R - measure[0].dR / 3600.0;
+     measure[0].dR = (ra_fit - ra_obs) * 3600.0;
+     measure[0].dR = (ra_fit - average[0].R + measure[0].dR/3600) * 3600.0
+     measure[0].dR = (ra_fit - average[0].R) * 3600.0 + measure[0].dR;
+  */
+
+  /* the measure carries the instantaneous mean position at the epoch t */ 
+  measure[0].dR += (ra_fit - average[0].R) * 3600.0;
+
+  /* possible corrections to mean ra:
+
+  - proper-motion and parallax
+  - abberation
+  - precession and nutation, etc
+  - refraction
+  - DCR
+
+  */
+
+  return (TRUE);
+}
+
+int setMeanD (double dec_fit, Measure *measure, Average *average, SecFilt *secfilt) {
+
+  /* math to get from new fitted position to new measure offset
+     dec_obs = average[0].D - measure[0].dD / 3600.0;
+     measure[0].dD = (dec_fit - dec_obs) * 3600.0;
+     measure[0].dD = (dec_fit - average[0].D + measure[0].dD/3600) * 3600.0
+     measure[0].dD = (dec_fit - average[0].D) * 3600.0 + measure[0].dD;
+  */
+
+  /* the measure carries the instantaneous mean position at the epoch t */ 
+  measure[0].dD += (dec_fit - average[0].D) * 3600.0;
+
+  /* possible corrections to mean ra:
+
+  - proper-motion and parallax
+  - abberation
+  - precession and nutation, etc
+  - refraction
+  - DCR
+
+  */
+
+  return (TRUE);
+}
Index: trunk/Ohana/src/relastro/src/fitpoly.c
===================================================================
--- trunk/Ohana/src/relastro/src/fitpoly.c	(revision 12332)
+++ trunk/Ohana/src/relastro/src/fitpoly.c	(revision 12332)
@@ -0,0 +1,295 @@
+# include "relastro.h"
+
+/* these functions support simultaneous 2D fits to
+   x2 = \sum a_i,j x1^i y1^j 
+   y2 = \sum b_i,j x1^i y1^j 
+
+   the order of the fit (largest coefficient) is fixed to a single
+   value for both x1,y1 terms and for x2,y2 fits
+
+   the code is currently confusing because we limit to i+j <= order.
+   this could be cleaner if we used masks and allowed i <= order, j <= order
+*/
+
+double **array_init (int Nx, int Ny) {
+
+  int i;
+  double **array;
+
+  ALLOCATE (array, double *, Nx);
+  for (i = 0; i < Nx; i++) {
+    ALLOCATE (array[i], double, Ny);
+    memset (array[i], 0, Ny*sizeof(double));
+  }    
+  return (array);
+}
+
+void array_free (double **array, int Nx) {
+
+  int i;
+
+  for (i = 0; i < Nx; i++) {
+    free (array[i]);
+  }    
+  free (array);
+}
+
+// XXX define a fit structure and drop the file static variables?
+CoordFit *fit_init (int order) {
+
+  CoordFit *fit;
+
+  ALLOCATE (fit, CoordFit, 1);
+
+  fit[0].Npts   = 0;
+  fit[0].Norder = order;
+  fit[0].Nterms = order + 1;
+  fit[0].Nsums  = 2*order + 1;
+  fit[0].Nelems = SQ(order + 1);
+
+  /* summing arrays for fit solution */
+
+  // xsum[i][j] holds \sum (x2 wt x1^i y1^j)
+  // ysum[i][j] holds \sum (y2 wt x1^i y1^j)
+  fit[0].xsum = array_init (fit[0].Nterms, fit[0].Nterms);
+  fit[0].ysum = array_init (fit[0].Nterms, fit[0].Nterms);
+
+  // xfit[i][j] holds x2 coeff for x1^i y1^j
+  // yfit[i][j] holds y2 coeff for x1^i y1^j
+  fit[0].xfit = array_init (fit[0].Nterms, fit[0].Nterms);
+  fit[0].yfit = array_init (fit[0].Nterms, fit[0].Nterms);
+
+  // sum[i][j] holds \sum (wt x1^i y1^j)
+  fit[0].sum = array_init (fit[0].Nsums, fit[0].Nsums);
+
+  return (fit);
+}
+
+void fit_free (CoordFit *fit) {
+
+  array_free (fit[0].xfit, fit[0].Nterms);
+  array_free (fit[0].yfit, fit[0].Nterms);
+
+  array_free (fit[0].sum, fit[0].Nsums);
+  array_free (fit[0].xsum, fit[0].Nterms);
+  array_free (fit[0].ysum, fit[0].Nterms);
+  
+  free (fit);
+}
+  
+// XXX use implicit masks as below or explicit masks (with function to set?)
+// XXX eg, add a global mask to this file and 
+void fit_add (CoordFit *fit, double x1, double y1, double x2, double y2, double wt) {
+
+  int ix, iy;
+  double xterm, yterm, term;
+
+  xterm = 1;
+  for (ix = 0; ix < fit[0].Nsums; ix++) {
+    yterm = 1;
+    for (iy = 0; iy < fit[0].Nsums; iy++) {
+      term = xterm*yterm*wt;
+      fit[0].sum[ix][iy] += term;
+      if ((iy < fit[0].Nterms) && (ix < fit[0].Nterms)) {
+	fit[0].xsum[ix][iy] += x2*term;
+	fit[0].ysum[ix][iy] += y2*term;
+      }
+      yterm *= y1;
+    }
+    xterm *= x1;
+  }
+  fit[0].Npts ++;
+}
+
+/* convert the xsum,ysum,sum terms into vector,matrix and solve */
+void fit_eval (CoordFit *fit) {
+
+  int i, j, ix, iy, jx, jy;
+  double **matrix, **vector;
+
+  if (fit[0].Npts == 0) {
+    fprintf (stderr, "warning: no valid pts\n");
+  }
+
+  // matrix, vector hold the final linear system
+  matrix = array_init (fit[0].Nelems, fit[0].Nelems);
+  vector = array_init (fit[0].Nelems, 2);
+
+  /* remap the xsum,ysum terms into the vector */
+  for (i = 0; i < fit[0].Nelems; i++) {
+    ix = i % fit[0].Nterms;
+    iy = i / fit[0].Nterms;
+    vector[i][0] = fit[0].xsum[ix][iy];
+    vector[i][1] = fit[0].ysum[ix][iy];
+
+    for (j = 0; j < fit[0].Nelems; j++) {
+      jx = j % fit[0].Nterms;
+      jy = j / fit[0].Nterms;
+      matrix[i][j] = fit[0].sum[ix+jx][iy+jy];
+    }
+  }
+
+  dgaussj (matrix, fit[0].Nelems, vector, 2); 
+
+  for (i = 0; i < fit[0].Nelems; i++) {
+    ix = i % fit[0].Nterms;
+    iy = i / fit[0].Nterms;
+    fprintf (stderr, "x2 : x^%dy^%d: %10.4g    y2 : x^%dy^%d: %10.4g \n", 
+	     ix, iy, vector[i][0], ix, iy, vector[i][1]);
+  }	
+
+  /* remap the vector terms into xfit,yfit */
+  for (i = 0; i < fit[0].Nelems; i++) {
+    ix = i % fit[0].Nterms;
+    iy = i / fit[0].Nterms;
+    fit[0].xfit[ix][iy] = vector[i][0];
+    fit[0].yfit[ix][iy] = vector[i][1];
+  }
+
+  array_free (matrix, fit[0].Nelems);
+  array_free (vector, fit[0].Nelems);
+}
+
+// Nx, Ny is the number of terms in x and in y
+double **poly2d_dx (double **poly, int Nx, int Ny) {
+
+  int i, j, Nxout, Nyout;
+  double **dpoly;
+
+  Nxout = Nx - 1;
+  Nyout = Ny;
+
+  // poly[i][j] holds coeff for x1^i y1^j
+  dpoly = array_init (Nxout, Nyout);
+
+  for (i = 0; i < Nxout; i++) {
+    for (j = 0; j < Nyout; j++) {
+      dpoly[i][j] = poly[i+1][j] * (i+1);
+    }
+  }
+  return dpoly;
+}
+
+// Nx, Ny is the number of terms in x and in y
+double **poly2d_dy (double **poly, int Nx, int Ny) {
+
+  int i, j, Nxout, Nyout;
+  double **dpoly;
+
+  Nxout = Nx;
+  Nyout = Ny - 1;
+
+  // poly[i][j] holds coeff for x1^i y1^j
+  dpoly = array_init (Nxout, Nyout);
+
+  for (i = 0; i < Nxout; i++) {
+    for (j = 0; j < Nyout; j++) {
+      dpoly[i][j] = poly[i][j+1] * (j+1);
+    }
+  }
+  return dpoly;
+}
+
+// Nx, Ny is the number of terms in x and in y
+double **poly2d_copy (double **poly, int Nx, int Ny) {
+
+  int i, j;
+  double **out;
+
+  // poly[i][j] holds coeff for x1^i y1^j
+  out = array_init (Nx, Ny);
+
+  for (i = 0; i < Nx; i++) {
+    for (j = 0; j < Ny; j++) {
+      out[i][j] = poly[i][j+1];
+    }
+  }
+  return out;
+}
+
+double poly2d_eval (double **poly, int Nx, int Ny, double x, double y) {
+
+  int i, j;
+  double xsum, ysum, sum;
+
+  sum = 0;
+  xsum = ysum = 1.0;
+
+  for (i = 0; i < Nx; i++) {
+    ysum = xsum;
+    for (j = 0; j < Ny; j++) {
+      sum += poly[i][j] * ysum;
+      ysum *= y;
+    }
+    xsum *= x;
+  }
+  return (sum);
+}
+
+/* linear portion of fit : NORDER is 1 */
+/* this should only apply to the polynomial, not the projection terms */
+/* compare with psastro supporting code */
+void fit_apply_coords (CoordFit *fit, Coords *coords) {
+
+  double c11, c12;
+  double c21, c22;
+  double Xo, Yo, R;
+  CoordFit *modfit;
+
+  /* I have L,M = fit(X,Y). set corresponding terms for coords */
+
+  // L = a[0][0] + a[1][0]x^1 y^0 + a[0][1] x^0 y^1 + ...
+  // L = pc1_1*cd1*(x - cp1) + pc1_2*cd2*(y - cp2) + ...
+
+  CoordsGetCenter (fit, 0.001, &Xo, &Yo);
+  coords[0].crpix1 = Xo;
+  coords[0].crpix2 = Yo;
+  
+  // resulting fit should have zero constant terms
+  modfit = CoordsSetCenter (fit, Xo, Yo);
+
+  // set pc1_1, pc1_2, pc2_1, pc2_2 (cd1,cd2 = 1.0)
+  coords[0].pc1_1 = modfit[0].xfit[1][0];
+  coords[0].pc1_2 = modfit[0].xfit[0][1];
+  coords[0].pc2_1 = modfit[0].yfit[1][0];
+  coords[0].pc2_2 = modfit[0].yfit[0][1];
+
+  // set the polyterm elements 
+  coords[0].polyterms[0][0] = modfit[0].xfit[2][0];
+  coords[0].polyterms[1][0] = modfit[0].xfit[1][1];
+  coords[0].polyterms[2][0] = modfit[0].xfit[0][2];
+
+  coords[0].polyterms[0][1] = modfit[0].yfit[2][0];
+  coords[0].polyterms[1][1] = modfit[0].yfit[1][1];
+  coords[0].polyterms[2][1] = modfit[0].yfit[0][2];
+
+  // I need to validate Norder
+  coords[0].polyterms[3][0] = modfit[0].xfit[3][0];
+  coords[0].polyterms[4][0] = modfit[0].xfit[2][1];
+  coords[0].polyterms[5][0] = modfit[0].xfit[1][2];
+  coords[0].polyterms[6][0] = modfit[0].xfit[0][3];
+
+  coords[0].polyterms[3][1] = modfit[0].yfit[3][0];
+  coords[0].polyterms[4][1] = modfit[0].yfit[2][1];
+  coords[0].polyterms[5][1] = modfit[0].yfit[1][2];
+  coords[0].polyterms[6][1] = modfit[0].yfit[0][3];
+
+  /* we do not modify crval1,2: these are kept at the default values */
+
+  /* normalize pc11,etc */
+  c11 = coords[0].pc1_1;
+  c21 = coords[0].pc1_1;
+  c12 = coords[0].pc1_1;
+  c22 = coords[0].pc1_1;
+  coords[0].cdelt1 = coords[0].cdelt2 = sqrt(fabs(c11*c22 - c12*c21));
+  R = 1 / coords[0].cdelt1;
+
+  coords[0].pc1_1  = c11*R;
+  coords[0].pc2_1  = c21*R;
+  coords[0].pc1_2  = c12*R;
+  coords[0].pc2_2  = c22*R;
+
+  fit_free (modfit);
+  /* keep the order and type from initial values */
+}
+
Index: trunk/Ohana/src/relastro/src/global_stats.c
===================================================================
--- trunk/Ohana/src/relastro/src/global_stats.c	(revision 12067)
+++ trunk/Ohana/src/relastro/src/global_stats.c	(revision 12332)
@@ -1,3 +1,3 @@
-# include "relphot.h"
+# include "relastro.h"
 
 void global_stats (Catalog *catalog, int Ncatalog) {
Index: trunk/Ohana/src/relastro/src/initialize.c
===================================================================
--- trunk/Ohana/src/relastro/src/initialize.c	(revision 12067)
+++ trunk/Ohana/src/relastro/src/initialize.c	(revision 12332)
@@ -1,24 +1,35 @@
-# include "relphot.h"
+# include "relastro.h"
 
 void initialize (int argc, char **argv) {
 
-  int N;
+  int NPHOTCODES;
+  char *codename, *ptr;
 
   ConfigInit (&argc, argv);
   args (argc, argv);
 
-  N = UserPatchSelect ? 1 : 2;
-  if ((photcode = GetPhotcodebyName (argv[N])) == NULL) {
-    fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", argv[N]);
-    exit (1);
+  /* XXX we need to build a list of accepted photcodes. these will be used by bcatalog to accept or
+     reject loaded data */
+
+  Nphotcodes = 0;
+  photcodes = NULL;
+  if (PHOTCODE_LIST != NULL) {
+    NPHOTCODES = 10;
+    ALLOCATE (photcodes, PhotCode *, NPHOTCODES);
+
+    /* parse the comma-separated list of photcodes */
+    strtok_r (PHOTCODE_LIST, ",", &ptr);
+    while ((codename = strtok_r (NULL, ",", &ptr)) != NULL) {
+      fprintf (stderr, "PHOTCODE_LIST: %s\n", PHOTCODE_LIST);
+      fprintf (stderr, "codename: %s\n", codename);
+      if ((photcodes[Nphotcodes] = GetPhotcodebyName (codename)) == NULL) {
+	fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename);
+	exit (1);
+      }
+      Nphotcodes ++;
+      CHECK_REALLOCATE (photcodes, PhotCode *, NPHOTCODES, Nphotcodes, 10);
+    }
   }
-  if ((photcode[0].type != PHOT_PRI) && (photcode[0].type != PHOT_SEC)) {
-    fprintf (stderr, "photcode %s is not a primary or secondary filter\n", argv[N]);
-    exit (1);
-  }
-
-  PhotPrimary = (photcode[0].type == PHOT_PRI);
   PhotNsec = GetPhotcodeNsecfilt ();
-  PhotSec = GetPhotcodeNsec (photcode[0].code);
 
   initstats (STATMODE);
@@ -28,4 +39,5 @@
   MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR | ID_MEAS_SKIP | ID_MEAS_AREA;
 
+  /* XXX drop irrelevant entries */
   if (SHOW_PARAMS) {
     fprintf (stderr, "current parameter settings:\n");
Index: trunk/Ohana/src/relastro/src/liststats.c
===================================================================
--- trunk/Ohana/src/relastro/src/liststats.c	(revision 12067)
+++ trunk/Ohana/src/relastro/src/liststats.c	(revision 12332)
@@ -1,3 +1,3 @@
-# include "relphot.h"
+# include "relastro.h"
 
 enum {M_MEAN, M_MEDIAN, M_WT_MEAN, M_INNER_MEAN, 
Index: trunk/Ohana/src/relastro/src/load_catalogs.c
===================================================================
--- trunk/Ohana/src/relastro/src/load_catalogs.c	(revision 12067)
+++ trunk/Ohana/src/relastro/src/load_catalogs.c	(revision 12332)
@@ -1,8 +1,8 @@
-# include "relphot.h"
+# include "relastro.h"
 
-Catalog *load_catalogs (SkyList *skylist, int *Ncatalog) {
+Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int subselect) {
 
   int i, Nstar;
-  Catalog *catalog, tcatalog;
+  Catalog *catalog, *pcatalog, tcatalog;
 
   if (VERBOSE) fprintf (stderr, "loading catalog data\n");
@@ -13,20 +13,24 @@
   for (i = 0; i < skylist[0].Nregions; i++) {
 
+    pcatalog = subselect ? &tcatalog : &catalog[i];
+
     // set up the basic catalog info
-    tcatalog.filename = skylist[0].filename[i];
-    tcatalog.catformat = dvo_catalog_catformat (CATFORMAT);    // set the default catformat from config data
-    tcatalog.catmode   = dvo_catalog_catmode (CATMODE);        // set the default catmode from config data
-    tcatalog.catflags  = LOAD_AVES | LOAD_MEAS | LOAD_SECF; // don't need to load all data at this point
+    pcatalog[0].filename  = skylist[0].filename[i];
+    pcatalog[0].catformat = dvo_catalog_catformat (CATFORMAT);    // set the default catformat from config data
+    pcatalog[0].catmode   = dvo_catalog_catmode (CATMODE);        // set the default catmode from config data
+    pcatalog[0].catflags  = LOAD_AVES | LOAD_MEAS | LOAD_SECF; // don't need to load all data at this point
 
-    if (!dvo_catalog_open (&tcatalog, skylist[0].regions[i], VERBOSE, "r")) {
-      fprintf (stderr, "ERROR: failure reading catalog %s\n", tcatalog.filename);
+    if (!dvo_catalog_open (pcatalog, skylist[0].regions[i], VERBOSE, "w")) {
+      fprintf (stderr, "ERROR: failure reading catalog %s\n", pcatalog[0].filename);
       exit (1);
     }
-    if (VERBOSE && !tcatalog.Nave_disk) fprintf (stderr, "no data in %s, skipping\n", tcatalog.filename);
+    if (VERBOSE && !pcatalog[0].Nave_disk) fprintf (stderr, "no data in %s, skipping\n", pcatalog[0].filename);
 
     // select only the brighter stars
-    bcatalog (&catalog[i], &tcatalog);
-    dvo_catalog_unlock (&tcatalog);
-    dvo_catalog_free (&tcatalog);
+    if (subselect) {
+      bcatalog (&catalog[i], &tcatalog);
+      dvo_catalog_unlock (&tcatalog);
+      dvo_catalog_free (&tcatalog);
+    }
   }
 
@@ -37,6 +41,5 @@
   }
   if (Nstar < 2) { 
-    fprintf (stderr, "insufficient stars %d\n", Nstar);
-    exit (0);
+    fprintf (stderr, "warning: insufficient stars %d\n", Nstar);
   }
 
Index: trunk/Ohana/src/relastro/src/load_images.c
===================================================================
--- trunk/Ohana/src/relastro/src/load_images.c	(revision 12067)
+++ trunk/Ohana/src/relastro/src/load_images.c	(revision 12332)
@@ -1,8 +1,8 @@
-# include "relphot.h"
+# include "relastro.h"
 
-SkyList *load_images (FITS_DB *db, char *regionName, SkyRegion *region, int RegionSelect) {
+SkyList *load_images (FITS_DB *db, SkyRegion *region) {
 
   Image     *image, *subset;
-  int        Nimage, Nsubset, Nchar;
+  int        Nimage, Nsubset;
   int       *LineNumber;
 
@@ -15,11 +15,5 @@
   
   // determine the populated SkyRegions overlapping the requested area
-  if (RegionSelect) {
-    skylist = SkyListByPatch (sky, -1, region);
-  } else {
-    Nchar = strlen(regionName);
-    if (!strcmp (&regionName[Nchar-4], ".cpt")) regionName[Nchar-4] = 0;
-    skylist = SkyListByName (sky, regionName);
-  }
+  skylist = SkyListByPatch (sky, -1, region);
 
   // convert database table to internal structure
@@ -34,4 +28,7 @@
   initMosaics (subset, Nsubset);
   
+  /* unlock, if we can (else, unlocked below) */
+  if (!UPDATE) dvo_image_unlock (db); 
+
   return (skylist);
 }
Index: trunk/Ohana/src/relastro/src/mkpolyterm.c
===================================================================
--- trunk/Ohana/src/relastro/src/mkpolyterm.c	(revision 12332)
+++ trunk/Ohana/src/relastro/src/mkpolyterm.c	(revision 12332)
@@ -0,0 +1,159 @@
+# include "relastro.h"
+
+double factorial (int N) {
+
+  int i;
+  double F;
+
+  F = 1;
+  for (i = N; i > 1; i--) {
+    F *= i;
+  }
+  return F;
+}
+
+// XXX use a separate structure for the map (x2,y2 = f,g(x1,y1)) ?
+
+// given a 2D transformation -- L(x,y),M(x,y) -- find the coordinates x,y
+// for which L,M = 0,0. tol is the allowed error on x,y.
+int CoordsGetCenter (CoordFit *fit, double tol, double *xo, double *yo) {
+
+  int i, Nx, Ny;
+  double R, Xo, Yo, dPos;
+  double **XdX, **XdY, **YdX, **YdY, **alpha, **beta;
+  double **xfit, **yfit;
+
+  xfit = fit[0].xfit;
+  yfit = fit[0].yfit;
+
+  // solve for Xo,Yo where L,M(Xo,Yo) = 0,0
+  // start with linear solution for Xo,Yo:
+  R  = (xfit[1][0]*yfit[0][1] - xfit[0][1]*yfit[1][0]);
+  Xo = (yfit[0][0]*xfit[0][1] - xfit[0][0]*yfit[0][1])/R;
+  Yo = (xfit[0][0]*yfit[1][0] - yfit[0][0]*xfit[1][0])/R;
+
+  Nx = fit[0].Nterms;
+  Ny = fit[0].Nterms;
+
+  // iterate to actual solution: requires small non-linear terms
+  if (fit[0].Norder > 1) {
+    XdX = poly2d_dx (xfit, Nx, Ny);
+    XdY = poly2d_dy (xfit, Nx, Ny);
+    YdX = poly2d_dx (yfit, Nx, Ny);
+    YdY = poly2d_dy (yfit, Nx, Ny);
+
+    alpha = array_init (2, 2);
+    beta  = array_init (2, 1);
+
+    /* this loop uses the Newton-Raphson method to solve for Xo,Yo
+     * it needs the high order terms to be small 
+     * Xo,Yo are in pixels;
+     */
+    dPos = tol + 1;
+    for (i = 0; (dPos > tol) && (i < 20); i++) {
+      // NOTE: order for alpha is: [y][x]
+      alpha[0][0] = poly2d_eval (XdX, Nx-1, Ny,   Xo, Yo);
+      alpha[1][0] = poly2d_eval (XdY, Nx,   Ny-1, Xo, Yo);
+      alpha[0][1] = poly2d_eval (YdX, Nx-1, Ny,   Xo, Yo);
+      alpha[1][1] = poly2d_eval (YdY, Nx,   Ny-1, Xo, Yo);
+
+      beta[0][0] = poly2d_eval (xfit, Nx, Ny, Xo, Yo);
+      beta[1][0] = poly2d_eval (yfit, Nx, Ny, Xo, Yo);
+
+      dgaussj (alpha, 2, beta, 1);
+
+      Xo -= beta[0][0];
+      Yo -= beta[1][0];
+      dPos = hypot(beta[0][0], beta[1][0]);
+    }
+    array_free (alpha, 2);
+    array_free (beta, 2);
+    array_free (XdX, Nx - 1);
+    array_free (XdY, Nx);
+    array_free (YdX, Nx - 1);
+    array_free (YdY, Nx);
+  }
+  *xo = Xo;
+  *yo = Yo;
+  return TRUE;
+}
+
+// convert a transformation L(x,y) to L'(x-xo,y-yo)
+CoordFit *CoordsSetCenter (CoordFit *input, double Xo, double Yo) {
+
+  int i, j, Nx, Ny;
+  double **xPx, **yPx, **xPy, **yPy, **tmp;
+
+  CoordFit *output;
+  output = fit_init (input->Norder);
+
+  /* given two equivalent polynomial representations L(x,y) = \sum_i \sum_j A_{i,j} x^i y^j
+   * we can transform L(x,y) into L'(x-xo,y-yo) by taking the derivatives of both sides and 
+   * noting that the constant term in each is the coefficient in the case of L(x,y) and is the 
+   * value of L'(-xo,-yo) in the second case.
+   */
+
+  Nx = input->Nterms;
+  Ny = input->Nterms;
+
+  xPx = poly2d_copy (input->xfit, Nx, Ny);
+  yPx = poly2d_copy (input->yfit, Nx, Ny);
+
+  for (i = 0; i <= input->Nterms; i++) {
+    xPy = poly2d_copy (xPx, Nx, Ny);
+    yPy = poly2d_copy (yPx, Nx, Ny);
+    for (j = 0; j <= input->Nterms; j++) {
+      output->xfit[i][j] = poly2d_eval (xPy, Nx, Ny, Xo, Yo) / factorial(i) / factorial(j);
+      output->yfit[i][j] = poly2d_eval (yPy, Nx, Ny, Xo, Yo) / factorial(i) / factorial(j);
+
+      // take the next derivative wrt y, catch output (is NULL on last pass)
+      if (Ny > 0) {
+	tmp = poly2d_dy(xPy, Nx, Ny);
+	array_free (xPy, Nx);
+	xPy = tmp;
+
+	tmp = poly2d_dy(yPy, Nx, Ny);
+	array_free (yPy, Nx);
+	yPy = tmp;
+      } else {
+	array_free (xPy, Nx);
+	array_free (yPy, Nx);
+      }
+
+      Ny --;
+    }
+    Ny = input->Nterms;
+    // take the next derivative wrt x, catch output (is NULL on last pass)
+    if (Nx > 0) {
+      tmp = poly2d_dx(xPx, Nx, Ny);
+      array_free (xPx, Nx);
+      xPx = tmp;
+      tmp = poly2d_dx(yPx, Nx, Ny);
+      array_free (yPx, Nx);
+      yPx = tmp;
+      Nx --;
+    } else {
+      array_free (xPx, Nx);
+      array_free (yPx, Nx);
+    }
+  }
+  return output;
+}
+
+/*
+  Coords uses a rigid sequence for the coefficients:
+
+  x^0 y^0 : crpix1,2 (but note this is applied before higher order terms)
+
+  x^1 y^0 : pc1_1, pc2_1
+  x^0 y^1 : pc1_2, pc2_2
+
+  x^2 y^0 : polyterm[0][0,1]
+  x^1 y^1 : polyterm[1][0,1]
+  x^0 y^2 : polyterm[2][0,1]
+
+  x^3 y^0 : polyterm[3][0,1]
+  x^2 y^1 : polyterm[4][0,1]
+  x^1 y^2 : polyterm[5][0,1]
+  x^0 y^3 : polyterm[6][0,1]
+*/
Index: trunk/Ohana/src/relastro/src/plot_scatter.c
===================================================================
--- trunk/Ohana/src/relastro/src/plot_scatter.c	(revision 12067)
+++ trunk/Ohana/src/relastro/src/plot_scatter.c	(revision 12332)
@@ -1,6 +1,7 @@
-# include "relphot.h"
+# include "relastro.h"
    
 void plot_scatter (Catalog *catalog, int Ncatalog) {
 
+# if (0)   
   int i, j, k, m, N, Ntot;
   float Mrel, Mcal, Mmos, Mgrid;
@@ -28,9 +29,6 @@
       for (k = 0; k < catalog[i].average[j].Nm; k++, m++) {
 	if (catalog[i].measure[m].flags & MEAS_BAD) continue;
-	if ((Mcal  = getMcal  (m, i)) == NO_MAG) continue;
-	if ((Mmos  = getMmos  (m, i)) == NO_MAG) continue;
-	if ((Mgrid = getMgrid (m, i)) == NO_MAG) continue;
 
-	Mrel = PhotPrimary ? catalog[i].average[j].M : catalog[i].secfilt[PhotNsec*j+PhotSec].M_PS;
+	Mrel = catalog[i].secfilt[PhotNsec*j+PhotSec].M;
 	xlist[N] = Mrel;
 	ylist[N] = PhotSys  (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*PhotNsec]) - Mcal - Mmos - Mgrid - Mrel;
@@ -55,3 +53,7 @@
   free (ylist);
   free (ilist);
+
+# endif
 }
+
+/* XXX this should become astrometrically relevant */
Index: trunk/Ohana/src/relastro/src/plotstuff.c
===================================================================
--- trunk/Ohana/src/relastro/src/plotstuff.c	(revision 12067)
+++ trunk/Ohana/src/relastro/src/plotstuff.c	(revision 12332)
@@ -1,3 +1,3 @@
-# include "relphot.h"
+# include "relastro.h"
 # include <signal.h>
 
@@ -49,5 +49,5 @@
   if (Xgraph[N] == 0) return;
 
-  KiiPS (Xgraph[N], TRUE, KAPA_PS_NEWPLOT, filename);
+  KiiPS (Xgraph[N], TRUE, KAPA_PS_NEWPLOT, filename, "default");
   return;
 }
Index: trunk/Ohana/src/relastro/src/relastro.c
===================================================================
--- trunk/Ohana/src/relastro/src/relastro.c	(revision 12067)
+++ trunk/Ohana/src/relastro/src/relastro.c	(revision 12332)
@@ -3,5 +3,5 @@
 int main (int argc, char **argv) {
 
-  int i, status, Ncatalog;
+  int status, Ncatalog;
   Catalog *catalog;
   FITS_DB db;
@@ -22,38 +22,40 @@
 
   /* load regions and images based on specified sky patch */
-  // XXX need to mimic old-style load by passing patch name
-  // XXX need to reduce number of global variables in use.
-  // XXX this is fairly lame: argv[1] is photcode if UserPatchSelect is true
-  skylist = load_images (&db, argv[1], &UserPatch, UserPatchSelect);
+  skylist = load_images (&db, &UserPatch);
 
-  /* unlock, if we can (else, unlocked below) */
-  if (!UPDATE) dvo_image_unlock (&db); 
-
-  /* load catalog data from region files */
-  catalog = load_catalogs (skylist, &Ncatalog);
+  /* load catalog data from region files : subselect only if we are not doing the objects */
+  catalog = load_catalogs (skylist, &Ncatalog, (FIT_TARGET != TARGET_OBJECTS));
 
   /* match measurements with images, mosaics */
   initImageBins  (catalog, Ncatalog);
-  initMosaicBins (catalog, Ncatalog);
-  initGridBins   (catalog, Ncatalog);
-  initMrel (catalog, Ncatalog);
+  // initMosaicBins (catalog, Ncatalog);
 
   findImages (catalog, Ncatalog);
-  findMosaics (catalog, Ncatalog);  /* also sets Grid values */
-  SAVEPLOT = FALSE;
-
-  setExclusions (catalog, Ncatalog);
+  // findMosaics (catalog, Ncatalog);  /* also sets Grid values */
 
   if (PLOTSTUFF) {
-    plot_star_coords (catalog, Ncatalog);
+    // plot_star_coords (catalog, Ncatalog);
     plot_mosaic_fields (catalog);
   }
 
   /* major modes */
-  if (UpdateObjects) {
-    UpdateObjects (catalog, Ncatalog);
-  }
-  if (UpdateImages) {
-    UpdateImages (catalog, Ncatalog);
+  switch (FIT_TARGET) {
+    case TARGET_OBJECTS:
+      UpdateObjects (catalog, Ncatalog);
+      break;
+    case TARGET_SIMPLE:
+      UpdateSimple (catalog, Ncatalog);
+      break;
+
+    case TARGET_CHIPS:
+      UpdateChips (catalog, Ncatalog);
+      break;
+
+    case TARGET_MOSAICS:
+      UpdateMosaic (catalog, Ncatalog);
+      break;
+    default:
+      fprintf (stderr, "programming error at %s:%s", __FILE__, __LINE__);
+      exit (2);
   }
 
@@ -69,7 +71,10 @@
   if (!UPDATE) exit (0);
 
-  /* need to figure out how to update images, etc */
-  dvo_image_update (&db, VERBOSE);
-  dvo_image_unlock (&db); 
+  if (FIT_TARGET == TARGET_OBJECTS) {
+    save_catalogs (catalog, Ncatalog);
+  } else {
+    dvo_image_update (&db, VERBOSE);
+    dvo_image_unlock (&db); 
+  }
 
   exit (0);
Index: trunk/Ohana/src/relastro/src/reload_catalogs.c
===================================================================
--- trunk/Ohana/src/relastro/src/reload_catalogs.c	(revision 12067)
+++ trunk/Ohana/src/relastro/src/reload_catalogs.c	(revision 12332)
@@ -1,3 +1,3 @@
-# include "relphot.h"
+# include "relastro.h"
 
 void reload_catalogs (SkyList *skylist) {
@@ -38,10 +38,8 @@
     initImageBins  (&catalog, 1);
     initMosaicBins (&catalog, 1);
-    initGridBins   (&catalog, 1);
 
     findImages (&catalog, 1);
     findMosaics (&catalog, 1);
 
-    setMrelFinal (&catalog);
     dvo_catalog_save (&catalog, VERBOSE); 
     dvo_catalog_unlock (&catalog);
@@ -50,5 +48,4 @@
     freeImageBins (1);
     freeMosaicBins (1);
-    freeGridBins (1);
   }
 }
Index: trunk/Ohana/src/relastro/src/save_catalogs.c
===================================================================
--- trunk/Ohana/src/relastro/src/save_catalogs.c	(revision 12332)
+++ trunk/Ohana/src/relastro/src/save_catalogs.c	(revision 12332)
@@ -0,0 +1,14 @@
+# include "relastro.h"
+
+void save_catalogs (Catalog *catalog, int Ncatalog) {
+
+  int i;
+
+  /* load data from each region file */
+  for (i = 0; i < Ncatalog; i++) {
+
+    dvo_catalog_save (&catalog[i], VERBOSE); 
+    dvo_catalog_unlock (&catalog[i]);
+    dvo_catalog_free (&catalog[i]);
+  }
+}
Index: trunk/Ohana/src/relastro/src/select_images.c
===================================================================
--- trunk/Ohana/src/relastro/src/select_images.c	(revision 12067)
+++ trunk/Ohana/src/relastro/src/select_images.c	(revision 12332)
@@ -1,3 +1,3 @@
-# include "relphot.h"
+# include "relastro.h"
 
 /* this function returns a list of all images which overlap the given SkyList (set of
@@ -17,5 +17,5 @@
   Image *image;
   int i, j, k, m, found, nimage, NIMAGE;
-  int InRange, ecode;
+  int InRange;
   double Ri[5], Di[5], Xi[5], Yi[5], dx, dy;
   int *line_number;
@@ -75,7 +75,12 @@
   for (i = 0; i < Ntimage; i++) {
       
-    /* exclude images by photcode */
-    ecode = GetPhotcodeEquivCodebyCode (timage[i].source);
-    if (ecode != photcode[0].code) continue;
+    /* exclude images by photcode, if selected */
+    if (Nphotcodes > 0) {
+      found = FALSE;
+      for (j = 0; (j < Nphotcodes) && !found; j++) {
+	if (photcodes[j][0].code == timage[i].photcode) found = TRUE;
+      }
+      if (!found) continue;
+    }
 
     /* exclude images by time */
@@ -135,5 +140,5 @@
       if (RESET) {
 	assignMcal (&image[nimage], (double *) NULL, -1);
-	image[nimage].dMcal_PS = NO_MAG;
+	image[nimage].dMcal = NO_MAG;
 	image[nimage].code &= ~ID_IMAGE_POOR;
       }
Index: trunk/Ohana/src/relastro/src/setExclusions.c
===================================================================
--- trunk/Ohana/src/relastro/src/setExclusions.c	(revision 12067)
+++ trunk/Ohana/src/relastro/src/setExclusions.c	(revision 12332)
@@ -1,7 +1,8 @@
-# include "relphot.h"
+# include "relastro.h"
 
+/* XXX I think all of these, except for X,Y selection, are done / can be done in bcatalog */
 int setExclusions (Catalog *catalog, int Ncatalog) {
 
-  int i, j, k, m, Narea, Nnocal, ecode;
+  int i, j, k, m, Narea, Nnocal;
   Coords *coords;
   double r, d, x, y;
@@ -14,6 +15,7 @@
 
 	/* select measurements by photcode */
-	ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[m].source);
-	if (ecode != photcode[0].code) goto mark_nocal;
+	// XXXX is this done when loaded, or is it needed for reload?
+	// ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[m].photcode);
+	// if (ecode != photcode[0].code) goto mark_nocal;
 	
 	/* select measurements by time */
@@ -25,6 +27,6 @@
 	/* select measurements by mag limit */
 	if (AreaSelect) {
-	  r = catalog[i].average[j].R + catalog[i].measure[m].dR_PS / 3600.0;
-	  d = catalog[i].average[j].D + catalog[i].measure[m].dD_PS / 3600.0;
+	  r = catalog[i].average[j].R + catalog[i].measure[m].dR / 3600.0;
+	  d = catalog[i].average[j].D + catalog[i].measure[m].dD / 3600.0;
 	  if ((coords = getCoords (m, i)) == NULL) goto markbad;
 	  RD_to_XY (&x, &y, r, d, coords);
Index: trunk/Ohana/src/relastro/src/setMrelFinal.c
===================================================================
--- trunk/Ohana/src/relastro/src/setMrelFinal.c	(revision 12067)
+++ trunk/Ohana/src/relastro/src/setMrelFinal.c	(revision 12332)
@@ -1,3 +1,3 @@
-# include "relphot.h"
+# include "relastro.h"
 
 void setMrelFinal (Catalog *catalog) {
@@ -11,7 +11,7 @@
 
     for (i = 0; i < catalog[0].Naverage; i++) {
-      p = (PhotPrimary) ? &catalog[0].average[i].M  : &catalog[0].secfilt[PhotNsec*i+PhotSec].M_PS;
+      p = (PhotPrimary) ? &catalog[0].average[i].M  : &catalog[0].secfilt[PhotNsec*i+PhotSec].M;
       *p = NO_MAG;
-      p = (PhotPrimary) ? &catalog[0].average[i].dM : &catalog[0].secfilt[PhotNsec*i+PhotSec].dM_PS;
+      p = (PhotPrimary) ? &catalog[0].average[i].dM : &catalog[0].secfilt[PhotNsec*i+PhotSec].dM;
       *p = NO_MAG;
       q = (PhotPrimary) ? &catalog[0].average[i].Xm : &catalog[0].secfilt[PhotNsec*i+PhotSec].Xm;
@@ -22,5 +22,5 @@
 	
 	/* select measurements by photcode */
-	ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].source);
+	ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].photcode);
 	if (ecode != photcode[0].code) continue;
 	
@@ -31,5 +31,5 @@
 	}
 	
-	catalog[0].measure[m].Mcal_PS = 0;
+	catalog[0].measure[m].Mcal = 0;
 	catalog[0].measure[m].flags &= 0xff00;
 	catalog[0].measure[m].flags &= ~ID_MEAS_POOR;
@@ -97,5 +97,5 @@
       /** never use these measurements (wrong photcode, bad time range) */
       /* skipped via NOCAL, don't mark as skipped */
-      ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].source);
+      ecode = GetPhotcodeEquivCodebyCode (catalog[0].measure[m].photcode);
       if (ecode != photcode[0].code) continue;
 
Index: trunk/Ohana/src/relastro/src/sort.c
===================================================================
--- trunk/Ohana/src/relastro/src/sort.c	(revision 12067)
+++ trunk/Ohana/src/relastro/src/sort.c	(revision 12332)
@@ -1,3 +1,3 @@
-# include "relphot.h"
+# include "relastro.h"
 
 void sortA (double *X, int N) {
Index: trunk/Ohana/src/relastro/src/write_coords.c
===================================================================
--- trunk/Ohana/src/relastro/src/write_coords.c	(revision 12067)
+++ trunk/Ohana/src/relastro/src/write_coords.c	(revision 12332)
@@ -1,3 +1,3 @@
-# include "relphot.h"
+# include "relastro.h"
 
 # define CD_COORDS 1