Index: /trunk/Ohana/src/mosastro/Makefile
===================================================================
--- /trunk/Ohana/src/mosastro/Makefile	(revision 3294)
+++ /trunk/Ohana/src/mosastro/Makefile	(revision 3294)
@@ -0,0 +1,114 @@
+include ../../Configure
+HOME    =       $(ROOT)/src/mosastro.v3
+PROGRAM =       mosastro
+default: $(PROGRAM)
+
+BIN	=	$(HOME)/bin
+INC	= 	$(HOME)/include
+SRC	=	$(HOME)/src
+MAN	=	$(HOME)/doc
+DESTBIN	=	$(LBIN)
+DESTLIB	=	$(LLIB)
+DESTINC	=	$(LINC)
+DESTMAN	=	$(LMAN)
+
+#
+INCS	= 	-I$(INC) -I$(LINC) -I$(XINC)
+LIBS	= 	-L$(LLIB) -lohana -lFITS -lm 
+CFLAGS	=	$(INCS)
+LFLAGS	=	$(LIBS)
+
+MOS = \
+$(SRC)/mosastro.$(ARCH).o \
+$(SRC)/args.$(ARCH).o \
+$(SRC)/LoadStars.$(ARCH).o \
+$(SRC)/match.$(ARCH).o \
+$(SRC)/mkheader.$(ARCH).o \
+$(SRC)/fakefield.$(ARCH).o \
+$(SRC)/greference.$(ARCH).o \
+$(SRC)/project.$(ARCH).o \
+$(SRC)/transforms.$(ARCH).o \
+$(SRC)/dump.$(ARCH).o \
+$(SRC)/sort.$(ARCH).o \
+$(SRC)/ConfigInit.$(ARCH).o \
+$(SRC)/get2mass.$(ARCH).o \
+$(SRC)/getgsc.$(ARCH).o \
+$(SRC)/getptolemy.$(ARCH).o \
+$(SRC)/gptolemy.$(ARCH).o \
+$(SRC)/gregions.$(ARCH).o \
+$(SRC)/gcatalog.$(ARCH).o \
+$(SRC)/getusno.$(ARCH).o \
+$(SRC)/wstars.$(ARCH).o \
+$(SRC)/field.$(ARCH).o \
+$(SRC)/chips.$(ARCH).o \
+$(SRC)/parse_time.$(ARCH).o
+
+mosastro: $(BIN)/mosastro.$(ARCH)
+$(BIN)/mosastro.$(ARCH) : $(MOS)
+
+STD = \
+$(SRC)/mkstandards.$(ARCH).o \
+$(SRC)/mkheader.$(ARCH).o \
+$(SRC)/wstars.$(ARCH).o
+
+mkstandards: $(BIN)/mkstandards.$(ARCH)
+$(BIN)/mkstandards.$(ARCH) : $(STD)
+
+OBS = \
+$(SRC)/mkobs.$(ARCH).o \
+$(SRC)/mkheader.$(ARCH).o \
+$(SRC)/fakefield.$(ARCH).o \
+$(SRC)/greference.$(ARCH).o \
+$(SRC)/project.$(ARCH).o \
+$(SRC)/transforms.$(ARCH).o \
+$(SRC)/dump.$(ARCH).o \
+$(SRC)/ConfigInit.$(ARCH).o \
+$(SRC)/get2mass.$(ARCH).o \
+$(SRC)/getgsc.$(ARCH).o \
+$(SRC)/getptolemy.$(ARCH).o \
+$(SRC)/gptolemy.$(ARCH).o \
+$(SRC)/gregions.$(ARCH).o \
+$(SRC)/gcatalog.$(ARCH).o \
+$(SRC)/getusno.$(ARCH).o \
+$(SRC)/wstars.$(ARCH).o
+
+mkobs: $(BIN)/mkobs.$(ARCH)
+$(BIN)/mkobs.$(ARCH) : $(OBS)
+
+warptest: $(BIN)/warptest.(ARCH)
+$(BIN)/warptest.(ARCH) : $(SRC)/warptest.$(ARCH).o
+
+# dependancy rules for binary code #########################
+.PRECIOUS: %.$(ARCH).o
+.PRECIOUS: $(BIN)/%.$(ARCH)
+
+%.$(ARCH).o : %.c
+	$(CC) $(CFLAGS) -c $< -o $@
+
+$(BIN)/%.$(ARCH) : $(SRC)/%.$(ARCH).o
+	@if [ ! -d $(BIN) ]; then mkdir -p $(BIN); fi
+	$(CC) $^ -o $@ $(LFLAGS)
+
+$(DESTBIN)/%: $(BIN)/%.$(ARCH)
+	@if [ ! -d $(DESTBIN) ]; then mkdir -p $(DESTBIN); fi
+	rm -f $(DESTBIN)/$*
+	cp $(BIN)/$*.$(ARCH) $(DESTBIN)/$*
+
+$(INSTALL) $(DEVEL): % : $(BIN)/%.$(ARCH)
+
+%.clean :
+	rm -f $(BIN)/$*.$(ARCH)
+
+%.install:
+	make $(DESTBIN)/$*
+
+# utilities #################################################
+
+install:
+	for i in $(INSTALL); do make $$i.install; done
+
+clean:	
+	rm -f $(BIN)/*.$(ARCH)
+	rm -f `find . -name "*.o"`
+	rm -f `find . -name "*~"`
+	rm -f `find . -name "#*"`
Index: /trunk/Ohana/src/mosastro/doc/notes.txt
===================================================================
--- /trunk/Ohana/src/mosastro/doc/notes.txt	(revision 3294)
+++ /trunk/Ohana/src/mosastro/doc/notes.txt	(revision 3294)
@@ -0,0 +1,59 @@
+
+- need to update header astrometry 
+  - determine higher-order polynomial terms from fit?
+  - insert all distortion & warping terms into header?
+
+- better input parameters:
+  - use camera to define chips we expect to find
+  - camera -> default warp
+  - add dispersion stats to header in output (CERROR, CPRECISE)
+
+- allow chips to have full linear terms (after warp)
+
+- load 2mass, USNO SA data
+
+-----------------------------------------
+- load *.smp / *.cmp : 
+  chip[i].pt[j].x
+  chip[i].pt[j].y
+  chip[i].coords
+  (global Ro, Do = RA_DEG, DEC_DEG)
+
+- load usno in region
+  stars[i].r, stars[i].d
+
+- match usno:smp:
+  chip[i].pt[j].r, chip[i].pt[j].d
+
+- project ra,dec to p,q (arcsec on sky relative to Ro, Do)
+  chip[i].pt[j].p, chip[i].pt[j].q
+
+- init chip paramters:
+  chip[i].Px, chip[i].Py - parities
+  chip[i].Po, chip[i].Qo, chip[i].To
+
+  * p' = Px*(p - Po), q' = Py*(q - Qo)
+  * x = cos(To) p' + sin(To) q'
+  * y = cos(To) p' + sin(To) q'
+
+- other parameter:
+  * field.Ro, field.Do - not fitted
+  field.A[0] - field.A[5]
+  field.D[0] - field.D[5]
+
+
+  (r,d) -> field.coords -> (p,q) [arcsec from field center]
+  (p,q) -> field.A,D    -> (L,M) [pixels from field center]
+  (L,M) -> Px,Py,Po,Qo  -> (P,Q) [pixels from chip corner]
+  (P,Q) -> To		-> (x,y) [pixels on chip]
+
+  L = A[0] + A[1]*p + A[2]*q + A[3]*p^2 + A[4]*q^2 + A[5]*p*q
+  M = D[0] + D[1]*p + D[2]*q + D[3]*p^2 + D[4]*q^2 + D[5]*p*q
+
+  p' = Px*(L - Po)
+  q' = Py*(M - Qo)
+
+  x = cos(To) p' + sin(To) q'
+  y = cos(To) p' + sin(To) q'
+
+  
Index: /trunk/Ohana/src/mosastro/include/mosastro.h
===================================================================
--- /trunk/Ohana/src/mosastro/include/mosastro.h	(revision 3294)
+++ /trunk/Ohana/src/mosastro/include/mosastro.h	(revision 3294)
@@ -0,0 +1,104 @@
+# include <ohana.h>
+# include <loneos.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 type;
+} StarData;
+
+typedef struct {
+  double RA[2], DEC[2];
+  double Area;
+  char *name;
+} CatStats;
+
+typedef struct {
+  char *file;
+  int Nstars;
+  StarData *stars;
+
+  int Nmatch;
+  StarData *raw, *ref;
+
+  Coords coords;
+  Coords map;
+} Chip;
+
+typedef struct {
+  double Ro, Do, PSx, PSy, To;
+  double Rmin, Rmax, Dmin, Dmax;
+  int    Norder;
+  double **A;
+  double **D;
+  Coords project;
+  Coords distort;
+  int    fit;
+} Field;
+
+typedef struct {
+  double *dPdL;
+  double *dPdM;
+  double *dQdL;
+  double *dQdM;
+  double *Lo;
+  double *Mo;
+  int Npts;
+} Gradients;
+
+typedef struct {
+  double X;
+  double Y;
+  double M, dM;
+  char   dophot;
+  double sky;
+  double fx, fy, df;
+  double Mgal, Map;
+} Stars;
+
+int  Nchip;
+Chip *chip;
+Field field;
+double Year;  /** carried for precession - probably put this in chip data **/
+double RADIUS; /** raw / ref matching radius (units?) **/
+
+char CDROM[256];
+char REFCAT[256];
+char CATDIR[256];
+char DateKeyword[256];
+char DateMode[256];
+char UTKeyword[256];
+char MJDKeyword[256];
+char JDKeyword[256];
+char TWO_MASS_DIR[256];
+int VERBOSE;
+
+char GSCFILE[256];
+char GSC_DIR[256];
+
+StarData *getusno (CatStats *catstats, int *Nstars);
+StarData *getgsc (CatStats *catstats, int *NSTARS);
+StarData *get2mass (CatStats *catstats, int *NSTARS);
+StarData *getptolemy (CatStats *catstats, int *NSTARS);
+void GetConfig (char *config, char *field, char *format, int N, void *ptr);
+StarData *greference (int *Nrefcat);
+
+Header *mkheader (int Nx, int Ny, int Nstars, Coords *coords);
+
+StarData *gcatalog (char *filename, int *Nstars);
+
+CatStats *gregions (CatStats *patch, int *nregion);
+int find_dec_bands (CatStats *area);
+int load_ra_blocks (int Ndec, CatStats *area);
+int parse_GSC_line (CatStats *tregion, char *line);
+void init_regions ();
+void add_to_regions (CatStats *area);
+void area_of_region (CatStats *region);
+void set_catalog (char *catdir);
+
+StarData *gptolemy (char *filename, int *NSTARS);
+int fake_field (double RA, double DEC);
+Gradients *GetGradients ();
Index: /trunk/Ohana/src/mosastro/src/ConfigInit.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/ConfigInit.c	(revision 3294)
+++ /trunk/Ohana/src/mosastro/src/ConfigInit.c	(revision 3294)
@@ -0,0 +1,50 @@
+# include "mosastro.h"
+
+void ConfigInit (int *argc, char **argv) {
+  
+  char *config, *file;
+
+  VERBOSE = TRUE;
+
+  /*** load configuration info ***/
+  file = SelectConfigFile (argc, argv, "ptolemy");
+  config = LoadConfigFile (file);
+  if (config == (char *) NULL) {
+    fprintf (stderr, "ERROR: can't find configuration file %s\n", file);
+    if (file != (char *) NULL) free (file);
+    exit (0);
+  }
+  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, "ASTRO_REFCAT",      "%s",  0, REFCAT);
+  ScanConfig (config, "CATDIR",            "%s",  0, CATDIR);
+  ScanConfig (config, "GSCFILE",           "%s",  0, GSCFILE);
+  ScanConfig (config, "GSCDIR",            "%s",  0, GSC_DIR);
+
+  GetConfig (config, "RADIUS",            "%lf", 0, &RADIUS);
+
+  free (config);
+  free (file);
+
+  return;
+}
+
+void GetConfig (char *config, char *field, char *format, int N, void *ptr) {
+
+  char *status;
+
+  status = ScanConfig (config, field, format, N, ptr);
+  if (status == NULL) {
+    fprintf (stderr, "error in config, cannot find %s\n", field);
+    exit (1);
+  }
+  return;
+}
Index: /trunk/Ohana/src/mosastro/src/LoadStars.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/LoadStars.c	(revision 3294)
+++ /trunk/Ohana/src/mosastro/src/LoadStars.c	(revision 3294)
@@ -0,0 +1,91 @@
+# include "mosastro.h"
+# define BYTES_STAR 66
+
+int LoadStars (int Nfile, char **file) {
+
+  Header header;
+  int i, j, NCHIP, itmp, Nbytes, nbytes;
+  FILE *f;
+  char *buffer;
+  time_t tsval;
+  struct tm *tmval;
+
+  ALLOCATE (chip, Chip, Nfile);
+
+  /* Nchip is number of valid chips */
+
+  Nchip = 0;
+  for (i = 0; i < Nfile; i++) {
+    chip[Nchip].file = strcreate (file[i]);
+
+    /* load header */
+    if (!fits_read_header (chip[Nchip].file, &header)) {
+      fprintf (stderr, "ERROR: can't read header for %s\n", file[i]);
+      free (chip[Nchip].file);
+      continue;
+    }
+    f = fopen (chip[Nchip].file, "r");
+    if (f == NULL) {
+      fprintf (stderr, "ERROR: can't read data from %s\n", file[i]);
+      free (chip[Nchip].file);
+      continue;
+    }
+    fseek (f, header.size, SEEK_SET); 
+
+    /* get astrometry information */
+    if (!GetCoords (&chip[Nchip].coords, &header)) {
+      fprintf (stderr, "ERROR: no astrometric solution in header\n");
+      free (chip[Nchip].file);
+      continue;
+    }
+    while (chip[Nchip].coords.crval1 < 0) chip[Nchip].coords.crval1 += 360.0;
+    while (chip[Nchip].coords.crval1 > 360.0) chip[Nchip].coords.crval1 -= 360.0;
+
+    itmp = 0;
+    fits_scan (&header, "NASTRO",   "%d", 1, &itmp);
+    if (itmp == 0) {
+      fprintf (stderr, "ERROR: bad astrometric solution in header %s\n", file[i]);
+      free (chip[Nchip].file);
+      continue;
+    }
+
+    /* get time info */
+    tsval = parse_time (&header);
+    tmval = gmtime (&tsval);
+    Year = tmval[0].tm_year;
+
+    /* find expected number of stars */
+    fits_scan (&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]);
+      free (chip[Nchip].file);
+      continue;
+    }
+    ALLOCATE (chip[Nchip].stars, StarData, chip[Nchip].Nstars);
+    
+    Nbytes = BYTES_STAR * chip[Nchip].Nstars;
+    ALLOCATE (buffer, char, Nbytes + 1);
+    nbytes = Fread (buffer, 1, Nbytes, f, "char");
+    if (nbytes != Nbytes) { exit (1); }
+
+    /** need to save all input lines so we can write on output? (gastro does not) **/
+    for (j = 0; j < chip[Nchip].Nstars; j++) {
+      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]);
+    }
+    /** 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);
+
+    Nchip ++;
+  }
+ 
+  if (Nchip == 0) {
+    fprintf (stderr, "no valid data in chip files, exiting\n");
+    exit (1);
+  }
+
+  return (TRUE);
+}
Index: /trunk/Ohana/src/mosastro/src/args.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/args.c	(revision 3294)
+++ /trunk/Ohana/src/mosastro/src/args.c	(revision 3294)
@@ -0,0 +1,61 @@
+# include "mosastro.h"
+
+void print_help () {
+
+  fprintf (stderr, "mosastro -- mosaic astrometry\n");
+  fprintf (stderr, "\n"); 
+  exit (0);
+
+}
+
+void args (int *argc, char **argv) {
+  
+  int N;
+  char line[500];
+
+  if (get_argument (*argc, argv, "-help") ||
+      get_argument (*argc, argv, "-h")) {
+    print_help ();
+  }
+
+  VERBOSE = FALSE;
+  if ((N = get_argument (*argc, argv, "-v"))) {
+    VERBOSE = TRUE;
+    remove_argument (N, argc, argv);
+  }
+
+  /*
+  DUMP_REF = FALSE;
+  if ((N = get_argument (*argc, argv, "-dumpref"))) {
+    DUMP_REF = TRUE;
+    remove_argument (N, argc, argv);
+  }
+
+  DUMP_RAW = FALSE;
+  if ((N = get_argument (*argc, argv, "-dumpraw"))) {
+    DUMP_RAW = TRUE;
+    remove_argument (N, argc, argv);
+  }
+
+  CHIPS = (char *) NULL;
+  if ((N = get_argument (*argc, argv, "-chips"))) {
+    remove_argument (N, argc, argv);
+    CHIPS = strcreate (argv[N]);
+    remove_argument (N, argc, argv);
+  }
+
+  FIELD = (char *) NULL;
+  if ((N = get_argument (*argc, argv, "-field"))) {
+    remove_argument (N, argc, argv);
+    FIELD = strcreate (argv[N]);
+    remove_argument (N, argc, argv);
+  }
+  */
+
+  field.Norder = 0;
+  if ((N = get_argument (*argc, argv, "-order"))) {
+    remove_argument (N, argc, argv);
+    field.Norder = atoi (argv[N]);
+    remove_argument (N, argc, argv);
+  }
+}
Index: /trunk/Ohana/src/mosastro/src/chips.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/chips.c	(revision 3294)
+++ /trunk/Ohana/src/mosastro/src/chips.c	(revision 3294)
@@ -0,0 +1,45 @@
+# include "mosastro.h"
+
+/* set chip model based on initial header info */
+int init_chips () {
+
+  int i, j;
+  double R, D, P, Q, L, M;
+
+# 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");
+    chip[i].map.crpix1 = 0.0;
+    chip[i].map.crpix2 = 0.0;
+    chip[i].map.cdelt1 = 1.0;
+    chip[i].map.cdelt2 = 1.0;
+
+    /* find (L,M) coords of reference pixel (0,0) */
+    XY_to_RD (&R, &D, 0.0, 0.0, &chip[i].coords);
+    RD_to_XY (&P, &Q, R, D, &field.project);
+    RD_to_XY (&L, &M, P, Q, &field.distort);
+    chip[i].map.crval1 = L;
+    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;
+
+    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;
+    }
+  }
+  return (1);
+}
Index: /trunk/Ohana/src/mosastro/src/dump.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/dump.c	(revision 3294)
+++ /trunk/Ohana/src/mosastro/src/dump.c	(revision 3294)
@@ -0,0 +1,17 @@
+# include "mosastro.h"
+
+int dump_stars (FILE *f, StarData *stars, int Nstars) {
+
+  int i;
+
+  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", 
+	     i, 
+	     stars[i].R, stars[i].D,
+	     stars[i].P, stars[i].Q,
+	     stars[i].L, stars[i].M,
+	     stars[i].X, stars[i].Y,
+	     stars[i].Mag, stars[i].dMag);
+  }
+  return (1);
+}
Index: /trunk/Ohana/src/mosastro/src/fakefield.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/fakefield.c	(revision 3294)
+++ /trunk/Ohana/src/mosastro/src/fakefield.c	(revision 3294)
@@ -0,0 +1,59 @@
+# include "mosastro.h"
+
+/* determine an initial guess to field parameters from data */ 
+int fake_field (double RA, double DEC) {
+
+  int i;
+
+  field.Rmin = RA  - 0.5;
+  field.Rmax = RA  + 0.5;
+  field.Dmin = DEC - 0.5;
+  field.Dmax = DEC + 0.5;
+
+  /* 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);
+  field.project.crpix1 = 0;
+  field.project.crpix2 = 0;
+  
+  /* measure average plate scale - would be better using parabolic min... */
+  field.project.cdelt1 = 1.0/3600.0;
+  field.project.cdelt2 = 1.0/3600.0;
+
+  /** allow guess at field rotation?? **/
+  field.project.pc1_1  = 1;
+  field.project.pc2_2  = 1;
+  field.project.pc1_2  = 0;
+  field.project.pc2_1  = 0;
+  field.project.Npolyterms = 1;
+
+  /* bore site center guess */
+  strcpy (field.distort.ctype, "RA---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 = 3;
+  for (i = 0; i < 7; i++) {
+    field.distort.polyterms[0][i] = 0;
+    field.distort.polyterms[1][i] = 0;
+  }
+
+  field.distort.polyterms[0][3] = 1e-10;
+  field.distort.polyterms[1][6] = 1e-10;
+
+  return (1);
+}
+
Index: /trunk/Ohana/src/mosastro/src/field.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/field.c	(revision 3294)
+++ /trunk/Ohana/src/mosastro/src/field.c	(revision 3294)
@@ -0,0 +1,88 @@
+# include "mosastro.h"
+
+/* determine an initial guess to field parameters from data */ 
+int init_field (int Norder) {
+
+  int i, Nterm;
+  double PS;
+
+# if (0)
+  if (FIELD != (char *) NULL) {
+    load_field (FIELD);
+    return (1);
+  }
+# endif
+
+  /* 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);
+  field.project.crpix1 = 0;
+  field.project.crpix2 = 0;
+  
+  /* measure average plate scale - would be better using parabolic min... */
+  field.project.cdelt1 = 0;
+  field.project.cdelt2 = 0;
+  for (i = 0; i < Nchip; i++) {
+    field.project.cdelt1 += chip[i].coords.cdelt1;
+    field.project.cdelt2 += chip[i].coords.cdelt2;
+  }
+  field.project.cdelt1 /= Nchip;
+  field.project.cdelt2 /= Nchip;
+  /* force starting guess to have equal x & y plate scales? */
+
+  /** allow guess at field rotation?? **/
+  field.project.pc1_1  = 1;
+  field.project.pc2_2  = 1;
+  field.project.pc1_2  = 0;
+  field.project.pc2_1  = 0;
+  field.project.Npolyterms = 1;
+
+  /* bore site center guess */
+  strcpy (field.distort.ctype, "RA---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 = 1;
+  for (i = 0; i < 7; i++) {
+    field.distort.polyterms[0][i] = 0;
+    field.distort.polyterms[1][i] = 0;
+  }
+  return (1);
+}
+
+void field_stats () {
+  
+  int i, j;
+  double Rmin, Rmax, Dmin, Dmax;
+
+  /* find range for single image */
+  Rmin = Dmin = 360.0;
+  Rmax = Dmax = -90.0;
+
+  for (i = 0; i < Nchip; i++) {
+    for (j = 0; j < chip[i].Nstars; j++) {
+      Rmin = MIN (Rmin, chip[i].stars[j].R);
+      Dmin = MIN (Dmin, chip[i].stars[j].D);
+      Rmax = MAX (Rmax, chip[i].stars[j].R);
+      Dmax = MAX (Dmax, chip[i].stars[j].D);
+    }
+  }
+  field.Rmin = Rmin;
+  field.Rmax = Rmax;
+  field.Dmin = Dmin;
+  field.Dmax = Dmax;
+
+}
Index: /trunk/Ohana/src/mosastro/src/gcatalog.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/gcatalog.c	(revision 3294)
+++ /trunk/Ohana/src/mosastro/src/gcatalog.c	(revision 3294)
@@ -0,0 +1,47 @@
+# include "mosastro.h"
+
+# define BYTES_STAR 23
+# define BLOCK 1000
+
+StarData *gcatalog (char *filename, int *Nstars) {
+  
+  StarData *stars;
+  int i, NSTAR, nstar, Nbytes, nbytes;
+  char *buffer;
+  FILE *f;
+
+  f = fopen (filename, "r");
+  if (f == NULL) {
+    fprintf (stderr, "ERROR: can't find catalog file %s\n", filename);
+    return (NULL);
+  }
+  
+  nstar = 0;
+  NSTAR = 1000;
+  ALLOCATE (stars, StarData, NSTAR);
+
+  Nbytes = BLOCK*BYTES_STAR;
+  ALLOCATE (buffer, char, Nbytes);
+  while ((nbytes = fread (buffer, 1, Nbytes, f)) > 0) {
+    for (i = 0; i < nbytes / BYTES_STAR; i++) {
+      dparse (&stars[nstar].R,   1, &buffer[i*BYTES_STAR]);
+      dparse (&stars[nstar].D,   2, &buffer[i*BYTES_STAR]);
+      dparse (&stars[nstar].Mag, 3, &buffer[i*BYTES_STAR]);
+      nstar++;
+      if (nstar == NSTAR) {
+	NSTAR += 1000;
+	REALLOCATE (stars, StarData, NSTAR);
+      }
+    }
+  }
+  free (buffer);
+
+  *Nstars = nstar;
+  return (stars);
+}
+
+/** this function reads a file in the format of my extracted GSC files
+    this is a fairly weak way of storing photometry / astrometry data 
+    eventually, this should be moved to ptolemy / DVO or something else
+    or else, the existing 2mass and HST-GSC files should be parsed
+**/
Index: /trunk/Ohana/src/mosastro/src/get2mass.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/get2mass.c	(revision 3294)
+++ /trunk/Ohana/src/mosastro/src/get2mass.c	(revision 3294)
@@ -0,0 +1,35 @@
+# include "mosastro.h"
+
+StarData *get2mass (CatStats *catstats, int *NSTARS) {
+  
+  int i, j, k, Ns, Ngsc, Nregions, Nstars; 
+  StarData *gsc;
+  StarData *stars;
+  CatStats *regions;
+
+  Nstars = 0;
+  ALLOCATE (stars, StarData, 1);
+
+  set_catalog (TWO_MASS_DIR);
+  regions = gregions (catstats, &Nregions);
+  
+  for (i = 0; i < Nregions; i++) {
+    gsc = gcatalog (regions[i].name, &Ngsc);
+    if (gsc == NULL) continue;
+
+    Ns = Nstars;
+    Nstars += Ngsc;
+
+    REALLOCATE (stars, StarData, MAX (1, Nstars));
+    for (k = Ns, j = 0; j < Ngsc; k++, j++) {
+      stars[k].R = gsc[j].R;
+      stars[k].D = gsc[j].D;
+      stars[k].M = gsc[j].M;
+    }      
+    free (gsc);
+  }
+  
+  if (VERBOSE) fprintf (stderr, "%d stars from 2MASS\n", Nstars);
+  *NSTARS = Nstars;
+  return (stars);
+}  
Index: /trunk/Ohana/src/mosastro/src/getgsc.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/getgsc.c	(revision 3294)
+++ /trunk/Ohana/src/mosastro/src/getgsc.c	(revision 3294)
@@ -0,0 +1,36 @@
+# include "mosastro.h"
+
+StarData *getgsc (CatStats *catstats, int *NSTARS) {
+  
+  int i, j, k, Ns, Ngsc, Nregions, Nstars; 
+  StarData *gsc;
+  StarData *stars;
+  CatStats *regions;
+
+  Nstars = 0;
+  ALLOCATE (stars, StarData, 1);
+
+  set_catalog (GSC_DIR);
+  regions = gregions (catstats, &Nregions);
+  
+  for (i = 0; i < Nregions; i++) {
+    gsc = gcatalog (regions[i].name, &Ngsc);
+    if (gsc == NULL) continue;
+
+    Ns = Nstars;
+    Nstars += Ngsc;
+
+    REALLOCATE (stars, StarData, MAX (1, Nstars));
+    for (k = Ns, j = 0; j < Ngsc; k++, j++) {
+      stars[k].R   = gsc[j].R;
+      stars[k].D   = gsc[j].D;
+      stars[k].Mag = gsc[j].Mag;
+    }      
+    free (gsc);
+  }
+  
+  if (VERBOSE) fprintf (stderr, "%d stars from HST GSC\n", Nstars);
+  *NSTARS = Nstars;
+  return (stars);
+}  
+
Index: /trunk/Ohana/src/mosastro/src/getptolemy.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/getptolemy.c	(revision 3294)
+++ /trunk/Ohana/src/mosastro/src/getptolemy.c	(revision 3294)
@@ -0,0 +1,32 @@
+# include "mosastro.h"
+
+StarData *getptolemy (CatStats *catstats, int *NSTARS) {
+  
+  int i, j, k, Ns, Nptolemy, Nregions, Nstars; 
+  CatStats *regions;
+  StarData *stars, *ptolemy;
+
+  Nstars = 0;
+  ALLOCATE (stars, StarData, 1);
+
+  set_catalog (CATDIR);
+  regions = gregions (catstats, &Nregions);
+  
+  for (i = 0; i < Nregions; i++) {
+    ptolemy = gptolemy (regions[i].name, &Nptolemy);
+    if (ptolemy == NULL) continue;
+
+    Ns = Nstars;
+    Nstars += Nptolemy;
+
+    REALLOCATE (stars, StarData, MAX (1, Nstars));
+    for (k = Ns, j = 0; j < Nptolemy; k++, j++) {
+      stars[k] = ptolemy[j];
+    }
+    free (ptolemy);
+  }
+
+  if (VERBOSE) fprintf (stderr, "%d stars from PTOLEMY\n", Nstars);
+  *NSTARS = Nstars;
+  return (stars);
+}  
Index: /trunk/Ohana/src/mosastro/src/getusno.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/getusno.c	(revision 3294)
+++ /trunk/Ohana/src/mosastro/src/getusno.c	(revision 3294)
@@ -0,0 +1,136 @@
+# include "mosastro.h"
+# define NZONE 24
+
+int SPDzone[] = {
+  0, 75, 450, 375, 1500, 1650, 300, 1425, 1725, 525, 1275, 225, 
+  675, 150, 600, 1575, 750, 975, 900, 1050, 1125, 1200, 825, 1350};
+
+int USNOdisk[] = {
+  1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4, 5, 5, 6, 6, 6, 7, 7, 8, 8, 9, 9, 10, 10};
+
+StarData *getusno (CatStats *catstats, int *Nstars) {
+
+  long int offset;
+  int i, bin, first, last, nitems, Nitems, Nbins;
+  float hours[100];
+  int start[100], number[100], *buffer, *buf;
+  char filename[128], c;
+  FILE *f;
+  int iRA0, iRA1, iDEC0, iDEC1;
+  double RA0, RA1, DEC0, DEC1, dec;
+  int spd, spd_start, spd_end, disk;
+  int NUSNO, Nusno;
+  StarData *stars;
+
+  RA0  = catstats[0].RA[0]; 
+  DEC0 = catstats[0].RA[1]; 
+  RA1  = catstats[0].DEC[0];
+  DEC1 = catstats[0].DEC[1];
+
+  /* identify ra & dec range of interest */
+  iRA0 = RA0 * 360000.0;
+  iRA1 = RA1 * 360000.0;
+  iDEC0 = (DEC0 + 90.0) * 360000.0;
+  iDEC1 = (DEC1 + 90.0) * 360000.0;
+  
+  /* data is organized in south-pole distance zones */
+  spd_start = (int)((DEC0 + 90) / 7.5) * 75.0;
+  dec = (DEC1 + 90) / 7.5;
+  if (dec > (int)(dec)) {
+    spd_end =   (int)(1 + (DEC1 + 90) / 7.5) * 75.0;
+  } else {
+    spd_end =   (int)(0 + (DEC1 + 90) / 7.5) * 75.0;
+  }
+
+  Nusno = 0;
+  NUSNO = 5000;
+  ALLOCATE (stars, StarData, NUSNO);
+
+  for (spd = spd_start; spd < spd_end; spd += 75) {
+    disk = -1;
+    for (i = 0; i < NZONE; i++) {
+      if (spd == SPDzone[i]) 
+	disk = USNOdisk[i];
+    }
+    if (disk < 0) {
+      fprintf (stderr, "ERROR: can't find cdrom for spd %d\n",  spd);
+      exit (0);
+    }
+    
+    /* load accelerator file */
+    sprintf (filename, "%s/zone%04d.acc", CDROM, spd); 
+    if (VERBOSE) fprintf (stderr, "reading from %s\n", filename);
+    f = fopen (filename, "r");
+    if (f == (FILE *) NULL) {
+      fprintf (stderr, "can't open file %s, is cdrom %d in drive?\n", filename, disk);
+      fprintf (stderr, "press return when ready to continue: ");
+      fscanf (stdin, "%c", &c);
+      fprintf (stderr, "trying again...\n");
+      f = fopen (filename, "r");
+      if (f == (FILE *) NULL) {
+	fprintf (stderr, "ERROR: can't open file %s, is cdrom %d in drive?\n", filename, disk);
+	exit (1);  
+      }
+    }
+    for (i = 0; fscanf (f, "%f %d %d", &hours[i], &start[i], &number[i]) != EOF; i++);
+    Nbins = i;
+    fclose (f);
+    
+    first = RA0 / 3.75;
+    if ((RA1 / 3.75) == (int) (RA1 / 3.75)) 
+      last  = RA1 / 3.75;
+    else 
+      last  = 1 + RA1 / 3.75;
+
+    if ((first > Nbins) || (last > Nbins)) {
+      fprintf (stderr, "ERROR: RA out of range\n");
+      exit (1);
+    }
+    
+    /* open data file */
+    sprintf (filename, "%s/zone%04d.cat", CDROM, spd);
+    if (VERBOSE) fprintf (stderr, "reading from %s\n", filename);
+    f = fopen (filename, "r");
+    if (f == (FILE *) NULL) {
+      fprintf (stderr, "ERROR: can't open file %s\n", filename);
+      exit (1);
+    }
+    /* advance file pointer to first slice */
+    offset = 3*sizeof(int)*(start[first] - 1);
+    fseek (f, offset, SEEK_SET);
+    /* on each loop, load data from an RA slice of the catalog */
+    for (bin = first; bin < last; bin++) {
+      Nitems = 3*number[bin];
+      ALLOCATE (buffer, int, Nitems);
+      nitems = Fread (buffer, sizeof(int), Nitems, f, "int");
+      if (nitems != Nitems) {
+	fprintf (stderr, "ERROR: failure reading data from file %s\n", filename);
+	exit (1);
+      }
+      buf = buffer;
+      /* print out data from slice within RA and DEC range */
+      for (i = 0; i < number[bin]; i++, buf+=3) {
+	if ((buf[0] > iRA0) && (buf[0] < iRA1) &&
+	    (buf[1] > iDEC0) && (buf[1] < iDEC1)) {
+	  stars[Nusno].R = buf[0]/360000.0;
+	  stars[Nusno].D = buf[1]/360000.0 - 90.0;
+	  stars[Nusno].Mag = fabs (0.1*(buf[2] - 1000*((int)(buf[2]/1000))));
+	  /* b = 0.1*((int)(buf[2] - 1000000*((int)(buf[2]/1000000))) / 1000); */
+	  Nusno ++;
+	  if (Nusno == NUSNO) {
+	    NUSNO += 5000;
+	    REALLOCATE (stars, StarData, NUSNO);
+	  }	  
+	}
+      }
+      free (buffer);
+    }
+    fclose (f);
+  }
+
+  *Nstars = Nusno;
+  if (VERBOSE) fprintf (stderr, "%d stars from USNO 1.0\n", Nusno);
+  return (stars);
+}
+
+
Index: /trunk/Ohana/src/mosastro/src/gptolemy.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/gptolemy.c	(revision 3294)
+++ /trunk/Ohana/src/mosastro/src/gptolemy.c	(revision 3294)
@@ -0,0 +1,40 @@
+# include "mosastro.h"
+
+StarData *gptolemy (char *filename, int *NSTARS) {
+  
+  int i, Nstars;
+  Catalog catalog;
+  StarData *stars;
+
+  Nstars = 0;
+  catalog.filename = filename;
+  switch (lock_catalog (&catalog, LCK_SOFT)) {
+  case 0:
+  case 2:
+    fprintf (stderr, "can't lock catalog data %s\n", filename);
+    *NSTARS = Nstars;
+    return (NULL);
+  case 1:
+    break;
+  }
+  if (!load_catalog (&catalog, LOAD_AVES | LOAD_MEAS, FALSE)) {
+    fprintf (stderr, "can't load catalog data %s\n", filename);
+    *NSTARS = Nstars;
+    return (NULL);
+  }
+  unlock_catalog (&catalog);
+
+  Nstars = catalog.Naverage;
+  ALLOCATE (stars, StarData, MAX (1, Nstars));
+  for (i = 0; i < catalog.Naverage; i++) {
+    bzero (&stars[i], sizeof(StarData));
+    stars[i].R 	 = catalog.average[i].R;
+    stars[i].D 	 = catalog.average[i].D;
+    stars[i].Mag = 0.001*catalog.measure[catalog.average[i].offset].M;
+  }
+  free (catalog.average);
+  free (catalog.measure);
+
+  *NSTARS = Nstars;
+  return (stars);
+}
Index: /trunk/Ohana/src/mosastro/src/greference.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/greference.c	(revision 3294)
+++ /trunk/Ohana/src/mosastro/src/greference.c	(revision 3294)
@@ -0,0 +1,46 @@
+# include "mosastro.h"
+
+StarData *greference (int *Nrefcat) {
+
+  int Nstars;
+  StarData *stars;
+  CatStats catstats;
+
+  if (VERBOSE) fprintf (stderr, "loading astrometric reference data from %s\n", REFCAT); 
+
+  Nstars = 0;
+  catstats.RA[0]  = field.Rmin;
+  catstats.RA[1]  = field.Rmax;
+  catstats.DEC[0] = field.Dmin;
+  catstats.DEC[1] = field.Dmax;
+
+  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 USNO catalog for the given region */
+  if (!strcmp (REFCAT, "USNO")) {
+    stars = getusno (&catstats, &Nstars);
+  }
+
+  /* get stars from the HST GSC catalog for the given region */
+  if (!strcmp (REFCAT, "GSC")) {
+    stars = getgsc (&catstats, &Nstars);
+  }
+  
+  /* get stars from the HST GSC catalog for the given region */
+  if (!strcmp (REFCAT, "2MASS")) {
+    stars = get2mass (&catstats, &Nstars);
+  }
+  
+  /* get stars from the HST GSC catalog for the given region */
+  if (!strcmp (REFCAT, "PTOLEMY")) {
+    stars = getptolemy (&catstats, &Nstars);
+  }
+  
+  if (Nstars == 0) {
+    fprintf (stderr, "no ref objs: %s\n", REFCAT);
+    exit (1);
+  }
+  *Nrefcat = Nstars;
+  return (stars);
+
+}
Index: /trunk/Ohana/src/mosastro/src/gregions.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/gregions.c	(revision 3294)
+++ /trunk/Ohana/src/mosastro/src/gregions.c	(revision 3294)
@@ -0,0 +1,255 @@
+# include "mosastro.h"
+# define MAX_NAME 256
+
+static int Nregion, NREGION;
+static CatStats *region;
+static char *CatDir;
+
+/* data which defines the layout of the regions */
+# define NDECBANDS 12
+double DecBand[] = {0.0, +7.5, +15.0, +22.5, +30.0, +37.5, +45.0, +52.5, +60.0, +67.5, +75.0, +82.5, +90.0,
+		    0.0, -7.5, -15.0, -22.5, -30.0, -37.5, -45.0, -52.5, -60.0, -67.5, -75.0, -82.5, -90.0};
+
+int NDecLines[] =  {593, 584, 551, 530, 522, 465, 406, 362, 280, 198, 123, 24, 0, 
+		    597, 578, 574, 577, 534, 499, 442, 376, 294, 212, 144, 48, 0};
+
+int NRaBlocks [] = {48, 47, 45, 43, 40, 36, 32, 28, 21, 15, 9, 3, 0, 
+		    48, 47, 45, 43, 40, 36, 32, 28, 21, 15, 9, 3, 0};
+
+char *Dec2Sections[] = {"n0000", "n0730", "n1500", "n2230", "n3000", "n3730", "n4500", 
+			"n5230", "n6000", "n6730", "n7500", "n8230", "weirdness", 
+			"s0000", "s0730", "s1500", "s2230", "s3000", "s3730", "s4500", 
+			"s5230", "s6000", "s6730", "s7500", "s8230", "weirdness"};
+
+/* given patch on sky with RA[0] - RA[1] & DEC[0] - DEC[1], return all that fall in range */
+CatStats *gregions (CatStats *patch, int *nregion) {
+  
+  CatStats area;
+  double dtmp;
+  char file[MAX_NAME];
+
+  init_regions ();
+
+  /* force RA[i] to be in range 0 - 360 , DEC[0] < DEC[1] */
+  while (patch[0].RA[0] > 360.0) { patch[0].RA[0] -= 360.0; }
+  while (patch[0].RA[0] < 0.0)   { patch[0].RA[0] += 360.0; }
+  while (patch[0].RA[1] > 360.0) { patch[0].RA[1] -= 360.0; }
+  while (patch[0].RA[1] < 0.0)   { patch[0].RA[1] += 360.0; }
+  if (patch[0].DEC[0] > patch[0].DEC[1]) {
+    dtmp = patch[0].DEC[0];
+    patch[0].DEC[0] = patch[0].DEC[1];
+    patch[0].DEC[1] = dtmp;
+  }
+
+  /* include the pole if dec > 86.25 is included (ADD SOUTH POLE!)  */
+  if (patch[0].DEC[1] > 86.25) {
+    area.RA[0]  =   0.00;
+    area.RA[1]  = 360.00;
+    area.DEC[0] =  86.25;
+    area.DEC[1] =  90.00;
+    sprintf (file, "n8230/pole");
+    area.name = strcreate (file);
+    add_to_regions (&area);
+    free (area.name);
+  }
+  
+  /* if RA[0] > RA[1], split in two passes */
+  if (patch[0].RA[0] > patch[0].RA[1]) {
+    area = patch[0];
+    area.RA[0]  = 0.0;
+    find_dec_bands (&area);
+
+    area = patch[0];
+    area.RA[1]  = 360.0;
+    find_dec_bands (&area);
+  } else {
+    area = patch[0];
+    find_dec_bands (&area);
+  }
+
+  *nregion = Nregion;
+  return (region);
+}
+
+/* given an area [RA[0] - RA[1], DEC[0] - DEC[1],
+   where 0.0 < RA[0] < RA[1] < 360.0 
+   and -90.0 < (DEC[0], DEC[1]) < 90.0 
+   and DEC[0] < DEC[1]
+   search through the dec bands to identify ones which overlap the area */
+
+int find_dec_bands (CatStats *area) {
+
+  int i;
+
+  /* for NORTH region (DecBand[i] < DecBand[i+1]) */
+  for (i = 0; i < NDECBANDS; i++) {
+    if ((area[0].DEC[1] > DecBand[i]) || (area[0].DEC[0] < DecBand[i+1])) {
+      load_ra_blocks (i, area);
+    }
+  }
+
+  /* for SOUTH region (DecBand[i] > DecBand[i+1]) */
+  for (i = NDECBANDS+1; i < 2*NDECBANDS; i++) {
+    if ((area[0].DEC[0] < DecBand[i]) || (area[0].DEC[1] > DecBand[i+1])) {
+      load_ra_blocks (i, area);
+    }
+  }
+  return (TRUE);
+}
+
+/* given the number of a DecBlock and a defined area, return all area names in that block */
+/* area is defined so that RA[0] < RA[1], DEC[0] < DEC[1] */
+int load_ra_blocks (int Ndec, CatStats *area) {
+
+  CatStats tregion;
+  int i, Nlines, Nskip;
+  char *buffer, file[MAX_NAME];
+  FILE *f;
+
+  /* count lines before Dec Band, and in Dec Band */
+  Nskip = 0;
+  for (i = 0; i < Ndec; i++) Nskip += NDecLines[i];
+  Nlines = NDecLines[Ndec];
+
+  ALLOCATE (buffer, char, 48*Nlines);
+  ALLOCATE (tregion.name, char, MAX_NAME);
+
+  /* get data from GSC file */ 
+  f = fopen (GSCFILE, "r");
+  if (f == NULL) {
+    fprintf (stderr, "ERROR: can't find GSC regions file %s\n", GSCFILE);
+    exit (1);
+  }
+  fseek (f, 5*2880 + 48*Nskip, SEEK_SET);
+  fread (buffer, 48*Nlines, 1, f);
+  fclose (f);
+
+  /* find overlapping regions */
+  for (i = 0; i < Nlines; i++) {
+    parse_GSC_line (&tregion, &buffer[i*48]);
+    if (area[0].RA[0] > tregion.RA[1]) continue;
+    if (area[0].RA[1] < tregion.RA[0]) continue;
+    if (area[0].DEC[0] > tregion.DEC[1]) continue;
+    if (area[0].DEC[1] < tregion.DEC[0]) continue;
+    sprintf (file, "%s/%s", Dec2Sections[Ndec], tregion.name);
+    strcpy (tregion.name, file);
+    add_to_regions (&tregion);
+  }
+
+  free (buffer);
+  free (tregion.name);
+  return (TRUE);
+}
+
+/* parse the GSC line, placing values in 'region'. force RA[0] < RA[1], DEC[0] < DEC[1] */
+int parse_GSC_line (CatStats *tregion, char *line) {
+  
+  double v0, s, dtmp;
+  char tmp[10];
+
+  /* copy name (4 chars) */
+  bzero (tregion[0].name, 10);
+  strncpy (tregion[0].name, &line[1], 4);
+
+  /* get RA[0] */
+  bzero (tmp, 10);
+  strncpy (tmp, &line[7], 2);
+  v0 = atof (tmp);
+  
+  bzero (tmp, 10);
+  strncpy (tmp, &line[10], 2);
+  v0 += atof (tmp) / 60.0;
+  
+  bzero (tmp, 10);
+  strncpy (tmp, &line[13], 5);
+  v0 += atof (tmp) / 3600.0;
+  tregion[0].RA[0] = v0 * 15.0;
+
+  /* get RA[1] */
+  bzero (tmp, 10);
+  strncpy (tmp, &line[19], 2);
+  v0 = atof (tmp);
+  
+  bzero (tmp, 10);
+  strncpy (tmp, &line[22], 2);
+  v0 += atof (tmp) / 60.0;
+  
+  bzero (tmp, 10);
+  strncpy (tmp, &line[25], 5);
+  v0 += atof (tmp) / 3600.0;
+  tregion[0].RA[1] = v0 * 15.0;
+
+  /* this only happens at the 0,360 boundary */
+  if (tregion[0].RA[1] < tregion[0].RA[0]) tregion[0].RA[1] += 360.0;
+
+  /* get DEC[0] */
+  if (line[31] == '-') s = -1;
+  else s = +1;
+
+  bzero (tmp, 10);
+  strncpy (tmp, &line[32], 2);
+  v0 = atof (tmp);
+  
+  bzero (tmp, 10);
+  strncpy (tmp, &line[35], 2);
+  v0 += atof (tmp) / 60.0;
+  tregion[0].DEC[0] = s * v0;
+
+  /* get DEC[1] */
+  if (line[40] == '-') s = -1;
+  else s = +1;
+
+  bzero (tmp, 10);
+  strncpy (tmp, &line[41], 2);
+  v0 = atof (tmp);
+  
+  bzero (tmp, 10);
+  strncpy (tmp, &line[44], 2);
+  v0 += atof (tmp) / 60.0;
+  tregion[0].DEC[1] = s * v0;
+
+  if (tregion[0].DEC[0] > tregion[0].DEC[1]) {
+    dtmp = tregion[0].DEC[1];
+    tregion[0].DEC[1] = tregion[0].DEC[0];
+    tregion[0].DEC[0] = dtmp;
+  }
+  return (TRUE);
+}
+
+void init_regions () {
+
+  Nregion = 0;
+  NREGION = 100;
+
+  ALLOCATE (region, CatStats, NREGION);
+
+}
+
+void add_to_regions (CatStats *area) {
+
+  region[Nregion] = area[0];
+  ALLOCATE (region[Nregion].name, char, MAX_NAME);
+  sprintf (region[Nregion].name, "%s/%s.%s", CatDir, area[0].name, "cpt");
+
+  area_of_region (&region[Nregion]);
+
+  Nregion ++;
+
+  if (Nregion == NREGION) {
+    NREGION += 100;
+    REALLOCATE (region, CatStats, NREGION);
+  }
+
+}
+
+void area_of_region (CatStats *region) {
+  
+  double area;
+
+  area = DEG_RAD*(region[0].RA[1] - region[0].RA[0])*(sin(region[0].DEC[1]*RAD_DEG) - sin(region[0].DEC[0]*RAD_DEG));
+  region[0].Area = area;
+}
+
+void set_catalog (char *catdir) {
+  CatDir = catdir;
+}
Index: /trunk/Ohana/src/mosastro/src/match.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/match.c	(revision 3294)
+++ /trunk/Ohana/src/mosastro/src/match.c	(revision 3294)
@@ -0,0 +1,95 @@
+# include "mosastro.h"
+/* match is done on tangent plane : should be focal plane in some cases? */
+
+int match (StarData *refcat, int Nrefcat) {
+
+  int i, j, k, N, K, Ntotal;
+  int Nmatch, NMATCH;
+  double dp, dq, radius;
+  double *p, *q, *P, *Q;
+  int *U, *u;
+
+  /* sort the REFCAT data by P in tangent plane */ 
+  ALLOCATE (P, double, Nrefcat);
+  ALLOCATE (Q, double, Nrefcat);
+  ALLOCATE (U, int, Nrefcat);
+  for (i = 0; i < Nrefcat; i++) {
+    U[i] = i;
+    P[i] = refcat[i].P;
+    Q[i] = refcat[i].Q;
+  }
+  sortthree (P, Q, U, Nrefcat);
+
+
+  for (i = 0; i < Nchip; i++) {
+
+    Nmatch = 0;
+    NMATCH = 1000;
+    ALLOCATE (chip[i].raw, StarData, NMATCH);
+    ALLOCATE (chip[i].ref, StarData, NMATCH);
+
+    /* sort the star data by P in tangent plane */ 
+    ALLOCATE (p, double, chip[i].Nstars);
+    ALLOCATE (q, double, chip[i].Nstars);
+    ALLOCATE (u, int, chip[i].Nstars);
+    for (j = 0; j < chip[i].Nstars; j++) {
+      u[j] = j;
+      p[j] = chip[i].stars[j].P;
+      q[j] = chip[i].stars[j].Q;
+    }
+    sortthree (p, q, u, chip[i].Nstars);
+
+    /* find star matches in the tangent plane coord system */
+    for (j = k = 0; (j < chip[i].Nstars) && (k < Nrefcat); ) {
+      dp = p[j] - P[k];
+	
+      if (dp <= -2*RADIUS) {
+	j++;
+	continue;
+      }
+      if (dp >= 2*RADIUS) {
+	k++;
+	continue;
+      }
+
+      K = 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) {
+	  chip[i].raw[Nmatch] = chip[i].stars[u[j]];
+	  chip[i].ref[Nmatch] = refcat[U[k]];
+	  Nmatch ++;
+	  if (Nmatch == NMATCH) {
+	    NMATCH += 1000;
+	    REALLOCATE (chip[i].raw, StarData, NMATCH);
+	    REALLOCATE (chip[i].ref, StarData, NMATCH);
+	  }
+	  goto done;
+	}
+      }
+    done:
+      k = K;
+      j++;
+    }
+    REALLOCATE (chip[i].raw, StarData, MAX (1, Nmatch));
+    REALLOCATE (chip[i].ref, StarData, MAX (1, Nmatch));
+    chip[i].Nmatch = Nmatch;
+    free (p);
+    free (q);
+    free (u);
+  }
+  free (P);
+  free (Q);
+  free (U);
+  
+  Ntotal = 0;
+  for (i = 0; i < Nchip; i++) {
+    Ntotal += chip[i].Nmatch;
+    if (VERBOSE) fprintf (stderr, "chip %d: %d of %d stars\n", 
+			  i, chip[i].Nmatch, chip[i].Nstars);
+  }
+  fprintf (stderr, "Nchips: %d  Nmatch:  %d\n", Nchip, Ntotal);
+  return (1);
+}
Index: /trunk/Ohana/src/mosastro/src/mkheader.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/mkheader.c	(revision 3294)
+++ /trunk/Ohana/src/mosastro/src/mkheader.c	(revision 3294)
@@ -0,0 +1,68 @@
+# include "mosastro.h"
+
+Header *mkheader (int Nx, int Ny, int Nstars, Coords *coords) {
+
+  Header *header;
+
+  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] = Nx;
+  header[0].Naxis[1] = Ny;
+
+  fits_create_header (header);
+
+  fits_modify (header, "NSTARS",   "%d", 1, Nstars);
+  fits_modify (header, "PHOTCODE", "%s", 1, "STD.R");
+  fits_modify (header, "DATE-OBS", "%s", 1, "2004-04-22");
+  fits_modify (header, "UTC-OBS",  "%s", 1, "14:27:45.30");
+  fits_modify (header, "ZERO_PT", "%lf", 1, 25.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");
+  }    
+  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/mkobs.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/mkobs.c	(revision 3294)
+++ /trunk/Ohana/src/mosastro/src/mkobs.c	(revision 3294)
@@ -0,0 +1,115 @@
+# include "mosastro.h"
+
+int main (int argc, char **argv) {
+
+  int i, j, N, Nrefcat, Nchip, Nstars, NSTARS;
+  double RA, DEC;
+  double Ro, Do, Po, Qo, Lo, Mo;
+  double dX, dY;
+  char filename[64];
+  Coords map, coords;
+  StarData *refcat;
+  Stars *stars;
+  Header *header;
+
+  ConfigInit (&argc, argv);
+
+  if (argc != 4) {
+    fprintf (stderr, "USAGE: mkobs (RA) (DEC) (output)\n");
+    exit (2);
+  }
+
+  RA  = atof(argv[1]);
+  DEC = atof(argv[2]);
+  dX = dY = 3600*0.1;
+
+  fake_field (RA, DEC);
+
+  refcat = greference (&Nrefcat);
+  project_refcat (refcat, Nrefcat);
+
+  Nchip = 0;
+  for (i = -2; i < 3; i++) {
+    for (j = -2; j < 3; j++) {
+
+      /* find chip center in TP, FP and Sky */
+      Po = 0.2*i*3600.0;
+      Qo = 0.2*j*3600.0;
+      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");
+      map.crval1 = Lo;
+      map.crval2 = Mo;
+      map.crpix1 = dX;
+      map.crpix2 = dX;
+  
+      map.cdelt1 = 1.0;
+      map.cdelt2 = 1.0;
+
+      map.pc1_1  = 1;
+      map.pc2_2  = 1;
+      map.pc1_2  = 0;
+      map.pc2_1  = 0;
+      map.Npolyterms = 1;
+
+      /* project catalog stars to chip */
+      FPtoChip (refcat, Nrefcat, &map);
+
+      Nstars = 0;
+      NSTARS = 1000;
+      ALLOCATE (stars, Stars, NSTARS);
+
+      /* find catalog stars on chip */
+      for (N = 0; N < Nrefcat; N++) {
+	if (refcat[N].X < 0) continue;
+	if (refcat[N].Y < 0) continue;
+	if (refcat[N].X >= 2*dX) continue;
+	if (refcat[N].Y >= 2*dY) continue;
+
+	stars[Nstars].X = refcat[N].X;
+	stars[Nstars].Y = refcat[N].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;
+
+	Nstars ++;
+	if (Nstars >= NSTARS) {
+	  NSTARS += 1000;
+	  REALLOCATE (stars, Stars, NSTARS);
+	}
+      }
+
+      /* Chip-Sky terms */
+      strcpy (coords.ctype, "RA---WRP");
+      coords.crval1 = Ro;
+      coords.crval2 = Do;
+      coords.crpix1 = dX;
+      coords.crpix2 = dX;
+  
+      coords.cdelt1 = 1/3600.0;
+      coords.cdelt2 = 1/3600.0;
+
+      coords.pc1_1  = 1;
+      coords.pc2_2  = 1;
+      coords.pc1_2  = 0;
+      coords.pc2_1  = 0;
+      coords.Npolyterms = 1;
+
+      header = mkheader (2*dX, 2*dY, Nstars, &coords);
+      sprintf (filename, "%s.%02d.fits", argv[3], Nchip);
+      wstars (filename, stars, Nstars, header);
+      fits_free_header (header);
+      free (stars);
+      Nchip ++;
+    }
+  }
+  exit (0);
+}
Index: /trunk/Ohana/src/mosastro/src/mosastro.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/mosastro.c	(revision 3294)
+++ /trunk/Ohana/src/mosastro/src/mosastro.c	(revision 3294)
@@ -0,0 +1,85 @@
+# include "mosastro.h"
+# include <glob.h>
+
+int main (int argc, char **argv) {
+
+  int Nrefcat;
+  glob_t pglob;
+  StarData *refcat;
+  Gradients *grad;
+
+  ConfigInit (&argc, argv);
+  args (&argc, argv); 
+
+  if (argc != 3) { 
+    fprintf (stderr, "USAGE: mosastro (inglob) (ext)\n");
+    exit (2);
+  }
+
+  pglob.gl_offs = 0;
+  glob (argv[1], 0, NULL, &pglob);
+
+  LoadStars (pglob.gl_pathc, pglob.gl_pathv);
+
+  deproject_stars ();
+  field_stats ();
+  init_field ();
+  init_chips ();
+
+  refcat = greference (&Nrefcat);
+  project_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 ();
+
+  grad = GetGradients ();
+  FitGradients (grad);
+
+  deproject_raw ();
+  project_ref ();
+
+  FitChips (1);
+
+  deproject_raw ();
+  project_ref ();
+
+  ClipOnFP (3.0);
+
+  FitField (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/parse_time.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/parse_time.c	(revision 3294)
+++ /trunk/Ohana/src/mosastro/src/parse_time.c	(revision 3294)
@@ -0,0 +1,126 @@
+# include "mosastro.h"
+
+parse_time (Header *header) {
+
+  double jd;
+  int Ny, Nf, mode;
+  int Nsec, hour, min, sec, year, month, day;
+  char *py, *pm, *pd, *c;
+  char line[256];
+
+  /* we want to find JD or MJD to get Nsec (seconds since 01/01/1970) */
+
+  /* try JD first */
+  if (strcasecmp (JDKeyword, "NONE")) {
+    uppercase (JDKeyword);
+    fits_scan (header, JDKeyword, "%lf", 1, &jd);
+    Nsec = (jd - 2440587.5)*86400;
+    return (Nsec);
+  }
+
+  /* try MJD next */
+  if (strcasecmp (MJDKeyword, "NONE")) {
+    uppercase (MJDKeyword);
+    fits_scan (header, MJDKeyword, "%lf", 1, &jd);
+    Nsec = (jd - 40587.0)*86400;
+    return (Nsec);
+  }
+    
+  /* get UT and DATE */
+  uppercase (UTKeyword);
+  fits_scan (header, UTKeyword, "%s", 1, line);
+  /* remove ':' characters */
+  for (c = strchr (line, 0x3a); c != (char *) NULL; c = strchr (line, 0x3a)) { *c = ' '; }
+  sscanf (line, "%d %d %d", &hour, &min, &sec);
+
+  /* parse mode line */
+  uppercase (DateMode);
+  for (Ny = 0, c = strchr (DateMode, 'Y'); c != (char ) NULL; c = strchr (c + 1, 'Y'), Ny++);
+  if ((Ny != 2) && (Ny != 4)) {
+    fprintf (stderr, "error in DATE-MODE format: %s\n", DateMode);
+    exit (1);
+  }
+  py = strchr (DateMode, 'Y');
+  pm = strchr (DateMode, 'M');
+  pd = strchr (DateMode, 'D');
+  if ((py == (char *) NULL) || (pm == (char *) NULL) || (pd == (char *) NULL)) {
+    fprintf (stderr, "error in DATE-MODE format: %s\n", DateMode);
+    exit (1);
+  }
+  if ((py > pm) && (py < pd)) {
+    fprintf (stderr, "error in DATE-MODE format: %s\n", DateMode);
+    exit (1);
+  }
+  if ((py > pd) && (py < pm)) {
+    fprintf (stderr, "error in DATE-MODE format: %s\n", DateMode);
+    exit (1);
+  }
+  mode = 0;
+  if ((py < pm) && (pm < pd)) { mode = 1; }  /* yyyy-mm-dd */
+  if ((py < pm) && (pm > pd)) { mode = 2; }  /* yyyy-dd-mm */
+  if ((py > pm) && (pm < pd)) { mode = 3; }  /* mm-dd-yyyy */
+  if ((py > pm) && (pm > pd)) { mode = 4; }  /* dd-mm-yyyy */
+  if (!mode) {
+    fprintf (stderr, "error in DATE-MODE format: %s\n", DateMode);
+    exit (1);
+  }
+
+  /* parse date entry */
+  uppercase (DateKeyword);
+  fits_scan (header, DateKeyword, "%s",  1, line);
+  /* remove possible separators: ':', '/' '.', '-' */
+  for (c = strchr (line, 0x3a); c != (char *) NULL; c = strchr (line, 0x3a)) { *c = ' '; }
+  for (c = strchr (line, 0x2f); c != (char *) NULL; c = strchr (line, 0x2f)) { *c = ' '; }
+  for (c = strchr (line, 0x2e); c != (char *) NULL; c = strchr (line, 0x2e)) { *c = ' '; }
+  for (c = strchr (line, 0x2d); c != (char *) NULL; c = strchr (line, 0x2d)) { *c = ' '; }
+
+  switch (mode) {
+  case 1:
+    Nf = sscanf (line, "%d %d %d", &year, &month, &day);
+    break;
+  case 2:
+    Nf = sscanf (line, "%d %d %d", &year, &day, &month);
+    break;
+  case 3:
+    Nf = sscanf (line, "%d %d %d", &month, &day, &year);
+    break;
+  case 4:
+    Nf = sscanf (line, "%d %d %d", &day, &month, &year);
+    break;
+  }
+  if (Nf != 3) {
+    fprintf (stderr, "error in date entry (%s) or DATE-MODE format (%s)\n", line, DateMode);
+    exit (1);
+  }
+
+  if (year > 1000) {
+    if (Ny == 2) {
+      fprintf (stderr, "warning: mode line claims 2 digit year, but 4 digit year found\n");
+    }
+  } else {
+    if (Ny == 4) {
+      fprintf (stderr, "warning: mode line claims 4 digit year, but 2 digit year found\n");
+    }
+    if (year < 50) year += 100;
+    year += 1900;
+  }    
+
+  /* convert yy.mm.dd hh.mm.ss to Nsec since 1970 (jd = 2440587.5) */
+  /* note that in this section, tm_mon has range 1-12, unlike for gmtime () */
+  jd = day - 32075 + (int)(1461*(year + 4800 + (int)(((month)-14)/12))/4)
+    + (int)(367*((month) - 2 - (int)(((month) - 14)/12)*12)/12)
+    - (int)(3*(int)((year + 4900 + (int)(((month) - 14)/12))/100)/4) - 0.5;
+  /* jd is the julian day of the whole day only not the time */
+  Nsec = (jd - 2440587.5)*86400 + 3600.0*hour + min*60.0 + sec;
+  
+  return (Nsec);
+
+}
+
+uppercase (char *string) {
+
+  int i;
+    
+  for (i = 0; i < strlen (string); i++) string[i] = toupper (string[i]);
+
+}
Index: /trunk/Ohana/src/mosastro/src/project.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/project.c	(revision 3294)
+++ /trunk/Ohana/src/mosastro/src/project.c	(revision 3294)
@@ -0,0 +1,46 @@
+# include "mosastro.h"
+
+int project_refcat (StarData *refcat, int Nrefcat) {
+
+  int i;
+
+  SkyToTP (refcat, Nrefcat, &field.project);
+  TPtoFP  (refcat, Nrefcat, &field.distort);
+
+  return (1);
+}
+
+int deproject_stars () {
+
+  int i, j;
+
+  for (i = 0; i < Nchip; i++) {
+    ChipToSky (chip[i].stars, chip[i].Nstars, &chip[i].coords);
+    SkyToTP   (chip[i].stars, chip[i].Nstars, &field.project);
+  }
+  return (1);
+}
+
+int deproject_raw () {
+
+  int i, j;
+
+  for (i = 0; i < Nchip; i++) {
+    ChipToFP (chip[i].raw, chip[i].Nmatch, &chip[i].map);
+    FPtoTP   (chip[i].raw, chip[i].Nmatch, &field.distort);
+    TPtoSky  (chip[i].raw, chip[i].Nmatch, &field.project);
+  }
+  return (1);
+}
+
+int project_ref () {
+
+  int i, j;
+
+  for (i = 0; i < Nchip; i++) {
+    SkyToTP  (chip[i].ref, chip[i].Nmatch, &field.project);
+    TPtoFP   (chip[i].ref, chip[i].Nmatch, &field.distort);
+    FPtoChip (chip[i].ref, chip[i].Nmatch, &chip[i].map);
+  }
+  return (1);
+}
Index: /trunk/Ohana/src/mosastro/src/sort.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/sort.c	(revision 3294)
+++ /trunk/Ohana/src/mosastro/src/sort.c	(revision 3294)
@@ -0,0 +1,50 @@
+# include "mosastro.h"
+
+int sortthree (double *X, double *Y, int *Z, int N) {
+
+  int l,j,ir,i;
+  double tX, tY;
+  int tZ;
+  
+  if (N < 2) return (1);
+  l = N >> 1;
+  ir = N - 1;
+  for (;;) {
+    if (l > 0) {
+      l--;
+      tX = X[l];
+      tY = Y[l];
+      tZ = Z[l];
+    }
+    else {
+      tX = X[ir];
+      X[ir] = X[0];
+      tY = Y[ir];
+      Y[ir] = Y[0];
+      tZ = Z[ir];
+      Z[ir] = Z[0];
+      if (--ir == 0) {
+        X[0] = tX;
+        Y[0] = tY;
+        Z[0] = tZ;
+        return (1);
+      }
+    }
+    i = l;
+    j = (l << 1) + 1;
+    while (j <= ir) {
+      if (j < ir && X[j] < X[j+1]) j++;
+      if (tX < X[j]) {
+        X[i] = X[j];
+        Y[i] = Y[j];
+        Z[i] = Z[j];
+        j += (i=j) + 1;
+      }
+      else j = ir + 1;
+    }
+    X[i] = tX;
+    Y[i] = tY;
+    Z[i] = tZ;
+  }
+  return (1);
+}
Index: /trunk/Ohana/src/mosastro/src/transforms.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/transforms.c	(revision 3294)
+++ /trunk/Ohana/src/mosastro/src/transforms.c	(revision 3294)
@@ -0,0 +1,54 @@
+# include "mosastro.h"
+
+ChipToSky (StarData *stars, int Nstars, Coords *coords) {
+  int i;
+  for (i = 0; i < Nstars; i++) {
+    XY_to_RD (&stars[i].R, &stars[i].D, stars[i].X, stars[i].Y, coords);
+    while (stars[i].R < 0.0)    stars[i].R += 360.0;
+    while (stars[i].R >= 360.0) stars[i].R -= 360.0;
+  }
+}
+
+ChipToFP (StarData *stars, int Nstars, Coords *coords) {
+  int i;
+  for (i = 0; i < Nstars; i++) {
+    XY_to_RD (&stars[i].L, &stars[i].M, stars[i].X, stars[i].Y, coords);
+  }
+}
+
+FPtoTP (StarData *stars, int Nstars, Coords *coords) {
+  int i;
+  for (i = 0; i < Nstars; i++) {
+    XY_to_RD (&stars[i].P, &stars[i].Q, stars[i].L, stars[i].M, coords);
+  }
+}
+
+TPtoSky (StarData *stars, int Nstars, Coords *coords) {
+  int i;
+  for (i = 0; i < Nstars; i++) {
+    XY_to_RD (&stars[i].R, &stars[i].D, stars[i].P, stars[i].Q, coords);
+    while (stars[i].R < 0.0)    stars[i].R += 360.0;
+    while (stars[i].R >= 360.0) stars[i].R -= 360.0;
+  }
+}
+
+SkyToTP (StarData *stars, int Nstars, Coords *coords) {
+  int i;
+  for (i = 0; i < Nstars; i++) {
+    RD_to_XY (&stars[i].P, &stars[i].Q, stars[i].R, stars[i].D, coords);
+  }
+}
+
+TPtoFP (StarData *stars, int Nstars, Coords *coords) {
+  int i;
+  for (i = 0; i < Nstars; i++) {
+    RD_to_XY (&stars[i].L, &stars[i].M, stars[i].P, stars[i].Q, coords);
+  }
+}
+
+FPtoChip (StarData *stars, int Nstars, Coords *coords) {
+  int i;
+  for (i = 0; i < Nstars; i++) {
+    RD_to_XY (&stars[i].X, &stars[i].Y, stars[i].L, stars[i].M, coords);
+  }
+}
Index: /trunk/Ohana/src/mosastro/src/wstars.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/wstars.c	(revision 3294)
+++ /trunk/Ohana/src/mosastro/src/wstars.c	(revision 3294)
@@ -0,0 +1,43 @@
+# include "mosastro.h"
+# define NCHAR 66 /* 65 char EXCLUDING return */
+
+void wstars (char *filename, Stars *stars, int Nstars, Header *header) {
+  
+  int i, Nchar;
+  FILE *g;
+  char line[NCHAR + 3];
+
+  g = fopen (filename, "w");
+  if (g == (FILE *) NULL) {
+    fprintf (stderr, "ERROR: can't create output file %s\n", filename);
+    exit (1);
+  }
+
+  fwrite (header[0].buffer, 1, header[0].size, g);
+
+  for (i = 0; i < Nstars; i++) {
+    if (stars[i].dophot == 0) continue;
+    Nchar = snprintf (line, NCHAR, "%6.1f %6.1f %6.3f %03d %2d %3.1f %6.3f %6.3f %6.2f %6.2f %5.1f", 
+		      stars[i].X, stars[i].Y, stars[i].M, 
+		      (int)(1000*stars[i].dM), stars[i].dophot, stars[i].sky, 
+		      stars[i].Mgal, stars[i].Map, stars[i].fx, stars[i].fy, stars[i].df);
+    
+    /* this is just a little funny.  NCHAR (in) includes the trailing NULL, Nchar (out) excludes it */
+    if (Nchar != NCHAR - 1) {
+      fprintf (stderr, "funny line %d (%d)\n%s\n", i, Nchar, line);
+    } else {
+      fprintf (g, "%s\n", line);
+    }
+  }
+
+  fclose (g);
+
+}
+
+/*
+
+  63.6 2869.5 17.568 157 17.568 17.568 25.01 25.00 360.0 7 2.9
+
+  63.6 2869.5 17.568 157 7 2.9
+
+*/