Index: /trunk/Ohana/src/mosastro/Makefile
===================================================================
--- /trunk/Ohana/src/mosastro/Makefile	(revision 3322)
+++ /trunk/Ohana/src/mosastro/Makefile	(revision 3323)
@@ -1,7 +1,9 @@
+default: mosastro
+help:
+	@echo "make options: mosastro (default) mkstandards mkobs warptest"
+
 include ../../Configure
+
 HOME    =       $(ROOT)/src/mosastro.v3
-PROGRAM =       mosastro
-default: $(PROGRAM)
-
 BIN	=	$(HOME)/bin
 INC	= 	$(HOME)/include
@@ -32,4 +34,5 @@
 $(SRC)/sort.$(ARCH).o \
 $(SRC)/ConfigInit.$(ARCH).o \
+$(SRC)/getstone.$(ARCH).o \
 $(SRC)/get2mass.$(ARCH).o \
 $(SRC)/getgsc.$(ARCH).o \
@@ -42,4 +45,15 @@
 $(SRC)/field.$(ARCH).o \
 $(SRC)/chips.$(ARCH).o \
+$(SRC)/GetGradients.$(ARCH).o \
+$(SRC)/FitGradients.$(ARCH).o \
+$(SRC)/FitChips.$(ARCH).o \
+$(SRC)/fitpoly.$(ARCH).o \
+$(SRC)/output.$(ARCH).o \
+$(SRC)/wstars.$(ARCH).o \
+$(SRC)/gaussj.$(ARCH).o \
+$(SRC)/mkpolyterm.$(ARCH).o \
+$(SRC)/mkmosaic.$(ARCH).o \
+$(SRC)/GetScatter.$(ARCH).o \
+$(SRC)/clip.$(ARCH).o \
 $(SRC)/parse_time.$(ARCH).o
 
@@ -50,4 +64,5 @@
 $(SRC)/mkstandards.$(ARCH).o \
 $(SRC)/mkheader.$(ARCH).o \
+$(SRC)/random.$(ARCH).o \
 $(SRC)/wstars.$(ARCH).o
 
@@ -56,4 +71,5 @@
 
 OBS = \
+$(SRC)/args_obs.$(ARCH).o \
 $(SRC)/mkobs.$(ARCH).o \
 $(SRC)/mkheader.$(ARCH).o \
@@ -67,8 +83,11 @@
 $(SRC)/getgsc.$(ARCH).o \
 $(SRC)/getptolemy.$(ARCH).o \
+$(SRC)/getstone.$(ARCH).o \
 $(SRC)/gptolemy.$(ARCH).o \
 $(SRC)/gregions.$(ARCH).o \
 $(SRC)/gcatalog.$(ARCH).o \
 $(SRC)/getusno.$(ARCH).o \
+$(SRC)/mkpolyterm.$(ARCH).o \
+$(SRC)/random.$(ARCH).o \
 $(SRC)/wstars.$(ARCH).o
 
@@ -78,4 +97,10 @@
 warptest: $(BIN)/warptest.(ARCH)
 $(BIN)/warptest.(ARCH) : $(SRC)/warptest.$(ARCH).o
+
+$(MOS) : $(INC)/mosastro.h
+$(OBS) : $(INC)/mosastro.h
+$(STD) : $(INC)/mosastro.h
+
+INSTALL = mosastro mkstandards mkobs
 
 # dependancy rules for binary code #########################
Index: /trunk/Ohana/src/mosastro/include/mosastro.h
===================================================================
--- /trunk/Ohana/src/mosastro/include/mosastro.h	(revision 3322)
+++ /trunk/Ohana/src/mosastro/include/mosastro.h	(revision 3323)
@@ -8,5 +8,5 @@
   double X, Y;  /* Chip Coords   - pixels  */
   double Mag, dMag;
-  int type;
+  int mask;
 } StarData;
 
@@ -20,5 +20,9 @@
   char *file;
   int Nstars;
+  int NX, NY;
   StarData *stars;
+  Header header;
+  char *buffer;
+  int Nbuffer;
 
   int Nmatch;
@@ -60,13 +64,18 @@
 } Stars;
 
+int  ChipOrder;
 int  Nchip;
 Chip *chip;
 Field field;
 double Year;  /** carried for precession - probably put this in chip data **/
-double RADIUS; /** raw / ref matching radius (units?) **/
+double RADIUS; /** raw / ref matching radius (pixels on Focal Plane) **/
+double SIGMA;
+double SIGMA_LIM;
 
+char StoneRegions[256];
 char CDROM[256];
 char REFCAT[256];
 char CATDIR[256];
+char ExptimeKeyword[256];
 char DateKeyword[256];
 char DateMode[256];
@@ -76,7 +85,14 @@
 char TWO_MASS_DIR[256];
 int VERBOSE;
+int NO_CHIPS;
 
 char GSCFILE[256];
 char GSC_DIR[256];
+
+char *FIELD;
+char *CHIPS;
+char *OUTPUT;
+char *DUMP;
+char *FOCAL_PLANE;
 
 StarData *getusno (CatStats *catstats, int *Nstars);
@@ -103,2 +119,7 @@
 int fake_field (double RA, double DEC);
 Gradients *GetGradients ();
+void FitChips (int Norder);
+void FitChipLinear (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
+void FitChipResid (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
+double GetScatter ();
+Header *mkmosaic (int Nx, int Ny, int Nstars, Coords *coords);
Index: /trunk/Ohana/src/mosastro/src/ConfigInit.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/ConfigInit.c	(revision 3322)
+++ /trunk/Ohana/src/mosastro/src/ConfigInit.c	(revision 3323)
@@ -17,18 +17,28 @@
   if (VERBOSE) fprintf (stderr, "loaded config file: %s\n", file);
 
-  GetConfig (config, "USNO_CDROM",        "%s", 0, CDROM);
-  GetConfig (config, "DATE-KEYWORD",      "%s", 0, DateKeyword);
-  GetConfig (config, "DATE-MODE",         "%s", 0, DateMode);
-  GetConfig (config, "UT-KEYWORD",        "%s", 0, UTKeyword);
-  GetConfig (config, "MJD-KEYWORD",       "%s", 0, MJDKeyword);
-  GetConfig (config, "JD-KEYWORD",        "%s", 0, JDKeyword);
-  ScanConfig (config, "2MASS_DIR",         "%s",  0, TWO_MASS_DIR);
+  GetConfig (config, "USNO_CDROM",        "%s",  0, CDROM);
+  GetConfig (config, "DATE-KEYWORD",      "%s",  0, DateKeyword);
+  GetConfig (config, "DATE-MODE",         "%s",  0, DateMode);
+  GetConfig (config, "UT-KEYWORD",        "%s",  0, UTKeyword);
+  GetConfig (config, "MJD-KEYWORD",       "%s",  0, MJDKeyword);
+  GetConfig (config, "JD-KEYWORD",        "%s",  0, JDKeyword);
+  ScanConfig (config, "2MASS_DIR",        "%s",  0, TWO_MASS_DIR);
 
   GetConfig (config, "ASTRO_REFCAT",      "%s",  0, REFCAT);
-  ScanConfig (config, "CATDIR",            "%s",  0, CATDIR);
-  ScanConfig (config, "GSCFILE",           "%s",  0, GSCFILE);
-  ScanConfig (config, "GSCDIR",            "%s",  0, GSC_DIR);
+  ScanConfig (config, "CATDIR",           "%s",  0, CATDIR);
+  ScanConfig (config, "GSCFILE",          "%s",  0, GSCFILE);
+  ScanConfig (config, "GSCDIR",           "%s",  0, GSC_DIR);
+  ScanConfig (config, "STONE_DIR",        "%s",  0, StoneRegions);
+
+  /* abstracted header keywords */
+  ScanConfig (config, "DATE-KEYWORD",     "%s",  0, DateKeyword);
+  ScanConfig (config, "DATE-MODE",        "%s",  0, DateMode);
+  ScanConfig (config, "UT-KEYWORD",       "%s",  0, UTKeyword);
+  ScanConfig (config, "MJD-KEYWORD",      "%s",  0, MJDKeyword);
+  ScanConfig (config, "JD-KEYWORD",       "%s",  0, JDKeyword);
+  ScanConfig (config, "EXPTIME-KEYWORD",  "%s",  0, ExptimeKeyword);
 
   GetConfig (config, "RADIUS",            "%lf", 0, &RADIUS);
+  GetConfig (config, "SIGMA_LIM",         "%lf", 0, &SIGMA_LIM);
 
   free (config);
Index: /trunk/Ohana/src/mosastro/src/FitChips.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/FitChips.c	(revision 3323)
+++ /trunk/Ohana/src/mosastro/src/FitChips.c	(revision 3323)
@@ -0,0 +1,25 @@
+# include "mosastro.h"
+
+void FitChips (int Norder) {
+
+  int i;
+
+  for (i = 0; i < Nchip; i++) {
+    chip[i].map.Npolyterms = Norder;
+    FitChip (chip[i].raw, chip[i].ref, chip[i].Nmatch, &chip[i].map);
+  }
+  deproject_raw ();
+  project_ref ();
+}
+
+void FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords) {
+
+  int i;
+
+  fit_init (coords[0].Npolyterms);
+  for (i = 0; i < Nmatch; i++) {
+    fit_add (raw[i].X, raw[i].Y, ref[i].L, ref[i].M);
+  }
+  fit_eval ();
+  fit_apply_coords (coords);
+}
Index: /trunk/Ohana/src/mosastro/src/FitGradients.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/FitGradients.c	(revision 3323)
+++ /trunk/Ohana/src/mosastro/src/FitGradients.c	(revision 3323)
@@ -0,0 +1,48 @@
+# include "mosastro.h"
+
+/* fit dPdL, etc vs L,M */
+void FitGradients (Gradients *grad) {
+
+  int i, j, Nterm, Norder;
+  double **Pf, **Qf, chisq;
+  Gradients grfix;
+
+  ALLOCATE (grfix.dPdL, double, grad[0].Npts);
+  ALLOCATE (grfix.dPdM, double, grad[0].Npts);
+  ALLOCATE (grfix.dQdL, double, grad[0].Npts);
+  ALLOCATE (grfix.dQdM, double, grad[0].Npts);
+
+  ALLOCATE (grfix.Lo,   double, grad[0].Npts);
+  ALLOCATE (grfix.Mo,   double, grad[0].Npts);
+  grfix.Npts = grad[0].Npts;
+
+  /* where do we set field.distort.Npolyterms? */
+  Norder = 3;
+
+  fit_init (Norder - 1);
+  for (i = 0; i < grad[0].Npts; i++) {
+    fit_add (grad[0].Lo[i], grad[0].Mo[i], grad[0].dPdL[i], grad[0].dPdM[i]);
+  }
+  fit_eval ();
+  fit_apply_grads (&field.distort, &field.project, 0);
+  fit_correct_grads (grad, &grfix, 0);
+  fit_free ();
+
+  fit_init (Norder - 1);
+  for (i = 0; i < grad[0].Npts; i++) {
+    fit_add (grad[0].Lo[i], grad[0].Mo[i], grad[0].dQdL[i], grad[0].dQdM[i]);
+  }
+  fit_eval ();
+  fit_apply_grads (&field.distort, &field.project, 1);
+  fit_correct_grads (grad, &grfix, 1);
+  fit_free ();
+
+  /* use new model of field to get TP & FP coords */
+  deproject_raw ();
+  project_ref ();
+
+  if ((DUMP != NULL) && !strcmp (DUMP, "grads")) {
+    dump_grads (&grfix, "gradfix.dat");
+    exit (0);
+  }
+}
Index: /trunk/Ohana/src/mosastro/src/GetGradients.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/GetGradients.c	(revision 3323)
+++ /trunk/Ohana/src/mosastro/src/GetGradients.c	(revision 3323)
@@ -0,0 +1,108 @@
+# include "mosastro.h"
+
+/* measure local gradient vs field position : dP, dQ = f(L,M) 
+   gradient is measured in small boxes across the field */
+Gradients *GetGradients () {
+
+  int i, j, nx, ny, Nval, Nx, Ny, Npts, DX, DY;
+  double **a, **b;
+  double L, M, dP, dQ;
+  double Xmin, Xmax, Ymin, Ymax;
+
+  Gradients *grad;
+
+  ALLOCATE (grad, Gradients, 1);
+
+  /** this should not be hard-wired **/
+  Nx = 2;
+  Ny = 2;
+
+  ALLOCATE (grad[0].dPdL, double, Nchip*Nx*Ny);
+  ALLOCATE (grad[0].dPdM, double, Nchip*Nx*Ny);
+  ALLOCATE (grad[0].dQdL, double, Nchip*Nx*Ny);
+  ALLOCATE (grad[0].dQdM, double, Nchip*Nx*Ny);
+
+  ALLOCATE (grad[0].Lo, double, Nchip*Nx*Ny);
+  ALLOCATE (grad[0].Mo, double, Nchip*Nx*Ny);
+  
+  ALLOCATE (a, double *, 3);
+  ALLOCATE (b, double *, 3);
+  for (i = 0; i < 3; i++) {
+    ALLOCATE (a[i], double, 3);
+    ALLOCATE (b[i], double, 2);
+  }
+
+  Nval = 0;
+  for (i = 0; i < Nchip; i++) {
+    DX = chip[i].NX / Nx;
+    DY = chip[i].NY / Ny;
+    for (nx = 0; nx < Nx; nx++) {
+      for (ny = 0; ny < Ny; ny++) {
+	Xmin = nx*DX;
+	Xmax = Xmin + DX;
+	Ymin = ny*DY;
+	Ymax = Ymin + DY;
+
+	for (j = 0; j < 3; j++) {
+	  bzero (a[j], 3*sizeof(double));
+	  bzero (b[j], 2*sizeof(double));
+	}
+  
+	/* find local gradient in limited box */
+	for (j = 0; j < chip[i].Nmatch; j++) {
+	  if (chip[i].raw[j].X < Xmin) continue;
+	  if (chip[i].raw[j].X > Xmax) continue;
+	  if (chip[i].raw[j].Y < Ymin) continue;
+	  if (chip[i].raw[j].Y > Ymax) continue;
+	  
+	  L   = chip[i].ref[j].L;
+	  M   = chip[i].ref[j].M;
+	  dP  = chip[i].ref[j].P - chip[i].raw[j].P;
+	  dQ  = chip[i].ref[j].Q - chip[i].raw[j].Q;
+	  
+	  a[0][0] += 1;
+	  a[1][0] += L;
+	  a[2][0] += M;
+	  a[1][1] += L*L;
+	  a[2][1] += L*M;
+	  a[2][2] += M*M;
+
+	  b[0][0] += dP;
+	  b[1][0] += dP*L;
+	  b[2][0] += dP*M;
+
+	  b[0][1] += dQ;
+	  b[1][1] += dQ*L;
+	  b[2][1] += dQ*M;
+
+	}
+	/* fitting dP = a[0][0] + a[1][0]*L + a[2][0]*M
+	           dQ = a[0][1] + a[1][1]*L + a[2][1]*M 
+	*/
+	a[0][1] = a[1][0];
+	a[0][2] = a[2][0];
+	a[1][2] = a[2][1];
+
+	Npts = a[0][0];
+	grad[0].Lo[Nval] = a[1][0] / a[0][0];
+	grad[0].Mo[Nval] = a[2][0] / a[0][0];
+	/* point-weighted average coordinate */
+
+	if (Npts < 5) continue;
+	if (!gaussj (a, 3, b, 2)) continue;
+
+	/* we only care about the slopes, not the offsets */
+	grad[0].dPdL[Nval] = b[1][0];
+	grad[0].dPdM[Nval] = b[2][0];
+	grad[0].dQdL[Nval] = b[1][1];
+	grad[0].dQdM[Nval] = b[2][1];
+
+	Nval ++;
+
+      }
+    }
+  }
+  grad[0].Npts = Nval;
+  if ((DUMP != NULL) && !strcmp (DUMP, "grads")) dump_grads (grad, "grads.dat");
+  return (grad);
+}
Index: /trunk/Ohana/src/mosastro/src/GetScatter.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/GetScatter.c	(revision 3323)
+++ /trunk/Ohana/src/mosastro/src/GetScatter.c	(revision 3323)
@@ -0,0 +1,26 @@
+# include "mosastro.h"
+
+/* measure scatter on the focal plane */
+double GetScatter (int *Nscatter) {
+
+  int i, j, Ntotal;
+  double dL, dM, dR2, dR;
+  StarData *raw, *ref;
+
+  Ntotal = dR2 = 0.0;
+
+  for (i = 0; i < Nchip; i++) {
+    raw = chip[i].raw;
+    ref = chip[i].ref;
+    for (j = 0; j < chip[i].Nmatch; j++) {
+      if (raw[j].mask) continue;
+      dL = raw[j].L - ref[j].L;
+      dM = raw[j].M - ref[j].M;
+      dR2 += dL*dL + dM*dM;
+      Ntotal ++;
+    }
+  }
+  dR = sqrt (dR2 / Ntotal);
+  *Nscatter = Ntotal;
+  return (dR);
+}
Index: /trunk/Ohana/src/mosastro/src/LoadStars.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/LoadStars.c	(revision 3322)
+++ /trunk/Ohana/src/mosastro/src/LoadStars.c	(revision 3323)
@@ -4,8 +4,7 @@
 int LoadStars (int Nfile, char **file) {
 
-  Header header;
+  Header *header;
   int i, j, NCHIP, itmp, Nbytes, nbytes;
   FILE *f;
-  char *buffer;
   time_t tsval;
   struct tm *tmval;
@@ -20,5 +19,5 @@
 
     /* load header */
-    if (!fits_read_header (chip[Nchip].file, &header)) {
+    if (!fits_read_header (chip[Nchip].file, &chip[Nchip].header)) {
       fprintf (stderr, "ERROR: can't read header for %s\n", file[i]);
       free (chip[Nchip].file);
@@ -31,8 +30,12 @@
       continue;
     }
-    fseek (f, header.size, SEEK_SET); 
+    fseek (f, chip[Nchip].header.size, SEEK_SET); 
+
+    /* get image dimensions */
+    fits_scan (&chip[Nchip].header, "NAXIS1",   "%d", 1, &chip[Nchip].NX);
+    fits_scan (&chip[Nchip].header, "NAXIS2",   "%d", 1, &chip[Nchip].NY);
 
     /* get astrometry information */
-    if (!GetCoords (&chip[Nchip].coords, &header)) {
+    if (!GetCoords (&chip[Nchip].coords, &chip[Nchip].header)) {
       fprintf (stderr, "ERROR: no astrometric solution in header\n");
       free (chip[Nchip].file);
@@ -43,5 +46,5 @@
 
     itmp = 0;
-    fits_scan (&header, "NASTRO",   "%d", 1, &itmp);
+    fits_scan (&chip[Nchip].header, "NASTRO",   "%d", 1, &itmp);
     if (itmp == 0) {
       fprintf (stderr, "ERROR: bad astrometric solution in header %s\n", file[i]);
@@ -51,10 +54,10 @@
 
     /* get time info */
-    tsval = parse_time (&header);
+    tsval = parse_time (&chip[Nchip].header);
     tmval = gmtime (&tsval);
     Year = tmval[0].tm_year;
 
     /* find expected number of stars */
-    fits_scan (&header, "NSTARS", "%d", 1, &chip[Nchip].Nstars);
+    fits_scan (&chip[Nchip].header, "NSTARS", "%d", 1, &chip[Nchip].Nstars);
     if (chip[Nchip].Nstars == 0) {
       fprintf (stderr, "ERROR: can't get NSTARS from header for %s\n", file[i]);
@@ -65,17 +68,17 @@
     
     Nbytes = BYTES_STAR * chip[Nchip].Nstars;
-    ALLOCATE (buffer, char, Nbytes + 1);
-    nbytes = Fread (buffer, 1, Nbytes, f, "char");
+    ALLOCATE (chip[Nchip].buffer, char, Nbytes + 1);
+    nbytes = Fread (chip[Nchip].buffer, 1, Nbytes, f, "char");
     if (nbytes != Nbytes) { exit (1); }
+    chip[Nchip].Nbuffer = Nbytes;
 
-    /** need to save all input lines so we can write on output? (gastro does not) **/
     for (j = 0; j < chip[Nchip].Nstars; j++) {
       bzero (&chip[Nchip].stars[j], sizeof(StarData));
-      dparse (&chip[Nchip].stars[j].X,   1, &buffer[j*BYTES_STAR]);
-      dparse (&chip[Nchip].stars[j].Y,   2, &buffer[j*BYTES_STAR]);
-      dparse (&chip[Nchip].stars[j].Mag, 3, &buffer[j*BYTES_STAR]);
+      dparse (&chip[Nchip].stars[j].X,    1, &chip[Nchip].buffer[j*BYTES_STAR]);
+      dparse (&chip[Nchip].stars[j].Y,    2, &chip[Nchip].buffer[j*BYTES_STAR]);
+      dparse (&chip[Nchip].stars[j].Mag,  3, &chip[Nchip].buffer[j*BYTES_STAR]);
+      dparse (&chip[Nchip].stars[j].dMag, 4, &chip[Nchip].buffer[j*BYTES_STAR]);
+      chip[Nchip].stars[j].dMag *= 0.001;  /* millimag errors stored in file */
     }
-    /** free this or keep with chip data? **/
-    free (buffer);
 
     if (VERBOSE) fprintf (stderr, "loaded %d stars from %s\n", chip[Nchip].Nstars, chip[Nchip].file);
Index: /trunk/Ohana/src/mosastro/src/args.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/args.c	(revision 3322)
+++ /trunk/Ohana/src/mosastro/src/args.c	(revision 3323)
@@ -25,14 +25,8 @@
   }
 
-  /*
-  DUMP_REF = FALSE;
-  if ((N = get_argument (*argc, argv, "-dumpref"))) {
-    DUMP_REF = TRUE;
+  DUMP = NULL;
+  if ((N = get_argument (*argc, argv, "-dump"))) {
     remove_argument (N, argc, argv);
-  }
-
-  DUMP_RAW = FALSE;
-  if ((N = get_argument (*argc, argv, "-dumpraw"))) {
-    DUMP_RAW = TRUE;
+    DUMP = strcreate(argv[N]);
     remove_argument (N, argc, argv);
   }
@@ -51,6 +45,6 @@
     remove_argument (N, argc, argv);
   }
-  */
 
+  /** currently unused **/
   field.Norder = 0;
   if ((N = get_argument (*argc, argv, "-order"))) {
@@ -59,3 +53,11 @@
     remove_argument (N, argc, argv);
   }
+
+  /** currently unused **/
+  ChipOrder = 1;
+  if ((N = get_argument (*argc, argv, "-chiporder"))) {
+    remove_argument (N, argc, argv);
+    ChipOrder = atoi (argv[N]);
+    remove_argument (N, argc, argv);
+  }
 }
Index: /trunk/Ohana/src/mosastro/src/args_obs.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/args_obs.c	(revision 3323)
+++ /trunk/Ohana/src/mosastro/src/args_obs.c	(revision 3323)
@@ -0,0 +1,91 @@
+# include "mosastro.h"
+
+void print_help () {
+  fprintf (stderr, "USAGE: mkobs (RA) (DEC) (output) [-p param value]\n");
+  exit (1);
+}
+
+void args (int *argc, char **argv) {
+  
+  int N, No, Np, Nx, Ny;
+  double theta;
+  char line[500];
+
+  if (get_argument (*argc, argv, "--help")) print_help ();
+  if (get_argument (*argc, argv, "-h")) print_help ();
+
+  FOCAL_PLANE = NULL;
+  if ((N = get_argument (*argc, argv, "-fp"))) {
+    remove_argument (N, argc, argv);
+    FOCAL_PLANE = strcreate (argv[N]);
+    remove_argument (N, argc, argv);
+  }
+
+  NO_CHIPS = FALSE;
+  if ((N = get_argument (*argc, argv, "-nochips"))) {
+    remove_argument (N, argc, argv);
+    NO_CHIPS = TRUE;
+  }
+
+  SIGMA = 0;
+  if ((N = get_argument (*argc, argv, "-sigma"))) {
+    remove_argument (N, argc, argv);
+    SIGMA = atof (argv[N]);
+    remove_argument (N, argc, argv);
+  }
+
+  fake_field_defaults ();
+
+  while ((N = get_argument (*argc, argv, "-p"))) {
+    remove_argument (N, argc, argv);
+    
+    if (!strcmp (argv[N], "help")) {
+      fprintf (stderr, "valid parameters:\n");
+      fprintf (stderr, "cdelt (arcsec/pix)\n");
+      fprintf (stderr, "crpix (Xo) (Yo)\n");
+      fprintf (stderr, "theta (angle)\n");
+      exit (2);
+    }
+
+    if (!strcmp (argv[N], "cdelt")) {
+      remove_argument (N, argc, argv);
+      field.project.cdelt2 = field.project.cdelt1 = atof(argv[N])/3600.0;
+      remove_argument (N, argc, argv);
+      continue;
+    }
+    if (!strcmp (argv[N], "crpix")) {
+      remove_argument (N, argc, argv);
+      field.project.crpix1 = atof(argv[N]);
+      remove_argument (N, argc, argv);
+      field.project.crpix2 = atof(argv[N]);
+      remove_argument (N, argc, argv);
+      continue;
+    }
+    if (!strcmp (argv[N], "theta")) {
+      remove_argument (N, argc, argv);
+      theta = atof(argv[N]);
+      remove_argument (N, argc, argv);
+      field.project.pc1_1 = +cos (RAD_DEG*theta);
+      field.project.pc1_2 = -sin (RAD_DEG*theta);
+      field.project.pc2_1 = +sin (RAD_DEG*theta);
+      field.project.pc2_2 = +cos (RAD_DEG*theta);
+      continue;
+    }
+    if (!strncmp (argv[N], "pca", 3)) {
+      No = argv[N][3] - '0';
+      Nx = argv[N][5] - '0';
+      Ny = argv[N][7] - '0';
+      if ((Nx + Ny > 3) || (No > 1)) {
+	fprintf (stderr, "PCA out of range\n");
+	exit (1);
+      }
+      Np = mkpolyterm (Nx, Ny);
+      remove_argument (N, argc, argv);
+      field.distort.polyterms[Np][No] = atof(argv[N]);
+      remove_argument (N, argc, argv);
+      continue;
+    }
+  }
+
+  if (*argc != 4) print_help ();
+}
Index: /trunk/Ohana/src/mosastro/src/chips.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/chips.c	(revision 3322)
+++ /trunk/Ohana/src/mosastro/src/chips.c	(revision 3323)
@@ -5,17 +5,15 @@
 
   int i, j;
-  double R, D, P, Q, L, M;
+  double R, D, P, Q, L, M, Scale, Det;
 
-# if (0)
   if (CHIPS != (char *) NULL) {
     load_chips (CHIPS);
     return (1);
   }
-# endif
 
   for (i = 0; i < Nchip; i++) {
 
     /* bore site center guess */
-    strcpy (chip[i].map.ctype, "RA---WRP");
+    strcpy (chip[i].map.ctype, "DEC--WRP");
     chip[i].map.crpix1 = 0.0;
     chip[i].map.crpix2 = 0.0;
@@ -30,16 +28,27 @@
     chip[i].map.crval2 = M;
 
-    /** allow guess at field rotation?? add correct handling for negative cdelt1,2 **/
-    chip[i].map.pc1_1  = chip[i].coords.pc1_1;
-    chip[i].map.pc2_2  = chip[i].coords.pc2_2;
-    chip[i].map.pc1_2  = chip[i].coords.pc1_2;
-    chip[i].map.pc2_1  = chip[i].coords.pc2_1;
+    /** we preserve the rotation and parity of coords.pc_ij, but renormalize to unity scale **/
+    Det = chip[i].coords.pc1_1*chip[i].coords.pc2_2 - chip[i].coords.pc1_2*chip[i].coords.pc2_1;
+    Scale = 1.0 / sqrt(fabs(chip[i].coords.cdelt1*chip[i].coords.cdelt2*Det));
+
+    /** test for NaN Scale **/
+
+    chip[i].map.pc1_1  = Scale * chip[i].coords.pc1_1 * chip[i].coords.cdelt1;
+    chip[i].map.pc2_2  = Scale * chip[i].coords.pc2_2 * chip[i].coords.cdelt2;
+    chip[i].map.pc1_2  = Scale * chip[i].coords.pc1_2 * chip[i].coords.cdelt2;
+    chip[i].map.pc2_1  = Scale * chip[i].coords.pc2_1 * chip[i].coords.cdelt1;
 
     chip[i].map.Npolyterms = 1;
     for (j = 0; j < 7; j++) {
-      chip[i].map.polyterms[0][j] = 0;
-      chip[i].map.polyterms[1][j] = 0;
+      chip[i].map.polyterms[j][0] = 0;
+      chip[i].map.polyterms[j][1] = 0;
     }
   }
   return (1);
 }
+
+int load_chips (char *filename) {
+
+  fprintf (stderr, "not ready yet\n");
+  exit (1);
+}
Index: /trunk/Ohana/src/mosastro/src/clip.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/clip.c	(revision 3323)
+++ /trunk/Ohana/src/mosastro/src/clip.c	(revision 3323)
@@ -0,0 +1,28 @@
+# include "mosastro.h"
+
+int ClipOnFP (double Nsigma) {
+
+  int i, j, Nscatter;
+  double dL, dM, dV;
+  double sigma;
+  StarData *raw, *ref;
+
+  sigma = GetScatter (&Nscatter);
+
+  for (i = 0; i < Nchip; i++) {
+    raw = chip[i].raw;
+    ref = chip[i].ref;
+    for (j = 0; j < chip[i].Nstars; j++) {
+      dL = raw[j].L - ref[j].L;
+      dM = raw[j].M - ref[j].M;
+      dV = hypot(dL, dM);
+      if (dV > Nsigma*sigma) {
+	raw[j].mask = TRUE;
+      } else {
+	raw[j].mask = FALSE;	
+      }
+    }
+  }
+  return (TRUE);
+}
+
Index: /trunk/Ohana/src/mosastro/src/dump.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/dump.c	(revision 3322)
+++ /trunk/Ohana/src/mosastro/src/dump.c	(revision 3323)
@@ -6,5 +6,5 @@
 
   for (i = 0; i < Nstars; i++) {
-    fprintf (f, "%4d  %10.6f %10.6f  %8.2f %8.2f  %8.2f %8.2f  %7.2f %7.2f\n", 
+    fprintf (f, "%4d  %10.6f %10.6f  %8.2f %8.2f  %8.2f %8.2f  %7.2f %7.2f  %7.2f %7.2f\n", 
 	     i, 
 	     stars[i].R, stars[i].D,
@@ -16,2 +16,59 @@
   return (1);
 }
+
+int dump_grads (Gradients *grad, char *filename) {
+
+  int i;
+  FILE *f;
+
+  fprintf (stderr, "printing to file %s\n", filename);
+
+  f = fopen (filename, "w");
+
+  for (i = 0; i < grad[0].Npts; i++) {
+    fprintf (f, "%4d  %10.6f %10.6f  %10.6f %10.6f   %10.6f %10.6f\n", 
+	     i, 
+	     grad[0].dPdL[i], grad[0].dPdM[i], 
+	     grad[0].dQdL[i], grad[0].dQdM[i], 
+	     grad[0].Lo[i], grad[0].Mo[i]);
+  }
+  fclose (f);
+  return (1);
+}
+
+int dump_match () {
+
+  int i;
+  FILE *f, *g;
+  f = fopen ("raw.dat", "w");
+  g = fopen ("ref.dat", "w");
+  for (i = 0; i < Nchip; i++) {
+    dump_stars (f, chip[i].raw, chip[i].Nmatch);
+    dump_stars (g, chip[i].ref, chip[i].Nmatch);
+  }
+  fclose (f);
+  fclose (g);
+  exit (1);
+}
+
+int dump_rawstars () {
+
+  int i;
+  FILE *f;
+  f = fopen ("stars.dat", "w");
+  for (i = 0; i < Nchip; i++) {
+    dump_stars (f, chip[i].stars, chip[i].Nstars);
+  }
+  fclose (f);
+  exit (1);
+}
+
+int dump_refcat (StarData *refcat, int Nrefcat) {
+
+  int i;
+  FILE *f;
+  f = fopen ("refcat.dat", "w");
+  dump_stars (f, refcat, Nrefcat);
+  fclose (f);
+  exit (1);
+}
Index: /trunk/Ohana/src/mosastro/src/fakefield.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/fakefield.c	(revision 3322)
+++ /trunk/Ohana/src/mosastro/src/fakefield.c	(revision 3323)
@@ -2,7 +2,5 @@
 
 /* determine an initial guess to field parameters from data */ 
-int fake_field (double RA, double DEC) {
-
-  int i;
+int fake_field_center (double RA, double DEC) {
 
   field.Rmin = RA  - 0.5;
@@ -12,11 +10,21 @@
 
   /* bore site center guess */
-  strcpy (field.project.ctype, "RA---TAN");
   field.project.crval1 = 0.5*(field.Rmin + field.Rmax);
   field.project.crval2 = 0.5*(field.Dmin + field.Dmax);
+  
+  return (1);
+}
+
+/* set default field parameters (overridden in args_obs) */ 
+int fake_field_defaults () {
+
+  int i;
+
+  /* bore site center guess */
+  strcpy (field.project.ctype, "DEC--TAN");
   field.project.crpix1 = 0;
   field.project.crpix2 = 0;
   
-  /* measure average plate scale - would be better using parabolic min... */
+  /* set TP plate scale */
   field.project.cdelt1 = 1.0/3600.0;
   field.project.cdelt2 = 1.0/3600.0;
@@ -29,6 +37,6 @@
   field.project.Npolyterms = 1;
 
-  /* bore site center guess */
-  strcpy (field.distort.ctype, "RA---WRP");
+  /** distort only has power in polyterms **/
+  strcpy (field.distort.ctype, "DEC--WRP");
   field.distort.crval1 = 0.0;
   field.distort.crval2 = 0.0;
@@ -38,5 +46,4 @@
   field.distort.cdelt2 = 1.0;
   
-  /** allow guess at field rotation?? **/
   field.distort.pc1_1  = 1;
   field.distort.pc2_2  = 1;
@@ -44,23 +51,10 @@
   field.distort.pc2_1  = 0;
 
-  /* allow 2nd and 3rd order? */
   /* how do we handle renormalization? (fixed at 1000 pixels??) */
   field.distort.Npolyterms = 3;
   for (i = 0; i < 7; i++) {
-    field.distort.polyterms[0][i] = 0;
-    field.distort.polyterms[1][i] = 0;
+    field.distort.polyterms[i][0] = 0;
+    field.distort.polyterms[i][1] = 0;
   }
-
-  /* megacam distortion : 50 arcsec @ 1800 arcsec */
-  field.distort.polyterms[0][3] = 1e-8;
-  field.distort.polyterms[1][6] = 1e-8;
-
-  /* tiny distortion    : 0.5 arcsec @ 1800 arcsec */
-  field.distort.polyterms[0][3] = 1e-10;
-  field.distort.polyterms[1][6] = 1e-10;
-
-  /* minor distortion    : 5 arcsec @ 1800 arcsec */
-  field.distort.polyterms[0][3] = 1e-9;
-  field.distort.polyterms[1][6] = 1e-9;
 
   return (1);
Index: /trunk/Ohana/src/mosastro/src/field.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/field.c	(revision 3322)
+++ /trunk/Ohana/src/mosastro/src/field.c	(revision 3323)
@@ -7,13 +7,11 @@
   double PS;
 
-# if (0)
   if (FIELD != (char *) NULL) {
     load_field (FIELD);
-    return (1);
+    return (TRUE);
   }
-# endif
 
   /* bore site center guess */
-  strcpy (field.project.ctype, "RA---TAN");
+  strcpy (field.project.ctype, "DEC--TAN");
   field.project.crval1 = 0.5*(field.Rmin + field.Rmax);
   field.project.crval2 = 0.5*(field.Dmin + field.Dmax);
@@ -40,5 +38,5 @@
 
   /* bore site center guess */
-  strcpy (field.distort.ctype, "RA---WRP");
+  strcpy (field.distort.ctype, "DEC--WRP");
   field.distort.crval1 = 0.0;
   field.distort.crval2 = 0.0;
@@ -58,8 +56,8 @@
   field.distort.Npolyterms = 1;
   for (i = 0; i < 7; i++) {
-    field.distort.polyterms[0][i] = 0;
-    field.distort.polyterms[1][i] = 0;
+    field.distort.polyterms[i][0] = 0;
+    field.distort.polyterms[i][1] = 0;
   }
-  return (1);
+  return (TRUE);
 }
 
@@ -87,2 +85,76 @@
 
 }
+
+void field_combine () {
+
+  int i;
+
+  strcpy (field.project.ctype, "DEC--WRP");
+
+  field.project.Npolyterms = field.distort.Npolyterms;
+  for (i = 0; i < 7; i++) {
+    field.project.polyterms[i][0] = field.distort.polyterms[i][0];
+    field.project.polyterms[i][1] = field.distort.polyterms[i][1];
+  }
+}
+
+int load_field (char *filename) {
+
+  int i;
+  Coords coords;
+  Header header;
+
+  /* load header */
+  if (!fits_read_header (filename, &header)) {
+    fprintf (stderr, "ERROR: can't read header for %s\n", filename);
+    exit (1);
+  }
+  /* get astrometry information */
+  if (!GetCoords (&coords, &header)) {
+    fprintf (stderr, "ERROR: no astrometric solution in field %s\n", filename);
+    exit (1);
+  }
+ 
+  /* bore site center guess */
+  strcpy (field.project.ctype, "DEC--TAN");
+  field.project.crval1 = 0.5*(field.Rmin + field.Rmax);
+  field.project.crval2 = 0.5*(field.Dmin + field.Dmax);
+  field.project.crpix1 = 0;
+  field.project.crpix2 = 0;
+  
+  /* measure average plate scale - would be better using parabolic min... */
+  field.project.cdelt1 = coords.cdelt1;
+  field.project.cdelt2 = coords.cdelt2;
+
+  /** allow guess at field rotation?? **/
+  field.project.pc1_1  = coords.pc1_1;
+  field.project.pc2_2  = coords.pc2_2;
+  field.project.pc1_2  = coords.pc1_2;
+  field.project.pc2_1  = coords.pc2_1;
+  field.project.Npolyterms = 1;
+
+  /* bore site center guess */
+  strcpy (field.distort.ctype, "DEC--WRP");
+  field.distort.crval1 = 0.0;
+  field.distort.crval2 = 0.0;
+  field.distort.crpix1 = 0.0;
+  field.distort.crpix2 = 0.0;
+  field.distort.cdelt1 = 1.0;
+  field.distort.cdelt2 = 1.0;
+  
+  /** allow guess at field rotation?? **/
+  field.distort.pc1_1  = 1;
+  field.distort.pc2_2  = 1;
+  field.distort.pc1_2  = 0;
+  field.distort.pc2_1  = 0;
+
+  /* allow 2nd and 3rd order? */
+  /* how do we handle renormalization? (fixed at 1000 pixels??) */
+  field.distort.Npolyterms = coords.Npolyterms;
+  for (i = 0; i < 7; i++) {
+    field.distort.polyterms[i][0] = coords.polyterms[i][0];
+    field.distort.polyterms[i][1] = coords.polyterms[i][1];
+  }
+  return (TRUE);
+}
+
Index: /trunk/Ohana/src/mosastro/src/fitpoly.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/fitpoly.c	(revision 3323)
+++ /trunk/Ohana/src/mosastro/src/fitpoly.c	(revision 3323)
@@ -0,0 +1,281 @@
+# include "mosastro.h"
+
+static int NTERM, NPOWER, NPARS, NORDER, Npts;
+static double **sum, **xsum, **ysum;
+static double **matrix, **vector;
+
+void fit_init (int order) {
+
+  int i;
+
+  Npts  = 0;
+  NORDER = order;
+  NPOWER = NORDER + 1;
+  NTERM = 2*NORDER + 1;
+  NPARS = (NORDER + 1)*(NORDER + 2) / 2;
+
+  /* allocate arrays for fit solution */
+  ALLOCATE (sum, double *, NTERM);
+  ALLOCATE (xsum, double *, NTERM);
+  ALLOCATE (ysum, double *, NTERM);
+  for (i = 0; i < NTERM; i++) {
+    ALLOCATE (sum[i], double, NTERM);
+    bzero (sum[i], NTERM*sizeof(double));
+    ALLOCATE (xsum[i], double, NTERM);
+    bzero (xsum[i], NTERM*sizeof(double));
+    ALLOCATE (ysum[i], double, NTERM);
+    bzero (ysum[i], NTERM*sizeof(double));
+  }
+  ALLOCATE (matrix, double *, NPARS);
+  ALLOCATE (vector, double *, NPARS);
+  for (i = 0; i < NPARS; i++) {
+    ALLOCATE (matrix[i], double, NPARS);
+    ALLOCATE (vector[i], double, 2);
+    bzero (vector[i], 2*sizeof(double));
+    bzero (matrix[i], NPARS*sizeof(double));
+  }
+
+}
+
+# define SCALE 1.0
+void fit_add (double x1, double y1, double x2, double y2) {
+
+  int n, m;
+  double xterm, yterm, term;
+
+  xterm = 1;
+  for (n = 0; n < NTERM; n++) {
+    yterm = 1;
+    for (m = 0; m < NTERM; m++) {
+      term = xterm*yterm;
+      if (n+m < NTERM) {
+	sum[n][m] += term;
+      }
+      if (n+m < NPOWER) {
+	xsum[n][m] += x2*term;
+	ysum[n][m] += y2*term;
+      }
+      yterm *= y1/SCALE;
+    }
+    xterm *= x1/SCALE;
+  }
+  Npts ++;
+}
+
+/** I am renormalizing here by the max pivots to keep gaussj sane **
+ ** would not be needed if the fit used scaled ind. variables **/
+void fit_eval () {
+
+  int i, j, n, m, M, N;
+  double max;
+
+  i = 0;
+  for (m = 0; m < NPOWER; m++) {
+    for (n = 0; n < NPOWER - m; n++, i++) {
+      vector[i][0] = xsum[n][m];
+      vector[i][1] = ysum[n][m];
+    }	
+  }
+  j = 0;
+  for (M = 0; M < NPOWER; M++) {
+    for (N = 0; N < NPOWER - M; N++, j++) {
+      i = 0;
+      for (m = 0; m < NPOWER; m++) {
+	for (n = 0; n < NPOWER - m; n++, i++) {
+	  matrix[i][j] = sum[n+N][m+M];
+	}	
+      }
+    }
+  }       
+# if (0)
+  max = 0.0;
+  for (i = 0; i < NPARS; i++) {
+    for (j = 0; j < NPARS; j++) {
+      max = MAX (max, fabs(matrix[i][j]));
+    }
+    max = MAX (max, fabs(vector[i][0]));
+    max = MAX (max, fabs(vector[i][1]));
+  }
+  for (i = 0; i < NPARS; i++) {
+    for (j = 0; j < NPARS; j++) {
+      matrix[i][j] /= max;
+    }
+    vector[i][0] /= max;
+    vector[i][1] /= max;
+  }
+# endif
+
+  gaussj (matrix, NPARS, vector, 2); 
+
+# if (0)
+  i = 0;
+  for (m = 0; m < NPOWER; m++) {
+    for (n = 0; n < NPOWER - m; n++, i++) {
+      fprintf (stderr, "RA x^%dy^%d: %10.4g    DEC x^%dy^%d: %10.4g \n", 
+	       n, m, vector[i][0], n, m, vector[i][1]);
+    }	
+  }
+# endif
+}
+
+/* linear portion of fit : NORDER is 1 */
+void fit_apply_coords (Coords *coords) {
+
+  int i, j, Np, Nv, N;
+  double c10, c11, c12;
+  double c20, c21, c22;
+  double R, D;
+
+  /* update the higher order terms */
+  if (NORDER > 1) {
+    for (i = 0; i < NPOWER; i++) {
+      for (j = 0; j < (NPOWER - i); j++) {
+	if (i + j < 2) continue;
+	Np = mkpolyterm (i, j);
+	Nv = mkvector (i, j, NORDER);
+	coords[0].polyterms[Np][0] = vector[Nv][0];
+	coords[0].polyterms[Np][1] = vector[Nv][1];
+      }
+    }
+  }
+
+  /* get the correct vector entries for the linear terms */
+  N = mkvector (0, 0, NORDER);
+  coords[0].crval1 = vector[N][0];  
+  coords[0].crval2 = vector[N][1];
+
+  N = mkvector (1, 0, NORDER);
+  c11 = vector[N][0];  
+  c21 = vector[N][1];
+  N = mkvector (0, 1, NORDER);
+  c12 = vector[N][0];  
+  c22 = vector[N][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;
+
+  coords[0].crpix1 = 0;
+  coords[0].crpix2 = 0;
+
+  coords[0].Npolyterms = NORDER;
+  strcpy (coords[0].ctype, "DEC--WRP");
+}
+
+/*
+  if we have just linear terms, the following holds for crpix1,2:
+  D = R / (coords[0].pc1_1*coords[0].pc2_2 - coords[0].pc1_2*coords[0].pc2_1);
+  coords[0].crpix1 = D * (coords[0].pc1_2*c20 - coords[0].pc2_2*c10);
+  coords[0].crpix2 = D * (coords[0].pc2_1*c10 - coords[0].pc1_1*c20);
+*/
+
+/* NORDER is order of gradient fit : Npolyterms is Norder + 1 */
+fit_apply_grads (Coords *distort, Coords *project, int term) {
+
+  int i, j, Np, Nv1, Nv2;
+
+  distort[0].Npolyterms = NORDER + 1;
+  if (distort[0].Npolyterms > 1) strcpy (distort[0].ctype, "DEC--WRP");
+
+  for (i = 0; i < NPOWER + 1; i++) {
+    for (j = 0; j < (NPOWER + 1 - i); j++) {
+      
+      if (i + j < 2) continue;
+      Np  = mkpolyterm (i, j);
+      Nv1 = mkvector (i-1, j, NORDER);
+      Nv2 = mkvector (i, j-1, NORDER);
+
+      /** why do we have the negative sign? **/
+      if (j == 0) {
+	distort[0].polyterms[Np][term] = vector[Nv1][0] / i;
+      }
+      if (i == 0) {
+	distort[0].polyterms[Np][term] = vector[Nv2][1] / j;
+      }
+      if ((i > 0) && (j > 0)) {
+	distort[0].polyterms[Np][term] = 0.5*(vector[Nv1][0] / i + vector[Nv2][1] / j);
+      }
+    }
+  }
+
+  Nv1 = mkvector (0, 0, NORDER);
+  if (term == 0) {
+    project[0].pc1_1 = project[0].pc1_1 * (1 + vector[Nv1][0]);
+    project[0].pc1_2 = project[0].pc1_2 * (1 + vector[Nv1][1]);
+  } else {
+    project[0].pc2_1 = project[0].pc2_1 * (1 + vector[Nv1][0]);
+    project[0].pc2_2 = project[0].pc2_2 * (1 + vector[Nv1][1]);
+  }
+}
+
+fit_correct_grads (Gradients *in, Gradients *out, int term) {
+
+  int i, k, m, n;
+  double x, y, dx, dy, dz1, dz2;
+
+  for (i = 0; i < in[0].Npts; i++) {
+    
+    dx = in[0].Lo[i];
+    dy = in[0].Mo[i];
+    dz1 = dz2 = 0.0;
+
+    k = 0;
+    x = y = 1;
+    for (m = 0; m < NPOWER; m++) {
+      x = y;
+      for (n = 0; n < NPOWER - m; n++, k++) {
+	dz1 += vector[k][0]*x;
+	dz2 += vector[k][1]*x;
+	x = x * dx / SCALE;
+      }
+      y = y * dy / SCALE;
+    }
+
+    out[0].Lo[i] = dx;
+    out[0].Mo[i] = dy;
+    if (term == 0) {
+      out[0].dPdL[i] = in[0].dPdL[i] - dz1;
+      out[0].dPdM[i] = in[0].dPdM[i] - dz2;
+    } else {
+      out[0].dQdL[i] = in[0].dQdL[i] - dz1;
+      out[0].dQdM[i] = in[0].dQdM[i] - dz2;
+    }
+  }
+}
+
+int mkvector (int n, int m, int norder) {
+  
+  int i, N;
+  
+  N = 0;
+  for (i = 0; i < m; i++) {
+    N += (norder - i + 1);
+  }
+  N += n;
+  return (N);
+}
+
+void fit_free () {
+
+  int i;
+
+  for (i = 0; i < NTERM; i++) {
+    free (sum[i]);
+    free (xsum[i]);
+    free (ysum[i]);
+  }
+  free (sum);
+  free (xsum);
+  free (ysum);
+
+  for (i = 0; i < NPARS; i++) {
+    free (matrix[i]);
+    free (vector[i]);
+  }
+  free (matrix);
+  free (vector);
+}
+  
Index: /trunk/Ohana/src/mosastro/src/gaussj.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/gaussj.c	(revision 3323)
+++ /trunk/Ohana/src/mosastro/src/gaussj.c	(revision 3323)
@@ -0,0 +1,75 @@
+# include "mosastro.h"
+
+int gaussj (double **a, int n, double **b, int m) {
+
+  int *indxc,*indxr,*ipiv;
+  int i, icol, irow, j, k, l, ll;
+  double big,dum,pivinv;
+  
+  ALLOCATE (indxc, int, n);
+  ALLOCATE (indxr, int, n);
+  ALLOCATE (ipiv, int, n);
+  for (j = 0; j < n; j++) 
+    ipiv[j] = 0;
+
+  for (i = 0; i < n; i++) {
+    big = 0.0;
+    for (j = 0; j < n; j++) {
+      if (ipiv[j] != 1) {
+	for (k = 0; k < n; k++) {
+	  if (ipiv[k] == 0) {
+	    if (fabs (a[j][k]) >= big) {
+	      big  = fabs (a[j][k]);
+	      irow = j;
+	      icol = k;
+	    }
+	  } 
+	  else 
+	    if (ipiv[k] > 1) {
+	      /* fprintf (stderr, "GAUSSJ: Singular Matrix! (1)\n"); */
+	      return (0);
+	    }
+	}
+      }
+    }
+    ipiv[icol]++;
+    if (irow != icol) {
+      for (l = 0; l < n; l++) 
+	SWAP (a[irow][l], a[icol][l]);
+      for (l = 0; l < m; l++) 
+	SWAP (b[irow][l], b[icol][l]);
+    }
+    indxr[i] = irow;
+    indxc[i] = icol;
+    if (a[icol][icol] == 0.0) {
+      /* fprintf (stderr, "GAUSSJ: Singular Matrix! (2)\n"); */
+      return (0);
+    }
+    pivinv = 1.0 / a[icol][icol];
+    a[icol][icol] = 1.0;
+    for (l = 0; l < n; l++) 
+      a[icol][l] *= pivinv;
+    for (l = 0; l < m; l++) 
+      b[icol][l] *= pivinv;
+    for (ll = 0; ll < n; ll++) {
+      if (ll != icol) {
+	dum = a[ll][icol];
+	a[ll][icol] = 0.0;
+	for (l = 0; l < n; l++) 
+	  a[ll][l] -= a[icol][l]*dum;
+	for (l = 0; l < m; l++) 
+	  b[ll][l] -= b[icol][l]*dum;
+      }
+    }
+  }
+
+  for (l = n - 1; l >= 0; l--) {
+    if (indxr[l] != indxc[l])
+      for (k = 0; k < n; k++)
+	SWAP (a[k][indxr[l]], a[k][indxc[l]]);
+  }
+  free (ipiv);
+  free (indxr);
+  free (indxc);
+  return (1);
+}
Index: /trunk/Ohana/src/mosastro/src/getstone.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/getstone.c	(revision 3323)
+++ /trunk/Ohana/src/mosastro/src/getstone.c	(revision 3323)
@@ -0,0 +1,89 @@
+# include "mosastro.h"
+
+StarData *getstone (CatStats *input, int *nstars) {
+
+  FILE *f;
+  int i, Nregion, NREGION, Nstars, NSTARS;
+  double r, d, R, D, M, Tr, Td, dRdT, dDdT;
+  char **regname, filename[1024], line[1024];
+  StarData *stars;
+  CatStats *regions;
+
+  sprintf (filename, "%s/Regions.dat", StoneRegions);
+  f = fopen (filename, "r");
+  if (f == (FILE *) NULL) {
+    fprintf (stderr, "ERROR: can't open stone %s\n", filename);
+    exit (1);
+  }
+
+  Nregion = 0;
+  NREGION = 20;
+  ALLOCATE (regions, CatStats, NREGION);
+  ALLOCATE (regname, char *, NREGION);
+
+  /* strip off first commented line */
+  scan_line (f, line);
+  while (scan_line (f, line) != EOF) {
+    sscanf (line, "%s %lf %lf %lf %lf", filename, 
+	    &regions[Nregion].RA[0], &regions[Nregion].RA[1], 
+	    &regions[Nregion].DEC[0], &regions[Nregion].DEC[1]);
+    sprintf (line, "%s/%s", StoneRegions, filename);
+    regname[Nregion] = strcreate (line);
+    Nregion ++;
+    if (Nregion == NREGION) {
+      NREGION += 50;
+      REALLOCATE (regions, CatStats, NREGION);
+      REALLOCATE (regname, char *, NREGION);
+    }	  
+  }
+  fclose (f);
+
+  Nstars = 0;
+  NSTARS = 1000;
+  ALLOCATE (stars, StarData, NSTARS);
+
+  /* find the region(s) that overlap the input region */
+  for (i = 0; i < Nregion; i++) {
+    if (input[0].RA[0] > regions[i].RA[1]) continue;
+    if (input[0].RA[1] < regions[i].RA[0]) continue;
+    if (input[0].DEC[0] > regions[i].DEC[1]) continue;
+    if (input[0].DEC[1] < regions[i].DEC[0]) continue;
+    /* this will fail on 0,360 boundary */
+    fprintf (stderr, "loading data from %s\n", regname[i]);
+
+    f = fopen (regname[i], "r");
+    if (f == (FILE *) NULL) {
+      fprintf (stderr, "ERROR: can't open stone %s\n", regname[i]);
+      exit (1);
+    }
+
+    while (scan_line (f, line) != EOF) {
+    
+      dparse (&R, 3, line);
+      dparse (&D, 4, line);
+      dparse (&M, 5, line);
+      
+      dparse (&Tr, 12, line);
+      dparse (&Td, 13, line);
+      dparse (&dRdT, 14, line);
+      dparse (&dDdT, 15, line);
+      
+      r = R + dRdT*(Year - Tr)/(100.0*3600.0);
+      d = D + dDdT*(Year - Td)/(100.0*3600.0);
+      
+      stars[Nstars].M = M;
+      stars[Nstars].R = r*15.0;
+      stars[Nstars].D = d;
+      Nstars ++;
+      if (Nstars == NSTARS) {
+	NSTARS += 5000;
+	REALLOCATE (stars, StarData, NSTARS);
+      }	  
+    }
+    fclose (f);
+  }
+
+  *nstars = Nstars;
+  return (stars);
+}
+
Index: /trunk/Ohana/src/mosastro/src/getusno.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/getusno.c	(revision 3322)
+++ /trunk/Ohana/src/mosastro/src/getusno.c	(revision 3323)
@@ -24,6 +24,6 @@
 
   RA0  = catstats[0].RA[0]; 
-  DEC0 = catstats[0].RA[1]; 
-  RA1  = catstats[0].DEC[0];
+  RA1  = catstats[0].RA[1]; 
+  DEC0 = catstats[0].DEC[0];
   DEC1 = catstats[0].DEC[1];
 
Index: /trunk/Ohana/src/mosastro/src/greference.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/greference.c	(revision 3322)
+++ /trunk/Ohana/src/mosastro/src/greference.c	(revision 3323)
@@ -16,4 +16,9 @@
 
   if (VERBOSE) fprintf (stderr, "full region: %f - %f, %f - %f\n", catstats.RA[0], catstats.RA[1], catstats.DEC[0], catstats.DEC[1]);
+
+  /* get stars from the Stone et al catalog for the given region */
+  if (!strcmp (REFCAT, "STONE")) {
+    stars = getstone (&catstats, &Nstars);
+  }
 
   /* get stars from the USNO catalog for the given region */
Index: /trunk/Ohana/src/mosastro/src/match.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/match.c	(revision 3322)
+++ /trunk/Ohana/src/mosastro/src/match.c	(revision 3323)
@@ -6,7 +6,10 @@
   int i, j, k, N, K, Ntotal;
   int Nmatch, NMATCH;
-  double dp, dq, radius;
+  double dp, dq, radius, Radius;
   double *p, *q, *P, *Q;
   int *U, *u;
+
+  /* requested radius is in arcsec ; internally radius is in pixels */
+  Radius = RADIUS / (3600.0 * field.project.cdelt2);
 
   /* sort the REFCAT data by P in tangent plane */ 
@@ -42,11 +45,19 @@
     /* find star matches in the tangent plane coord system */
     for (j = k = 0; (j < chip[i].Nstars) && (k < Nrefcat); ) {
+      /* skip anything with dMag too large */
+      if (SIGMA_LIM > 0.0) {
+	if (chip[i].stars[u[j]].dMag > SIGMA_LIM) {
+	  j++;
+	  continue;
+	}
+      }
+	
       dp = p[j] - P[k];
 	
-      if (dp <= -2*RADIUS) {
+      if (dp <= -2*Radius) {
 	j++;
 	continue;
       }
-      if (dp >= 2*RADIUS) {
+      if (dp >= 2*Radius) {
 	k++;
 	continue;
@@ -54,9 +65,9 @@
 
       K = k;
-      for (; (dp > -2*RADIUS) && (k < Nrefcat); k++) {
+      for (; (dp > -2*Radius) && (k < Nrefcat); k++) {
 	dp = p[j] - P[k];
 	dq = q[j] - Q[k];
 	radius = hypot (dp, dq);
-	if (radius < RADIUS) {
+	if (radius < Radius) {
 	  chip[i].raw[Nmatch] = chip[i].stars[u[j]];
 	  chip[i].ref[Nmatch] = refcat[U[k]];
Index: /trunk/Ohana/src/mosastro/src/mkheader.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/mkheader.c	(revision 3322)
+++ /trunk/Ohana/src/mosastro/src/mkheader.c	(revision 3323)
@@ -7,4 +7,5 @@
   ALLOCATE (header, Header, 1);
 
+  header[0].simple = TRUE;
   header[0].bitpix = -32;
   header[0].unsign = FALSE;
@@ -22,46 +23,9 @@
   fits_modify (header, "UTC-OBS",  "%s", 1, "14:27:45.30");
   fits_modify (header, "ZERO_PT", "%lf", 1, 25.0);
+  fits_modify (header, "EXPTIME", "%lf", 1, 2.0);
 
-  if (coords[0].Npolyterms > 1) {
-    fits_modify (header, "CTYPE1",   "%s",  1, "RA---PLY");
-    fits_modify (header, "CTYPE2",   "%s",  1, "DEC--PLY");
-  } else {
-    fits_modify (header, "CTYPE1",   "%s",  1, "RA---TAN");
-    fits_modify (header, "CTYPE2",   "%s",  1, "DEC--TAN");
-  }    
+  PutCoords (coords, header);
+
   fits_modify (header, "NASTRO",   "%d", 1, 1); 
-
-  fits_modify (header, "CDELT1",   "%le", 1, coords[0].cdelt1); 
-  fits_modify (header, "CDELT2",   "%le", 1, coords[0].cdelt2);
-  fits_modify (header, "CRVAL1",   "%lf", 1, coords[0].crval1);
-  fits_modify (header, "CRVAL2",   "%lf", 1, coords[0].crval2);  
-  fits_modify (header, "CRPIX1",   "%lf", 1, coords[0].crpix1);
-  fits_modify (header, "CRPIX2",   "%lf", 1, coords[0].crpix2);
-  fits_modify (header, "PC001001", "%le", 1, coords[0].pc1_1);
-  fits_modify (header, "PC001002", "%le", 1, coords[0].pc1_2);
-  fits_modify (header, "PC002001", "%le", 1, coords[0].pc2_1);
-  fits_modify (header, "PC002002", "%le", 1, coords[0].pc2_2);
-  fits_modify (header, "NPLYTERM", "%d", 1, coords[0].Npolyterms);
-  if (coords[0].Npolyterms > 1) {
-    fits_modify (header, "PCA1X2Y0", "%le", 1, coords[0].polyterms[0][0]);
-    fits_modify (header, "PCA1X1Y1", "%le", 1, coords[0].polyterms[1][0]);
-    fits_modify (header, "PCA1X0Y2", "%le", 1, coords[0].polyterms[2][0]);
-
-    fits_modify (header, "PCA2X2Y0", "%le", 1, coords[0].polyterms[0][1]);
-    fits_modify (header, "PCA2X1Y1", "%le", 1, coords[0].polyterms[1][1]);
-    fits_modify (header, "PCA2X0Y2", "%le", 1, coords[0].polyterms[2][1]);
-  }
-      
-  if (coords[0].Npolyterms > 2) {
-    fits_modify (header, "PCA1X3Y0", "%le", 1, coords[0].polyterms[3][0]);
-    fits_modify (header, "PCA1X2Y1", "%le", 1, coords[0].polyterms[4][0]);
-    fits_modify (header, "PCA1X1Y2", "%le", 1, coords[0].polyterms[5][0]);
-    fits_modify (header, "PCA1X0Y3", "%le", 1, coords[0].polyterms[6][0]);
-
-    fits_modify (header, "PCA2X3Y0", "%le", 1, coords[0].polyterms[3][1]);
-    fits_modify (header, "PCA2X2Y1", "%le", 1, coords[0].polyterms[4][1]);
-    fits_modify (header, "PCA2X1Y2", "%le", 1, coords[0].polyterms[5][1]);
-    fits_modify (header, "PCA2X0Y3", "%le", 1, coords[0].polyterms[6][1]);
-  }
 
   return (header);
Index: /trunk/Ohana/src/mosastro/src/mkmosaic.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/mkmosaic.c	(revision 3323)
+++ /trunk/Ohana/src/mosastro/src/mkmosaic.c	(revision 3323)
@@ -0,0 +1,78 @@
+# include "mosastro.h"
+
+Header *mkmosaic (int Nx, int Ny, int Nstars, Coords *coords) {
+
+  int i, Nastro;
+  double Xmin, Xmax, Ymin, Ymax, tmp, Cerror, Cprecise;
+  Header *header;
+  char line[80];
+
+  /* mosaic header has info to define the outline of the mosaic */
+  RD_to_XY (&Xmin, &Ymin, field.Rmin, field.Dmin, coords);
+  RD_to_XY (&Xmax, &Ymax, field.Rmax, field.Dmax, coords);
+
+  ALLOCATE (header, Header, 1);
+
+  header[0].bitpix = -32;
+  header[0].unsign = FALSE;
+  header[0].bscale = 1.0;
+  header[0].bzero  = 0.0;
+  header[0].Naxes = 2;
+  header[0].Naxis[0] = Xmax - Xmin;
+  header[0].Naxis[1] = Ymax - Ymin;
+
+  fits_create_header (header);
+
+  /* calculate image-wide astrometry scatter */
+  Cerror   = GetScatter (&Nastro);
+  Cprecise = Cerror / sqrt (Nastro);
+  fits_modify (header, "CERROR",   "%lf", 1, Cerror);
+  fits_modify (header, "CPRECISE", "%lf", 1, Cprecise);
+  fits_modify (header, "NSTARS",   "%d", 1, 1);  /*** modify addstar to allow NSTARS = 0 ***/
+  fits_modify (header, "NASTRO",   "%d", 1, Nastro);
+
+  /* make a better selection for the mosaic photcode? */
+  fits_scan (&chip[0].header, "PHOTCODE", "%s", 1, line);
+  fits_modify (header, "PHOTCODE", "%s", 1, line);
+
+  fits_scan (&chip[0].header, "ZERO_PT", "%lf", 1, &tmp);
+  fits_modify (header, "ZERO_PT", "%lf", 1, tmp);
+
+  fits_scan (&chip[0].header, ExptimeKeyword,  "%lf", 1, &tmp);
+  fits_modify (header, ExptimeKeyword,  "%lf", 1, tmp);
+
+  /* insert time data from chip[0] */
+  /* try JD first */
+  if (strcasecmp (JDKeyword, "NONE")) {
+    uppercase (JDKeyword);
+    fits_scan (&chip[0].header, JDKeyword, "%lf", 1, &tmp);
+    fits_modify (header, JDKeyword, "%lf", 1, tmp);
+    goto got_time;
+  }
+
+  /* try MJD next */
+  if (strcasecmp (MJDKeyword, "NONE")) {
+    uppercase (MJDKeyword);
+    fits_scan (&chip[0].header, MJDKeyword, "%lf", 1, &tmp);
+    fits_modify (header, MJDKeyword, "%lf", 1, tmp);
+    goto got_time;
+  }
+    
+  /* get UT and DATE */
+  if (strcasecmp (UTKeyword, "NONE") && strcasecmp (DateKeyword, "NONE")) {
+    uppercase (UTKeyword);
+    fits_scan (&chip[0].header, UTKeyword, "%s", 1, line);
+    fits_modify (header, UTKeyword, "%s", 1, line);
+    uppercase (DateKeyword);
+    fits_scan (&chip[0].header, DateKeyword, "%s",  1, line);
+    fits_modify (header, DateKeyword, "%s", 1, line);
+    goto got_time;
+  }
+  fprintf (stderr, "ERROR: missing time abstraction in config\n");
+  exit (1);
+
+got_time:
+
+  PutCoords (coords, header);
+  return (header);
+}
Index: /trunk/Ohana/src/mosastro/src/mkobs.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/mkobs.c	(revision 3322)
+++ /trunk/Ohana/src/mosastro/src/mkobs.c	(revision 3323)
@@ -4,7 +4,7 @@
 
   int i, j, N, Nrefcat, Nchip, Nstars, NSTARS;
+  double Ro, Do, Po, Qo, Lo, Mo;
+  double dX, dY, Dist;
   double RA, DEC;
-  double Ro, Do, Po, Qo, Lo, Mo;
-  double dX, dY;
   char filename[64];
   Coords map, coords;
@@ -14,18 +14,25 @@
 
   ConfigInit (&argc, argv);
-
-  if (argc != 4) {
-    fprintf (stderr, "USAGE: mkobs (RA) (DEC) (output)\n");
-    exit (2);
-  }
+  args (&argc, argv);
+  init_random ();
 
   RA  = atof(argv[1]);
   DEC = atof(argv[2]);
-  dX = dY = 3600*0.1;
+  OUTPUT = argv[3];
+  fake_field_center (RA, DEC);
 
-  fake_field (RA, DEC);
+  /* chip size in pixels (output is 5x5 grid of chips */
+  dX = 0.1 / field.project.cdelt1;
+  dY = 0.1 / field.project.cdelt2;
 
   refcat = greference (&Nrefcat);
   project_refcat (refcat, Nrefcat);
+  if (FOCAL_PLANE != NULL) {
+    FILE *f;
+    f = fopen (FOCAL_PLANE, "w");
+    dump_stars (f, refcat, Nrefcat);
+    fclose (f);
+  }
+  if (NO_CHIPS) exit (0);
 
   Nchip = 0;
@@ -34,11 +41,11 @@
 
       /* find chip center in TP, FP and Sky */
-      Po = 0.2*i*3600.0;
-      Qo = 0.2*j*3600.0;
+      Po = 2*dX*i;
+      Qo = 2*dY*j;
       RD_to_XY (&Lo, &Mo, Po, Qo, &field.distort);
       XY_to_RD (&Ro, &Do, Po, Qo, &field.project);
 
       /* FP-Chip terms */
-      strcpy (map.ctype, "RA---WRP");
+      strcpy (map.ctype, "DEC--WRP");
       map.crval1 = Lo;
       map.crval2 = Mo;
@@ -53,5 +60,11 @@
       map.pc1_2  = 0;
       map.pc2_1  = 0;
-      map.Npolyterms = 1;
+
+      map.Npolyterms = 2;
+      for (N = 0; N < 7; N++) {
+	map.polyterms[N][0] = 0;
+	map.polyterms[N][1] = 0;
+      }
+      map.polyterms[2][1] = 1e-5;
 
       /* project catalog stars to chip */
@@ -69,6 +82,7 @@
 	if (refcat[N].Y >= 2*dY) continue;
 
-	stars[Nstars].X = refcat[N].X;
-	stars[Nstars].Y = refcat[N].Y;
+	/* add random noise - not gaussian noise */
+	stars[Nstars].X = refcat[N].X + SIGMA*(drand48() - 0.5);
+	stars[Nstars].Y = refcat[N].Y + SIGMA*(drand48() - 0.5);
 	stars[Nstars].M = 16.0;
 	stars[Nstars].dM = 0.02;
@@ -88,13 +102,25 @@
       }
 
+
       /* Chip-Sky terms */
-      strcpy (coords.ctype, "RA---WRP");
+      strcpy (coords.ctype, "DEC--WRP");
       coords.crval1 = Ro;
       coords.crval2 = Do;
       coords.crpix1 = dX;
-      coords.crpix2 = dX;
+      coords.crpix2 = dY;
   
-      coords.cdelt1 = 1/3600.0;
-      coords.cdelt2 = 1/3600.0;
+      { 
+	double dP, dQ, dL, dM;
+	double Mx, Lx, scale;
+
+	dP = 10;
+	dQ = 10;
+	RD_to_XY (&Lx, &Mx, Po+dP, Qo+dQ, &field.distort);
+	dL = Lx - Lo;
+	dM = Mx - Mo;
+	scale = hypot (10.0, 10.0) / hypot (dL, dM);
+	coords.cdelt1 = scale/3600.0;
+	coords.cdelt2 = scale/3600.0;
+      }
 
       coords.pc1_1  = 1;
@@ -105,8 +131,6 @@
 
       header = mkheader (2*dX, 2*dY, Nstars, &coords);
-      fits_modify (header, "C_L0", "%lf", 1, Lo);
-      fits_modify (header, "C_M0", "%lf", 1, Mo);
 
-      sprintf (filename, "%s.%02d.fits", argv[3], Nchip);
+      sprintf (filename, "%s.%02d.fits", OUTPUT, Nchip);
       wstars (filename, stars, Nstars, header);
       fits_free_header (header);
Index: /trunk/Ohana/src/mosastro/src/mkpolyterm.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/mkpolyterm.c	(revision 3323)
+++ /trunk/Ohana/src/mosastro/src/mkpolyterm.c	(revision 3323)
@@ -0,0 +1,14 @@
+# include "mosastro.h"
+
+int mkpolyterm (int n, int m) {
+  
+  int i, nt, N;
+  
+  N = 0;
+  nt = n + m;
+  for (i = 2; i < nt; i++) {
+    N += i + 1;
+  }
+  N += m;
+  return (N);
+}
Index: /trunk/Ohana/src/mosastro/src/mkstandards.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/mkstandards.c	(revision 3322)
+++ /trunk/Ohana/src/mosastro/src/mkstandards.c	(revision 3323)
@@ -1,17 +1,27 @@
 # include "mosastro.h"
-# define DX 100.0
-# define DY 100.0
+# define DX 30.0
+# define DY 30.0
+
+extern double drand48();
 
 /* build a grid of reference stars */
 int main (int argc, char **argv) {
 
-  int Nstars, NSTARS;
-  double x, y, dX, dY;
+  int i, N, Nstars, NSTARS, Random;
+  double x, y, m, dX, dY;
   Header *header;
   Coords coords;
   Stars *stars;
 
+  Random = FALSE;
+  if ((N = get_argument (argc, argv, "-random"))) {
+    Random = TRUE;
+    remove_argument (N, &argc, argv);
+    Nstars = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
   if (argc != 5) {
-    fprintf (stderr, "USAGE: mkstandards (RA) (DEC) (half-width) (output)\n");
+    fprintf (stderr, "USAGE: mkstandards (RA) (DEC) (half-width) (output) [-random Nstars]\n");
     exit (2);
   }
@@ -21,5 +31,5 @@
 
   /* bore site center guess */
-  strcpy (coords.ctype, "RA---TAN");
+  strcpy (coords.ctype, "DEC--TAN");
   coords.crval1 = atof(argv[1]);
   coords.crval2 = atof(argv[2]);
@@ -36,28 +46,47 @@
   coords.Npolyterms = 1;
 
-  Nstars = 0;
-  NSTARS = 1000;
-  ALLOCATE (stars, Stars, NSTARS);
+  if (Random) {
+    
+    init_random ();
 
-  for (x = 0; x < 2*dX; x += DX) {
-    for (y = 0; y < 2*dY; y += DY) {
-      stars[Nstars].X = x;
-      stars[Nstars].Y = y;
-      stars[Nstars].M = 16.0;
-      stars[Nstars].dM = 0.02;
-      stars[Nstars].dophot = 1;
-      stars[Nstars].sky = 1.0;
-      stars[Nstars].Mgal = 16.0;
-      stars[Nstars].Map = 16.0;
-      stars[Nstars].fx = 1.0;
-      stars[Nstars].fy = 1.0;
-      stars[Nstars].df = 0.0;
+    ALLOCATE (stars, Stars, Nstars);
+    for (i = 0; i < Nstars; i++) {
+      x = 2*dX*drand48();
+      y = 2*dY*drand48();
+      m = 14.0 + 4.0*drand48();
+      stars[i].X = x;
+      stars[i].Y = y;
+      stars[i].M = m;
+      stars[i].dM = 0.02;
+    }
+  } else {
 
-      Nstars ++;
-      if (Nstars >= NSTARS) {
-	NSTARS += 1000;
-	REALLOCATE (stars, Stars, NSTARS);
+    Nstars = 0;
+    NSTARS = 1000;
+    ALLOCATE (stars, Stars, NSTARS);
+
+    for (x = 0; x < 2*dX; x += DX) {
+      for (y = 0; y < 2*dY; y += DY) {
+	stars[Nstars].X = x;
+	stars[Nstars].Y = y;
+	stars[Nstars].M = 16.0;
+
+	Nstars ++;
+	if (Nstars >= NSTARS) {
+	  NSTARS += 1000;
+	  REALLOCATE (stars, Stars, NSTARS);
+	}
       }
     }
+  }
+
+  for (i = 0; i < Nstars; i++) {
+    stars[i].dophot = 1;
+    stars[i].sky = 1.0;
+    stars[i].Mgal = 16.0;
+    stars[i].Map = 16.0;
+    stars[i].fx = 1.0;
+    stars[i].fy = 1.0;
+    stars[i].df = 0.0;
   }
 
Index: /trunk/Ohana/src/mosastro/src/mosastro.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/mosastro.c	(revision 3322)
+++ /trunk/Ohana/src/mosastro/src/mosastro.c	(revision 3323)
@@ -4,5 +4,6 @@
 int main (int argc, char **argv) {
 
-  int Nrefcat;
+  int Nrefcat, Nastro;
+  double Cerror;
   glob_t pglob;
   StarData *refcat;
@@ -12,6 +13,6 @@
   args (&argc, argv); 
 
-  if (argc != 3) { 
-    fprintf (stderr, "USAGE: mosastro (inglob) (ext)\n");
+  if (argc != 4) { 
+    fprintf (stderr, "USAGE: mosastro (inglob) (ext) (mosaic.phu)\n");
     exit (2);
   }
@@ -30,60 +31,44 @@
   /* use field model to get TP & FP coords */
   project_stars ();
+  if ((DUMP != NULL) && !strcmp (DUMP, "rawstars")) dump_rawstars();
 
   refcat = greference (&Nrefcat);
   project_refcat (refcat, Nrefcat);
+  if ((DUMP != NULL) && !strcmp (DUMP, "refcat")) dump_refcat(refcat, Nrefcat);
 
   match (refcat, Nrefcat);
-
-  if (1) { 
-    int i;
-    FILE *f, *g;
-    f = fopen ("ref.dat", "w");
-    g = fopen ("raw.dat", "w");
-    for (i = 0; i < Nchip; i++) {
-      dump_stars (f, chip[i].raw, chip[i].Nmatch);
-      dump_stars (g, chip[i].ref, chip[i].Nmatch);
-    }
-    fclose (f);
-    fclose (g);
-    exit (1);
-  }
-
   deproject_raw ();
   project_ref ();
+  Cerror   = GetScatter (&Nastro);
+  fprintf (stderr, "scatter (raw)  : %f for %d stars\n", Cerror, Nastro);
+  if ((DUMP != NULL) && !strcmp (DUMP, "rawmatch")) dump_match();
 
   grad = GetGradients ();
   FitGradients (grad);
+  Cerror   = GetScatter (&Nastro);
+  fprintf (stderr, "scatter (grad) : %f for %d stars\n", Cerror, Nastro);
+  if ((DUMP != NULL) && !strcmp (DUMP, "fitgrads")) dump_match();
 
-  deproject_raw ();
-  project_ref ();
-
-  FitChips (1);
-
-  deproject_raw ();
-  project_ref ();
+  FitChips (ChipOrder);
+  Cerror   = GetScatter (&Nastro);
+  fprintf (stderr, "scatter (chip) : %f for %d stars\n", Cerror, Nastro);
+  if ((DUMP != NULL) && !strcmp (DUMP, "fitchips_unclip")) dump_match();
 
   ClipOnFP (3.0);
+  Cerror   = GetScatter (&Nastro);
+  fprintf (stderr, "scatter (clip) : %f for %d stars\n", Cerror, Nastro);
 
+  FitChips (ChipOrder);
+  Cerror   = GetScatter (&Nastro);
+  fprintf (stderr, "scatter (chip) : %f for %d stars\n", Cerror, Nastro);
+  if ((DUMP != NULL) && !strcmp (DUMP, "fitchips")) dump_match();
+
+  output (argv[2], argv[3]);
+  exit (0);
+
+/*
   FitField (3);
+  FitChips (3);
+*/
 
-  deproject_raw ();
-  project_ref ();
-
-  FitChips (3);
-
-  deproject_raw ();
-  project_ref ();
-
-  wstars ();
-  exit (0);
 }
-
-/** bore site terms : field.coords.crval1, field.coords.crval2, field.coords.pcN_M **/
-
-
-  /* on the first pass, we only have the rough R,D for each input star */
-  /* on the second pass, we have the reference R,D for each input star */
-
-  /* pt[0].R,D here are from the chip coords */
-  /* pt[0].R,D here are USNO coords */
Index: /trunk/Ohana/src/mosastro/src/output.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/output.c	(revision 3323)
+++ /trunk/Ohana/src/mosastro/src/output.c	(revision 3323)
@@ -0,0 +1,43 @@
+# include "mosastro.h"
+
+void output (char *ext, char *phu) {
+
+  int i, N;
+  char *p;
+  char outname[256];
+  
+  for (i = 0; i < Nchip; i++) {
+    
+    /* convert chip[i].file to output name */
+    
+    p = strrchr (chip[i].file, '.');
+    N = p - chip[i].file + 1;
+    if (p == NULL) N = strlen (chip[i].file);
+
+    bzero (outname, 256);
+    strncpy (outname, chip[i].file, N);
+
+    strcat (outname, ext);
+
+    PutCoords (&chip[i].map, &chip[i].header);
+    wchip (outname, &chip[i]);
+  }
+
+  { 
+
+    Header *header;
+    FILE *f;
+    
+    field_combine ();
+
+    f = fopen (phu, "w");
+    if (f == (FILE *) NULL) {
+      fprintf (stderr, "ERROR: can't create output file %s\n", phu);
+      exit (1);
+    }
+
+    header = mkmosaic (1, 1, 0, &field.project);
+    fwrite (header[0].buffer, 1, header[0].size, f);
+    fclose (f);
+  }
+}
Index: /trunk/Ohana/src/mosastro/src/random.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/random.c	(revision 3323)
+++ /trunk/Ohana/src/mosastro/src/random.c	(revision 3323)
@@ -0,0 +1,10 @@
+# include <stdio.h>
+
+init_random () {
+
+  long A, B;
+
+  A = time(NULL);
+  for (B = 0; A == time(NULL); B++);
+  srand48(B);
+}
Index: /trunk/Ohana/src/mosastro/src/warptest.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/warptest.c	(revision 3322)
+++ /trunk/Ohana/src/mosastro/src/warptest.c	(revision 3323)
@@ -38,6 +38,6 @@
   coords.Npolyterms = 3;
   for (i = 0; i < 7; i++) {
-    coords.polyterms[0][i] = 0;
-    coords.polyterms[1][i] = 0;
+    coords.polyterms[i][0] = 0;
+    coords.polyterms[i][1] = 0;
   }
   coords.polyterms[3][0] = 1e-10;
Index: /trunk/Ohana/src/mosastro/src/wstars.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/wstars.c	(revision 3322)
+++ /trunk/Ohana/src/mosastro/src/wstars.c	(revision 3323)
@@ -1,4 +1,20 @@
 # include "mosastro.h"
 # define NCHAR 66 /* 65 char EXCLUDING return */
+
+void wchip (char *filename, Chip *data) {
+  
+  FILE *g;
+
+  g = fopen (filename, "w");
+  if (g == (FILE *) NULL) {
+    fprintf (stderr, "ERROR: can't create output file %s\n", filename);
+    exit (1);
+  }
+
+  fwrite (data[0].header.buffer, 1, data[0].header.size, g);
+  fwrite (data[0].buffer, 1, data[0].Nbuffer, g);
+  fclose (g);
+
+}
 
 void wstars (char *filename, Stars *stars, int Nstars, Header *header) {