Index: /branches/eam_branches/relastro.20100326/Makefile
===================================================================
--- /branches/eam_branches/relastro.20100326/Makefile	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/Makefile	(revision 27487)
@@ -0,0 +1,61 @@
+default: relastro
+help:
+	@echo "make options: relastro (default)"
+
+include ../../Makefile.System
+HOME    =       $(ROOT)/src/relastro
+BIN	=	$(HOME)/bin
+LIB	=	$(HOME)/lib
+SRC	=	$(HOME)/src
+MAN	=	$(HOME)/doc
+INC	= 	$(HOME)/include
+include ../../Makefile.Common
+
+# programs may add their own internal requirements here
+FULL_CFLAGS   = $(BASE_CFLAGS)
+FULL_CPPFLAGS = $(BASE_CPPFLAGS)
+FULL_LDFLAGS  = -lkapa -ldvo -lFITS -lohana $(BASE_LDFLAGS)
+
+relastro: $(BIN)/relastro.$(ARCH)
+install: $(DESTBIN)/relastro
+
+RELASTRO = \
+$(SRC)/ConfigInit.$(ARCH).o	     \
+$(SRC)/FitChip.$(ARCH).o             \
+$(SRC)/FitMosaic.$(ARCH).o           \
+$(SRC)/FitPM.$(ARCH).o               \
+$(SRC)/FitPar.$(ARCH).o              \
+$(SRC)/FitPMandPar.$(ARCH).o         \
+$(SRC)/FitSimple.$(ARCH).o           \
+$(SRC)/ImageOps.$(ARCH).o	     \
+$(SRC)/MosaicOps.$(ARCH).o	     \
+$(SRC)/ParFactor.$(ARCH).o           \
+$(SRC)/SetSignals.$(ARCH).o 	     \
+$(SRC)/Shutdown.$(ARCH).o 	     \
+$(SRC)/UpdateChips.$(ARCH).o         \
+$(SRC)/UpdateMosaic.$(ARCH).o        \
+$(SRC)/UpdateObjects.$(ARCH).o       \
+$(SRC)/UpdateObjectOffsets.$(ARCH).o \
+$(SRC)/UpdateSimple.$(ARCH).o        \
+$(SRC)/UpdateMeasures.$(ARCH).o      \
+$(SRC)/GetAstromError.$(ARCH).o      \
+$(SRC)/args.$(ARCH).o		     \
+$(SRC)/bcatalog.$(ARCH).o	     \
+$(SRC)/dvo_astrom_ops.$(ARCH).o      \
+$(SRC)/fitpoly.$(ARCH).o             \
+$(SRC)/initialize.$(ARCH).o	     \
+$(SRC)/liststats.$(ARCH).o	     \
+$(SRC)/load_catalogs.$(ARCH).o	     \
+$(SRC)/load_images.$(ARCH).o	     \
+$(SRC)/mkpolyterm.$(ARCH).o          \
+$(SRC)/plot_scatter.$(ARCH).o	     \
+$(SRC)/plotstuff.$(ARCH).o	     \
+$(SRC)/select_images.$(ARCH).o	     \
+$(SRC)/relastro.$(ARCH).o	     \
+$(SRC)/relastro_objects.$(ARCH).o    \
+$(SRC)/save_catalogs.$(ARCH).o       \
+$(SRC)/write_coords.$(ARCH).o        \
+$(SRC)/relastroVisual.$(ARCH).o
+
+$(RELASTRO): $(INC)/relastro.h $(KAPA_INCS)
+$(BIN)/relastro.$(ARCH): $(RELASTRO) $(KAPA_LIBS)
Index: /branches/eam_branches/relastro.20100326/doc/notes.txt
===================================================================
--- /branches/eam_branches/relastro.20100326/doc/notes.txt	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/doc/notes.txt	(revision 27487)
@@ -0,0 +1,183 @@
+
+2010.03.24
+
+  I have finished only part of the optimizations listed below: I have
+  made the detection->image and image->mosaic lookups use bracket
+  searches.  I have NOT converted the image->mosaic lookups to use the
+  parentID, so these matches are still based on time ranges.  I have
+  also added a bit of positional constraint to the image match
+  function in 'select_images' (comparing Rmax of the sky regions to
+  Rmin of the images).  These have had a huge improvement.  Here is
+  the execution time for the SAS data (~450 images) with the old vs
+  the new code on pikake:
+
+  old : 39min
+  new :  7min
+
+  The bulk of the time for the new code is now disk I/O.
+
+2010.03.23
+
+  I've been working on optimization.  I've updated ImageOps to make
+  the measure->image relationship generation use the image IDs along
+  with a bracket search -- this goes much faster as a result.  More
+  work is needed:
+
+  * select_images needs to apply astrometry to all image corners: 
+    -> cache the image center and radius and use this to narrow down
+       the searches.
+
+  * MosaicOps uses the time to match mosaics.  I've added sorting to
+    speed this up, but this should be done based on the image and
+    parent index.
+
+    -> make a tool to create the parent ID for existing DBs.
+    -> make sure addstar is populating the parent IDs
+    -> use the parent IDs to make the link.
+
+load image data: 0.006347 sec
+  setup sky: 0.091475 sec
+  convert image table: 0.091498 sec
+  select images: 2.925978 sec
+  init images: 2.926957 sec
+  init mosaics: 2.927319 sec
+load images: 2.933762 sec
+load catalog data: 9.772259 sec
+make image bins: 9.780618 sec
+set up image indexes: 10.196519 sec
+
+load image data: 0.006746 sec
+  setup sky: 0.091440 sec
+  convert image table: 0.091467 sec
+  select images: 0.374193 sec
+  init images: 0.375355 sec
+  init mosaics: 0.375753 sec
+load images: 0.382606 sec
+load catalog data: 7.260094 sec
+make image bins: 7.268564 sec
+set up image indexes: 7.687485 sec
+
+
+2008.03.01
+
+  relastro (and relphot for that matter) needs to have a somewhat
+  flexible error model for the detections.  Some data have an accurate
+  astrometric error; some have a reported error, but it is either
+  insufficient (eg, only from poisson error), or inaccurate.  A
+  complete error model would be something like this:
+
+  dR_total^2 =  dR_sys^2 + AS * dR_obs^2 + MS * dM_obs^2
+
+  dR_sys : systematicAstrometryError
+  AS     : astrometryErrorScale
+  MS     : astrometryErrorMagScale
+
+2007.11.12
+
+  relastro is working for the SIMPLE (single-chip) and CHIP
+  (mosaic-chip) modes.  it is not yet working for the MOSAIC mode.
+  The image table contains one mosaic entry and N chip entries for a
+  mosaic with N chips.  The mosaic entry has coord.ctype = DIS, while
+  the corresponding chips have coord.ctype = WRP (and matching
+  time/photcodes).
+
+  (this is fairly weak. we should just define the relationship when we
+  register the chips...)
+
+2007.10.06
+
+  relastro currently loads the set of detections, averages, and image,
+  and sets up arrays of pointers to define the relationships beteen
+  the detections, objects, and images.  This is probably an acceptable
+  scheme.  
+
+  However, the generation of those index arrays is probably very
+  inefficient.  Currently it is written for the old DVO tables in
+  which there was no real link between the detection and image.  We
+  are trying to move to a situation where each detection has a
+  reference to the image (imageID).  If this is unique, then we should
+  be able to make these links directly.
+
+  Outstanding issues:
+
+  * am I getting the correct set of matches for the mosaic images? 
+    ie, do getImageRaw & getImageRef return the correct list when the
+    image is a DIS mosaic image?
+
+  * is the operation in fit_apply_coords doing the correct thing?
+    compare this to psastro code.  This should not be a problem since
+    we are only fitting in the linear frames.  we are just fitting
+    simple polynomials.  We are not changing the reference coordinate
+    for the projection.
+
+  * measure the errors?
+
+  * apply the magnitude errors as weight?
+
+  * improve matching process with the index information available?
+
+  * add the ability to down-weight groups of detections by photcode?
+
+2007.02.11
+
+  relastro major modes:
+
+  - update the astrometry of objects in the images
+    - load objects within a region
+    - foreach object
+      - calculate average R,D
+      - update db tables (dR, dD)
+    - include external trends
+      - parallax factor
+      - atm trends
+
+  - update the astrometric parameters of images
+    - use the average R,D for objects, update image terms
+    * include external refs with adjustable weighting
+
+  - image parameter smoothing
+
+  - simultaneous fit of objects and images
+    - identify parameters to constrain
+    - limited number of internal objects?
+
+
+2006.04.08 : relastro contemplations
+
+relastro will perform the astrometry equivalent to relphot.  For
+region, it loads in the average objects, detections, and image
+parameters.  The goal is to determine improved astrometric positions
+of objects in the images and to determine improved astrometric
+parameters for the images.  Just like relphot, this is an iterative
+process in which the object postions are improved in one step, the
+image parameters are improved in the next, then back to the images.
+
+Relastro needs to deal with outlier objects in a robust way.  There
+are a few reasons why objects should have poor astrometric solutions,
+and these can be flagged up front: saturation and too many bad pixels.
+For detections which come from psphot, we can use the quality
+parameter to exclude objects which have excessively high bad pixel
+counts.  For detections from sextractor or other analysis tools,
+this information is missing.  In this case, we can use the distance
+from the image boundaries as a trigger.
+
+An open question to be addressed is the issue of external calibration
+sources.  Unlike photometry, external astrometric measurements do not
+need to have an unknown filter transformation applied.  They can thus
+be included as part of the solution automatically.  It is not clear,
+however, how to include *proper-motion* data in the analysis.  If I
+introduce, for example, USNO objects with proper motions, the
+information provided consists in practice of a velocity vector and
+errors.  How does this get folded into the fit?  One answer is to
+construct additional artificial datapoints for each entry with the
+proper errors and epochs.  Seems like something of a hack.  
+
+We may want to allow for different errors / weightings for the
+different external sources.
+
+After the initial iterations, we will can solve for proper-motion and
+parallax for (all?) objects
+
+- calculate the parallax normal vector for each image
+- calculate the atmospheric / orbital motion effects for each image
+
Index: /branches/eam_branches/relastro.20100326/doc/parallax.pro
===================================================================
--- /branches/eam_branches/relastro.20100326/doc/parallax.pro	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/doc/parallax.pro	(revision 27487)
@@ -0,0 +1,348 @@
+
+macro sample
+  PMsim 15 30 20 10 $ANGLE
+  set dX = zero(X) + 10
+  set dY = zero(Y) + 10
+
+  style -c black
+  plot -x 2 -pt 2 -sz 2.0 X Y -dx dX -dy dY
+
+  # set x = (1000.0/$dS)*cos(t*$TP) + $vX*t
+  # set y = (1000.0/$dS)*sin(t*$TP)*dsin($ANGLE) + $vY*t
+  # plot -c grey60 -x 2 -sz 2.0 -pt 7 x y
+
+  create t 0 1.5 0.01
+  set x = (1000.0/$dS)*cos(t*$TP) + $vX*t
+  set y = (1000.0/$dS)*sin(t*$TP)*dsin($ANGLE) + $vY*t
+  plot -c grey60 -x 0 x y
+
+  # label -x "offset east (milliarcsec)"
+  # label -y "offset north (milliarcsec)"
+end
+
+macro PMsim
+  if ($0 != 6)
+    echo "USAGE: PMsim Do vX vY dP dT"
+    echo "Do : distance in parsec"
+    echo "vX, vY : velocity (milliarcsec/year)"
+    echo "dP : astrometric scatter (arcsec)"
+    echo "dT : ecliptic lattitude (degrees)"
+    break
+  end
+
+  $TP = 2*3.1416
+
+  $dS = $1
+  $vX = $2
+  $vY = $3
+  $dP = $4
+  $dT = $5
+
+  mkbase
+
+  set x = (1000.0/$dS)*cos(t*$TP) + $vX*t
+  set y = (1000.0/$dS)*sin(t*$TP)*dsin($dT) + $vY*t
+
+  gaussdev dx x[] 0.0 $dP
+  gaussdev dy y[] 0.0 $dP 
+
+  set X = x + dx
+  set Y = y + dy
+
+  if ($PLOT)
+    lim -200 200 -200 200; clear; box;
+    style -c black
+    plot -x 2 -pt 2 -sz 2.0 X Y
+    plot -x 0 X Y
+  end
+
+  PMfit X Y t $dT
+
+  if ($PLOT) 
+    label -x "offset east (milliarcsec)"
+    label -y "offset north (milliarcsec)"
+    section a 0 0 1 1
+    lim 0 1 0 1
+    sprintf line "input: D: %5.1f (pc), V: %5.1f, %5.1f (mas/yr), dS: %4.1f" $1 $2 $3 $4
+    textline -fn courier 14 0.15 0.95 "$line"
+    sprintf line "fit:    D: %5.1f (pc), V: %5.1f, %5.1f (mas/yr)" $dS $vX $vY
+    textline -fn courier 14 0.15 0.90 "$line"
+    section default
+  end
+end
+
+macro PMfit
+  if ($0 != 5)
+    echo "USAGE: PMfit X Y t (theta)"
+    echo "X,Y are offsets in milliarcsec"
+    echo "t is time in years"
+    break
+  end
+
+  $TP = 2*3.1416
+  set x = $1
+  set y = $2
+  set t = $3
+  set cs = cos(t*$TP)
+  set snx = sin(t*$TP)*dsin($4)
+
+  mcreate A 3 3
+  create B 0 3
+
+  # define A values:
+  set tmp = cs^2 + snx^2
+  vstat -q tmp
+  zap A 0 0 1 1 -v $TOTAL
+
+  set tmp = t*cs
+  vstat -q tmp
+  zap A 1 0 1 1 -v $TOTAL
+  zap A 0 1 1 1 -v $TOTAL
+
+  set tmp = t*snx
+  vstat -q tmp
+  zap A 2 0 1 1 -v $TOTAL
+  zap A 0 2 1 1 -v $TOTAL
+
+  set tmp = t^2
+  vstat -q tmp
+  zap A 1 1 1 1 -v $TOTAL
+  zap A 2 2 1 1 -v $TOTAL
+
+  zap A 2 1 1 1 -v 0
+  zap A 1 2 1 1 -v 0
+
+  # define B values:
+  set tmp = x*cs + y*snx
+  vstat -q tmp
+  B[0] = $TOTAL
+
+  set tmp = x*t
+  vstat -q tmp
+  B[1] = $TOTAL
+
+  set tmp = y*t
+  vstat -q tmp
+  B[2] = $TOTAL
+
+  gaussj A B
+
+  $dS = 1000/B[0]
+  $vX = B[1]
+  $vY = B[2]
+
+  if ($PLOT)
+    echo "Do: $dS"
+    echo "vX: $vX"
+    echo "vY: $vY"
+
+    set x = (1000.0/$dS)*cos(t*$TP) + $vX*t
+    set y = (1000.0/$dS)*sin(t*$TP)*dsin($4) + $vY*t
+    plot -c red -x 2 -sz 2.0 -pt 7 x y
+  end
+end
+
+macro PMstats
+  
+  if ($0 != 7)
+    echo "USAGE: PMsim Do vX vY dP dT (Niter)"
+    echo "Do : distance in parsec"
+    echo "vX, vY : velocity (milliarcsec/year)"
+    echo "dP : astrometric scatter (arcsec)"
+    echo "dT : ecliptic lattitude (degrees)"
+    break
+  end
+
+  delete -q vx vy ds
+  $PLOT = 0
+
+  for i 0 $6
+    PMsim $1 $2 $3 $4 $5
+    concat $dS ds
+    concat $vX vx
+    concat $vY vy
+  end
+
+  vstat ds
+  vstat vx
+  vstat vy
+
+  set gm = 1000/ds
+  vstat gm
+
+  clear -s -n 1
+  label -fn courier 14
+  # section a 0.0 0.0 0.5 1.0
+  lim {$2-30} {$2+30} {$3-30} {$3+30}; box
+  plot -x 2 -pt 2 -sz 0.5 vx vy
+  label -x "p.m. (east, mas/yr)"
+  label -y "p.m. (east, mas/yr)"
+
+  clear -s -n 2
+  # section b 0.5 0.0 0.5 1.0
+  histogram gm Ng {1000/$1 - 30.0} {1000/$1 + 30.0} 1.0
+  create dg {1000/$1 - 30.0} {1000/$1 + 30.0} 1.0
+  lim dg Ng; box
+  plot -x 1 dg Ng
+
+  label -x "parallax (mas)"
+  label -y "# of tests"
+
+end
+
+macro PMtrend
+
+  PMstats 10 5.0 5.0 10 90 100
+  echo 10
+  PMstats 10 5.0 5.0 {750/2.38/25} 90 100
+  echo {750/2.38/25} 
+  PMstats 10 5.0 5.0 {750/2.38/10} 90 100
+  echo {750/2.38/10} 
+  PMstats 10 5.0 5.0 {750/2.38/5} 90 100
+  echo {750/2.38/5} 
+end
+
+list dist
+  5
+  10
+  15
+  20
+  30
+  40
+  50
+  60
+  70
+  80
+end
+
+macro ParDetect
+
+  local i
+  delete -q xp yp dyp dym
+  
+  $ANGLE = 45
+
+  for i 0 $dist:n
+    PMstats $dist:$i 5.0 5.0 10 $ANGLE 100
+    concat $dist:$i xp
+    concat {1000/$MEAN} yp
+    concat {1000/($MEAN - $SIGMA) - 1000/$MEAN} dyp
+    concat {1000/$MEAN - 1000/($MEAN + $SIGMA)} dym
+  end
+
+  section a 0.0 0.0 0.5 1.0
+
+  lim xp -1.0 101; clear; box; 
+  plot -c black -x 2 xp yp -dy dym +dy dyp
+  
+  delete -q xp yp 
+  concat 0 xp
+  concat 100 xp
+  set yp = xp
+  plot -x 0 -c blue xp yp
+
+  set yp = zero(xp) + 62.5
+  plot -x 0 -c red -lt 1 xp yp
+  style -lt 0  
+  label -x "input dist (pc)" -y "output dist" -fn courier 14
+
+  $PLOT = 0
+  section b 0.5 0.0 0.5 1.0
+  lim -120 120 -120 120; box -labels 1001;
+
+  for i 0 5
+   sample
+  end
+
+  label -x "offset east (milliarcsec)"
+  label +y "offset north"
+  ps -name ParallaxDist.ps
+end
+
+macro fields
+
+  $Dx = 0.25*4800*8/3600
+  $Fo = 1.09
+
+  $dD = 2.65/1.09
+
+  $Nf = 0
+  for i {-30+$dD/2} 90 $dD
+    $Np = int(dcos($i)*1.09*360/2.65) + 1
+    $Nf = $Nf + $Np
+    fprintf "%5.1f %5.1f %4d %4d" $i {dcos($i)*1.09*360/2.65} $Np $Nf 
+  end
+end
+
+macro SNsim
+  $w  = 0.76
+  $m1 = 24.8
+  $mu = 20.4
+
+  lim 13.5 26.5 0.0 4.0; clear; box
+
+  $t = 1
+  SN
+
+  $t = 5
+  SN
+
+  $t = 30
+  SN
+
+  $t = 100
+  SN
+
+  set sn = 5 + zero (m)
+  set lsn = log(sn)
+  plot -lt 2 m lsn
+
+  set sn = 25 + zero (m)
+  set lsn = log(sn)
+  plot -lt 2 m lsn
+
+  set sn = 100 + zero (m)
+  set lsn = log(sn)
+  plot -lt 2 m lsn
+  style -lt 0
+end
+
+macro SN
+  create m 14 26 0.1
+  set sn = ten (-0.2*(2*m - $mu - $m1)) * sqrt($t/$w^2/3.1416)
+  set lsn = log(sn)
+  plot m lsn
+end
+
+# r: 30s @ 5sig = 22.65  (0.6")
+# i: 30s @ 5sig = 22.42  (0.6")
+# i: 30s @ 5sig = 22.27  (for 0.3 mag lower zp)
+
+macro mkbase
+  delete -q t
+  concat 0.0 t
+  concat {10.0/365.0} t
+  concat {20.0/365.0} t
+  concat { 6.0/52.0} t
+
+  concat { 3.0/12.0} t
+  concat { 6.0/12.0} t
+  concat {12.0/12.0} t
+  concat {18.0/12.0} t
+end
+
+
+macro parfactor
+
+  $RA = $1
+  $Dec = $2
+  $e = 23 + 27/60
+  
+  create s 0 360
+  set pR = dcos($e)*d
+
+  set pR =  +(dcos($e)*dsin(s)*dcos($RA) - dcos(s)*dsin($RA))
+  set pD =  -(dcos($e)*dsin(s)*dsin($RA) + dcos(s)*dcos($RA))*dsin($Dec) + dsin($e)*dsin(s)*dcos($Dec)
+
+  lim -1.1 1.1 -1.1 1.1; clear; box; plot pR pD
+end
Index: /branches/eam_branches/relastro.20100326/doc/relastro.txt
===================================================================
--- /branches/eam_branches/relastro.20100326/doc/relastro.txt	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/doc/relastro.txt	(revision 27487)
@@ -0,0 +1,18 @@
+
+Config (.ptolemyrc) variables used by relastro:
+
+  RELASTRO_SIGMA_LIM   : exclude measurements with dMag > SIGMA_LIM for image calibrations
+  PM_DT_MIN            : minimum temporal baseline needed to calculate a proper motion
+  PM_TOOFEW            : minimum number of valid detections to calculate a proper motion
+  POS_TOOFEW           : minimum number of valid detections to calculate an average position
+
+  CATDIR               : database directory
+
+  GSCFILE              : template sky layout definition table (used for CATDIRs created before this feature)
+  SKY_DEPTH            : define the output sky catalog size
+  SKY_TABLE            : specify an alternate name for the sky catalog table
+
+  PHOTCODE_FILE        : template photcode table (used for CATDIRs created before the photcode table was added to catdir)
+
+  CATMODE              : define the output catalog mode (mef,split,etc) -- this is ignored
+  CATFORMAT            : define the output catalog format (eg, PS1_DEV_1) -- this is ignored
Index: /branches/eam_branches/relastro.20100326/include/relastro.h
===================================================================
--- /branches/eam_branches/relastro.20100326/include/relastro.h	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/include/relastro.h	(revision 27487)
@@ -0,0 +1,304 @@
+# include <ohana.h>
+# include <dvo.h>
+# include <kapa.h>
+# include <signal.h>
+
+typedef enum {
+  MODE_SIMPLE,
+  MODE_CHIP,
+  MODE_MOSAIC,
+} CoordMode;
+
+typedef enum {ERROR_MODE_RA, ERROR_MODE_DEC, ERROR_MODE_POS} ErrorMode;
+
+typedef enum {FIT_NONE, FIT_AVERAGE, FIT_PM_ONLY, FIT_PAR_ONLY, FIT_PM_AND_PAR} FitMode;
+
+typedef enum {TARGET_NONE, TARGET_OBJECTS, TARGET_SIMPLE, TARGET_CHIPS, TARGET_MOSAICS} FitTarget;
+
+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;
+  double dMag;
+  double dPos;
+  int mask;
+} StarData;
+
+// structure to hold coordinate fitting terms
+typedef struct {
+    int Npts;
+    int Nterms;
+    int Norder;
+    int Nsums;
+    int Nelems;
+    double **sum;
+    double **xsum;
+    double **ysum;
+    double **xfit;
+    double **yfit;
+} CoordFit;
+
+typedef struct {
+  double Ro, dRo;
+  double Do, dDo;
+
+  double uR, duR;
+  double uD, duD;
+
+  double p, dp;
+
+  double chisq;
+  int Nfit;
+} PMFit;
+
+typedef struct {
+  unsigned int start;
+  unsigned int stop;
+  float Mcal;
+  float dMcal;
+  short Xm;
+  float secz;
+  char flags;
+  Coords coords;
+} Mosaic;
+
+typedef struct {
+  double median;
+  double mean;
+  double sigma;
+  double error;
+  double chisq;
+  double min;
+  double max;
+  double total;
+  int    Nmeas;
+} StatType;
+
+/* global variables set in parameter file */
+char   ImageCat[256];
+char   GSCFILE[256];
+char   CATDIR[256];
+char   CATMODE[16];    /* raw, mef, split, mysql */
+char   CATFORMAT[16];  /* internal, elixir, loneos, panstarrs */
+char   SKY_TABLE[256];
+int    SKY_DEPTH;  /** XXX EAM : depth of catalog tables, fix usage */
+
+double SIGMA_LIM;
+int SRC_MEAS_TOOFEW; //catalog objects wich fewer detections then this are ignored
+double MIN_ERROR;
+
+int    VERBOSE;
+
+int    RESET;
+int    UPDATE;
+int    PLOTSTUFF;
+int    SAVEPLOT;
+int    SHOW_PARAMS;
+char   STATMODE[32];
+int    POS_TOOFEW;
+int    PM_TOOFEW;
+double PM_DT_MIN;
+int    PLOTDELAY;
+int    CHIPORDER;
+
+char          *PHOTCODE_KEEP_LIST, *PHOTCODE_SKIP_LIST;
+int           NphotcodesKeep,      NphotcodesSkip;
+PhotCode     **photcodesKeep,     **photcodesSkip;
+
+int AreaSelect;
+double AreaXmin, AreaXmax, AreaYmin, AreaYmax;
+
+int ImagSelect, ImagMin, ImagMax;
+
+double  PlotMmin, PlotMmax, PlotdMmin, PlotdMmax;
+
+int DophotSelect, DophotValue;
+int PhotFlagSelect, PhotFlagPoor, PhotFlagBad;
+
+int TimeSelect;
+time_t TSTART, TSTOP;
+
+int FlagOutlier;
+int    CLIP_THRESH;
+
+FitMode FIT_MODE;
+
+FitTarget FIT_TARGET;
+
+SkyRegion UserPatch;
+
+int DoUpdateObjects;
+int DoUpdateSimple;
+int DoUpdateChips;
+int DoUpdateMosaics;
+
+/*** relphot prototypes ***/
+void          ConfigInit          PROTO((int *argc, char **argv));
+void          GetConfig           PROTO((char *config, char *field, char *format, int N, void *ptr));
+char         *GetPhotnamebyCode   PROTO((PhotCodeData *photcodes, int code));
+void          InterpolateGrid     PROTO((float *buffer, int Nx, int Ny, Coords *ccd, Coords *gcoords));
+off_t          *SelectRefMosaic     PROTO((Mosaic **refmosaic, off_t *Nimage));
+int           args                PROTO((int argc, char **argv));
+int           bcatalog            PROTO((Catalog *subcatalog, Catalog *catalog));
+void          clean_images        PROTO((void));
+void          clean_measures      PROTO((Catalog *catalog, int Ncatalog, int final));
+void          clean_mosaics       PROTO((void));
+void          clean_stars         PROTO((Catalog *catalog, int Ncatalog));
+int           corner_check        PROTO((double *x1, double *y1, double *x2, double *y2));
+void          dumpGrid            PROTO((void));
+void          dump_grid           PROTO((void));
+int           edge_check          PROTO((double *x1, double *y1, double *x2, double *y2));
+void          findImages          PROTO((Catalog *catalog, int Ncatalog));
+int           findMosaics         PROTO((Catalog *catalog, int Ncatalog));
+Image        *find_images         PROTO((FITS_DB *db, GSCRegion *region, off_t Nregion, off_t *Nimage, off_t **LineNum));
+void set_db (FITS_DB *in);
+int Shutdown (char *format, ...) OHANA_FORMAT(printf, 1, 2);
+void TrapSignal (int sig);
+void SetProtect (int mode);
+int SetSignals (void);
+
+GSCRegion    *find_regions        PROTO((Image *image, off_t Nimage, int *Nregions, GSCRegion *fullregion));
+void          freeGridBins        PROTO((int Ncatalog));
+void          freeImageBins       PROTO((int Ncatalog));
+void          freeMosaicBins      PROTO((int Ncatalog));
+void          free_catalogs       PROTO((Catalog *catalog, int Ncatalog));
+int           gcatalog            PROTO((Catalog *catalog, int FINAL));
+Coords       *getCoords           PROTO((off_t meas, int cat));
+float         getMcal             PROTO((off_t meas, int cat));
+float         getMgrid            PROTO((off_t meas, int cat));
+float         getMmos             PROTO((off_t meas, int cat));
+float         getMrel             PROTO((Catalog *catalog, off_t meas, int cat));
+GSCRegion    *get_regions         PROTO((double minRa, double maxRa, double minDec, double maxDec, off_t *Nregions));
+void          getfullregion       PROTO((Image *image, off_t Nimage, GSCRegion *fullregion));
+Image        *getimage            PROTO((off_t N));
+Image        *getimages           PROTO((off_t *N));
+void          global_stats        PROTO((Catalog *catalog, int Ncatalog));
+void          initGrid            PROTO((int dX, int dY));
+void          initGridBins        PROTO((Catalog *catalog, int Ncatalog));
+void          initImageBins       PROTO((Catalog *catalog, int Ncatalog));
+void          initImages          PROTO((Image *input, off_t N));
+void          initMosaicBins      PROTO((Catalog *catalog, int Ncatalog));
+void          initMosaicGrid      PROTO((Image *image, off_t Nimage));
+void          initMosaics         PROTO((Image *image, off_t Nimage));
+void          initMrel            PROTO((Catalog *catalog, int Ncatalog));
+void          initialize          PROTO((int argc, char **argv));
+void          initstats           PROTO((char *mode));
+int           liststats           PROTO((double *value, double *dvalue, int N, StatType *stats));
+Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int subselect));
+SkyList      *load_images         PROTO((FITS_DB *db, SkyRegion *region));
+Image        *select_images       PROTO((SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage));
+
+void check_permissions (char *basefile);
+void lock_image_db (FITS_DB *db, char *filename);
+void unlock_image_db (FITS_DB *db);
+void create_image_db (FITS_DB *db);
+void save_catalogs (Catalog *catalog, int Ncatalog);
+
+int           main                PROTO((int argc, char **argv));
+void          mark_images         PROTO((Image *image, off_t Nimage, Image *timage, off_t Ntimage));
+void          matchImage          PROTO((Catalog *catalog, off_t meas, int cat));
+void          matchMosaics        PROTO((Catalog *catalog, off_t meas, int cat));
+GSCRegion    *name_region         PROTO((char *name, off_t *Nregions));
+double        opening_angle       PROTO((double x1, double y1, double x2, double y2, double x3, double y3));
+void          plot_chisq          PROTO((Catalog *catalog, int Ncatalog));
+void          plot_defaults       PROTO((Graphdata *graphdata));
+void          plot_grid           PROTO((Catalog *catalog));
+void          plot_images         PROTO((void));
+void          plot_list           PROTO((Graphdata *graphdata, double *xlist, double *ylist, int N, char *label, char *file));
+void          plot_mosaic_fields  PROTO((Catalog *catalog));
+void          plot_mosaics        PROTO((void));
+void          plot_scatter        PROTO((Catalog *catalog, int Ncatalog));
+void          plot_star_coords    PROTO((Catalog *catalog, int Ncatalog));
+void          plot_stars          PROTO((Catalog *catalog, int Ncatalog));
+void          reload_catalogs     PROTO((SkyList *skylist));
+int           setExclusions       PROTO((Catalog *catalog, int Ncatalog));
+void          setMcal             PROTO((Catalog *catalog, int Poor));
+void          setMcalFinal        PROTO((void));
+int           setMcalOutput       PROTO((Catalog *catalog, int Ncatalog));
+void          setMgrid            PROTO((Catalog *catalog));
+int           setMmos             PROTO((Catalog *catalog, int Poor));
+int           setMrel             PROTO((Catalog *catalog, int Ncatalog));
+void          setMrelFinal        PROTO((Catalog *catalog));
+int           setMrelOutput       PROTO((Catalog *catalog, int Ncatalog, int mark));
+void          set_ZP              PROTO((double ZERO));
+int           setrefcode          PROTO((Image *image, off_t Nimage));
+void          skip_measurements   PROTO((Catalog *catalog, int pass));
+void          sortA               PROTO((double *X, int N));
+void          sortB               PROTO((double *X, double *Y, int N));
+void          sortC               PROTO((double *X, double *Y, double *F1, double *F2, int N));
+void          sortD               PROTO((double *X, double *Y, double *Z, int N));
+StatType      statsImageM         PROTO((Catalog *catalog));
+StatType      statsImageN         PROTO((Catalog *catalog));
+StatType      statsImageX         PROTO((Catalog *catalog));
+StatType      statsImagedM        PROTO((Catalog *catalog));
+StatType      statsMosaicM        PROTO((Catalog *catalog));
+StatType      statsMosaicN        PROTO((Catalog *catalog));
+StatType      statsMosaicX        PROTO((Catalog *catalog));
+StatType      statsMosaicdM       PROTO((Catalog *catalog));
+StatType      statsStarN          PROTO((Catalog *catalog, int Ncatalog));
+StatType      statsStarS          PROTO((Catalog *catalog, int Ncatalog));
+StatType      statsStarX          PROTO((Catalog *catalog, int Ncatalog));
+void          wcatalog            PROTO((Catalog *catalog));
+void          wimages             PROTO((void));
+void          write_coords        PROTO((Header *header, Coords *coords));
+
+double **array_init (int Nx, int Ny);
+void array_free (double **array, int Nx);
+CoordFit *fit_init (int order);
+void fit_free (CoordFit *fit);
+void fit_add (CoordFit *fit, double x1, double y1, double x2, double y2, double wt);
+void fit_eval (CoordFit *fit);
+void fit_apply (CoordFit *fit, double *x2, double *y2, double x1, double y1);
+double **poly2d_dx (double **poly, int Nx, int Ny);
+double **poly2d_dy (double **poly, int Nx, int Ny);
+double **poly2d_copy (double **poly, int Nx, int Ny);
+double poly2d_eval (double **poly, int Nx, int Ny, double x, double y);
+CoordFit *fit_apply_coords (CoordFit *fit, Coords *coords);
+int CoordsGetCenter (CoordFit *fit, double tol, double *xo, double *yo);
+CoordFit *CoordsSetCenter (CoordFit *input, double Xo, double Yo);
+void FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
+void FitMosaic (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
+void FitSimple (StarData *raw, StarData *ref, int Nmatch, Coords *coords);
+void initObjectData (Catalog *catalog, int Ncatalog);
+int UpdateObjects (Catalog *catalog, int Ncatalog);
+int UpdateSimple (Catalog *catalog, int Ncatalog);
+int UpdateChips (Catalog *catalog, int Ncatalog);
+int UpdateMosaic (Catalog *catalog, int Ncatalog);
+int UpdateMeasures (Catalog *catalog, int Ncatalog);
+void fixImageRaw (Catalog *catalog, int Ncatalog, off_t im);
+void FlagOutliers(Catalog *catalog);
+int MeasFilterTest(Measure *measure);
+
+int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon);
+int ParFactor (double *pR, double *pD, double R, double D, time_t T);
+int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts);
+int FitPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *pR, double *pD, int Npts);
+int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts);
+
+Mosaic *getMosaicForImage (off_t N);
+
+StarData *getImageRef (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode);
+StarData *getImageRaw (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode);
+
+Mosaic *getmosaics (off_t *N);
+void initMosaics (Image *image, off_t Nimage);
+StarData *getMosaicRaw (Catalog *catalog, int Ncatalog, off_t mos, off_t *Nstars);
+StarData *getMosaicRef (Catalog *catalog, int Ncatalog, off_t mos, off_t *Nstars);
+Mosaic *getMosaicForImage (off_t im);
+
+double getMeanR (Measure *measure, Average *average, SecFilt *secfilt);
+double getMeanD (Measure *measure, Average *average, SecFilt *secfilt);
+int setMeanR (double ra_fit, Measure *measure, Average *average, SecFilt *secfilt);
+int setMeanD (double dec_fit, Measure *measure, Average *average, SecFilt *secfilt);
+
+float GetAstromError (Measure *measure, int mode);
+int relastro_objects (void);
+int UpdateObjectOffsets (SkyList *skylist);
+
+int relastroVisualPlotRawRef(StarData *raw, StarData *ref, double dRmax, int numObj);
+int relastroVisualPlotScatter(double values[], double thresh, int npts);
+int relastroVisualPlotOutliers(Catalog *catalog, int offset, int Nmeasure, 
+			       StatType statsR, StatType statsD, double thresh);
+
Index: /branches/eam_branches/relastro.20100326/src/ConfigInit.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/ConfigInit.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/ConfigInit.c	(revision 27487)
@@ -0,0 +1,83 @@
+# include "relastro.h"
+
+void ConfigInit (int *argc, char **argv) {
+
+  char  *config, *file;
+  char CatdirPhotcodeFile[256];
+  char MasterPhotcodeFile[256];
+  struct stat filestat;
+  int status;
+
+  /*** 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, "RELASTRO_SIGMA_LIM",     "%lf", 0, &SIGMA_LIM);
+  GetConfig (config, "RELASTRO_SRC_MEAS_TOOFEW",   "%d",  0, &SRC_MEAS_TOOFEW);
+
+  // XXX these are used in relphot to identify poor stars / images -- define
+  // an equivalent concept here?
+  // GetConfig (config, "STAR_SCATTER",           "%lf", 0, &STAR_SCATTER);
+  // GetConfig (config, "IMAGE_SCATTER",          "%lf", 0, &IMAGE_SCATTER);
+  // GetConfig (config, "IMAGE_OFFSET",           "%lf", 0, &IMAGE_OFFSET);
+  // GetConfig (config, "STAR_CHISQ",             "%lf", 0, &STAR_CHISQ);
+
+  GetConfig (config, "PM_DT_MIN",              "%lf", 0, &PM_DT_MIN);
+  GetConfig (config, "PM_TOOFEW",              "%d",  0, &PM_TOOFEW);
+  GetConfig (config, "POS_TOOFEW",             "%d",  0, &POS_TOOFEW);
+
+  GetConfig (config, "GSCFILE",                "%s",  0, GSCFILE);
+  GetConfig (config, "CATDIR",                 "%s",  0, CATDIR);
+  ScanConfig(config, "CATMODE",                "%s",  0, CATMODE);
+  ScanConfig(config, "CATFORMAT",              "%s",  0, CATFORMAT);
+  ScanConfig(config, "PHOTCODE_FILE",          "%s",  0, MasterPhotcodeFile);
+
+  sprintf (ImageCat, "%s/Images.dat", CATDIR);
+
+  if (!ScanConfig (config, "SKY_DEPTH",         "%d",  0, &SKY_DEPTH)) {
+    SKY_DEPTH = 2;
+  }
+  if (!ScanConfig (config, "SKY_TABLE",         "%s",  0, SKY_TABLE)) {
+    SKY_TABLE[0] = 0;
+  }
+
+  if (*CATMODE == 0) strcpy (CATMODE, "RAW");
+  if (*CATFORMAT == 0) strcpy (CATFORMAT, "ELIXIR");
+
+  // check for existence of CATDIR
+  status = stat (CATDIR, &filestat);
+  if (status == -1) {
+    fprintf (stderr, "directory %s does not exist, giving up\n", CATDIR);
+    exit (1);
+  }
+
+  /* update master photcode table if not defined */
+  sprintf (CatdirPhotcodeFile, "%s/Photcodes.dat", CATDIR);
+  if (!LoadPhotcodes (CatdirPhotcodeFile, MasterPhotcodeFile, TRUE)) {
+    fprintf (stderr, "error loading photcode table %s or master file %s\n", CatdirPhotcodeFile, MasterPhotcodeFile);
+    exit (1);
+  }
+  SetZeroPoint (25.0);
+
+  free (config);
+  free (file);
+
+}
+
+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: /branches/eam_branches/relastro.20100326/src/FitChip.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/FitChip.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/FitChip.c	(revision 27487)
@@ -0,0 +1,151 @@
+# include "relastro.h"
+# include "relastroVisual.h"
+
+// XXX make these user parameters
+# define FIT_CHIP_MAX_ERROR 0.05
+# define FIT_CHIP_NITER     3
+# define FIT_CHIP_NSIGMA    3.0
+
+
+// XXX we should test if the fit is sufficiently constrained across the chip, or if the
+// new positions deviate too much from the old positions
+
+// XXX add visualization tools: per chip residual plots
+
+// XXX save the fit[0].Npts value in the image table?
+
+// XXX save measurements of the fit quality (scatter, chisq) in the image table
+
+void FitChip (StarData *raw, StarData *ref, int Nmatch, Coords *coords) {
+
+  int i, Nscatter, Niter, skip;
+  CoordFit *fit;
+  double dL, dM, dR, dRmax, *values;
+
+  ALLOCATE (values, double, Nmatch);
+  for (Niter = 0; Niter < FIT_CHIP_NITER; Niter ++) {
+
+    // measure the scatter distribution (use only the bright end detections)
+    for (i = Nscatter = 0; i < Nmatch; i++) {
+      if (raw[i].mask) continue;
+      if (isnan(raw[i].dMag) || raw[i].dMag > FIT_CHIP_MAX_ERROR) continue;
+
+      dL = raw[i].L - ref[i].L;
+      dM = raw[i].M - ref[i].M;
+      dR = hypot (dL, dM);
+
+      values[Nscatter] = dR;
+      Nscatter++;
+    }
+
+    // for a 2D Gaussian, 40% of the points are within 1 sigma; dRmax is ~ 3 sigma
+    dsort (values, Nscatter);
+    dRmax = FIT_CHIP_NSIGMA*values[(int)(0.40*Nscatter)];
+    relastroVisualPlotScatter(values, dRmax, Nscatter);
+    relastroVisualPlotRawRef(raw, ref, dRmax, Nmatch);
+
+    // fit the requested order polynomial
+    if (CHIPORDER > 0) {
+      coords[0].Npolyterms = CHIPORDER;
+    }
+    fit = fit_init (coords[0].Npolyterms);
+
+    // generate the fit matches
+    for (i = 0; i < Nmatch; i++) {
+      if (raw[i].mask) {
+	continue;
+      }
+      if (isnan(raw[i].dMag) || raw[i].dMag > FIT_CHIP_MAX_ERROR) {
+	continue;
+      }
+
+      // only keep objects within dRmax
+      dL = raw[i].L - ref[i].L;
+      dM = raw[i].M - ref[i].M;
+      dR = hypot (dL, dM);
+      if (dR > dRmax) continue;
+
+      fit_add (fit, raw[i].X, raw[i].Y, ref[i].L, ref[i].M, raw[i].dPos);
+    }
+
+    // check if the fit has enough data points for the polynomial order
+    skip = FALSE;
+    switch (coords[0].Npolyterms) {
+      case 0:
+      case 1:
+        skip = (fit[0].Npts < 8);
+        break;
+      case 2:
+        skip = (fit[0].Npts < 11);
+        break;
+      case 3:
+        skip = (fit[0].Npts < 15);
+        break;
+      default:
+        fprintf (stderr, "invalid chip order %d\n", coords[0].Npolyterms);
+        abort ();
+    }
+    if (skip) {
+      fprintf (stderr, "insufficient measurements (%d) for requested order (%d)\n", fit[0].Npts, coords[0].Npolyterms);
+      fit_free (fit);
+      free (values);
+      return;
+    }
+
+    fprintf (stderr, "scatter limit: %f based on %d detections; using %d of %d for fit\n", dRmax, Nscatter, fit[0].Npts, Nmatch);
+
+    // measure the fit, update the coords & object coordinates
+    fit_eval (fit);
+    fit_apply_coords (fit, coords);
+    fit_free (fit);
+
+    for (i = 0; i < Nmatch; i++) {
+      XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, coords);
+    }
+
+  }
+
+  free (values);
+  return;
+}
+
+/* in the mosaic case, we have four coord systems of interest:
+   R,D : the sky
+   P,Q : the tangent plane
+   L,M : the focal plane
+   X,Y : the chip
+
+   R,D -> P,Q (projection)
+   P,Q -> L,M (polynomial transformation : DIS)
+   L,M -> X,Y (polynomial transformation : WRP)
+
+Some details about this function:
+
+- the initial value of the clipping radius is set based on the distribution of the dR
+  values for the bright sources.
+
+- NITER clipping passes are performed
+
+- the per-star astrometric errors are included in the fit.  The photcode table controls
+  whether only the reported positional errors are used, or if the magnitudes errors are
+  scaled to determine the error, and if there is a systematic floor for all sources.
+
+- input detections are pre-filtered on the basis of the photFlags, etc, in bcatalog
+
+- outlying detections are clipped in the iterative passes, but the clipped detections are
+  not recorded
+
+*/
+
+/* example using fit_apply() :
+
+  f = fopen ("test3.dat", "w");
+
+  // apply new coords to raw (X,Y -> L,M)
+  for (i = 0; i < Nmatch; i++) {
+    fprintf (f, "%f %f  %f %f  ", raw[i].X, raw[i].Y, raw[i].L, raw[i].M);
+    fit_apply (newfit, &L1, &M1, raw[i].X - coords[0].crpix1, raw[i].Y - coords[0].crpix2);
+    fprintf (f, "%f %f\n", L1, M1);
+  }
+  fclose (f);
+*/
Index: /branches/eam_branches/relastro.20100326/src/FitMosaic.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/FitMosaic.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/FitMosaic.c	(revision 27487)
@@ -0,0 +1,48 @@
+# include "relastro.h"
+
+void FitMosaic (StarData *raw, StarData *ref, int Nmatch, Coords *coords) {
+
+  int i;
+  CoordFit *fit;
+  double dP, dQ, dR;
+
+  fit = fit_init (coords[0].Npolyterms);
+  for (i = 0; i < Nmatch; i++) {
+    if (raw[i].mask) continue;
+
+    // require radius of XXX arcsec
+    dP = raw[i].P - ref[i].P;
+    dQ = raw[i].Q - ref[i].Q;
+    dR = 3600.0 * hypot (dP, dQ);
+
+    // XXX the value needs to be set in a more intelligent way
+    if (dR > 0.15) continue;
+    
+    fit_add (fit, raw[i].L, raw[i].M, ref[i].P, ref[i].Q, 1.0);
+  }
+  if (fit[0].Npts == 0) {
+    fit_free (fit);
+    return;
+  }
+  fit_eval (fit);
+  fit_apply_coords (fit, coords);
+  fit_free (fit);
+
+  // apply new coords to raw (X,Y -> L,M)
+  for (i = 0; i < Nmatch; i++) {
+    XY_to_LM (&raw[i].P, &raw[i].Q, raw[i].L, raw[i].M, coords);
+  }
+}
+
+/* in the mosaic case, we have four coord systems of interest:
+   R,D : the sky
+   P,Q : the tangent plane
+   L,M : the focal plane
+   X,Y : the chip
+
+   R,D -> P,Q (projection)
+   P,Q -> L,M (polynomial transformation : DIS)
+   L,M -> X,Y (polynomial transformation : WRP)
+*/
+
+/* XXX I'm not using the errors at all : this could at least be done with the dMag values */
Index: /branches/eam_branches/relastro.20100326/src/FitPM.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/FitPM.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/FitPM.c	(revision 27487)
@@ -0,0 +1,87 @@
+# include "relastro.h"
+
+/* do we want an init function which does the alloc and a clear function to free? */
+int FitPM (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, int Npts) {
+
+  int i;
+
+  double **A, **B;
+  double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT;
+  double chisq, Xf, Yf;
+
+  /* do I need to do this as 2 2x2 matrix equations? */
+  A = array_init (4, 4);
+  B = array_init (4, 1);
+
+  Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0;
+  for (i = 0; i < Npts; i++) {
+    /* handle case where dX or dY = 0.0 */
+    wx = 1.0 / SQ(dX[i]);
+    wy = 1.0 / SQ(dY[i]);
+
+    Wx += wx;
+    Wy += wy;
+
+    Tx += T[i]*wx;
+    Ty += T[i]*wy;
+    
+    Tx2 += SQ(T[i])*wx;
+    Ty2 += SQ(T[i])*wy;
+    
+    Xs += X[i]*wx;
+    Ys += Y[i]*wy;
+
+    XT += X[i]*T[i]*wx;
+    YT += Y[i]*T[i]*wy;
+  }
+
+  A[0][0] = Wx;
+  A[0][1] = Tx;
+
+  A[1][0] = Tx;
+  A[1][1] = Tx2;
+
+  A[2][2] = Wy;
+  A[2][3] = Ty;
+
+  A[3][2] = Ty;
+  A[3][3] = Ty2;
+
+  B[0][0] = Xs;
+  B[1][0] = XT;
+  B[2][0] = Ys;
+  B[3][0] = YT;
+
+  dgaussjordan (A, B, 4, 1);
+
+  fit[0].Ro = B[0][0];
+  fit[0].uR = B[1][0];
+  fit[0].Do = B[2][0];
+  fit[0].uD = B[3][0];
+  fit[0].p  = 0.0;
+  
+  fit[0].dRo = sqrt(A[0][0]);
+  fit[0].duR = sqrt(A[1][1]);
+  fit[0].dDo = sqrt(A[2][2]);
+  fit[0].duD = sqrt(A[3][3]);
+  fit[0].dp  = 0.0;
+  
+  array_free (A, 4);
+  array_free (B, 4);
+
+  // add up the chi square for the fit
+  chisq = 0.0;
+  for (i = 0; i < Npts; i++) {
+      Xf = fit[0].Ro + fit[0].uR*T[i];
+      Yf = fit[0].Do + fit[0].uD*T[i];
+      chisq += SQ(X[i] - Xf) / SQ(dX[i]);
+      chisq += SQ(Y[i] - Yf) / SQ(dY[i]);
+  }
+  fit[0].Nfit = Npts;
+
+  // the reduced chisq is divided by (Ndof = 2*Npts - 4)
+  fit[0].chisq = chisq / (2.0*Npts - 4.0);
+  return (TRUE);
+}
+
+// XXX this function should (optionally?) iterate and clip outlier detections
Index: /branches/eam_branches/relastro.20100326/src/FitPMandPar.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/FitPMandPar.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/FitPMandPar.c	(revision 27487)
@@ -0,0 +1,98 @@
+# include "relastro.h"
+
+/* do we want an init function which does the alloc and a clear function to free? */
+int FitPMandPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *T, double *pR, double *pD, int Npts) {
+
+  int i;
+
+  double **A, **B;
+  double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT;
+  double PR, PD, PRT, PDT, PRX, PDY, PR2, PD2;
+
+  A = array_init (5, 5);
+  B = array_init (5, 1);
+
+  PR = PD = PRT = PDT = PRX = PDY = PR2 = PD2 = 0.0;
+  Wx = Wy = Tx = Ty = Tx2 = Ty2 = Xs = Ys = XT = YT = 0.0;
+  for (i = 0; i < Npts; i++) {
+    /* handle case where dX or dY = 0.0 */
+    wx = 1.0 / SQ(dX[i]);
+    wy = 1.0 / SQ(dY[i]);
+
+    Wx += wx;
+    Wy += wy;
+
+    Tx += T[i]*wx;
+    Ty += T[i]*wy;
+    
+    Tx2 += SQ(T[i])*wx;
+    Ty2 += SQ(T[i])*wy;
+    
+    PR += pR[i]*wx;
+    PD += pD[i]*wy;
+    
+    PRT += pR[i]*T[i]*wx;
+    PDT += pD[i]*T[i]*wy;
+    
+    PRX += pR[i]*X[i]*wx;
+    PDY += pD[i]*Y[i]*wy;
+    
+    PR2 += SQ(pR[i])*wx;
+    PD2 += SQ(pD[i])*wy;
+
+    Xs += X[i]*wx;
+    Ys += Y[i]*wy;
+
+    XT += X[i]*T[i]*wx;
+    YT += Y[i]*T[i]*wy;
+  }
+
+  A[0][0] = Wx;
+  A[0][1] = Tx;
+  A[0][4] = PR;
+
+  A[1][0] = Tx;
+  A[1][1] = Tx2;
+  A[1][4] = PRT;
+
+  A[2][2] = Wy;
+  A[2][3] = Ty;
+  A[2][4] = PD;
+
+  A[3][2] = Ty;
+  A[3][3] = Ty2;
+  A[3][4] = PDT;
+
+  A[4][0] = PR;
+  A[4][1] = PRT;
+  A[4][2] = PD;
+  A[4][3] = PDT;
+  A[4][4] = PR2 + PD2;
+
+  B[0][0] = Xs;
+  B[1][0] = XT;
+  B[2][0] = Ys;
+  B[3][0] = YT;
+  B[4][0] = PRX + PDY;
+
+  dgaussjordan (A, B, 5, 1);
+
+  fit[0].Ro = B[0][0];
+  fit[0].uR = B[1][0];
+  fit[0].Do = B[2][0];
+  fit[0].uD = B[3][0];
+  fit[0].p  = B[4][0];
+  
+  fit[0].dRo = sqrt(A[0][0]);
+  fit[0].duR = sqrt(A[1][1]);
+  fit[0].dDo = sqrt(A[2][2]);
+  fit[0].duD = sqrt(A[3][3]);
+  fit[0].dp  = sqrt(A[4][4]);
+  
+  array_free (A, 5);
+  array_free (B, 5);
+
+  /* get the chisq from the matrix values */
+
+  return (TRUE);
+}
Index: /branches/eam_branches/relastro.20100326/src/FitPar.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/FitPar.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/FitPar.c	(revision 27487)
@@ -0,0 +1,74 @@
+# include "relastro.h"
+
+/* do we want an init function which does the alloc and a clear function to free? */
+int FitPar (PMFit *fit, double *X, double *dX, double *Y, double *dY, double *pR, double *pD, int Npts) {
+
+  int i;
+
+  double **A, **B;
+  double wx, wy, Wx, Wy, Xs, Ys;
+  double PR, PD, PRX, PDY, PR2, PD2;
+
+  A = array_init (3, 3);
+  B = array_init (3, 1);
+
+  Wx = Wy = Xs = Ys = 0.0;
+  PR = PD = PRX = PDY = PR2 = PD2 = 0.0;
+  for (i = 0; i < Npts; i++) {
+    /* handle case where dX or dY = 0.0 */
+    wx = 1.0 / SQ(dX[i]);
+    wy = 1.0 / SQ(dY[i]);
+
+    Wx += wx;
+    Wy += wy;
+
+    PR += pR[i]*wx;
+    PD += pD[i]*wy;
+    
+    PRX += pR[i]*X[i]*wx;
+    PDY += pD[i]*Y[i]*wy;
+    
+    PR2 += SQ(pR[i])*wx;
+    PD2 += SQ(pD[i])*wy;
+
+    Xs += X[i]*wx;
+    Ys += Y[i]*wy;
+  }
+
+  A[0][0] = Wx;
+  A[0][2] = PR;
+
+  A[1][1] = Wy;
+  A[1][2] = PD;
+
+  A[2][0] = PR;
+  A[2][1] = PD;
+  A[2][2] = PR2 + PD2;
+
+  B[0][0] = Xs;
+  B[1][0] = Ys;
+  B[2][0] = PRX + PDY;
+
+  dgaussjordan (A, B, 3, 1);
+
+  fit[0].Ro = B[0][0];
+  fit[0].Do = B[1][0];
+  fit[0].p  = B[2][0];
+
+  fit[0].uR = 0.0;
+  fit[0].uD = 0.0;
+  
+  fit[0].dRo = sqrt(A[0][0]);
+  fit[0].dDo = sqrt(A[2][2]);
+  fit[0].dp  = sqrt(A[4][4]);
+  
+  fit[0].duR = 0.0;
+  fit[0].duD = 0.0;
+
+  array_free (A, 3);
+  array_free (B, 3);
+
+  /* get the chisq from the matrix values */
+
+  return (TRUE);
+}
Index: /branches/eam_branches/relastro.20100326/src/FitSimple.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/FitSimple.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/FitSimple.c	(revision 27487)
@@ -0,0 +1,52 @@
+# include "relastro.h"
+
+void FitSimple (StarData *raw, StarData *ref, int Nmatch, Coords *coords) {
+
+  int i;
+  CoordFit *fit;
+  double dP, dQ, dR;
+
+  fit = fit_init (coords[0].Npolyterms);
+  for (i = 0; i < Nmatch; i++) {
+    if (raw[i].mask) continue;
+
+    // require radius of XXX arcsec
+    dP = raw[i].P - ref[i].P;
+    dQ = raw[i].Q - ref[i].Q;
+    dR = 3600.0 * hypot (dP, dQ);
+
+    // XXX the value needs to be set in a more intelligent way
+    if (dR > 0.50) continue;
+    
+    fit_add (fit, raw[i].X, raw[i].Y, ref[i].P, ref[i].Q, 1.0);
+  }
+  if (fit[0].Npts == 0) {
+    fit_free (fit);
+    return;
+  }
+  fit_eval (fit);
+  fit_apply_coords (fit, coords);
+  fit_free (fit);
+
+  // apply new coords to raw (X,Y -> P,Q)
+  for (i = 0; i < Nmatch; i++) {
+    XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, coords);
+    raw[i].P = raw[i].L;
+    raw[i].Q = raw[i].M;
+  }
+}
+
+/* in the simple case, we only have three coord systems of interest:
+   R,D : the sky
+   P,Q : the tangent plane
+   X,Y : the chip
+
+   R,D -> P,Q (projection)
+   P,Q -> X,Y (polynomial transformation)
+
+   L,M is maintained, but is identical to P,Q
+*/
+
+/* XXX I'm not using the errors at all : this could at least be done with the dMag values */
+
+/* XXX See notes in FitChips.c */
Index: /branches/eam_branches/relastro.20100326/src/GetAstromError.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/GetAstromError.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/GetAstromError.c	(revision 27487)
@@ -0,0 +1,40 @@
+# include "relastro.h"
+
+float GetAstromError (Measure *measure, int mode) {
+
+  PhotCode *code;
+  float dPobs, dPsys, dPtotal, dM, AS, MS;
+
+  switch (mode) {
+    case ERROR_MODE_RA:
+      dPobs = measure[0].dXccd;  // need to redefine this as RAerr
+      break;
+    case ERROR_MODE_DEC:
+      dPobs = measure[0].dYccd;  // need to redefine this as RAerr
+      break;
+    case ERROR_MODE_POS:
+      dPobs = hypot (measure[0].dXccd, measure[0].dYccd);  // need to redefine this as RAerr
+      break;
+    default:
+      abort();
+  }
+
+  /* the astrometric errors are not being carried yet (but should be!) */
+  /* we use the photometric mag error as a weighting term */
+
+  code 	= GetPhotcodebyCode (measure[0].photcode);
+  AS   	= code[0].astromErrScale;
+  MS   	= code[0].astromErrMagScale;
+  dPsys = code[0].astromErrSys;
+  dM    = measure[0].dM;
+  
+  dPtotal = sqrt(SQ(dPsys) + AS*SQ(dPobs) + MS*SQ(dM));
+
+  //XXX dXccd, dYccd are now working correctly
+  //dPtotal = AS * dPobs;
+
+  dPtotal = MAX (dPtotal, MIN_ERROR);
+
+  return (dPtotal);
+}
+
Index: /branches/eam_branches/relastro.20100326/src/ImageOps.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/ImageOps.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/ImageOps.c	(revision 27487)
@@ -0,0 +1,767 @@
+# include "relastro.h"
+
+# define USE_IMAGE_ID 1
+
+static off_t       **bin;     // link from catalog,measure to image
+static int         **clist;   // catalog which supplied measurement on image
+static off_t       **mlist;   // measure reference for measurement on image
+static off_t        *Nlist;   // number of measurements on image
+static off_t        *NLIST;   // allocated number of measurements on image   
+
+static Image        *image;   // list of available images
+static off_t        Nimage;   // number of available images
+
+// if we search by image ID, we sort (imageIDs, imageIdx) by imageIDs to get a sorted
+// index
+
+# if USE_IMAGE_ID
+static off_t        *imageIDs; // list of all image IDs
+static off_t        *imageIdx; // list of index for image IDs 
+# else
+static unsigned int *start;
+static unsigned int *stop;
+# endif
+
+Image *getimages (off_t *N) {
+  *N = Nimage;
+  return (image);
+}
+
+Image *getimage (off_t N) {
+  return (&image[N]);
+}
+
+void initImages (Image *input, off_t N) {
+
+  off_t i;
+
+  image = input;
+  Nimage = N;
+
+# if USE_IMAGE_ID
+  ALLOCATE (imageIDs, off_t, Nimage);
+  ALLOCATE (imageIdx, off_t, Nimage);
+
+  for (i = 0; i < Nimage; i++) {
+    imageIdx[i] = i;
+    imageIDs[i] = image[i].imageID;
+  }
+  llsortpair (imageIDs, imageIdx, Nimage);
+# else
+  ALLOCATE (start, unsigned, Nimage);
+  ALLOCATE (stop, unsigned, Nimage);
+
+  for (i = 0; i < Nimage; i++) {
+    start[i] = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1);
+    stop[i]  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
+  }
+# endif
+}
+
+off_t getImageByID (off_t ID) {
+
+  // we have a pair of vectors (imageIDs, imageIdx) sorted by imageIDs
+  // use bisection to find the specified image ID
+
+# if USE_IMAGE_ID
+  off_t Nlo, Nhi, N;
+
+  Nlo = 0; Nhi = Nimage;
+  while (Nhi - Nlo > 10) {
+    N = 0.5*(Nlo + Nhi);
+    if (imageIDs[N] < ID) {
+      Nlo = MAX(N, 0);
+    } else {
+      Nhi = MIN(N + 1, Nimage);
+    }
+  }
+
+  for (N = Nlo; N < Nhi; N++) {
+    if (imageIDs[N] == ID)
+      return (imageIdx[N]);
+  }
+# endif
+
+  return (-1);
+}
+
+void initImageBins (Catalog *catalog, int Ncatalog) {
+
+  off_t i, j;
+
+  ALLOCATE (bin, off_t *, Ncatalog);
+  for (i = 0; i < Ncatalog; i++) {
+    ALLOCATE (bin[i], off_t, MAX (catalog[i].Nmeasure, 1));
+    for (j = 0; j < catalog[i].Nmeasure; j++) bin[i][j] = -1;
+  }
+
+  ALLOCATE (Nlist, off_t,   Nimage);
+  ALLOCATE (NLIST, off_t,   Nimage);
+  ALLOCATE (clist, int *,   Nimage);
+  ALLOCATE (mlist, off_t *, Nimage);
+
+  for (i = 0; i < Nimage; i++) {
+    Nlist[i] = 0;
+    NLIST[i] = 100;
+    ALLOCATE (clist[i], int, NLIST[i]);
+    ALLOCATE (mlist[i], off_t, NLIST[i]);
+  }
+}
+
+void freeImageBins (int Ncatalog) {
+
+  off_t i;
+
+  for (i = 0; i < Ncatalog; i++) {
+    free (bin[i]);
+  }
+  free (bin);
+  for (i = 0; i < Nimage; i++) {
+    free (clist[i]);
+    free (mlist[i]);
+  }
+  free (clist);
+  free (mlist);
+}
+
+/* match measurements to images */
+void findImages (Catalog *catalog, int Ncatalog) {
+
+  off_t i, j;
+  char *name;
+
+  for (i = 0; i < Ncatalog; i++) {
+    for (j = 0; j < catalog[i].Nmeasure; j++) {
+      // XXX : update to optionally restrict by photcode equiv a la relphot
+      // ecode = GetPhotcodeEquivCodebyCode (catalog[i].measure[j].photcode);
+      // if (photcode[0].code != ecode) continue;
+      matchImage (catalog, j, i);
+    }
+  }
+
+  for (i = 0; VERBOSE && (i < Nimage); i++) {
+    name = GetPhotcodeNamebyCode (image[i].photcode);
+    fprintf (stderr, "image %lld has %lld measures (%s, %s)\n", (long long) i, (long long) Nlist[i],
+             ohana_sec_to_date(image[i].tzero), name);
+  }
+}
+
+# if USE_IMAGE_ID
+// this is the imageID-based match
+void matchImage (Catalog *catalog, off_t meas, int cat) {
+
+  off_t idx, ID;
+  Measure *measure;
+
+  measure = &catalog[cat].measure[meas];
+
+  ID = measure[0].imageID;
+  idx = getImageByID (ID);
+  if (idx == -1) {
+    fprintf (stderr, "can't match detection to image?\n");
+    abort();
+  }
+
+  // index for (catalog, measure) -> image
+  bin[cat][meas] = idx;
+
+  // index for image, Nentry -> catalog
+  clist[idx][Nlist[idx]] = cat;
+
+  // index for image, Nentry -> measure
+  mlist[idx][Nlist[idx]] = meas;
+  Nlist[idx] ++;
+
+  if (Nlist[idx] == NLIST[idx]) {
+    NLIST[idx] += 100;
+    REALLOCATE (clist[idx], int,   NLIST[idx]);
+    REALLOCATE (mlist[idx], off_t, NLIST[idx]);
+  }
+  return;
+}
+
+# else
+// this is the time-based match
+void matchImage (Catalog *catalog, off_t meas, int cat) {
+
+  off_t i;
+  Measure *measure;
+
+  measure = &catalog[cat].measure[meas];
+
+  /* find the image that supplied this measurement */
+  for (i = 0; i < Nimage; i++) {
+    // let's try the very slow method first before adding a bisection search
+    // if (image[i].imageID != measure[0].imageID) continue;
+    if (image[0].photcode == -1) continue;
+    if (measure[0].photcode != image[i].photcode) continue;
+    if (measure[0].t < start[i]) continue;
+    if (measure[0].t > stop[i]) continue;
+
+    // index for (catalog, measure) -> image
+    bin[cat][meas] = i;
+
+    // index for image, Nentry -> catalog
+    clist[i][Nlist[i]] = cat;
+
+    // index for image, Nentry -> measure
+    mlist[i][Nlist[i]] = meas;
+    Nlist[i] ++;
+
+    if (Nlist[i] == NLIST[i]) {
+      NLIST[i] += 100;
+      REALLOCATE (clist[i], int,   NLIST[i]);
+      REALLOCATE (mlist[i], off_t, NLIST[i]);
+    }
+    return;
+  }
+  fprintf (stderr, "can't match detection to image?\n");
+  abort();
+}
+# endif
+
+Coords *getCoords (off_t meas, int cat) {
+
+  off_t i;
+
+  i = bin[cat][meas];
+  if (i == -1) return (NULL);
+  return (&image[i].coords);
+}
+
+void plot_images () {
+
+  off_t i, bin;
+  double *xlist, *Mlist, *dlist;
+  Graphdata graphdata;
+
+  ALLOCATE (xlist, double, Nimage);
+  ALLOCATE (dlist, double, Nimage);
+  ALLOCATE (Mlist, double, Nimage);
+
+  /**** dMcal vs airmass ****/
+  for (i = 0; i < Nimage; i++) {
+    Mlist[i] = image[i].Mcal;
+    dlist[i] = image[i].dMcal;
+    xlist[i] = image[i].secz;
+  }
+
+  plot_defaults (&graphdata);
+  graphdata.ymin = PlotdMmin;
+  graphdata.ymax = PlotdMmax;
+  plot_list (&graphdata, xlist, Mlist, Nimage, "airmass vs Mcal", "airmass.png");
+  plot_defaults (&graphdata);
+  plot_list (&graphdata, Mlist, dlist, Nimage, "Mcal vs dMcal", NULL);
+
+# define NBIN 20
+  REALLOCATE (xlist, double, NBIN);
+  REALLOCATE (Mlist, double, NBIN);
+
+  /**** dMcal histgram ****/
+  for (i = 0; i < NBIN; i++) xlist[i] = 0.00025*i;
+  bzero (Mlist, NBIN*sizeof(double));
+  for (i = 0; i < Nimage; i++) {
+    bin = image[i].dMcal / 0.00025;
+    bin = MAX (0, MIN (NBIN - 1, bin));
+    Mlist[bin] += 1.0;
+  }
+
+  plot_defaults (&graphdata);
+  graphdata.style = 1;
+  plot_list (&graphdata, xlist, Mlist, NBIN, "dMcal hist", "dMcalhist.png");
+
+  free (dlist);
+  free (xlist);
+  free (Mlist);
+}
+
+// return StarData values for detections in the specified image, converting coordinates from the
+// chip positions: X,Y -> L,M -> P,Q -> R,D
+void fixImageRaw (Catalog *catalog, int Ncatalog, off_t im) {
+
+  off_t i, m, c, n;
+  double X, Y, L, M, P, Q, R, D, dR, dD;
+
+  Mosaic *mosaic;
+  Coords *moscoords, *imcoords;
+
+  moscoords = NULL;
+  mosaic = getMosaicForImage (im);
+  if (mosaic != NULL) {
+    moscoords = &mosaic[0].coords;
+  }
+  imcoords = &image[im].coords;
+
+  for (i = 0; i < Nlist[im]; i++) {
+    m = mlist[im][i];
+    c = clist[im][i];
+
+    X = catalog[c].measure[m].Xccd;
+    Y = catalog[c].measure[m].Yccd;
+    n = catalog[c].measure[m].averef;
+
+    if (moscoords == NULL) {
+      // this is a Simple image (not a mosaic)
+      // note that for a Simple image, L,M = P,Q
+      XY_to_LM (&L, &M, X, Y, imcoords);
+      LM_to_RD (&R, &D, L, M, imcoords);
+    } else {
+      XY_to_LM (&L, &M, X, Y, imcoords);
+      XY_to_LM (&P, &Q, L, M, moscoords);
+      LM_to_RD (&R, &D, P, Q, moscoords);
+    }
+
+    // new dR, dD : test
+    dR = 3600.0*(catalog[c].average[n].R - R);
+    dD = 3600.0*(catalog[c].average[n].D - D);
+
+    if (fabs(catalog[c].measure[m].dR - dR) > 10.0) {
+      // XXXXX running into this still for last megacam exposure: wrong mosaic?
+      // ???? inconsistently hitting this????
+      fprintf (stderr, "!");
+      // abort ();
+    }
+    if (fabs(catalog[c].measure[m].dD - dD) > 10.0) {
+      fprintf (stderr, "*");
+      // abort ();
+    }
+
+    catalog[c].measure[m].dR = dR;
+    catalog[c].measure[m].dD = dD;
+
+    if (catalog[c].measure[m].dR > +180.0*3600.0) {
+      // average on high end of boundary, move star up
+      R += 360.0;
+      catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R);
+    }
+    if (catalog[c].measure[m].dR < -180.0*3600.0) {
+      // average on low end of boundary, move star down
+      R -= 360.0;
+      catalog[c].measure[m].dR = 3600.0*(catalog[c].average[n].R - R);
+    }
+  }
+  return;
+}
+
+// return StarData values for detections in the specified image, converting coordinates from the
+// chip positions: X,Y -> L,M -> P,Q -> R,D.  This function is used by the image fitting steps, for
+// which the detections have already been filtered when they were loaded (bcatalog)
+StarData *getImageRaw (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode) {
+
+  off_t i, m, c, n;
+
+  Mosaic *mosaic;
+  Coords *moscoords;
+  StarData *raw;
+
+  ALLOCATE (raw, StarData, Nlist[im]);
+
+  mosaic = NULL;
+  moscoords = NULL;
+  if (mode == MODE_MOSAIC) {
+    mosaic = getMosaicForImage (im);
+    if (mosaic == NULL) {
+      fprintf (stderr, "mosaic not found for image %s\n", image[im].name);
+      exit (1);
+    }
+    moscoords = &mosaic[0].coords;
+  }
+
+  for (i = 0; i < Nlist[im]; i++) {
+    m = mlist[im][i];
+    c = clist[im][i];
+
+    /* apply the current image transformation or use the current value of R+dR, D+dD? */
+    raw[i].X = catalog[c].measure[m].Xccd;
+    raw[i].Y = catalog[c].measure[m].Yccd;
+
+    raw[i].Mag  = catalog[c].measure[m].M;
+    raw[i].dMag = catalog[c].measure[m].dM;
+    raw[i].dPos = GetAstromError (&catalog[c].measure[m], ERROR_MODE_POS);
+
+    n = catalog[c].measure[m].averef;
+
+    // an object with only one detection provides no information about the image calibration
+    //XXX this is already taken care of in bcatalog
+    raw[i].mask = FALSE;
+    int mask = FALSE;
+    if (catalog[c].average[n].Nmeasure <= SRC_MEAS_TOOFEW) {
+      mask = TRUE;
+    }
+    if (!finite(catalog[c].measure[m].dR) || !finite(catalog[c].measure[m].dD)) {
+      mask = TRUE;
+    }
+
+    raw[i].mask = mask;
+
+
+    switch (mode) {
+      case MODE_SIMPLE:
+        /* note that for a Simple image, L,M = P,Q */
+        XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords);
+        raw[i].P = raw[i].L;
+        raw[i].Q = raw[i].M;
+        LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, &image[im].coords);
+        break;
+      case MODE_MOSAIC:
+        XY_to_LM (&raw[i].L, &raw[i].M, raw[i].X, raw[i].Y, &image[im].coords);
+        XY_to_LM (&raw[i].P, &raw[i].Q, raw[i].L, raw[i].M, moscoords);
+        LM_to_RD (&raw[i].R, &raw[i].D, raw[i].P, raw[i].Q, moscoords);
+        break;
+      default:
+	fprintf (stderr, "error: invalid mode in getImageRaw");
+	abort ();
+    }
+  }
+
+  *Nstars = Nlist[im];
+  return (raw);
+}
+
+// return StarData values for averages positions in the specified image, converting coordinates from
+// the sky positions: R,D -> P,Q -> L,M -> X,Y
+
+StarData *getImageRef (Catalog *catalog, int Ncatalog, off_t im, off_t *Nstars, CoordMode mode) {
+
+  off_t i, m, c, n;
+
+  Mosaic *mosaic;
+  Coords *moscoords;
+  StarData *ref;
+
+  ALLOCATE (ref, StarData, Nlist[im]);
+
+  mosaic = NULL;
+  moscoords = NULL;
+  if (mode == MODE_MOSAIC) {
+    mosaic = getMosaicForImage (im);
+    if (mosaic == NULL) {
+      fprintf (stderr, "mosaic not found for image %s\n", image[im].name);
+      exit (1);
+    }
+    moscoords = &mosaic[0].coords;
+  }
+
+  for (i = 0; i < Nlist[im]; i++) {
+    m = mlist[im][i];
+    c = clist[im][i];
+    n = catalog[c].measure[m].averef;
+
+    /* apply the current image transformation or use the current value of R+dR, D+dD? */
+    ref[i].R = catalog[c].average[n].R;
+    ref[i].D = catalog[c].average[n].D;
+
+    ref[i].Mag  = catalog[c].measure[m].M;
+    ref[i].dMag = catalog[c].measure[m].dM;
+    ref[i].dPos = GetAstromError (&catalog[c].measure[m], ERROR_MODE_POS);
+
+    ref[i].mask = FALSE;
+
+    /* note that for a Simple image, L,M = P,Q */
+    switch (mode) {
+      case MODE_SIMPLE:
+	RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, &image[im].coords);
+	ref[i].L = ref[i].P;
+	ref[i].M = ref[i].Q;
+	LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords);
+	break;
+      case MODE_MOSAIC:
+	RD_to_LM (&ref[i].P, &ref[i].Q, ref[i].R, ref[i].D, moscoords);
+	LM_to_XY (&ref[i].L, &ref[i].M, ref[i].P, ref[i].Q, moscoords);
+	LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords);
+	break;
+      default:
+        fprintf (stderr, "invalid case");
+        abort();
+    }
+  }
+
+  *Nstars = Nlist[im];
+  return (ref);
+}
+
+/** lifted from relphot/StarOps.clean_measures */
+void FlagOutliers2D(Catalog *catalog);
+
+void FlagOutliers (Catalog *catalog) {
+  // XXX FlagOutliers is just using FlagOutliers2D
+  FlagOutliers2D(catalog);
+  return;
+
+  int Ndel, Nave;
+  off_t i, j, k, m, N, Nmax, TOOFEW, Nsecfilt;
+  double Ns, theta, x, y;
+  double *R, *D, *dR, *dD;
+  StatType statsR, statsD;
+
+  Nsecfilt = GetPhotcodeNsecfilt();
+  assert(catalog[0].Nsecfilt == Nsecfilt);
+
+  if (VERBOSE) fprintf (stderr, "marking poor measures\n");
+  Nmax = 0;
+  for (i = 0; i < catalog[0].Naverage; i++) {
+    Nmax = MAX (Nmax, catalog[0].average[i].Nmeasure);
+  }
+
+  ALLOCATE (R, double, Nmax);
+  ALLOCATE (D, double, Nmax);
+  ALLOCATE (dR, double, Nmax);
+  ALLOCATE (dD, double, Nmax);
+
+  /* it makes no sense to mark 3-sigma outliers with <5 measurements */
+  TOOFEW = MAX (5, SRC_MEAS_TOOFEW);
+
+  Ns = CLIP_THRESH;
+  Ndel = Nave = 0;
+      
+  /* loop over each object in the catalog */
+  for (j = 0; j < catalog[0].Naverage; j++) {
+    
+    /* accumulate list of valid measurements */
+    m = catalog[0].average[j].measureOffset;
+    N = 0;
+    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
+      
+      // skip measurements based on user selected criteria
+      if (!MeasFilterTest(&catalog[0].measure[m])) continue;
+      R[N] = catalog[0].measure[m].dR;
+      D[N] = catalog[0].measure[m].dD;
+      dR[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_RA);
+      dD[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_DEC);
+      if (isnan(R[N]) || isnan(D[N])) continue;
+      N++;
+    }
+    if (N <= TOOFEW) continue;
+    
+    /* 3-sigma clip based on stats of inner 50% */
+    initstats ("MEAN");
+    liststats (R, dR, N, &statsR);
+    liststats (D, dD, N, &statsD);
+    
+    statsR.sigma = MAX (MIN_ERROR, statsR.sigma);
+    statsD.sigma = MAX (MIN_ERROR, statsD.sigma);
+    
+    /* compare per-object distance to this standard deviation, and flag outliers*/
+    m = catalog[0].average[j].measureOffset;
+    N = 0;
+    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
+      //reset flag on each invocation
+      catalog[0].measure[m].dbFlags &= ~ID_MEAS_POOR_ASTROM;
+
+      // skip measurements based on user selected criteria
+      if (!MeasFilterTest(&catalog[0].measure[m])) continue;
+      
+      x = catalog[0].measure[m].dR - statsR.median;
+      y = catalog[0].measure[m].dD - statsD.median;
+      theta = atan2(y,x);
+      if ((x*x + y*y) > (SQR(statsR.sigma * Ns * cos(theta)) + 
+			 SQR(statsD.sigma * Ns * sin(theta)))) {   
+	catalog[0].measure[m].dbFlags |= ID_MEAS_POOR_ASTROM;
+	Ndel++;
+      }
+      N++;
+      Nave ++;
+    }
+
+    //examine results
+    relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, 
+			       catalog[0].average[j].Nmeasure, 
+			       statsR, statsD, Ns);
+  }
+  
+  if (VERBOSE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave);
+  free (R);
+  free(dR);
+  free(D);
+  free(dD); 
+}
+
+
+/** an alternative outlier rejection scheme */
+void FlagOutliers2D (Catalog *catalog) {
+
+  int Ndel, Nave;
+  off_t i, j, k, m, N, Nmax, TOOFEW, Nsecfilt;
+  double *index;
+  double Ns, theta, x, y;
+  double *R, *D, *dR, *dD, *d2;
+  StatType statsR, statsD;
+
+  Nsecfilt = GetPhotcodeNsecfilt();
+  assert(catalog[0].Nsecfilt == Nsecfilt);
+
+  if (VERBOSE) fprintf (stderr, "marking poor measures\n");
+  Nmax = 0;
+  for (i = 0; i < catalog[0].Naverage; i++) {
+    Nmax = MAX (Nmax, catalog[0].average[i].Nmeasure);
+  }
+
+  ALLOCATE (R, double, Nmax);
+  ALLOCATE (D, double, Nmax);
+  ALLOCATE (dR, double, Nmax);
+  ALLOCATE (dD, double, Nmax);
+  ALLOCATE (d2, double, Nmax);
+  ALLOCATE (index, double, Nmax);
+
+  /* it makes no sense to mark 3-sigma outliers with <5 measurements */
+  TOOFEW = MAX (5, SRC_MEAS_TOOFEW);
+
+  Ns = CLIP_THRESH;
+  Ndel = Nave = 0;
+      
+  /* loop over each object in the catalog */
+  for (j = 0; j < catalog[0].Naverage; j++) {
+    
+    /* accumulate list of valid measurements */
+    m = catalog[0].average[j].measureOffset;
+    N = 0;
+    
+    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
+      
+      //reset flag on each invocation
+      catalog[0].measure[m].dbFlags &= ~ID_MEAS_POOR_ASTROM;
+      
+      // skip measurements based on user selected criteria
+      if (!MeasFilterTest(&catalog[0].measure[m])) continue;
+      R[N] = catalog[0].measure[m].dR;
+      D[N] = catalog[0].measure[m].dD;
+      dR[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_RA);
+      dD[N] = GetAstromError( &catalog[0].measure[m], ERROR_MODE_DEC);
+      if (isnan(R[N]) || isnan(D[N])) continue;
+      N++;
+    }
+    if (N <= TOOFEW) continue;
+    
+    /* calculate mean of all points*/
+    initstats ("MEAN");
+    liststats (R, dR, N, &statsR);
+    liststats (D, dD, N, &statsD);
+    statsR.sigma = MAX (MIN_ERROR, statsR.sigma);
+    statsD.sigma = MAX (MIN_ERROR, statsD.sigma);
+    
+    /* calculate deviations of all points*/
+    m = catalog[0].average[j].measureOffset;
+    N = 0;
+    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
+      //skip bad measurements
+      if (!MeasFilterTest(&catalog[0].measure[m])) continue;  
+      x = catalog[0].measure[m].dR - statsR.median;
+      y = catalog[0].measure[m].dD - statsD.median;
+      theta = atan2(y,x);
+      d2[N] = (x*x + y*y) / (SQR(statsR.sigma * Ns * cos(theta)) + 
+			     SQR(statsD.sigma * Ns * sin(theta)));      
+      index[N] = m;
+      N++;
+    }
+    
+    // sort d2
+    dsortpair(d2, index, N);
+    N = (N/2 > (N-1)) ? N/2 : N-1;
+
+    // recalculate image center, sigma based on closest 50% of points
+    for(k = 0;  k < N; k++) {
+      off_t ind = (off_t) index[k];
+      R[k] = catalog[0].measure[ind].dR;
+      D[k] = catalog[0].measure[ind].dD;
+      dR[k] = GetAstromError( &catalog[0].measure[ind], ERROR_MODE_RA);
+      dD[k] = GetAstromError( &catalog[0].measure[ind], ERROR_MODE_DEC);
+    }
+    liststats (R, dR, N, &statsR);
+    liststats (D, dD, N, &statsD);
+    statsR.sigma = MAX (MIN_ERROR, statsR.sigma);
+    statsD.sigma = MAX (MIN_ERROR, statsD.sigma);
+    
+    // use these new statistics to flag outliers 
+    m = catalog[0].average[j].measureOffset;
+    N = 0;
+    for (k = 0; k < catalog[0].average[j].Nmeasure; k++, m++) {
+      //skip bad measurements
+      if (!MeasFilterTest(&catalog[0].measure[m])) continue;  
+      x = catalog[0].measure[m].dR - statsR.median;
+      y = catalog[0].measure[m].dD - statsD.median;
+      theta = atan2(y,x);
+      d2[N] = (x*x + y*y) / (SQR(statsR.sigma * Ns * cos(theta)) + 
+			     SQR(statsD.sigma * Ns * sin(theta)));      
+      if ((d2[N]) > 1) {
+	catalog[0].measure[m].dbFlags |= ID_MEAS_POOR_ASTROM;
+	Ndel ++;
+      }
+      N++;
+      Nave++;
+    }  //done rejecting outliers
+
+    //examine results
+    relastroVisualPlotOutliers(catalog, catalog[0].average[j].measureOffset, 
+			       catalog[0].average[j].Nmeasure, 
+			       statsR, statsD, Ns);
+    
+  } //done looping over objects
+  
+  if (VERBOSE) fprintf (stderr, "%d measures marked poor, %d total\n", Ndel, Nave);
+  free (R);
+  free(dR);
+  free(D);
+  free(dD); 
+  free(d2);
+  free(index);
+}
+
+
+
+/** Determine whether a measurement should be included in the analysis, based on supplied filter criteria */ 
+int MeasFilterTest(Measure *measure) {
+  int found, k;
+  long mask;
+  PhotCode *code;
+  float mag;
+
+  if (!finite(measure[0].dR) || !finite(measure[0].dD)) return FALSE;
+  if (!finite(measure[0].M)) return FALSE; //XXX is this necessary for all relastro tasks?
+  
+  /* select measurements by photcode, or equiv photcode, if specified */
+  if (NphotcodesKeep > 0) {
+    found = FALSE;
+    for (k = 0; (k < NphotcodesKeep) && !found; k++) {
+      if (photcodesKeep[k][0].code == measure[0].photcode) found = TRUE;
+      if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
+    }
+    if (!found) return FALSE;
+  }
+  
+  if (NphotcodesSkip > 0) {
+    found = FALSE;
+    for (k = 0; (k < NphotcodesSkip) && !found; k++) {
+      if (photcodesSkip[k][0].code == measure[0].photcode) found = TRUE;
+      if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(measure[0].photcode)) found = TRUE;
+    }
+    if (found) return FALSE;
+  }  
+  
+  /* select measurements by time */
+  if (TimeSelect) {
+    if (measure[0].t < TSTART) return FALSE;
+    if (measure[0].t > TSTOP) return FALSE;
+  }
+  
+  /* select measurements by quality */
+  if (PhotFlagSelect) {
+    if (PhotFlagBad) {
+      mask = PhotFlagBad;
+    } else {
+      code = GetPhotcodebyCode (measure[0].photcode);
+      mask = code[0].astromBadMask;
+    }
+    if (mask & measure[0].photFlags) return FALSE;
+  }
+
+  /* select measurements by measurement error */
+  if ((SIGMA_LIM > 0) && (measure[0].dM > SIGMA_LIM)) return FALSE;
+  
+  /* select measurements by mag limit */
+  if (ImagSelect) {
+    mag = PhotInst (measure);
+    if (mag < ImagMin || mag > ImagMax) return FALSE;
+  }
+  
+  return TRUE;
+}
Index: /branches/eam_branches/relastro.20100326/src/MosaicOps.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/MosaicOps.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/MosaicOps.c	(revision 27487)
@@ -0,0 +1,246 @@
+# include "relastro.h"
+
+// array of mosaic definition structures
+static off_t   Nmosaic;
+static Mosaic *mosaic;
+
+// list of all images associated with a mosaic
+static off_t   *Nmosaic_own_images; // number of images for this mosaic
+static off_t   *Amosaic_own_images; // size of allocated array
+static off_t   **mosaic_own_images; // array of arrays: mosaic -> images
+
+// list of mosaic associated with each image  
+static off_t    Nmosaic_for_images; // number of images (for off_ternal checks)
+static off_t    *mosaic_for_images; // array of: image -> mosaic
+
+Mosaic *getmosaics (off_t *N) {
+  *N = Nmosaic;
+  return (mosaic);
+}
+
+off_t getMosaicByTimes (unsigned int start, unsigned int stop, unsigned int *startMos, unsigned int *stopMos, off_t *indexMos) {
+
+  // use bisection to find the overlapping mosaic
+
+  off_t Nlo, Nhi, N;
+
+  // find the last mosaic before start
+  Nlo = 0; Nhi = Nmosaic;
+  while (Nhi - Nlo > 10) {
+    N = 0.5*(Nlo + Nhi);
+    if (startMos[N] < start) {
+      Nlo = MAX(N, 0);
+    } else {
+      Nhi = MIN(N + 1, Nmosaic);
+    }
+  }
+
+  // check for the matched mosaic starting from Nlo 
+  // we may have to go much beyond Nlo since stop is not sorted
+  // can we use a sorted version of stop to check when we are beyond the valid range??
+  for (N = Nlo; N < Nmosaic; N++) { 
+    if (stop  < stopMos[N]) continue;
+    if (start > startMos[N])  continue;
+    if (stop  < startMos[N]) return (-1);
+    return (indexMos[N]);
+  }
+
+  return (-1);
+}
+
+// sort two times vectors and an index by first time vector
+void sort_mosaic_times (unsigned int *S, unsigned int *E, off_t *I, off_t N) {
+
+# define SWAPFUNC(A,B){ unsigned int tmp_t; off_t tmp_i; \
+  tmp_t = S[A]; S[A] = S[B]; S[B] = tmp_t; \
+  tmp_t = E[A]; E[A] = E[B]; E[B] = tmp_t; \
+  tmp_i = I[A]; I[A] = I[B]; I[B] = tmp_i; \
+}
+# define COMPARE(A,B)(S[A] < S[B])
+
+  OHANA_SORT (N, COMPARE, SWAPFUNC);
+
+# undef SWAPFUNC
+# undef COMPARE
+
+}
+
+// first, let's continue to use the time to make the match, but use bracketing to make it faster:
+
+// find mosaic frames (unique time periods & photcode name matches mosaic) 
+void initMosaics (Image *image, off_t Nimage) {
+
+  off_t i, Nmos, NMOSAIC;
+  unsigned int start, stop;
+
+  unsigned int *startMos, *stopMos;
+  off_t *indexMos;
+
+  Nmosaic = 0;
+  NMOSAIC = 10;
+  ALLOCATE (mosaic, Mosaic, NMOSAIC);
+
+  ALLOCATE (Nmosaic_own_images, off_t, NMOSAIC);
+  ALLOCATE (Amosaic_own_images, off_t, NMOSAIC);
+  ALLOCATE (mosaic_own_images, off_t *, NMOSAIC);
+  ALLOCATE (startMos, unsigned int, NMOSAIC);
+  ALLOCATE (stopMos, unsigned int, NMOSAIC);
+  ALLOCATE (indexMos, off_t, NMOSAIC);
+
+  /* find the mosaic images (coords.ctype = DIS); generate list of unique mosaics */
+  for (i = 0; i < Nimage; i++) {
+    if (strcmp(&image[i].coords.ctype[4], "-DIS")) continue;
+
+    /* set image time range */
+    start = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1);
+    stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
+
+    /* a new mosaic, define ranges */
+    mosaic[Nmosaic].start = start;
+    mosaic[Nmosaic].stop  = stop;
+    mosaic[Nmosaic].Mcal  = 0.0;
+    mosaic[Nmosaic].dMcal = 0.0;
+    mosaic[Nmosaic].Xm    = 0.0;
+    mosaic[Nmosaic].flags  = image[i].flags;
+    mosaic[Nmosaic].secz  = image[i].secz;
+    mosaic[Nmosaic].coords = image[i].coords;
+
+    // init the mosaic_own_images array data
+    Nmosaic_own_images[Nmosaic] = 0;
+    Amosaic_own_images[Nmosaic] = 10;
+    ALLOCATE (mosaic_own_images[Nmosaic], off_t, Amosaic_own_images[Nmosaic]);
+
+    startMos[Nmosaic] = start;
+    stopMos[Nmosaic] = stop;
+    indexMos[Nmosaic] = Nmosaic;
+
+    Nmosaic ++;
+    if (Nmosaic == NMOSAIC) {
+      NMOSAIC += 10;
+      REALLOCATE (mosaic, Mosaic, NMOSAIC);
+      REALLOCATE (mosaic_own_images, off_t *, NMOSAIC);
+      REALLOCATE (Nmosaic_own_images, off_t, NMOSAIC);
+      REALLOCATE (Amosaic_own_images, off_t, NMOSAIC);
+      REALLOCATE (startMos, unsigned int, NMOSAIC);
+      REALLOCATE (stopMos, unsigned int, NMOSAIC);
+      REALLOCATE (indexMos, off_t, NMOSAIC);
+    }
+  }
+
+  // sort the index, start, and stop by the start times:
+  sort_mosaic_times (startMos, stopMos, indexMos, Nmosaic);
+  
+  // array to store image->mosaic index
+  Nmosaic_for_images = Nimage;
+  ALLOCATE (mosaic_for_images, off_t, Nmosaic_for_images);
+
+  /* now assign the WRP images to these mosaics */
+  for (i = 0; i < Nimage; i++) {
+    mosaic_for_images[i] = -1; // default value for no mosaic found
+
+    if (strcmp(&image[i].coords.ctype[4], "-WRP")) continue;
+
+    /* set image time range */
+    start = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1);
+    stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
+
+    Nmos = getMosaicByTimes (start, stop, startMos, stopMos, indexMos);
+    if (Nmos == -1) continue;
+
+    // mosaic corresponding to this image
+    mosaic_for_images[i] = Nmos;
+
+    // add image to mosaic_own_image list 
+    mosaic_own_images[Nmos][Nmosaic_own_images[Nmos]] = i;
+    Nmosaic_own_images[Nmos] ++;
+    if (Nmosaic_own_images[Nmos] == Amosaic_own_images[Nmos]) {
+      Amosaic_own_images[Nmos] += 10;
+      REALLOCATE (mosaic_own_images[Nmos], off_t, Amosaic_own_images[Nmos]);
+    }
+    assert (Nmosaic_own_images[Nmos] < Amosaic_own_images[Nmos]);
+  }
+
+  free (startMos);
+  free (stopMos);
+  free (indexMos);
+  return;
+}
+
+// return StarData values for detections in the specified image, converting coordinates from the
+// chip positions: X,Y -> L,M -> P,Q -> R,D
+StarData *getMosaicRaw (Catalog *catalog, int Ncatalog, off_t mos, off_t *Nstars) {
+
+  off_t i, j, im, Nraw, Nnew;
+  StarData *raw, *new;
+
+  Nraw = 0;
+  ALLOCATE (raw, StarData, 1);
+
+  // loop over the images owned by this mosaic
+  for (i = 0; i < Nmosaic_own_images[mos]; i++) {
+
+    im = mosaic_own_images[mos][i];
+    
+    // retrieve stars for this chip, applying chip and mosaic astrometry
+    // this function does the reverse-lookup for the mosaic corresponding to this image
+    new = getImageRaw (catalog, Ncatalog, im, &Nnew, MODE_MOSAIC);
+    
+    // merge new and raw
+    REALLOCATE (raw, StarData, Nraw + Nnew);
+    for (j = 0; j < Nnew; j++) {
+      raw[Nraw+j] = new[j];
+    }
+    Nraw += Nnew;
+
+    free (new);
+  }
+
+  *Nstars = Nraw;
+  return (raw);
+}
+
+// return StarData values for averages positions in the specified image, converting coordinates from
+// the sky positions: R,D -> P,Q -> L,M -> X,Y
+StarData *getMosaicRef (Catalog *catalog, int Ncatalog, off_t mos, off_t *Nstars) {
+
+  off_t i, j, im, Nref, Nnew;
+  StarData *ref, *new;
+  
+  Nref = 0;
+  ALLOCATE (ref, StarData, 1);
+
+  for (i = 0; i < Nmosaic_own_images[mos]; i++) {
+
+    im = mosaic_own_images[mos][i];
+    
+    // retrieve stars for this chip, applying chip and mosaic astrometry
+    // this function does the reverse-lookup for the mosaic corresponding to this image
+    new = getImageRef (catalog, Ncatalog, im, &Nnew, MODE_MOSAIC);
+    
+    // merge new and ref
+    REALLOCATE (ref, StarData, Nref + Nnew);
+    for (j = 0; j < Nnew; j++) {
+      ref[Nref+j] = new[j];
+    }
+    Nref += Nnew;
+
+    free (new);
+  }
+
+  *Nstars = Nref;
+  return (ref);
+}
+
+Mosaic *getMosaicForImage (off_t im) {
+
+  off_t mos;
+
+  if (im < 0) abort();
+  if (im >= Nmosaic_for_images) abort();
+
+  // search for the mosaic that 
+  mos = mosaic_for_images[im];
+  if (mos < 0) return NULL;
+
+  return &mosaic[mos];
+}
Index: /branches/eam_branches/relastro.20100326/src/ParFactor.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/ParFactor.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/ParFactor.c	(revision 27487)
@@ -0,0 +1,62 @@
+# include "relastro.h"
+
+# if (0)
+/* Low precision formulae for the sun, from Almanac p. C24 (1990) */
+/* ra and dec are returned as decimal hours and decimal degrees. */
+void lpsun (double jd, double *ra, double *dec) {
+
+  double n, L, g, lambda,epsilon,alpha,delta,x,y,z;
+
+  n = jd - J2000;
+  L = 280.460 + 0.9856474 * n;
+  g = (357.528 + 0.9856003 * n)/DEG_IN_RADIAN;
+  lambda = (L + 1.915 * sin(g) + 0.020 * sin(2. * g))/DEG_IN_RADIAN;
+  epsilon = (23.439 - 0.0000004 * n)/DEG_IN_RADIAN;
+
+  // this is the conversion from ecliptic to celestial coords
+  x = cos(lambda);
+  y = cos(epsilon)*sin(lambda);
+  z = sin(epsilon)*sin(lambda);
+
+  *ra = (atan_circ(x,y))*HRS_IN_RADIAN;
+  *dec = (asin(z))*DEG_IN_RADIAN;
+}
+# endif
+
+/* code borrowed from Skycalc : fix this stuff XXX */
+/* Low precision formulae for the sun, from Almanac p. C24 (1990) */
+int sun_ecliptic (double jd, double *lambda, double *beta, double *epsilon) {
+
+  double n, L, g;
+
+# define J2000 2451545.       /* Julian date at standard epoch */
+
+  n = jd - J2000;
+  L = 280.460 + 0.9856474 * n;
+  g = (357.528 + 0.9856003 * n)*RAD_DEG;
+  *lambda = L + 1.915 * sin(g) + 0.020 * sin(2. * g); // longitude in degrees
+  *beta = 0.0;					  // approx latitude
+  *epsilon = (23.439 - 0.0000004 * n);		  // obliquity of ecliptic in degrees
+  return TRUE;
+}
+
+/* given RA, DEC, Time, calculate the parallax factor */
+int ParFactor (double *pR, double *pD, double R, double D, time_t T) {
+
+  double jd;
+  double L, B, E, e, s, r, d;
+
+  /* given a time T in UNIX seconds, determine the solar longitude S */
+
+  jd = ohana_sec_to_jd (T);
+  sun_ecliptic (jd, &L, &B, &E);
+
+  e = E * RAD_DEG;
+  s = L * RAD_DEG;
+  r = R * RAD_DEG;
+  d = D * RAD_DEG;
+  
+  *pR =  +(cos(e)*sin(s)*cos(r) - cos(s)*sin(r));
+  *pD =  -(cos(e)*sin(s)*sin(r) + cos(s)*cos(r))*sin(d) + sin(e)*sin(s)*cos(d);
+  return TRUE;
+}
Index: /branches/eam_branches/relastro.20100326/src/SetSignals.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/SetSignals.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/SetSignals.c	(revision 27487)
@@ -0,0 +1,122 @@
+# include "relastro.h"
+
+static int Protect = FALSE;
+static int Trapped = FALSE;
+
+void TrapSignal (int sig) {
+    fprintf (stderr, "trapped signal %d\n", sig);
+    if (sig == 11) {
+      fprintf (stderr, "seg fault\n");
+      exit (1);
+    }
+    if (Protect) {
+      Trapped = TRUE;
+      fprintf (stderr, "blocking until protected sections are clear\n");
+      return;
+    }
+    Shutdown ("halted by signal (trapped)");
+}    
+
+void SetProtect (int mode) {
+  Protect = mode;
+  if (Trapped && !Protect) Shutdown ("halted by signal (protect)");
+}
+
+int SetSignals () {
+
+  int i;
+
+  /* disable almost all signal interrupts */
+  for (i = 0; i < 36; i++) {
+    switch (i) {
+      /* can't redirect this signals */
+    case SIGKILL:    /* kill -9: cannot be caught or ignored */
+    case SIGSTOP:    /* SIGSTOP: cannot be caught or ignored */
+      /* ignore these signals */
+    case SIGCHLD:    /* child halted: ignore */
+    case SIGCONT:    /* continue - maintain this action */
+    case SIGTSTP:    /* stop signal sent from tty - why ignore? */
+    case SIGURG:     /* socket signal, ignore this */
+# ifdef SIGPWR
+    case SIGPWR:     /* power failure - why ignore this? (Sys V) */
+# endif
+# ifdef SIGWINCH
+    case SIGWINCH:   /* window resized (4.3BSD) */
+# endif
+      break;
+      
+    default:
+      signal (i, TrapSignal);
+    }
+  }
+  return (TRUE);
+}
+/*
+
+       Signal     Value     Action   Comment
+       -------------------------------------------------------------------------
+       SIGHUP        1        A      Hangup detected on controlling terminal
+                                     or death of controlling process
+       SIGINT        2        A      Interrupt from keyboard
+       SIGQUIT       3        A      Quit from keyboard
+       SIGILL        4        A      Illegal Instruction
+       SIGABRT       6        C      Abort signal from abort(3)
+       SIGFPE        8        C      Floating point exception
+       SIGKILL       9       AEF     Kill signal
+       SIGSEGV      11        C      Invalid memory reference
+       SIGPIPE      13        A      Broken pipe: write to pipe with no readers
+       SIGALRM      14        A      Timer signal from alarm(2)
+       SIGTERM      15        A      Termination signal
+       SIGUSR1   30,10,16     A      User-defined signal 1
+       SIGUSR2   31,12,17     A      User-defined signal 2
+       SIGCHLD   20,17,18     B      Child stopped or terminated
+       SIGCONT   19,18,25            Continue if stopped
+       SIGSTOP   17,19,23    DEF     Stop process
+       SIGTSTP   18,20,24     D      Stop typed at tty
+       SIGTTIN   21,21,26     D      tty input for background process
+       SIGTTOU   22,22,27     D      tty output for background process
+
+       Next various other signals.
+
+       Signal       Value     Action   Comment
+       ---------------------------------------------------------------------
+       SIGTRAP        5         CG     Trace/breakpoint trap
+       SIGIOT         6         CG     IOT trap. A synonym for SIGABRT
+       SIGEMT       7,-,7       G
+       SIGBUS      10,7,10      AG     Bus error
+       SIGSYS      12,-,12      G      Bad argument to routine (SVID)
+       SIGSTKFLT    -,16,-      AG     Stack fault on coprocessor
+       SIGURG      16,23,21     BG     Urgent condition on socket (4.2 BSD)
+       SIGIO       23,29,22     AG     I/O now possible (4.2 BSD)
+       SIGPOLL                  AG     A synonym for SIGIO (System V)
+       SIGCLD       -,-,18      G      A synonym for SIGCHLD
+       SIGXCPU     24,24,30     AG     CPU time limit exceeded (4.2 BSD)
+       SIGXFSZ     25,25,31     AG     File size limit exceeded (4.2 BSD)
+       SIGVTALRM   26,26,28     AG     Virtual alarm clock (4.2 BSD)
+       SIGPROF     27,27,29     AG     Profile alarm clock
+       SIGPWR      29,30,19     AG     Power failure (System V)
+       SIGINFO      29,-,-      G      A synonym for SIGPWR
+       SIGLOST      -,-,-       AG     File lock lost
+       SIGWINCH    28,28,20     BG     Window resize signal (4.3 BSD, Sun)
+       SIGUNUSED    -,31,-      AG     Unused signal
+       (Here - denotes that a signal is absent; there where three values are given, the first one is usually  valid  for  alpha  and
+       sparc,  the  middle  one  for i386 and ppc, the last one for mips. Signal 29 is SIGINFO / SIGPWR on an alpha but SIGLOST on a
+       sparc.)
+
+       The letters in the "Action" column have the following meanings:
+
+       A      Default action is to terminate the process.
+
+       B      Default action is to ignore the signal.
+
+       C      Default action is to dump core.
+
+       D      Default action is to stop the process.
+
+       E      Signal cannot be caught.
+
+       F      Signal cannot be ignored.
+
+       G      Not a POSIX.1 conformant signal.
+
+*/
Index: /branches/eam_branches/relastro.20100326/src/Shutdown.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/Shutdown.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/Shutdown.c	(revision 27487)
@@ -0,0 +1,30 @@
+# include "relastro.h"
+
+static FITS_DB *db;
+
+void set_db (FITS_DB *in) {
+  db = in;
+}
+
+/* clean up open / locked ImageCat before shutting down */
+int Shutdown (char *format, ...) {  
+  va_list argp;
+  char *formatplus;
+  
+  ALLOCATE (formatplus, char, strlen(format));
+  strcpy (formatplus, format);
+  strcat (formatplus, "\n");
+
+  va_start (argp, format);
+  vfprintf (stderr, formatplus, argp);
+  free (formatplus);
+  va_end (argp);
+
+  SetProtect (TRUE);
+  gfits_db_close (db);
+  fprintf (stderr, "ERROR: addstar halted\n");
+  exit (1);
+}
+
+
+/* XXX this is probably not needed anymore : just protect the write statements */
Index: /branches/eam_branches/relastro.20100326/src/UpdateChips.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/UpdateChips.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/UpdateChips.c	(revision 27487)
@@ -0,0 +1,36 @@
+# include "relastro.h"
+
+int UpdateChips (Catalog *catalog, int Ncatalog) {
+
+  /* we can measure new image parameters for each non-mosaic chip independently */
+  off_t i, Nimage, Nraw, Nref;
+  Image *image;
+  StarData *raw, *ref;
+
+  image = getimages (&Nimage);
+
+  for (i = 0; i < Nimage; i++) {
+
+    /* skip all except WRP images */
+    if (strcmp(&image[i].coords.ctype[4], "-WRP")) continue;
+
+    /* convert measure coordinates to raw entries */
+    raw = getImageRaw (catalog, Ncatalog, i, &Nraw, MODE_MOSAIC);
+
+    /* convert average coordinates to ref entries */
+    ref = getImageRef (catalog, Ncatalog, i, &Nref, MODE_MOSAIC);
+
+    // note that Nraw & Nref must be equal: if not, we made a programming error in one of these two functions.
+    assert (Nraw == Nref);
+
+    // FitChip does iterative, clipped fitting
+    fprintf (stderr, "image %lld : Nstars: %lld\n", (long long) i, (long long) Nraw);
+    FitChip (raw, ref, Nraw, &image[i].coords);
+
+    free (raw);
+    free (ref);
+  }
+
+  return (TRUE);
+}
+
Index: /branches/eam_branches/relastro.20100326/src/UpdateMeasures.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/UpdateMeasures.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/UpdateMeasures.c	(revision 27487)
@@ -0,0 +1,21 @@
+# include "relastro.h"
+
+int UpdateMeasures (Catalog *catalog, int Ncatalog) {
+
+  off_t i, Nimage;
+  Image *image;
+
+  image = getimages (&Nimage);
+
+  for (i = 0; i < Nimage; i++) {
+
+    /* skip DIS images (since they have no associated detections) */
+    if (!strcmp(&image[i].coords.ctype[4], "-DIS")) continue;
+
+    /* convert measure coordinates to raw entries */
+    fixImageRaw (catalog, Ncatalog, i);
+  }
+
+  return (TRUE);
+}
+
Index: /branches/eam_branches/relastro.20100326/src/UpdateMosaic.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/UpdateMosaic.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/UpdateMosaic.c	(revision 27487)
@@ -0,0 +1,29 @@
+# include "relastro.h"
+
+int UpdateMosaic (Catalog *catalog, int Ncatalog) {
+
+  /* we can measure new image parameters for each mosaic independently */
+  off_t i, Nmosaic, Nstars;
+  Mosaic *mosaic;
+  StarData *raw, *ref;
+
+  mosaic = getmosaics (&Nmosaic);
+
+  for (i = 0; i < Nmosaic; i++) {
+
+    /* convert measure coordinates to raw entries */
+    raw = getMosaicRaw (catalog, Ncatalog, i, &Nstars);
+
+    /* convert average coordinates to ref entries */
+    ref = getMosaicRef (catalog, Ncatalog, i, &Nstars);
+
+    // XXX : I'll need to supply these back to the image[] entry
+    FitMosaic (raw, ref, Nstars, &mosaic[i].coords);
+
+    free (raw);
+    free (ref);
+  }
+
+  return (TRUE);
+}
+
Index: /branches/eam_branches/relastro.20100326/src/UpdateObjectOffsets.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/UpdateObjectOffsets.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/UpdateObjectOffsets.c	(revision 27487)
@@ -0,0 +1,43 @@
+# include "relastro.h"
+
+int UpdateObjectOffsets (SkyList *skylist) {
+
+  int i;
+  Catalog catalog;
+
+  // load data from each region file, only use bright stars
+  for (i = 0; i < skylist[0].Nregions; i++) {
+
+    // set up the basic catalog info
+    catalog.filename  = skylist[0].filename[i];
+    catalog.catformat = dvo_catalog_catformat (CATFORMAT);    // set the default catformat from config data
+    catalog.catmode   = dvo_catalog_catmode (CATMODE);        // set the default catmode from config data
+    catalog.catflags  = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
+    catalog.Nsecfilt  = GetPhotcodeNsecfilt ();
+
+    if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) {
+      fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog.filename);
+      exit (1);
+    }
+    if (!catalog.Naves_disk) {
+      if (VERBOSE) fprintf (stderr, "no data in %s, skipping\n", catalog.filename);
+      dvo_catalog_unlock (&catalog);
+      dvo_catalog_free (&catalog);
+      continue;
+    }
+
+    // match measurements with images
+    initImageBins (&catalog, 1);
+    findImages (&catalog, 1);
+
+    // update the detection coordinates using the new image parameters
+    UpdateMeasures (&catalog, 1);
+
+    freeImageBins (1);
+
+    // write the updated detections to disk
+    save_catalogs (&catalog, 1);
+  }
+  
+  return (TRUE);
+}
Index: /branches/eam_branches/relastro.20100326/src/UpdateObjects.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/UpdateObjects.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/UpdateObjects.c	(revision 27487)
@@ -0,0 +1,303 @@
+# include "relastro.h"
+
+static off_t   Nmax;
+static double *X, *dX;
+static double *Y, *dY;
+static double *R, *dR;
+static double *D, *dD;
+static double *pX;
+static double *pY;
+static double *T;
+static double *dT;
+
+void initObjectData (Catalog *catalog, int Ncatalog) {
+
+  off_t i, j;
+  
+  Nmax = 0;
+  for (i = 0; i < Ncatalog; i++) {
+    for (j = 0; j < catalog[i].Naverage; j++) {
+      Nmax = MAX (Nmax, catalog[i].average[j].Nmeasure);
+    }
+  }
+
+  ALLOCATE (R, double, MAX (1, Nmax));
+  ALLOCATE (D, double, MAX (1, Nmax));
+  ALLOCATE (T, double, MAX (1, Nmax));
+  ALLOCATE (X, double, MAX (1, Nmax));
+  ALLOCATE (Y, double, MAX (1, Nmax));
+
+  ALLOCATE (dR, double, MAX (1, Nmax));
+  ALLOCATE (dD, double, MAX (1, Nmax));
+  ALLOCATE (dT, double, MAX (1, Nmax));
+  ALLOCATE (dX, double, MAX (1, Nmax));
+  ALLOCATE (dY, double, MAX (1, Nmax));
+
+  ALLOCATE (pX, double, MAX (1, Nmax));
+  ALLOCATE (pY, double, MAX (1, Nmax));
+}  
+
+int UpdateObjects (Catalog *catalog, int Ncatalog) {
+
+  int i, j, k, m, N, Nsecfilt;
+  StatType statsR, statsD;
+  Coords coords;
+  PMFit fit;
+  time_t To;
+  int mode, Nave, Npm, Npar, Nskip;
+  double Tmin, Tmax;
+
+  initObjectData (catalog, Ncatalog);
+
+  /* project coordinates to a plane centered on the object with units of arcsec */
+  coords.crval1 = 0;
+  coords.crval2 = 0;
+  coords.crpix1 = 0;
+  coords.crpix2 = 0;
+  coords.cdelt1 = coords.cdelt2 = 1.0 / 3600.0;
+  coords.pc1_1  = coords.pc2_2 = 1.0;
+  coords.pc1_2  = coords.pc2_1 = 0.0;
+  coords.Npolyterms = 1;
+  strcpy (coords.ctype, "RA---SIN");
+
+  // use J2000 as a reference time
+  To = ohana_date_to_sec ("2000/01/01");
+  Nave = Npar = Npm = 0;
+
+  // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
+  Nsecfilt = GetPhotcodeNsecfilt ();
+  assert (catalog[0].Nsecfilt == Nsecfilt);
+
+  for (i = 0; i < Ncatalog; i++) {
+
+    if (VERBOSE) fprintf (stderr, "astrometrize catalog %d : %lld ave, %lld meas\n", i, (long long) catalog[i].Naverage, (long long) catalog[i].Nmeasure);
+
+    Nskip = 0;
+    for (j = 0; j < catalog[i].Naverage; j++) {
+      /* calculate the average value of R,D for a single star */
+
+      // skip objects which are known to be problematic
+      // XXX include this code or not?
+      # if (0)
+      if (catalog[i].average[j].code & STAR_BAD) {
+	Nskip ++;
+	continue;  
+      }
+      # endif
+
+      N = 0;
+      m = catalog[i].average[j].measureOffset;
+
+      Tmin = Tmax = (catalog[i].measure[m].t - To) / (86400*365.25);
+      mode = FIT_MODE;
+
+      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
+
+	//does the measurement pass the supplied filtering constraints?
+	if (!MeasFilterTest(&catalog[i].measure[m])) {
+	  catalog[i].measure[m].dbFlags &= ~ID_MEAS_USED_OBJ;
+	  continue;
+	}
+
+	//outlier rejection
+	if (FlagOutlier && (catalog[i].measure[m].dbFlags & ID_MEAS_POOR_ASTROM)) {
+	  catalog[i].measure[m].dbFlags &= ~ID_MEAS_USED_OBJ;
+	  continue;
+	}
+
+	// exclude measurements by previous outlier detection
+	// XXX include this code or not?
+	# if (0)
+	if (catalog[i].measure[m].dbFlags & MEAS_BAD) { 
+	  catalog[i].measure[m].dbFlags |= ID_MEAS_SKIP_ASTROM;
+	  continue;
+	}
+	# endif
+
+	catalog[i].measure[m].dbFlags |= ID_MEAS_USED_OBJ;
+
+	R[N] = getMeanR (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+	D[N] = getMeanD (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+	T[N] = (catalog[i].measure[m].t - To) / (86400*365.25) ; // time relative to J2000 in years
+
+	Tmin = MIN(Tmin, T[N]);
+	Tmax = MAX(Tmax, T[N]);
+
+	dR[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_RA);
+	dD[N] = GetAstromError (&catalog[i].measure[m], ERROR_MODE_DEC);
+	dT[N] = catalog[i].measure[m].dt;
+
+	N++;
+      }
+
+      // if we have too few good detections for the desired fit, or too limited a baseline, use a
+      // fit with fewer parameters.  XXX if we have too few parameters for even the average
+      // position, consider including the lower-quality detections
+
+      catalog[i].average[j].flags &= ~ID_STAR_FEW;
+
+      // XXX add the parallax factor range as a criterion as well
+      if ((Tmax - Tmin) < PM_DT_MIN) mode = FIT_AVERAGE;
+      if ((mode == FIT_PM_ONLY) && (N < PM_TOOFEW)) mode = FIT_AVERAGE;
+
+      // too few measurements for average position (require 2 values)
+      if (N < SRC_MEAS_TOOFEW) {
+	// XXX need to define PHOTOM and ASTROM object flags
+	catalog[i].average[j].flags |= ID_STAR_FEW;
+	if (N < 2) continue;
+      }
+
+      /* we need to do the fit in a locally linear space; choose a ref coordinate */
+      coords.crval1 = R[0];
+      coords.crval2 = D[0];
+      
+      /* project all of the R,D coordinates to a plane centered on this coordinate */
+      for (k = 0; k < N; k++) {
+	RD_to_XY (&X[k], &Y[k], R[k], D[k], &coords);
+	dX[k] =  dR[k];
+	dY[k] =  dD[k];
+	// fprintf (stderr, "%d %f %f %f  %f %f\n", k, T[k], R[k], D[k], X[k], Y[k]);
+      }	  
+
+      /* fit the model components as needed */
+      switch (mode) {
+	case FIT_AVERAGE:
+	  liststats (R, dR, N, &statsR);
+	  liststats (D, dD, N, &statsD);
+
+	  fit.Ro = statsR.mean;
+	  fit.dRo = 3600.0*statsR.sigma;
+
+	  fit.Do = statsD.mean;
+	  fit.dDo = 3600.0*statsD.sigma;
+
+	  fit.chisq = 0.5*(statsR.chisq + statsD.chisq);
+	  fit.Nfit = N;
+
+	  fit.uR = fit.duR = 0.0;
+	  fit.uD = fit.duD = 0.0;
+	  fit.p  = fit.dp  = 0.0;
+
+	  Nave ++;
+	  break;
+
+	case FIT_PM_ONLY:
+	  FitPM (&fit, X, dX, Y, dY, T, N);
+	  // fprintf (stderr, "fitted:  %f - %f : %f %f : %f %f : %f\n", Tmin, Tmax, fit.Ro, fit.Do, fit.uR, fit.uD, fit.p);
+	  // project Ro, Do back to RA,DEC
+	  XY_to_RD (&fit.Ro, &fit.Do, fit.Ro, fit.Do, &coords);
+	  // fprintf (stderr, "project: %f %f : %f %f : %f\n", fit.Ro, fit.Do, fit.uR, fit.uD, fit.p);
+	  // continue;
+
+	  fit.p  = fit.dp  = 0.0;
+
+	  Npm ++;
+	  break;
+
+	case FIT_PAR_ONLY:
+	  fprintf (stderr, "programming error at %s, %d", __FILE__, __LINE__);
+	  exit (2);
+
+	  for (k = 0; k < N; k++) {
+	    ParFactor (&pX[k], &pY[k], R[k], D[k], T[k]);
+	  }
+	  FitPar (&fit, X, dX, Y, dY, pX, pY, N);
+
+	  // project Ro, Do back to RA,DEC
+	  XY_to_RD (&fit.Ro, &fit.Do, fit.Ro, fit.Do, &coords);
+
+	  fit.uR = fit.duR = 0.0;
+	  fit.uD = fit.duD = 0.0;
+
+	  Npar ++;
+	  break;
+
+	case FIT_PM_AND_PAR:
+	  fprintf (stderr, "programming error at %s, %d", __FILE__, __LINE__);
+	  exit (2);
+
+	  for (k = 0; k < N; k++) {
+	    ParFactor (&pX[k], &pY[k], R[k], D[k], T[k]);
+	  }
+	  FitPMandPar (&fit, X, dX, Y, dY, T, pX, pY, N);
+	  XY_to_RD (&fit.Ro, &fit.Do, fit.Ro, fit.Do, &coords);
+	  Npar ++;
+	  break;
+
+	default:
+	  fprintf (stderr, "programming error at %s, %d", __FILE__, __LINE__);
+	  exit (2);
+      }	  
+
+      if (0 && (j < 100)) {
+	fprintf (stderr, "%f %f -> %f %f (%f,%f)\n",
+		 catalog[i].average[j].R, 
+		 catalog[i].average[j].D, 
+		 fit.Ro, fit.Do, 
+		 3600*(catalog[i].average[j].R - fit.Ro), 
+		 3600*(catalog[i].average[j].D - fit.Do));
+      }
+
+      //make sure that the fit succeeded
+      assert(finite(fit.Ro) && finite(fit.Do) && 
+	     finite(fit.dRo) && finite(fit.dDo) &&
+	     finite(fit.uR) && finite(fit.uD) &&
+	     finite(fit.duR) && finite(fit.duD) &&
+	     finite(fit.p) && finite(fit.dp));
+
+      // the measure fields must be updated before the average fields
+      m = catalog[i].average[j].measureOffset;
+      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
+	// XXX why was this here?? if (catalog[i].measure[m].dbFlags & MEAS_BAD) continue;
+	setMeanR (fit.Ro, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+	setMeanD (fit.Do, &catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]);
+      }      
+
+      catalog[i].average[j].R  	= fit.Ro; // RA in degrees
+      catalog[i].average[j].D  	= fit.Do; // DEC in degrees
+      catalog[i].average[j].dR 	= fit.dRo; // RA scatter in arcsec
+      catalog[i].average[j].dD 	= fit.dDo; // DEC scatter in arcsec
+
+      catalog[i].average[j].uR  = fit.uR; // RA proper motion in arcsec/year
+      catalog[i].average[j].uD  = fit.uD; // DEC proper motion in arcsec/year
+      catalog[i].average[j].duR = fit.duR; // RA proper motion error in arcsec/year
+      catalog[i].average[j].duD = fit.duD; // DEC proper motion error in arcsec/year
+
+      catalog[i].average[j].P   = fit.p; // parallax in arcsec
+      catalog[i].average[j].dP  = fit.dp; // parallax error in arcsec
+
+      catalog[i].average[j].Xp  = (fit.Nfit > 1) ? 100.0*log10(fit.chisq) : NAN_S_SHORT;
+    }
+
+    if (VERBOSE) fprintf (stderr, "catalog %d : %d ave, %d pm, %d par : Nskip % d\n", i, Nave, Npm, Npar, Nskip);
+  }
+
+  if (VERBOSE) fprintf (stderr, "fitted %d objects (%d ave, %d pm, %d par)\n", Nave + Npm + Npar, Nave, Npm, Npar);
+  return (TRUE);
+}
+
+/* fitting proper-motion and parallax:
+
+given a source at position r,d, at a time t, we need to calculate a vector (pr,pd)
+
+let x,y be the coordinate in the linearized frame with y parallel to DEC lines
+
+L,B are the ecliptic longitude and latitude of the object, 
+dL and dB are the offsets in the L and B directions
+
+dL = sin(t - topp)
+dB = cos(t - topp)*sin(B)
+
+these need to be rotated to the R,D frame to yield pR,pD.  Then, the equation of motion
+for the source in the x,y frame is:
+
+x = Ro + uR * (t - to) + p * pR 
+y = Do + uD * (t - to) + p * pD
+
+the unknowns in these equations are Ro, uR, Do, uD, and p
+
+XXX think through the concepts for the pole a bit better.  all objects near the pole 
+move the same way with the same phase.  choose a projection center and define dL,dB relative 
+to that center point coordinate system?
+
+*/
Index: /branches/eam_branches/relastro.20100326/src/UpdateSimple.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/UpdateSimple.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/UpdateSimple.c	(revision 27487)
@@ -0,0 +1,32 @@
+# include "relastro.h"
+
+int UpdateSimple (Catalog *catalog, int Ncatalog) {
+
+  /* we can measure new image parameters for each non-mosaic chip independently */
+  off_t i, Nimage, Nstars;
+  Image *image;
+  StarData *raw, *ref;
+
+  image = getimages (&Nimage);
+
+  for (i = 0; i < Nimage; i++) {
+
+    /* skip WRP and DIS images */
+    if (!strcmp(&image[i].coords.ctype[4], "-WRP")) continue;
+    if (!strcmp(&image[i].coords.ctype[4], "-DIS")) continue;
+
+    /* convert measure coordinates to raw entries */
+    raw = getImageRaw (catalog, Ncatalog, i, &Nstars, MODE_SIMPLE);
+
+    /* convert average coordinates to ref entries */
+    ref = getImageRef (catalog, Ncatalog, i, &Nstars, MODE_SIMPLE);
+
+    FitSimple (raw, ref, Nstars, &image[i].coords);
+
+    free (raw);
+    free (ref);
+  }
+
+  return (TRUE);
+}
+
Index: /branches/eam_branches/relastro.20100326/src/args.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/args.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/args.c	(revision 27487)
@@ -0,0 +1,279 @@
+# include "relastro.h"
+void usage (void);
+
+int args (int argc, char **argv) {
+
+  int N;
+  double trange;
+
+  /* possible operations */
+  FIT_TARGET = TARGET_NONE;
+  FIT_MODE = FIT_AVERAGE;
+  if ((N = get_argument (argc, argv, "-update-objects"))) {
+    remove_argument (N, &argc, argv);
+    FIT_TARGET = TARGET_OBJECTS;
+
+    // check for object fitting modes (not valid for images)
+    if ((N = get_argument (argc, argv, "-pm"))) {
+	remove_argument (N, &argc, argv);
+	FIT_MODE = FIT_PM_ONLY;
+    }
+    if ((N = get_argument (argc, argv, "-par"))) {
+	remove_argument (N, &argc, argv);
+	FIT_MODE = FIT_PAR_ONLY;
+    }
+    if ((N = get_argument (argc, argv, "-pmpar"))) {
+	remove_argument (N, &argc, argv);
+	FIT_MODE = FIT_PM_AND_PAR;
+    }
+  }
+  if ((N = get_argument (argc, argv, "-update-simple"))) {
+    remove_argument (N, &argc, argv);
+    FIT_TARGET = TARGET_SIMPLE;
+  }
+  if ((N = get_argument (argc, argv, "-update-chips"))) {
+    remove_argument (N, &argc, argv);
+    FIT_TARGET = TARGET_CHIPS;
+  }
+  if ((N = get_argument (argc, argv, "-update-mosaics"))) {
+    remove_argument (N, &argc, argv);
+    FIT_TARGET = TARGET_MOSAICS;
+  }
+
+  FlagOutlier = FALSE;
+  if ((N = get_argument (argc, argv, "-clip"))) {
+    remove_argument (N, &argc, argv);
+    CLIP_THRESH = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    FlagOutlier = TRUE;
+  }
+
+  if (FIT_TARGET == TARGET_NONE) usage();
+
+  /* specify portion of the sky : allow default of all sky? */
+  UserPatch.Rmin = 0;
+  UserPatch.Rmax = 360;
+  UserPatch.Dmin = -90;
+  UserPatch.Dmax = +90;
+  if ((N = get_argument (argc, argv, "-region"))) {
+    remove_argument (N, &argc, argv);
+    UserPatch.Rmin = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    UserPatch.Rmax = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    UserPatch.Dmin = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    UserPatch.Dmax = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  } else {
+    usage ();
+  }
+
+  /* define time */
+  TimeSelect = FALSE;
+  if ((N = get_argument (argc, argv, "-time"))) {
+    remove_argument (N, &argc, argv);
+    if (!ohana_str_to_time (argv[N], &TSTART)) { 
+      fprintf (stderr, "ERROR: syntax error\n");
+      return (FALSE);
+    }
+    remove_argument (N, &argc, argv);
+    if (!ohana_str_to_dtime (argv[N], &trange)) { 
+      if (!ohana_str_to_time (argv[N], &TSTOP)) { 
+	fprintf (stderr, "ERROR: syntax error\n");
+	return (FALSE);
+      }
+    } else {
+      if (trange < 0) {
+	trange = fabs (trange);
+	TSTOP = TSTART;
+	TSTART -= trange;
+      } else {
+	TSTOP = TSTART + trange;
+      }
+    }
+    remove_argument (N, &argc, argv);
+    TimeSelect = TRUE;
+  }
+
+  PHOTCODE_KEEP_LIST = NULL;
+  if ((N = get_argument (argc, argv, "+photcode"))) {
+    remove_argument (N, &argc, argv);
+    PHOTCODE_KEEP_LIST = strcreate(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  PHOTCODE_SKIP_LIST = NULL;
+  if ((N = get_argument (argc, argv, "-photcode"))) {
+    remove_argument (N, &argc, argv);
+    PHOTCODE_SKIP_LIST = strcreate(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  VERBOSE = FALSE;
+  if ((N = get_argument (argc, argv, "-v"))) {
+    VERBOSE = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+
+  CHIPORDER = 0;
+  if ((N = get_argument (argc, argv, "-chiporder"))) {
+    remove_argument (N, &argc, argv);
+    CHIPORDER = atoi(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  SAVEPLOT = FALSE;
+  PLOTSTUFF = FALSE;
+  if ((N = get_argument (argc, argv, "-plot"))) {
+    PLOTSTUFF = TRUE;
+    remove_argument (N, &argc, argv);
+  }
+
+  PLOTDELAY = 500000;
+  if ((N = get_argument (argc, argv, "-plotdelay"))) {
+    remove_argument (N, &argc, argv);
+    PLOTDELAY = 1e6*atof(argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  // by default, require > 10pts to clip
+  strcpy (STATMODE, "CHI_INNER_80_WTMEAN");
+  if ((N = get_argument (argc, argv, "-statmode"))) {
+    remove_argument (N, &argc, argv);
+    strcpy (STATMODE, argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  RESET = FALSE;
+  if ((N = get_argument (argc, argv, "-reset"))) {
+    remove_argument (N, &argc, argv);
+    RESET = TRUE;
+  }
+
+  UPDATE = FALSE;
+  if ((N = get_argument (argc, argv, "-update"))) {
+    remove_argument (N, &argc, argv);
+    UPDATE = TRUE;
+  }
+
+  SHOW_PARAMS = FALSE;
+  if ((N = get_argument (argc, argv, "-params"))) {
+    remove_argument (N, &argc, argv);
+    SHOW_PARAMS = TRUE;
+  }
+
+  /* XXX update these for relevant plots */
+  PlotMmin = 10.0; PlotMmax = 20.0; PlotdMmin = -1.0; PlotdMmax = 1.0;
+  if ((N = get_argument (argc, argv, "-plrange"))) {
+    remove_argument (N, &argc, argv);
+    PlotMmin = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    PlotMmax = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    PlotdMmin = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    PlotdMmax = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  /* XXX update this */
+  MIN_ERROR = 0.001;
+  if ((N = get_argument (argc, argv, "-minerror"))) {
+    remove_argument (N, &argc, argv);
+    MIN_ERROR = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    /* require MIN_ERROR > 0 */
+  }  
+
+  AreaSelect = FALSE;
+  if ((N = get_argument (argc, argv, "-area"))) {
+    remove_argument (N, &argc, argv);
+    AreaXmin = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    AreaXmax = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    AreaYmin = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    AreaYmax = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    AreaSelect = TRUE;
+  }
+
+  ImagSelect = FALSE;
+  if ((N = get_argument (argc, argv, "-instmag"))) {
+    remove_argument (N, &argc, argv);
+    ImagMin = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    ImagMax = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    ImagSelect = TRUE;
+  }
+  
+  // for now, make the default to ignore the photflags
+  // XXX make it true by default instead?
+  PhotFlagSelect = FALSE;
+  if ((N = get_argument (argc, argv, "+photflags"))) {
+    remove_argument (N, &argc, argv);
+    PhotFlagSelect = TRUE;
+  }
+
+  PhotFlagBad = 0;
+  if ((N = get_argument (argc, argv, "-photflagbad"))) {
+    remove_argument (N, &argc, argv);
+    PhotFlagBad = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  PhotFlagPoor = 0;
+  if ((N = get_argument (argc, argv, "-photflagpoor"))) {
+    remove_argument (N, &argc, argv);
+    PhotFlagPoor = atoi (argv[N]);
+    remove_argument (N, &argc, argv);
+  }
+
+  /* XXX drop this? */
+  DophotSelect = FALSE;
+  if ((N = get_argument (argc, argv, "-dophot"))) {
+    remove_argument (N, &argc, argv);
+    DophotValue = atof (argv[N]);
+    remove_argument (N, &argc, argv);
+    DophotSelect = TRUE;
+  }
+
+  if (argc != 1) usage ();
+  return TRUE;
+}
+
+void usage () {
+  fprintf (stderr, "ERROR: USAGE: relastro -region RA RA DEC DEC\n");
+  fprintf (stderr, "  working options: \n");
+  fprintf (stderr, "  -update-objects\n");
+  fprintf (stderr, "    -pm\n");
+  fprintf (stderr, "    -par\n");
+  fprintf (stderr, "    -pmpar\n");
+  fprintf (stderr, "  -update-simple\n");
+  fprintf (stderr, "  -update-chips\n");
+  fprintf (stderr, "  -update-mosaics\n");
+  fprintf (stderr, "  -time (start)(stop)\n");
+  fprintf (stderr, "  +photcode (code)[,code,code...]\n");
+  fprintf (stderr, "  -photcode (code)[,code,code...]\n");
+  fprintf (stderr, "  -plot\n");
+  fprintf (stderr, "  -plotdelay (seconds)\n");
+  fprintf (stderr, "  -statmode (mode)\n");
+  fprintf (stderr, "  -reset");
+  fprintf (stderr, "  -update : apply new fit to database\n");
+  fprintf (stderr, "  -params\n");
+  fprintf (stderr, "  -plrange\n");
+  fprintf (stderr, "  -minerror\n");
+  fprintf (stderr, "  -area\n");
+  fprintf (stderr, "  -area Xmin Xmax Ymin Ymax\n");
+  fprintf (stderr, "  -instmag min max\n\n");
+  fprintf (stderr, "  +photflags\n");
+  fprintf (stderr, "  -photflagbad\n");
+  fprintf (stderr, "  -photflagpoor\n");
+  fprintf (stderr, "  -v\n");
+  fprintf (stderr, "  \n");
+  exit (2);
+} 
+
Index: /branches/eam_branches/relastro.20100326/src/bcatalog.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/bcatalog.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/bcatalog.c	(revision 27487)
@@ -0,0 +1,109 @@
+# include "relastro.h"
+
+int bcatalog (Catalog *subcatalog, Catalog *catalog) {
+
+  off_t i, j, offset;
+  off_t NAVERAGE, NMEASURE, Naverage, Nmeasure, Nm;
+  int Nsecfilt;
+
+  // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
+  Nsecfilt = GetPhotcodeNsecfilt ();
+  assert (catalog[0].Nsecfilt == Nsecfilt);
+
+  /* we are moving only the subset of measurements from catalog[0] to subcatalog[0] */
+  NAVERAGE = 50;
+  NMEASURE = 1000;
+  ALLOCATE (subcatalog[0].average, Average, NAVERAGE);
+  ALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt);
+  ALLOCATE (subcatalog[0].measure, Measure, NMEASURE);
+  Nmeasure = Naverage = 0;
+
+  /* exclude stars not in range or with too few measurements */
+  for (i = 0; i < catalog[0].Naverage; i++) {
+    if (catalog[0].average[i].Nmeasure <= SRC_MEAS_TOOFEW) {
+      continue;
+    }
+    
+    /* start with all stars good */
+    subcatalog[0].average[Naverage] = catalog[0].average[i];
+    subcatalog[0].average[Naverage].measureOffset = Nmeasure;
+    for (j = 0; j < Nsecfilt; j++) {
+      subcatalog[0].secfilt[Nsecfilt*Naverage+j] = catalog[0].secfilt[Nsecfilt*i+j];
+    }
+
+    if (RESET) {
+      // XXX reset the ra,dec coords?
+      // XXX define astrometry vs photometry average.dbFlags
+      // XXX put photometry flags in secfilt (ie, per average mag?)
+      subcatalog[0].average[Naverage].flags &= ~ID_STAR_FEW;
+      subcatalog[0].average[Naverage].flags &= ~ID_STAR_POOR;
+    }
+
+    Nm = 0;
+    for (j = 0; j < catalog[0].average[i].Nmeasure; j++) {
+
+      offset = catalog[0].average[i].measureOffset + j;
+      
+      //filter objects based on user supplied criteria
+      if (!MeasFilterTest(&catalog[0].measure[offset])) {
+	catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP;
+	continue;
+      }
+
+      //filter out outliers
+      if (FlagOutlier && (catalog[0].measure[offset].dbFlags & ID_MEAS_POOR_ASTROM)) {
+	catalog[0].measure[offset].dbFlags &= ~ID_MEAS_USED_CHIP;
+	continue;
+      }
+      catalog[0].measure[offset].dbFlags |= ID_MEAS_USED_CHIP;
+
+      // re-assess on each run of relastro if a measurement should be used
+
+      // NOCAL : this flag is used to mark measurements which must be ignored for this
+      // analysis.  an example would be time ranges or photcodes which are ignored.  NOCAL
+      // is internal only; it is not written out.
+
+      // SKIP : this flag is used to mark measurements currently ignored, but may be
+      // allowed before the analysis is complete.  for example, a star with measurements
+      // only near the edge of images will have this flag set for most of the analysis,
+      // but before the final average properties are calculated, these measurements may be
+      // allowed.
+
+      subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_SKIP_ASTROM;
+      subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_NOCAL;
+      subcatalog[0].measure[Nmeasure]          = catalog[0].measure[offset];
+      subcatalog[0].measure[Nmeasure].averef   = Naverage;
+      if (RESET) {
+        subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_POOR_ASTROM;
+        subcatalog[0].measure[Nmeasure].dbFlags &= ~ID_MEAS_AREA;
+      }
+     
+      Nmeasure ++;
+      Nm ++;
+      if (Nmeasure == NMEASURE) {
+        NMEASURE += 1000;
+        REALLOCATE (subcatalog[0].measure, Measure, NMEASURE);
+      }
+    }
+    subcatalog[0].average[Naverage].Nmeasure = Nm;
+    Naverage ++;
+    if (Naverage == NAVERAGE) {
+      NAVERAGE += 50;
+      REALLOCATE (subcatalog[0].average, Average, NAVERAGE);
+      REALLOCATE (subcatalog[0].secfilt, SecFilt, NAVERAGE*Nsecfilt);
+    }
+  }
+  REALLOCATE (subcatalog[0].average, Average, MAX (Naverage, 1));
+  REALLOCATE (subcatalog[0].measure, Measure, MAX (Nmeasure, 1));
+  REALLOCATE (subcatalog[0].secfilt, SecFilt, Nsecfilt*MAX (Naverage, 1));
+  subcatalog[0].Naverage = Naverage;
+  subcatalog[0].Nmeasure = Nmeasure;
+  subcatalog[0].Nsecfilt = catalog[0].Nsecfilt;
+  subcatalog[0].Nsecf_mem = Naverage * Nsecfilt;
+  assert (Nsecfilt == catalog[0].Nsecfilt);
+
+  if (VERBOSE) {
+    fprintf (stderr, "%lld: using %lld stars (%lld measures) for catalog\n", (long long) i, (long long) subcatalog[0].Naverage, (long long) subcatalog[0].Nmeasure);
+   }
+  return (TRUE);
+}
Index: /branches/eam_branches/relastro.20100326/src/dvo_astrom_ops.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/dvo_astrom_ops.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/dvo_astrom_ops.c	(revision 27487)
@@ -0,0 +1,91 @@
+# include "relastro.h"
+
+double getMeanR (Measure *measure, Average *average, SecFilt *secfilt) {
+
+  double ra;
+
+  /* the measure carries the instantaneous mean position at the epoch t */ 
+  ra = average[0].R - measure[0].dR / 3600.0;
+
+  /* possible corrections to mean ra:
+
+  - proper-motion and parallax
+  - abberation
+  - precession and nutation, etc
+  - refraction
+  - DCR
+
+  */
+
+  return (ra);
+}
+
+double getMeanD (Measure *measure, Average *average, SecFilt *secfilt) {
+
+  double dec;
+
+  /* the measure carries the instantaneous mean position at the epoch t */ 
+  dec = average[0].D - measure[0].dD / 3600.0;
+
+  /* possible corrections to mean ra:
+
+  - proper-motion and parallax
+  - abberation
+  - precession and nutation, etc
+  - refraction
+  - DCR
+
+  */
+
+  return (dec);
+}
+
+int setMeanR (double ra_fit, Measure *measure, Average *average, SecFilt *secfilt) {
+
+  /* math to get from new fitted position to new measure offset
+     ra_obs = average[0].R - measure[0].dR / 3600.0;
+     measure[0].dR = (ra_fit - ra_obs) * 3600.0;
+     measure[0].dR = (ra_fit - average[0].R + measure[0].dR/3600) * 3600.0
+     measure[0].dR = (ra_fit - average[0].R) * 3600.0 + measure[0].dR;
+  */
+
+  /* the measure carries the instantaneous mean position at the epoch t */ 
+  measure[0].dR += (ra_fit - average[0].R) * 3600.0;
+
+  /* possible corrections to mean ra:
+
+  - proper-motion and parallax
+  - abberation
+  - precession and nutation, etc
+  - refraction
+  - DCR
+
+  */
+
+  return (TRUE);
+}
+
+int setMeanD (double dec_fit, Measure *measure, Average *average, SecFilt *secfilt) {
+
+  /* math to get from new fitted position to new measure offset
+     dec_obs = average[0].D - measure[0].dD / 3600.0;
+     measure[0].dD = (dec_fit - dec_obs) * 3600.0;
+     measure[0].dD = (dec_fit - average[0].D + measure[0].dD/3600) * 3600.0
+     measure[0].dD = (dec_fit - average[0].D) * 3600.0 + measure[0].dD;
+  */
+
+  /* the measure carries the instantaneous mean position at the epoch t */ 
+  measure[0].dD += (dec_fit - average[0].D) * 3600.0;
+
+  /* possible corrections to mean ra:
+
+  - proper-motion and parallax
+  - abberation
+  - precession and nutation, etc
+  - refraction
+  - DCR
+
+  */
+
+  return (TRUE);
+}
Index: /branches/eam_branches/relastro.20100326/src/fitpoly.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/fitpoly.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/fitpoly.c	(revision 27487)
@@ -0,0 +1,334 @@
+# include "relastro.h"
+
+/* these functions support simultaneous 2D fits to
+   x2 = \sum a_i,j x1^i y1^j 
+   y2 = \sum b_i,j x1^i y1^j 
+
+   the order of the fit (largest coefficient) is fixed to a single
+   value for both x1,y1 terms and for x2,y2 fits
+
+   the code is currently confusing because we limit to i+j <= order.
+   this could be cleaner if we used masks and allowed i <= order, j <= order
+*/
+
+double **array_init (int Nx, int Ny) {
+
+  int i;
+  double **array;
+
+  ALLOCATE (array, double *, Nx);
+  for (i = 0; i < Nx; i++) {
+    ALLOCATE (array[i], double, Ny);
+    memset (array[i], 0, Ny*sizeof(double));
+  }    
+  return (array);
+}
+
+void array_free (double **array, int Nx) {
+
+  int i;
+
+  for (i = 0; i < Nx; i++) {
+    free (array[i]);
+  }    
+  free (array);
+}
+
+// XXX define a fit structure and drop the file static variables?
+CoordFit *fit_init (int order) {
+
+  CoordFit *fit;
+
+  ALLOCATE (fit, CoordFit, 1);
+
+  fit[0].Npts   = 0;
+  fit[0].Norder = order;
+  fit[0].Nterms = order + 1;
+  fit[0].Nsums  = 2*order + 1;
+  fit[0].Nelems = SQ(order + 1);
+
+  /* summing arrays for fit solution */
+
+  // xsum[i][j] holds \sum (x2 wt x1^i y1^j)
+  // ysum[i][j] holds \sum (y2 wt x1^i y1^j)
+  fit[0].xsum = array_init (fit[0].Nterms, fit[0].Nterms);
+  fit[0].ysum = array_init (fit[0].Nterms, fit[0].Nterms);
+
+  // xfit[i][j] holds x2 coeff for x1^i y1^j
+  // yfit[i][j] holds y2 coeff for x1^i y1^j
+  fit[0].xfit = array_init (fit[0].Nterms, fit[0].Nterms);
+  fit[0].yfit = array_init (fit[0].Nterms, fit[0].Nterms);
+
+  // sum[i][j] holds \sum (wt x1^i y1^j)
+  fit[0].sum = array_init (fit[0].Nsums, fit[0].Nsums);
+
+  return (fit);
+}
+
+void fit_free (CoordFit *fit) {
+
+  array_free (fit[0].xfit, fit[0].Nterms);
+  array_free (fit[0].yfit, fit[0].Nterms);
+
+  array_free (fit[0].sum, fit[0].Nsums);
+  array_free (fit[0].xsum, fit[0].Nterms);
+  array_free (fit[0].ysum, fit[0].Nterms);
+  
+  free (fit);
+}
+  
+// XXX use implicit masks as below or explicit masks (with function to set?)
+// XXX eg, add a global mask to this file and 
+void fit_add (CoordFit *fit, double x1, double y1, double x2, double y2, double wt) {
+
+  int ix, iy;
+  double xterm, yterm, term;
+
+  xterm = 1;
+  for (ix = 0; ix < fit[0].Nsums; ix++) {
+    yterm = 1;
+    for (iy = 0; iy < fit[0].Nsums; iy++) {
+      term = xterm*yterm*wt;
+      fit[0].sum[ix][iy] += term;
+      if ((iy < fit[0].Nterms) && (ix < fit[0].Nterms)) {
+	fit[0].xsum[ix][iy] += x2*term;
+	fit[0].ysum[ix][iy] += y2*term;
+      }
+      yterm *= y1;
+    }
+    xterm *= x1;
+  }
+  fit[0].Npts ++;
+}
+
+/* convert the xsum,ysum,sum terms into vector,matrix and solve */
+void fit_eval (CoordFit *fit) {
+
+  int i, j, ix, iy, jx, jy;
+  double **matrix, **vector;
+
+  if (fit[0].Npts == 0) {
+    fprintf (stderr, "warning: no valid pts\n");
+  }
+
+  // matrix, vector hold the final linear system
+  matrix = array_init (fit[0].Nelems, fit[0].Nelems);
+  vector = array_init (fit[0].Nelems, 2);
+
+  /* remap the xsum,ysum terms into the vector */
+  for (i = 0; i < fit[0].Nelems; i++) {
+    ix = i % fit[0].Nterms;
+    iy = i / fit[0].Nterms;
+    vector[i][0] = fit[0].xsum[ix][iy];
+    vector[i][1] = fit[0].ysum[ix][iy];
+
+    for (j = 0; j < fit[0].Nelems; j++) {
+      jx = j % fit[0].Nterms;
+      jy = j / fit[0].Nterms;
+      matrix[i][j] = fit[0].sum[ix+jx][iy+jy];
+    }
+
+    // mask the terms not represented by the Coords terms
+    if (ix + iy > fit[0].Norder) {
+      for (j = 0; j < fit[0].Nelems; j++) {
+	matrix[i][j] = 0.0;
+      }
+      vector[i][0] = 0.0;
+      vector[i][1] = 0.0;
+      matrix[i][i] = 1.0;
+    }      
+  }
+
+  for (i = 0; i < fit[0].Nelems; i++) {
+    ix = i % fit[0].Nterms;
+    iy = i / fit[0].Nterms;
+    // fprintf (stderr, "x2 : x^%dy^%d: %10.4g    y2 : x^%dy^%d: %10.4g \n", 
+    // ix, iy, vector[i][0], ix, iy, vector[i][1]);
+  }	
+
+  dgaussjordan (matrix, vector, fit[0].Nelems, 2); 
+
+  for (i = 0; i < fit[0].Nelems; i++) {
+    ix = i % fit[0].Nterms;
+    iy = i / fit[0].Nterms;
+    // fprintf (stderr, "x2 : x^%dy^%d: %10.4g    y2 : x^%dy^%d: %10.4g \n", 
+    // ix, iy, vector[i][0], ix, iy, vector[i][1]);
+  }	
+
+  /* remap the vector terms into xfit,yfit */
+  for (i = 0; i < fit[0].Nelems; i++) {
+    ix = i % fit[0].Nterms;
+    iy = i / fit[0].Nterms;
+    fit[0].xfit[ix][iy] = vector[i][0];
+    fit[0].yfit[ix][iy] = vector[i][1];
+  }
+
+  array_free (matrix, fit[0].Nelems);
+  array_free (vector, fit[0].Nelems);
+}
+
+void fit_apply (CoordFit *fit, double *x2, double *y2, double x1, double y1) {
+  
+  int ix, iy;
+  double xterm, yterm, term;
+  double x, y;
+
+  x = 0.0;
+  y = 0.0;
+
+  xterm = 1;
+  for (ix = 0; ix < fit[0].Nterms; ix++) {
+    yterm = 1;
+    for (iy = 0; iy < fit[0].Nterms; iy++) {
+      term = xterm*yterm;
+      x += fit[0].xfit[ix][iy]*term;
+      y += fit[0].yfit[ix][iy]*term;
+      yterm *= y1;
+    }
+    xterm *= x1;
+  }
+  *x2 = x;
+  *y2 = y;
+}
+
+// Nx, Ny is the number of terms in x and in y
+double **poly2d_dx (double **poly, int Nx, int Ny) {
+
+  int i, j, Nxout, Nyout;
+  double **dpoly;
+
+  Nxout = Nx - 1;
+  Nyout = Ny;
+
+  // poly[i][j] holds coeff for x1^i y1^j
+  dpoly = array_init (Nxout, Nyout);
+
+  for (i = 0; i < Nxout; i++) {
+    for (j = 0; j < Nyout; j++) {
+      dpoly[i][j] = poly[i+1][j] * (i+1);
+    }
+  }
+  return dpoly;
+}
+
+// Nx, Ny is the number of terms in x and in y
+double **poly2d_dy (double **poly, int Nx, int Ny) {
+
+  int i, j, Nxout, Nyout;
+  double **dpoly;
+
+  Nxout = Nx;
+  Nyout = Ny - 1;
+
+  // poly[i][j] holds coeff for x1^i y1^j
+  dpoly = array_init (Nxout, Nyout);
+
+  for (i = 0; i < Nxout; i++) {
+    for (j = 0; j < Nyout; j++) {
+      dpoly[i][j] = poly[i][j+1] * (j+1);
+    }
+  }
+  return dpoly;
+}
+
+// Nx, Ny is the number of terms in x and in y
+double **poly2d_copy (double **poly, int Nx, int Ny) {
+
+  int i, j;
+  double **out;
+
+  // poly[i][j] holds coeff for x1^i y1^j
+  out = array_init (Nx, Ny);
+
+  for (i = 0; i < Nx; i++) {
+    for (j = 0; j < Ny; j++) {
+      out[i][j] = poly[i][j];
+    }
+  }
+  return out;
+}
+
+double poly2d_eval (double **poly, int Nx, int Ny, double x, double y) {
+
+  int i, j;
+  double xsum, ysum, sum;
+
+  sum = 0;
+  xsum = ysum = 1.0;
+
+  for (i = 0; i < Nx; i++) {
+    ysum = xsum;
+    for (j = 0; j < Ny; j++) {
+      sum += poly[i][j] * ysum;
+      ysum *= y;
+    }
+    xsum *= x;
+  }
+  return (sum);
+}
+
+/* linear portion of fit : NORDER is 1 */
+/* this should only apply to the polynomial, not the projection terms */
+/* compare with psastro supporting code */
+CoordFit *fit_apply_coords (CoordFit *fit, Coords *coords) {
+
+  double Xo, Yo, R1, R2;
+  CoordFit *modfit;
+
+  /* I have L,M = fit(X,Y). set corresponding terms for coords */
+
+  // L = a[0][0] + a[1][0]x^1 y^0 + a[0][1] x^0 y^1 + ...
+  // L = pc1_1*cd1*(x - cp1) + pc1_2*cd2*(y - cp2) + ...
+
+  CoordsGetCenter (fit, 0.001, &Xo, &Yo);
+  coords[0].crpix1 = Xo;
+  coords[0].crpix2 = Yo;
+  
+  // resulting fit should have zero constant terms
+  modfit = CoordsSetCenter (fit, Xo, Yo);
+
+  /* we do not modify crval1,2: these are kept at the default values */
+
+  // set cdelt1, cdelt2
+  coords[0].cdelt1 = hypot (modfit[0].xfit[1][0], modfit[0].xfit[0][1]);
+  coords[0].cdelt2 = hypot (modfit[0].yfit[1][0], modfit[0].yfit[0][1]);
+  R1 = 1 / coords[0].cdelt1;
+  R2 = 1 / coords[0].cdelt2;
+
+  // set pc1_1, pc1_2, pc2_1, pc2_2 (cd1,cd2 = 1.0)
+  coords[0].pc1_1 = modfit[0].xfit[1][0] * R1;
+  coords[0].pc1_2 = modfit[0].xfit[0][1] * R2;
+  coords[0].pc2_1 = modfit[0].yfit[1][0] * R1;
+  coords[0].pc2_2 = modfit[0].yfit[0][1] * R2;
+
+  // set the polyterm elements 
+  if (coords->Npolyterms > 1) {
+    coords[0].polyterms[0][0] = modfit[0].xfit[2][0]*R1*R1;
+    coords[0].polyterms[1][0] = modfit[0].xfit[1][1]*R1*R2;
+    coords[0].polyterms[2][0] = modfit[0].xfit[0][2]*R2*R2;
+
+    coords[0].polyterms[0][1] = modfit[0].yfit[2][0]*R1*R1;
+    coords[0].polyterms[1][1] = modfit[0].yfit[1][1]*R1*R2;
+    coords[0].polyterms[2][1] = modfit[0].yfit[0][2]*R2*R2;
+  }
+
+  // I need to validate Norder
+  if (coords->Npolyterms > 2) {
+    coords[0].polyterms[3][0] = modfit[0].xfit[3][0]*R1*R1*R1;
+    coords[0].polyterms[4][0] = modfit[0].xfit[2][1]*R1*R1*R2;
+    coords[0].polyterms[5][0] = modfit[0].xfit[1][2]*R1*R2*R2;
+    coords[0].polyterms[6][0] = modfit[0].xfit[0][3]*R2*R2*R2;
+
+    coords[0].polyterms[3][1] = modfit[0].yfit[3][0]*R1*R1*R1;
+    coords[0].polyterms[4][1] = modfit[0].yfit[2][1]*R1*R1*R2;
+    coords[0].polyterms[5][1] = modfit[0].yfit[1][2]*R1*R2*R2;
+    coords[0].polyterms[6][1] = modfit[0].yfit[0][3]*R2*R2*R2;
+  }
+
+  fit_free (modfit);
+  /* keep the order and type from initial values */
+  
+  // XXX if desired in the future, return modfit (and free above)
+  return (NULL);
+}
+
Index: /branches/eam_branches/relastro.20100326/src/initialize.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/initialize.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/initialize.c	(revision 27487)
@@ -0,0 +1,78 @@
+# include "relastro.h"
+
+void initialize (int argc, char **argv) {
+
+  int NPHOTCODES;
+  char *codename, *ptr, *list;
+
+  ConfigInit (&argc, argv);
+  args (argc, argv);
+
+  /* build a list of accepted photcodes. these will be used by bcatalog to accept or
+     reject loaded data */
+
+  NphotcodesKeep = 0;
+  photcodesKeep = NULL;
+  if (PHOTCODE_KEEP_LIST != NULL) {
+    NPHOTCODES = 10;
+    ALLOCATE (photcodesKeep, PhotCode *, NPHOTCODES);
+
+    /* parse the comma-separated list of photcodesKeep */
+    list = PHOTCODE_KEEP_LIST;
+    while ((codename = strtok_r (list, ",", &ptr)) != NULL) {
+      list = NULL; // pass NULL on successive strtok_r calls
+      fprintf (stderr, "PHOTCODE_LIST: %s\n", PHOTCODE_KEEP_LIST);
+      fprintf (stderr, "codename: %s\n", codename);
+      if ((photcodesKeep[NphotcodesKeep] = GetPhotcodebyName (codename)) == NULL) {
+        fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename);
+        exit (1);
+      }
+      NphotcodesKeep ++;
+      CHECK_REALLOCATE (photcodesKeep, PhotCode *, NPHOTCODES, NphotcodesKeep, 10);
+    }
+  }
+
+  NphotcodesSkip = 0;
+  photcodesSkip = NULL;
+  if (PHOTCODE_SKIP_LIST != NULL) {
+    NPHOTCODES = 10;
+    ALLOCATE (photcodesSkip, PhotCode *, NPHOTCODES);
+
+    /* parse the comma-separated list of photcodesSkip */
+    list = PHOTCODE_SKIP_LIST;
+    while ((codename = strtok_r (list, ",", &ptr)) != NULL) {
+      list = NULL; // pass NULL on successive strtok_r calls
+      fprintf (stderr, "PHOTCODE_LIST: %s\n", PHOTCODE_SKIP_LIST);
+      fprintf (stderr, "codename: %s\n", codename);
+      if ((photcodesSkip[NphotcodesSkip] = GetPhotcodebyName (codename)) == NULL) {
+        fprintf (stderr, "ERROR: photcode %s not found in photcode table\n", codename);
+        exit (1);
+      }
+      NphotcodesSkip ++;
+      CHECK_REALLOCATE (photcodesSkip, PhotCode *, NPHOTCODES, NphotcodesSkip, 10);
+    }
+  }
+
+  initstats (STATMODE);
+
+  // IMAGE_BAD = ID_IMAGE_POOR | ID_IMAGE_FEW | ID_IMAGE_SKIP;
+  // STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
+  // MEAS_BAD  = ID_MEAS_NOCAL | ID_MEAS_POOR_ASTROM | ID_MEAS_SKIP_ASTROM | ID_MEAS_AREA;
+
+  /* XXX drop irrelevant entries */
+  if (SHOW_PARAMS) {
+    fprintf (stderr, "current parameter settings:\n");
+    if (TimeSelect) {
+      fprintf (stderr, "TimeSelect: TRUE (%s - %s)\n", ohana_sec_to_date (TSTART), ohana_sec_to_date (TSTOP));
+    } else {
+      fprintf (stderr, "TimeSelect: FALSE\n");
+    }
+    fprintf (stderr, "VERBOSE: %d, PLOTSTUFF: %d\n", VERBOSE, PLOTSTUFF);
+
+    fprintf (stderr, "IMAGE_CATALOG          %s\n",  ImageCat);
+    fprintf (stderr, "GSCFILE                %s\n",  GSCFILE);
+    fprintf (stderr, "CATDIR                 %s\n",  CATDIR);
+    exit (0);
+  }
+}
+
Index: /branches/eam_branches/relastro.20100326/src/liststats.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/liststats.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/liststats.c	(revision 27487)
@@ -0,0 +1,149 @@
+# include "relastro.h"
+
+enum {M_MEDIAN, 
+      M_MEAN, M_WTMEAN, 
+      M_INNER_MEAN, M_INNER_WTMEAN, 
+      M_INNER_80_MEAN, M_INNER_80_WTMEAN, 
+      M_CHI_INNER_MEAN, M_CHI_INNER_WTMEAN, 
+      M_CHI_INNER_80_MEAN, M_CHI_INNER_80_WTMEAN};
+
+static int statmode;
+
+void initstats (char *mode) {
+
+  statmode = -1;
+  if (!strcmp (mode, "MEAN")) statmode = M_MEAN;
+  if (!strcmp (mode, "MEDIAN")) statmode = M_MEDIAN;
+  if (!strcmp (mode, "WTMEAN")) statmode = M_WTMEAN;
+  if (!strcmp (mode, "INNER_MEAN")) statmode = M_INNER_MEAN;
+  if (!strcmp (mode, "INNER_WTMEAN")) statmode = M_INNER_WTMEAN;
+  if (!strcmp (mode, "INNER_80_MEAN")) statmode = M_INNER_80_MEAN;
+  if (!strcmp (mode, "INNER_80_WTMEAN")) statmode = M_INNER_80_WTMEAN;
+  if (!strcmp (mode, "CHI_INNER_MEAN")) statmode = M_CHI_INNER_MEAN;
+  if (!strcmp (mode, "CHI_INNER_WTMEAN")) statmode = M_CHI_INNER_WTMEAN;
+  if (!strcmp (mode, "CHI_INNER_80_MEAN")) statmode = M_CHI_INNER_80_MEAN;
+  if (!strcmp (mode, "CHI_INNER_80_WTMEAN")) statmode = M_CHI_INNER_80_WTMEAN;
+
+  if (statmode == -1) {
+    fprintf (stderr, "ERROR: invalid stats mode: %s\n", mode);
+    exit (1);
+  }
+}
+
+int liststats (double *value, double *dvalue, int N, StatType *stats) {
+  
+  int i, ks, ke, Nm;
+  double Mo, dMo, M, dM, X2, dS, *chi;
+
+  ke = ks = dMo = 0;
+
+  stats[0].Nmeas = N;
+  stats[0].mean  = 0;
+  stats[0].sigma = 0;
+  stats[0].error = 0;
+  stats[0].chisq = 0;
+  if (N < 1) return (FALSE);
+
+  dsortpair (value, dvalue, N);
+  if (N % 2) {
+    stats[0].median = value[(int)(0.5*N)];
+  } else {
+    stats[0].median = 0.5*(value[N/2] + value[N/2 - 1]);
+  }
+  stats[0].min    = value[0];
+  stats[0].max    = value[N-1];
+
+  switch (statmode) {
+  case M_MEDIAN:
+    ks = 0;
+    ke = N;
+    Mo = stats[0].median;
+    Nm = N;
+    goto chisq;
+    break;
+  case M_MEAN:
+  case M_WTMEAN:
+    ks = 0;
+    ke = N;
+    break;
+  case M_INNER_MEAN:
+  case M_INNER_WTMEAN:
+  case M_CHI_INNER_MEAN:
+  case M_CHI_INNER_WTMEAN:
+    ks = 0.25*N + 0.50;
+    ke = 0.75*N + 0.25;
+    if (N <= 3) {
+      ks = 0;
+      ke = N;
+    }
+    break;
+  case M_INNER_80_MEAN:
+  case M_INNER_80_WTMEAN:
+  case M_CHI_INNER_80_MEAN:
+  case M_CHI_INNER_80_WTMEAN:
+    ks = 0.1*N + 0.1;
+    ke = 0.8*N + 0.1;
+    if (N <= 10) {
+      ks = 0;
+      ke = N;
+    }
+    break;
+  }    
+
+  if ((statmode == M_CHI_INNER_MEAN) || 
+      (statmode == M_CHI_INNER_WTMEAN) || 
+      (statmode == M_CHI_INNER_80_MEAN) || 
+      (statmode == M_CHI_INNER_80_WTMEAN)) {
+    ALLOCATE (chi, double, N);
+    for (i = 0; i < N; i++) {
+      chi[i] = (value[i] - stats[0].median) / dvalue[i];
+    }
+    dsortthree (chi, value, dvalue, N);
+    free (chi);
+  }
+
+  /* calculating the per-star offset based on the weighted average */
+  M = dM = Nm = 0;
+  if ((statmode == M_WTMEAN) || 
+      (statmode == M_INNER_WTMEAN) || 
+      (statmode == M_INNER_80_WTMEAN) || 
+      (statmode == M_CHI_INNER_WTMEAN) || 
+      (statmode == M_CHI_INNER_80_WTMEAN)) {
+    for (i = ks; i < ke; i++) {
+      M   += value[i] / SQ (dvalue[i]);
+      dM  += 1.0 / SQ (dvalue[i]);
+      Nm  ++;  
+    }	
+    Mo = M / dM;
+    dMo = sqrt (1.0 / dM);
+  } else {
+    for (i = ks; i < ke; i++) {
+      M   += value[i];
+      dM  += SQ (dvalue[i]);
+      Nm  ++;  
+    }	
+    Mo = M / (double) Nm;
+    dMo = sqrt (dM / (double) Nm);
+  }
+
+ chisq:
+  /* find sigma and chisq */
+  X2 = dS = 0;
+  for (i = ks; i < ke; i++) {
+    M  = SQ (value[i] - Mo);
+    dM = SQ (dvalue[i]);
+    X2 += M / dM;
+    dS += M;
+  }
+  X2 = X2 / Nm;
+  dS = sqrt (dS / Nm);
+
+  stats[0].mean  = Mo;
+  stats[0].Nmeas = Nm;
+  stats[0].chisq = X2;
+  stats[0].sigma = dS;
+  stats[0].error = dMo;
+
+  return (TRUE);
+}
+
Index: /branches/eam_branches/relastro.20100326/src/load_catalogs.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/load_catalogs.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/load_catalogs.c	(revision 27487)
@@ -0,0 +1,73 @@
+# include "relastro.h"
+
+Catalog *load_catalogs (SkyList *skylist, int *Ncatalog, int subselect) {
+
+  int i, j, k, m, Nstar;
+  Catalog *catalog, *pcatalog, tcatalog;
+
+  if (VERBOSE) fprintf (stderr, "loading catalog data\n");
+
+  ALLOCATE (catalog, Catalog, skylist[0].Nregions);
+
+  // load data from each region file, only use bright stars
+  for (i = 0; i < skylist[0].Nregions; i++) {
+
+    // we only allow output if we do not use a subset.  in this case,
+    // the output parameters are correctly set for catalog[i] via pcatalog
+    pcatalog = subselect ? &tcatalog : &catalog[i];
+
+    // set up the basic catalog info
+    pcatalog[0].filename  = skylist[0].filename[i];
+    pcatalog[0].catformat = dvo_catalog_catformat (CATFORMAT);    // set the default catformat from config data
+    pcatalog[0].catmode   = dvo_catalog_catmode (CATMODE);        // set the default catmode from config data
+    pcatalog[0].catflags  = LOAD_AVES | LOAD_MEAS | LOAD_SECF; // don't need to load all data at this point
+    pcatalog[0].Nsecfilt  = GetPhotcodeNsecfilt ();
+
+    if (!dvo_catalog_open (pcatalog, skylist[0].regions[i], VERBOSE, "w")) {
+      fprintf (stderr, "ERROR: failure reading catalog %s\n", pcatalog[0].filename);
+      exit (1);
+    }
+    if (VERBOSE && !pcatalog[0].Naves_disk) fprintf (stderr, "no data in %s, skipping\n", pcatalog[0].filename);
+
+    //outlier rejection
+    if (FlagOutlier) {
+      FlagOutliers(&tcatalog);
+    }
+
+    // select only the brighter stars
+    if (subselect) {
+      bcatalog (&catalog[i], &tcatalog);
+      dvo_catalog_unlock (&tcatalog);
+      dvo_catalog_free (&tcatalog);
+    } else {
+      if (RESET) {
+	for (j = 0; j < catalog[i].Naverage; j++) {
+	  catalog[i].average[j].flags = 0;
+	  m = catalog[i].average[j].measureOffset;
+	  for (k = 0; k < catalog[i].average[j].Nmeasure; k++) {
+	    catalog[i].measure[m+k].dbFlags = 0;
+	  }
+	}
+      }
+    }
+  }
+
+  // XXX keep this test?
+  Nstar = 0;
+  for (i = 0; i < skylist[0].Nregions; i++) {
+    Nstar += catalog[i].Naverage;
+  }
+  if (Nstar < 2) { 
+    fprintf (stderr, "warning: insufficient stars %d\n", Nstar);
+  }
+
+  // XXX consider only returning the populated catalogs
+  *Ncatalog = skylist[0].Nregions;
+  return (catalog);
+}
+
+/* this function loads all relevant catalog files for the first pass.  it currently loads the data
+   read only (SOFT lock) since it assumes the image table has been locked. if we go to the new
+   addstar locking paradigm, in which the images and catalogs are updated independently, then we may
+   need to use an XCLD lock here.  
+*/
Index: /branches/eam_branches/relastro.20100326/src/load_images.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/load_images.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/load_images.c	(revision 27487)
@@ -0,0 +1,49 @@
+# include "relastro.h"
+
+# define MARKTIME(MSG,...) { \
+  float dtime; \
+  gettimeofday (&stop, (void *) NULL); \
+  dtime = DTIME (stop, start); \
+  fprintf (stderr, MSG, __VA_ARGS__); }
+
+SkyList *load_images (FITS_DB *db, SkyRegion *region) {
+
+  Image     *image, *subset;
+  off_t      Nimage, Nsubset;
+  off_t     *LineNumber;
+  struct timeval start, stop;
+
+  SkyTable *sky = NULL;
+  SkyList *skylist = NULL;
+
+  gettimeofday (&start, (void *) NULL);
+
+  // load the current sky table (layout of all SkyRegions) 
+  sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE);
+  SkyTableSetFilenames (sky, CATDIR, "cpt");
+  
+  // determine the populated SkyRegions overlapping the requested area
+  skylist = SkyListByPatch (sky, -1, region);
+  MARKTIME("  setup sky: %f sec\n", dtime);
+
+  // convert database table to internal structure
+  image = gfits_table_get_Image (&db[0].ftable, &Nimage, &db[0].swapped);
+  MARKTIME("  convert image table: %f sec\n", dtime);
+
+  // select the images which overlap the selected sky regions
+  subset = select_images (skylist, image, Nimage, &LineNumber, &Nsubset);
+  MARKTIME("  select images: %f sec\n", dtime);
+
+  gfits_vtable_from_ftable (&db[0].ftable, &db[0].vtable, LineNumber, Nsubset);
+
+  initImages (subset, Nsubset);
+  MARKTIME("  init images: %f sec\n", dtime);
+
+  initMosaics (subset, Nsubset);
+  MARKTIME("  init mosaics: %f sec\n", dtime);
+  
+  /* unlock, if we can (else, unlocked below) */
+  if (!UPDATE) dvo_image_unlock (db); 
+
+  return (skylist);
+}
Index: /branches/eam_branches/relastro.20100326/src/mkpolyterm.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/mkpolyterm.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/mkpolyterm.c	(revision 27487)
@@ -0,0 +1,166 @@
+# include "relastro.h"
+
+double factorial (int N) {
+
+  int i;
+  double F;
+
+  F = 1;
+  for (i = N; i > 1; i--) {
+    F *= i;
+  }
+  return F;
+}
+
+// XXX use a separate structure for the map (x2,y2 = f,g(x1,y1)) ?
+
+// given a 2D transformation -- L(x,y),M(x,y) -- find the coordinates x,y
+// for which L,M = 0,0. tol is the allowed error on x,y.
+int CoordsGetCenter (CoordFit *fit, double tol, double *xo, double *yo) {
+
+  int i, Nx, Ny;
+  double R, Xo, Yo, dPos, dPosRef;
+  double **XdX, **XdY, **YdX, **YdY, **alpha, **beta;
+  double **xfit, **yfit;
+
+  xfit = fit[0].xfit;
+  yfit = fit[0].yfit;
+
+  // solve for Xo,Yo where L,M(Xo,Yo) = 0,0
+  // start with linear solution for Xo,Yo:
+  R  = (xfit[1][0]*yfit[0][1] - xfit[0][1]*yfit[1][0]);
+  Xo = (yfit[0][0]*xfit[0][1] - xfit[0][0]*yfit[0][1])/R;
+  Yo = (xfit[0][0]*yfit[1][0] - yfit[0][0]*xfit[1][0])/R;
+
+  Nx = fit[0].Nterms;
+  Ny = fit[0].Nterms;
+
+  // iterate to actual solution: requires small non-linear terms
+  if (fit[0].Norder > 1) {
+    XdX = poly2d_dx (xfit, Nx, Ny);
+    XdY = poly2d_dy (xfit, Nx, Ny);
+    YdX = poly2d_dx (yfit, Nx, Ny);
+    YdY = poly2d_dy (yfit, Nx, Ny);
+
+    alpha = array_init (2, 2);
+    beta  = array_init (2, 1);
+
+    /* this loop uses the Newton-Raphson method to solve for Xo,Yo
+     * it needs the high order terms to be small 
+     * Xo,Yo are in pixels;
+     */
+    dPos = tol + 1;
+    dPosRef = dPos;
+    for (i = 0; (dPos > tol) && (i < 20); i++) {
+      // NOTE: order for alpha is: [y][x]
+      alpha[0][0] = poly2d_eval (XdX, Nx-1, Ny,   Xo, Yo);
+      alpha[1][0] = poly2d_eval (XdY, Nx,   Ny-1, Xo, Yo);
+      alpha[0][1] = poly2d_eval (YdX, Nx-1, Ny,   Xo, Yo);
+      alpha[1][1] = poly2d_eval (YdY, Nx,   Ny-1, Xo, Yo);
+
+      beta[0][0] = poly2d_eval (xfit, Nx, Ny, Xo, Yo);
+      beta[1][0] = poly2d_eval (yfit, Nx, Ny, Xo, Yo);
+
+      dgaussjordan (alpha, beta, 2, 1);
+
+      Xo -= beta[0][0];
+      Yo -= beta[1][0];
+      dPos = hypot(beta[0][0], beta[1][0]);
+      if (i == 0) {
+	dPosRef = dPos;
+      }
+    }
+    if (dPos > dPosRef) {
+      fprintf (stderr, "*** warning : non-convergence in model conversion (mkpolyterm.c;73) *** \n");
+    }
+    array_free (alpha, 2);
+    array_free (beta, 2);
+    array_free (XdX, Nx - 1);
+    array_free (XdY, Nx);
+    array_free (YdX, Nx - 1);
+    array_free (YdY, Nx);
+  }
+  *xo = Xo;
+  *yo = Yo;
+  return TRUE;
+}
+
+// convert a transformation L(x,y) to L'(x-xo,y-yo)
+CoordFit *CoordsSetCenter (CoordFit *input, double Xo, double Yo) {
+
+  int i, j, Nx, Ny;
+  double **xPx, **yPx, **xPy, **yPy, **tmp;
+
+  CoordFit *output;
+  output = fit_init (input->Norder);
+
+  /* given two equivalent polynomial representations L(x,y) = \sum_i \sum_j A_{i,j} x^i y^j
+   * we can transform L(x,y) into L'(x-xo,y-yo) by taking the derivatives of both sides and 
+   * noting that the constant term in each is the coefficient in the case of L(x,y) and is the 
+   * value of L'(-xo,-yo) in the second case.
+   */
+
+  Nx = input->Nterms;
+  Ny = input->Nterms;
+
+  xPx = poly2d_copy (input->xfit, Nx, Ny);
+  yPx = poly2d_copy (input->yfit, Nx, Ny);
+
+  for (i = 0; i < input->Nterms; i++) {
+    xPy = poly2d_copy (xPx, Nx, Ny);
+    yPy = poly2d_copy (yPx, Nx, Ny);
+    for (j = 0; j < input->Nterms; j++) {
+      output->xfit[i][j] = poly2d_eval (xPy, Nx, Ny, Xo, Yo) / factorial(i) / factorial(j);
+      output->yfit[i][j] = poly2d_eval (yPy, Nx, Ny, Xo, Yo) / factorial(i) / factorial(j);
+
+      // take the next derivative wrt y, catch output (is NULL on last pass)
+      if (Ny > 0) {
+	tmp = poly2d_dy(xPy, Nx, Ny);
+	array_free (xPy, Nx);
+	xPy = tmp;
+
+	tmp = poly2d_dy(yPy, Nx, Ny);
+	array_free (yPy, Nx);
+	yPy = tmp;
+      } else {
+	array_free (xPy, Nx);
+	array_free (yPy, Nx);
+      }
+
+      Ny --;
+    }
+    Ny = input->Nterms;
+    // take the next derivative wrt x, catch output (is NULL on last pass)
+    if (Nx > 0) {
+      tmp = poly2d_dx(xPx, Nx, Ny);
+      array_free (xPx, Nx);
+      xPx = tmp;
+      tmp = poly2d_dx(yPx, Nx, Ny);
+      array_free (yPx, Nx);
+      yPx = tmp;
+      Nx --;
+    } else {
+      array_free (xPx, Nx);
+      array_free (yPx, Nx);
+    }
+  }
+  return output;
+}
+
+/*
+  Coords uses a rigid sequence for the coefficients:
+
+  x^0 y^0 : crpix1,2 (but note this is applied before higher order terms)
+
+  x^1 y^0 : pc1_1, pc2_1
+  x^0 y^1 : pc1_2, pc2_2
+
+  x^2 y^0 : polyterm[0][0,1]
+  x^1 y^1 : polyterm[1][0,1]
+  x^0 y^2 : polyterm[2][0,1]
+
+  x^3 y^0 : polyterm[3][0,1]
+  x^2 y^1 : polyterm[4][0,1]
+  x^1 y^2 : polyterm[5][0,1]
+  x^0 y^3 : polyterm[6][0,1]
+*/
Index: /branches/eam_branches/relastro.20100326/src/plot_scatter.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/plot_scatter.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/plot_scatter.c	(revision 27487)
@@ -0,0 +1,63 @@
+# include "relastro.h"
+   
+void plot_scatter (Catalog *catalog, int Ncatalog) {
+
+# if (0)   
+  int i, j, k, m, N, Ntot, Nsecfilt;
+  float Mrel, Mcal, Mmos, Mgrid;
+  double *xlist, *ylist, *ilist;
+  Graphdata graphdata;
+
+  // XXX in the future, use catalog[0].Nsecfilt only?  allow catalogs to have variable Nsecfilt?
+  Nsecfilt = GetPhotcodeNsecfilt ();
+  assert (catalog[0].Nsecfilt == Nsecfilt);
+
+  Ntot = 0;
+  for (i = 0; i < Ncatalog; i++) {
+    for (j = 0; j < catalog[i].Naverage; j++) {
+      Ntot += catalog[i].average[j].Nmeasure;
+    }
+  }
+  ALLOCATE (xlist, double, Ntot);
+  ALLOCATE (ylist, double, Ntot);
+  ALLOCATE (ilist, double, Ntot);
+
+  N = 0;
+  for (i = 0; i < Ncatalog; i++) {
+    for (j = 0; j < catalog[i].Naverage; j++) {
+
+      /* calculate the average value for a single star */
+      if (catalog[i].average[j].code & STAR_BAD) continue;  
+      m = catalog[i].average[j].measureOffset;
+
+      for (k = 0; k < catalog[i].average[j].Nmeasure; k++, m++) {
+	if (catalog[i].measure[m].flags & MEAS_BAD) continue;
+
+	Mrel = catalog[i].secfilt[Nsecfilt*j+PhotSec].M;
+	xlist[N] = Mrel;
+	ylist[N] = PhotSys  (&catalog[i].measure[m], &catalog[i].average[j], &catalog[i].secfilt[j*Nsecfilt]) - Mcal - Mmos - Mgrid - Mrel;
+	ilist[N] = PhotInst (&catalog[i].measure[m]);
+	N++;
+      }
+    }
+  }
+
+  plot_defaults (&graphdata);
+  graphdata.xmin = PlotMmin;
+  graphdata.xmax = PlotMmax;
+  graphdata.ymin = PlotdMmin;
+  graphdata.ymax = PlotdMmax;
+  plot_list (&graphdata, xlist, ylist, N, "mag vs dmag", "Mag.png");
+
+  plot_defaults (&graphdata);
+  graphdata.ymin = PlotdMmin;
+  graphdata.ymax = PlotdMmax;
+  plot_list (&graphdata, ilist, ylist, N, "imag vs dmag", "iMag.png");
+  free (xlist);
+  free (ylist);
+  free (ilist);
+
+# endif
+}
+
+/* XXX this should become astrometrically relevant */
Index: /branches/eam_branches/relastro.20100326/src/plotstuff.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/plotstuff.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/plotstuff.c	(revision 27487)
@@ -0,0 +1,148 @@
+# include "relastro.h"
+# include <signal.h>
+
+static int Xgraph[5] = {0,0,0,0,0};
+static int active;
+
+enum {black, white, red, orange, yellow, green, blue, indigo, violet};
+
+/*
+static union { unsigned char c[4]; float f; } f_undef = { {0xff, 0xff, 0xff, 0xfe} };
+# define fUNDEF (f_undef.f)
+*/
+
+static union { unsigned char c[8]; float d; } d_undef = { {0xff, 0xff, 0xff, 0xfe, 0x00, 0x00, 0x00, 0x00} };
+# define dUNDEF (d_undef.d)
+
+void XDead () {
+  signal (SIGPIPE, XDead);
+  fprintf (stderr, "kapa is dead, must restart\n");
+  Xgraph[active] = -1;
+}
+
+int open_graph (int N) {
+
+  char name[100];
+
+  sprintf (name, "gastro [%d]", N);
+  Xgraph[N] = KapaOpen ("kapa", name);
+
+  if (Xgraph[N] < 0) {
+    fprintf (stderr, "error starting kapa\n");
+    return (FALSE);
+  }
+  return (TRUE);
+}
+
+void DonePlotting (Graphdata *graphmode, int N) {
+
+  if (Xgraph[N] == 0) return;
+
+  KapaBox (Xgraph[N], graphmode);
+  return;
+}
+
+void JpegPlot (Graphdata *graphmode, int N, char *filename) {
+
+  if (Xgraph[N] == 0) return;
+
+  KapaPNG (Xgraph[N], filename);
+  return;
+}
+
+void PSPlot (Graphdata *graphmode, int N, char *filename) {
+
+  if (Xgraph[N] == 0) return;
+
+  KiiPS (Xgraph[N], filename, TRUE, KAPA_PS_NEWPLOT, "default");
+  return;
+}
+
+void PrepPlotting (int Npts, Graphdata *graphmode, int N) {
+
+  if (Npts < 1) return;
+
+  if (Xgraph[N] < 1) if (!open_graph(N)) return;
+
+  KapaClearSections (Xgraph[N]);
+
+  KapaPrepPlot (Xgraph[N], Npts, graphmode);
+  return;
+}
+
+void PlotLabel (char *string, int N) {
+
+  if (Xgraph[N] == 0) return;
+
+  KapaSendLabel (Xgraph[N], string, 2);
+}
+
+void PlotVector (int Npts, double *vect, int mode, int N, char *type) {
+
+  float *values;
+  int i;
+
+  if (Npts < 1) return;
+
+  ALLOCATE (values, float, Npts);
+  for (i = 0; i < Npts; i++) {
+    values[i] = vect[i];
+  }
+
+  KapaPlotVector (Xgraph[N], Npts, values, type);
+  free (values);
+  return;
+}
+
+void plot_list (Graphdata *graphdata, double *xlist, double *ylist, int N, char *label, char *file) {
+
+  int i;
+  StatType stats;
+  
+  stats.min = stats.max = xlist[0];
+  for (i = 0; i < N; i++) {
+    stats.min = MIN (stats.min, xlist[i]);
+    stats.max = MAX (stats.max, xlist[i]);
+  }
+  if (graphdata[0].xmin == dUNDEF) graphdata[0].xmin = 1.05*stats.min - 0.05*stats.max;
+  if (graphdata[0].xmax == dUNDEF) graphdata[0].xmax = 1.05*stats.max - 0.05*stats.min;
+
+  stats.min = stats.max = ylist[0];
+  for (i = 0; i < N; i++) {
+    stats.min = MIN (stats.min, ylist[i]);
+    stats.max = MAX (stats.max, ylist[i]);
+  }
+  if (graphdata[0].ymin == dUNDEF) graphdata[0].ymin = 1.05*stats.min - 0.05*stats.max;
+  if (graphdata[0].ymax == dUNDEF) graphdata[0].ymax = 1.05*stats.max - 0.05*stats.min;
+
+  PrepPlotting (N, graphdata, 0);
+  PlotVector (N, xlist, 0, 0, "x");
+  PlotVector (N, ylist, 1, 0, "y");
+  if (label != NULL) PlotLabel (label, 0);
+  DonePlotting (graphdata, 0);
+
+  if ((file != NULL) && SAVEPLOT) JpegPlot (graphdata, 0, file);
+  if (PLOTDELAY > 500000) {
+    fprintf (stdout, "press return\n"); 
+    fscanf (stdin, "%*c");
+  } else {
+    usleep (PLOTDELAY);
+  }
+}
+
+void plot_defaults (Graphdata *graphdata) {
+
+  graphdata[0].style = 2;
+  graphdata[0].ptype = 2;
+  graphdata[0].ltype = 0;
+  graphdata[0].etype = 0;
+  graphdata[0].color = black;
+  graphdata[0].lweight = 0;
+  graphdata[0].size = 0.5;
+
+  graphdata[0].xmin = dUNDEF;
+  graphdata[0].xmax = dUNDEF;
+  graphdata[0].ymin = dUNDEF;
+  graphdata[0].ymax = dUNDEF;
+   
+}
Index: /branches/eam_branches/relastro.20100326/src/relastro.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/relastro.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/relastro.c	(revision 27487)
@@ -0,0 +1,90 @@
+# include "relastro.h"
+
+# define MARKTIME(MSG,...) { \
+  float dtime; \
+  gettimeofday (&stop, (void *) NULL); \
+  dtime = DTIME (stop, start); \
+  fprintf (stderr, MSG, __VA_ARGS__); }
+
+int main (int argc, char **argv) {
+
+  int status, Ncatalog;
+  Catalog *catalog;
+  FITS_DB db;
+  struct timeval start, stop;
+  SkyList *skylist = NULL;
+
+  gettimeofday (&start, (void *) NULL);
+
+  /* get configuration info, args */
+  initialize (argc, argv);
+
+  /* the object analysis is a separate process iterating over catalogs */
+  if (FIT_TARGET == TARGET_OBJECTS) {
+    relastro_objects ();
+    exit (0);
+  }
+
+  /* register database handle with shutdown procedure */
+  set_db (&db);
+
+  /* lock and load the image db table */
+  status = dvo_image_lock (&db, ImageCat, 60.0, (UPDATE ? LCK_XCLD : LCK_SOFT));
+  if (!status) Shutdown ("ERROR: failure to lock image catalog %s", db.filename);
+  if (db.dbstate == LCK_EMPTY) Shutdown ("ERROR: No images in catalog %s (1)", db.filename);
+  if (!dvo_image_load (&db, VERBOSE, FALSE)) Shutdown ("can't read image catalog %s", db.filename);
+  MARKTIME("load image data: %f sec\n", dtime);
+
+  /* load regions and images based on specified sky patch (default depth) */
+  skylist = load_images (&db, &UserPatch);
+  MARKTIME("load images: %f sec\n", dtime);
+
+  /* load catalog data from region files : subselect high-quality measurements */
+  catalog = load_catalogs (skylist, &Ncatalog, TRUE);
+  MARKTIME("load catalog data: %f sec\n", dtime);
+
+  /* match measurements with images */
+  initImageBins (catalog, Ncatalog);
+  MARKTIME("make image bins: %f sec\n", dtime);
+
+  findImages (catalog, Ncatalog);
+  MARKTIME("set up image indexes: %f sec\n", dtime);
+
+  if (PLOTSTUFF) {
+    // plot_star_coords (catalog, Ncatalog);
+    // plot_mosaic_fields (catalog);
+  }
+
+  /* major modes */
+  switch (FIT_TARGET) {
+    case TARGET_SIMPLE:
+      UpdateSimple (catalog, Ncatalog);
+      break;
+
+    case TARGET_CHIPS:
+      UpdateChips (catalog, Ncatalog);
+      break;
+
+    case TARGET_MOSAICS:
+      UpdateMosaic (catalog, Ncatalog);
+      break;
+
+    default:
+      fprintf (stderr, "programming error at %s:%d", __FILE__, __LINE__);
+      exit (2);
+  }
+
+  if (!UPDATE) exit (0);
+
+  // free the image / measurement pointers
+  freeImageBins (1);
+
+  // iterate over catalogs to make detection coordinates consistant
+  UpdateObjectOffsets (skylist);
+
+  // save the updated image parameters
+  dvo_image_update (&db, VERBOSE);
+  dvo_image_unlock (&db); 
+
+  exit (0);
+}
Index: /branches/eam_branches/relastro.20100326/src/relastroVisual.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/relastroVisual.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/relastroVisual.c	(revision 27487)
@@ -0,0 +1,512 @@
+/** Diagnostic plots for relastro
+ * @author Chris Beaumont, IfA
+ */
+
+#include "kapa.h"
+#include "relastro.h"
+
+#define TESTING
+
+#define KAPAX 700
+#define KAPAY 700
+
+static int kapa = -1;
+static int kapa2 = -1;
+static int kapa3 = -1;
+
+static int isVisual = FALSE;
+static int plotRawRef = FALSE;
+static int plotScatter = FALSE;
+static int plotResid = FALSE;
+static int plotVector = FALSE;
+static int plotOutliers = TRUE;
+
+/** Spawn a kapa window */
+static int initWindow(int *kapid) {
+  if (*kapid == -1) {
+    *kapid = KapaOpenNamedSocket("kapa", "relastro_plots");
+    if (*kapid == -1) {
+      fprintf(stderr, "Failure to open kapa.\n");
+      isVisual = 0;
+      return 0;
+    }
+    //    KapaResize (*kapid, KAPAX, KAPAY);
+  }
+  return 1;
+}
+
+/** Ask the user how to proceed */
+static int askUser(int *flag) {
+  char key[10];
+  fprintf (stdout, "[c]ontinue? [s]kip the rest of these plots? [a]bort all visual plots?");
+  if (!fgets(key, 8, stdin)) {
+    fprintf(stderr, "warning: Unable to read option");
+  }
+  if (key[0] == 's') {
+    *flag = 0;
+  }
+  if (key[0] == 'a') {
+    isVisual = 0;
+  }
+  return 1;
+}
+
+/** Size graphdata to encompass all points */
+static int scaleGraphdata(float x[], float y[], Graphdata *graphdata, int n) {
+    float xlo = FLT_MAX, xhi = -FLT_MAX, ylo = FLT_MAX, yhi = -FLT_MAX;
+    int i;
+    int goodData = 0;
+    for(i = 0; i < n; i++) {
+        goodData++;
+        if(!finite(x[i]) || !finite(y[i])) continue;
+        if(x[i] < xlo) xlo = x[i];
+        if(x[i] > xhi) xhi = x[i];
+        if(y[i] < ylo) ylo = y[i];
+        if(y[i] > yhi) yhi = y[i];
+    }
+    if (goodData < 2) return 0;
+    graphdata->xmin = xlo;
+    graphdata->ymin = ylo;
+    graphdata->xmax = xhi;
+    graphdata->ymax = yhi;
+    return 1;
+}
+
+static int residPlot(float x[], float y[],
+                     float xVec[], float yVec[],
+                     int npts, int *kapaID) {
+    if (!isVisual || !plotResid) return TRUE;
+    Graphdata graphdata;
+    KapaSection section;
+
+    if (!initWindow(kapaID)) return FALSE;
+
+    KapaInitGraph(&graphdata);
+    KapaClearPlots(*kapaID);
+    KapaSetFont(*kapaID, "helvetica", 14);
+
+    section.name = "0";
+    section.x = 0.0; section.y = 0.0;
+    section.dx = .45, section.dy = .45;
+    graphdata.ptype = 7;
+    graphdata.style = 2;
+
+    KapaSetSection(*kapaID, &section);
+    if(!scaleGraphdata(x, xVec, &graphdata, npts)) return 0;
+    KapaSetLimits(*kapaID, &graphdata);
+    KapaBox(*kapaID, &graphdata);
+    KapaSendLabel(*kapaID, "x", KAPA_LABEL_XM);
+    KapaSendLabel(*kapaID, "dx", KAPA_LABEL_YM);
+    KapaPrepPlot(*kapaID, npts, &graphdata);
+    KapaPlotVector(*kapaID, npts, x, "x");
+    KapaPlotVector(*kapaID, npts, xVec, "y");
+
+    section.x = .5; section.y = 0; section.name="1";
+    KapaSetSection(*kapaID, &section);
+    if(!scaleGraphdata(x, yVec, &graphdata, npts)) return 0;
+    KapaSetLimits(*kapaID, &graphdata);
+    KapaBox(*kapaID, &graphdata);
+    KapaSendLabel(*kapaID, "x", KAPA_LABEL_XM);
+    KapaSendLabel(*kapaID, "dy", KAPA_LABEL_YP);
+    KapaPrepPlot(*kapaID, npts, &graphdata);
+    KapaPlotVector(*kapaID, npts, x, "x");
+    KapaPlotVector(*kapaID, npts, yVec, "y");
+
+    section.x = .0; section.y = .5; section.name="2";
+    KapaSetSection(*kapaID, &section);
+    if(!scaleGraphdata(y, xVec, &graphdata, npts)) return 0;;
+    KapaSetLimits(*kapaID, &graphdata);
+    KapaBox(*kapaID, &graphdata);
+    KapaSendLabel(*kapaID, "y", KAPA_LABEL_XM);
+    KapaSendLabel(*kapaID, "dx", KAPA_LABEL_YM);
+    KapaPrepPlot(*kapaID, npts, &graphdata);
+    KapaPlotVector(*kapaID, npts, y, "x");
+    KapaPlotVector(*kapaID, npts, xVec, "y");
+
+    section.x = .5; section.y = .5; section.name="3";
+    KapaSetSection(*kapaID, &section);
+    if(!scaleGraphdata(y, yVec, &graphdata, npts)) return 0;
+    KapaSetLimits(*kapaID, &graphdata);
+    KapaBox(*kapaID, &graphdata);
+    KapaSendLabel(*kapaID, "y", KAPA_LABEL_XM);
+    KapaSendLabel(*kapaID, "dy", KAPA_LABEL_YP);
+    KapaPrepPlot(*kapaID, npts, &graphdata);
+    KapaPlotVector(*kapaID, npts, y, "x");
+    KapaPlotVector(*kapaID, npts, yVec, "y");
+
+    return 1;
+}
+
+
+/**Plot a vector field*/
+static int plotVectorField(float x[], float y[],
+                           float xVec[], float yVec[],
+                           int npts, int *kapaID, double maxVecLength) {
+
+    if(!plotVector) return 1;
+
+    Graphdata graphdata;
+    float singleX[2], singleY[2];
+    float vecScaleFactor;
+    float graphSize;
+    int i;
+    char plotTitle[50];
+    sprintf(plotTitle, "Maximum Vector Size = %5.1e", maxVecLength);
+
+    if (!initWindow(kapaID)) return 0;
+
+    KapaInitGraph(&graphdata);
+    KapaClearPlots(*kapaID);
+    if(!scaleGraphdata(x, y, &graphdata, npts)) return 0;
+
+    graphSize = graphdata.xmax - graphdata.xmin;
+    if ((graphdata.ymax - graphdata.ymin) > graphSize) {
+        graphSize = graphdata.ymax - graphdata.ymin;
+    }
+
+    vecScaleFactor = graphSize * 0.02 / (float)maxVecLength;
+#ifdef TESTING
+    fprintf(stderr, "GraphSize: %e\n maxVecLength: %e\n vecScaleFactor: %e\n",
+            graphSize, maxVecLength, vecScaleFactor);
+#endif
+
+
+    KapaSetFont (*kapaID, "helvetica", 14);
+    KapaSetLimits(*kapaID, &graphdata);
+    KapaBox(*kapaID, &graphdata);
+    KapaSendLabel(*kapaID, plotTitle, KAPA_LABEL_XP);
+
+    graphdata.ptype = 7;
+    graphdata.style = 2;
+    KapaPrepPlot(*kapaID, npts, &graphdata);
+    KapaPlotVector(*kapaID, npts, x, "x");
+    KapaPlotVector(*kapaID, npts, y, "y");
+
+    //plot each vector individually
+    graphdata.ptype = 0;
+    graphdata.style = 0;
+    graphdata.color = KapaColorByName("blue");
+    for(i = 0; i < npts; i++) {
+        singleX[0] = x[i];
+        singleY[0] = y[i];
+        singleX[1] = x[i] + xVec[i] * vecScaleFactor;
+        singleY[1] = y[i] + yVec[i] * vecScaleFactor;
+        KapaPrepPlot(*kapaID, 2, &graphdata);
+        KapaPlotVector(*kapaID, 2, singleX, "x");
+        KapaPlotVector(*kapaID, 2, singleY, "y");
+    }
+    return 1;
+}
+
+int relastroVisualPlotScatter(double values[], double thresh, int npts) {
+    float *x, *data;
+    int i;
+    float xline[2], yline[2];
+    if (!isVisual || !plotScatter) return 1;
+    if (!initWindow(&kapa2)) return 0;
+
+    ALLOCATE(x, float, npts);
+    ALLOCATE(data, float, npts);
+
+    for(i = 0; i < npts; i++) {
+        x[i] = i;
+        data[i] = (float) values[i];
+    }
+
+    Graphdata graphdata;
+    KapaSection section;
+    section.x = 0; section.y = 0; section.dx = 1; section.dy = 1;
+    section.name = "junk";
+
+    KapaInitGraph(&graphdata);
+    KapaClearPlots(kapa2);
+    KapaSetSection(kapa2, &section);
+
+    graphdata.ptype = 0;
+    graphdata.style = 0;
+    graphdata.xmin = 0;
+    graphdata.xmax = npts;
+    graphdata.ymin = 0;
+    graphdata.ymax = data[npts-1];
+
+    KapaSetFont(kapa2, "helvetica", 14);
+    KapaSetLimits(kapa2, &graphdata);
+    KapaBox(kapa2, &graphdata);
+    KapaSendLabel( kapa2, "Object", KAPA_LABEL_XM);
+    KapaSendLabel( kapa2, "Offset(pixels)", KAPA_LABEL_YM);
+    KapaSendLabel( kapa2, "Astrometric Offset with fit cutoff",
+                   KAPA_LABEL_XP);
+    KapaPrepPlot(kapa2, npts, &graphdata);
+    KapaPlotVector(kapa2, npts, x, "x");
+    KapaPlotVector(kapa2, npts, data, "y");
+
+    graphdata.color = KapaColorByName("red");
+    KapaPrepPlot(kapa2, 2, &graphdata);
+    yline[0] = (float) thresh;
+    yline[1] = (float) thresh;
+    xline[0] = graphdata.xmin;
+    xline[1] = graphdata.xmax;
+    KapaPlotVector(kapa2, 2, xline, "x");
+    KapaPlotVector(kapa2, 2, yline, "y");
+
+    askUser(&plotScatter);
+    return 1;
+}
+
+
+
+/** plot raw vs ref (L, M). Only those whose distance is < drMax are used in fit*/
+int relastroVisualPlotRawRef(StarData *raw, StarData *ref, double dRmax, int numObj) {
+
+  if( !isVisual || !plotRawRef) return 1;
+  if( !initWindow(&kapa)) return 0;
+
+  float *rawX, *rawY,  *refX, *refY;
+  float *rawXfit, *rawYfit, *refXfit, *refYfit;
+  float *magRaw, *magRef, *magRawfit, *magReffit;
+  float *xVec, *yVec;
+  int numFit = 0, numNoFit = 0;
+  double dL, dM, dR;
+
+  ALLOCATE(rawX,      float, numObj);
+  ALLOCATE(rawY,      float, numObj);
+  ALLOCATE(refX,      float, numObj);
+  ALLOCATE(refY,      float, numObj);
+  ALLOCATE(rawXfit,   float, numObj);
+  ALLOCATE(rawYfit,   float, numObj);
+  ALLOCATE(refXfit,   float, numObj);
+  ALLOCATE(refYfit,   float, numObj);
+  ALLOCATE(magRaw,    float, numObj);
+  ALLOCATE(magRef,    float, numObj);
+  ALLOCATE(magRawfit, float, numObj);
+  ALLOCATE(magReffit, float, numObj);
+
+  int i;
+  for(i = 0; i < numObj; i++) {
+    if  (raw[i].mask) continue;
+
+    dL = raw[i].L - ref[i].L;
+    dM = raw[i].M - ref[i].M;
+    dR = hypot (dL, dM);
+    if (dR > dRmax) {
+      rawX[numNoFit] = raw[i].X;
+      rawY[numNoFit] = raw[i].Y;
+      refX[numNoFit] = ref[i].X;
+      refY[numNoFit] = ref[i].Y;
+      magRaw[numNoFit] = raw[i].Mag;
+      magRef[numNoFit] = ref[i].Mag;
+      numNoFit++;
+    } else {
+      rawXfit[numFit] = raw[i].X;
+      rawYfit[numFit] = raw[i].Y;
+      refXfit[numFit] = ref[i].X;
+      refYfit[numFit] = ref[i].Y;
+      magRaw[numFit] = raw[i].Mag;
+      magRef[numFit] = ref[i].Mag;
+      numFit++;
+    }
+  }
+
+  if (numFit == 0) return 0;
+
+  Graphdata graphdata;
+
+  KapaInitGraph(&graphdata);
+  KapaClearPlots(kapa);
+
+  graphdata.ptype = 7;
+  graphdata.style = 2;
+
+  if(!scaleGraphdata(rawXfit, rawYfit, &graphdata, numFit)) {
+      fprintf(stderr, "Not enough finite points for plotting");
+      return 0;
+  }
+
+  KapaSetFont(kapa, "helvetica", 14);
+  KapaSetLimits(kapa, &graphdata);
+  KapaBox(kapa, &graphdata);
+  KapaSendLabel( kapa, "X", KAPA_LABEL_XM);
+  KapaSendLabel( kapa, "Y", KAPA_LABEL_YM);
+  KapaSendLabel( kapa, "orange, red, green, blue: (raw, ref), (nofit, fit)",
+                 KAPA_LABEL_XP);
+
+  graphdata.color = KapaColorByName("orange");
+  graphdata.size = 1;
+  KapaPrepPlot(kapa, numNoFit, &graphdata);
+  KapaPlotVector(kapa, numNoFit, rawX, "x");
+  KapaPlotVector(kapa, numNoFit, rawY, "y");
+
+  graphdata.color = KapaColorByName("red");
+  graphdata.size = 2;
+  KapaPrepPlot(kapa, numNoFit, &graphdata);
+  KapaPlotVector(kapa, numNoFit, refX, "x");
+  KapaPlotVector(kapa, numNoFit, refY, "y");
+
+  graphdata.color = KapaColorByName("green");
+  graphdata.size = 1;
+  KapaPrepPlot(kapa, numFit,  &graphdata);
+  KapaPlotVector(kapa, numFit, rawXfit, "x");
+  KapaPlotVector(kapa, numFit, rawYfit, "y");
+
+  graphdata.color = KapaColorByName("blue");
+  graphdata.size = 2;
+  KapaPrepPlot(kapa, numFit, &graphdata);
+  KapaPlotVector(kapa, numFit, refXfit, "x");
+  KapaPlotVector(kapa, numFit, refYfit, "y");
+
+  ALLOCATE(xVec, float, numFit);
+  ALLOCATE(yVec, float, numFit);
+
+  //plot the fitted objects as vectors
+  for(i = 0; i < numFit; i++) {
+      xVec[i] = rawXfit[i] - refXfit[i];
+      yVec[i] = rawYfit[i] - refYfit[i];
+  }
+
+  plotVectorField(rawXfit, rawYfit, xVec, yVec, numFit, &kapa3, dRmax);
+  if(!residPlot(rawXfit, rawYfit, xVec, yVec, numFit, &kapa2)) {
+      fprintf(stderr, "Unable to plot residuals");
+      return 0;
+  }
+
+  FREE(xVec);
+  FREE(yVec);
+
+  askUser(&plotRawRef);
+
+  FREE(rawX);
+  FREE(rawY);
+  FREE(refX);
+  FREE(refY);
+  FREE(rawXfit);
+  FREE(rawYfit);
+  FREE(refXfit);
+  FREE(refYfit);
+  FREE(magRaw);
+  FREE(magRef);
+  FREE(magRawfit);
+  FREE(magReffit);
+
+  return 1;
+}
+
+int relastroVisualPlotOutliers(Catalog *catalog, int offset, int Nmeasure,
+			       StatType statsR, StatType statsD, double thresh) {
+    
+  float *Din, *Rin, *Dout, *Rout;
+  double xmin, xmax, ymin, ymax, range;
+  float xCirc[100], yCirc[100];
+  int m, i;
+  int Nin, Nout;
+  Measure meas;
+  Graphdata graphdata;
+  KapaSection section;
+  
+  if (!isVisual || !plotOutliers) return 1;
+  if (!initWindow(&kapa)) return 0;
+  
+  // populate vectors
+  ALLOCATE(Din, float, Nmeasure);
+  ALLOCATE(Rin, float, Nmeasure);
+  ALLOCATE(Dout, float, Nmeasure);
+  ALLOCATE(Rout, float, Nmeasure);
+  
+  //create the threshhold ellipse
+  for(i = 0; i < 100; i++) {
+    xCirc[i] = statsR.median + thresh * statsR.sigma * cos(2 * 3.14 / 99. * i);
+    yCirc[i] = statsD.median + thresh * statsD.sigma * sin(2 * 3.14 / 99. * i);  
+  }
+
+  m = offset;
+  Nin = Nout = 0;
+  xmin = +FLT_MAX;
+  xmax = -FLT_MAX;
+  ymax = -FLT_MAX;
+  ymin = +FLT_MAX;
+  for(i = 0; i < Nmeasure; i++, m++) {
+    meas = catalog[0].measure[m];
+    if (!MeasFilterTest(&meas)) continue;
+    xmin = MIN(xmin, meas.dR);
+    xmax = MAX(xmax, meas.dR);
+    ymin = MIN(ymin, meas.dD);
+    ymax = MAX(ymax, meas.dD);
+    
+    if (meas.dbFlags & ID_MEAS_POOR_ASTROM) {
+      Rout[Nout] = (float)(meas.dR);
+      Dout[Nout] = (float)(meas.dD);
+      fprintf(stderr, "r: %f\td: %f\t outlier: 1\n", Rout[Nout], Dout[Nout]);
+      Nout++;
+    } else {
+      Rin[Nin] = (float)(meas.dR);
+      Din[Nin] = (float)(meas.dD);
+      fprintf(stderr, "r: %f\td: %f\t outlier: 0\n", Rin[Nin], Din[Nin]);
+      Nin++;
+    }
+  }
+ 
+  range = (xmax - xmin);
+  xmin -= .1 * range;
+  xmax += .1 * range;
+  range = (ymax - ymin);
+  ymax += .1 * range;
+  ymin -= .1 * range;
+
+  //temporary fix
+  xmin = -1; xmax = 1; ymin = -1; ymax = 1;
+
+
+  //initialize graph info
+  section.x = 0; section.y = 0; section.dx = 1; section.dy = 1;
+  section.name = "junk";
+  
+  KapaInitGraph(&graphdata);
+  KapaClearPlots(kapa);
+  KapaSetFont(kapa, "helvetica", 14);
+ 
+  graphdata.ptype = 7;
+  graphdata.style = 2;
+  graphdata.size = 3;
+  graphdata.xmin = xmin;
+  graphdata.xmax = xmax;
+  graphdata.ymin = ymin;
+  graphdata.ymax = ymax;
+
+  KapaSetSection(kapa, &section);
+  KapaSetLimits(kapa, &graphdata);
+  KapaBox(kapa, &graphdata);
+
+  KapaSendLabel( kapa, "RA (arcsec)", KAPA_LABEL_XM);
+  KapaSendLabel( kapa, "Dec (arcsec)", KAPA_LABEL_YM);
+  KapaSendLabel( kapa, "Points flagged as outliers (red)",
+		 KAPA_LABEL_XP);
+
+  graphdata.color = KapaColorByName("green");
+  KapaPrepPlot(kapa, Nin, &graphdata);
+  KapaPlotVector(kapa, Nin, Rin, "x");
+  KapaPlotVector(kapa, Nin, Din, "y");
+
+  graphdata.color = KapaColorByName("red");
+  KapaPrepPlot(kapa, Nout, &graphdata);
+  KapaPlotVector(kapa, Nout, Rout, "x");
+  KapaPlotVector(kapa, Nout, Dout, "y");
+
+  graphdata.color = KapaColorByName("black");
+  graphdata.ptype = 0;
+  graphdata.style = 0;
+  KapaPrepPlot(kapa, 100, &graphdata);
+  KapaPlotVector(kapa, 100, xCirc, "x");
+  KapaPlotVector(kapa, 100, yCirc, "y");
+
+
+  
+  askUser(&plotOutliers);
+
+  FREE(Rout);
+  FREE(Dout);
+  FREE(Rin);
+  FREE(Din);
+
+  return TRUE;
+}
+  
+  
Index: /branches/eam_branches/relastro.20100326/src/relastroVisual.h
===================================================================
--- /branches/eam_branches/relastro.20100326/src/relastroVisual.h	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/relastroVisual.h	(revision 27487)
@@ -0,0 +1,9 @@
+#ifndef RELASTRO_VISUAL_H
+#define RELASTROPVISUAL_H
+
+int relastroVisualPlotRawRef(StarData *raw, StarData *ref, double dRmax, int numObj);
+int relastroVisualPlotScatter(double values[], double thresh, int npts);
+int relastroVisualPlotOutliers(Catalog *catalog, int offset, int Nmeasure, 
+			       StatType statsR, StatType statsD, double thresh);
+
+#endif
Index: /branches/eam_branches/relastro.20100326/src/relastro_objects.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/relastro_objects.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/relastro_objects.c	(revision 27487)
@@ -0,0 +1,67 @@
+# include "relastro.h"
+
+int relastro_objects () {
+
+  int i, j, k, m;
+
+  SkyTable *sky = NULL;
+  SkyList *skylist = NULL;
+  Catalog catalog;
+
+
+  // load the current sky table (layout of all SkyRegions) 
+  sky = SkyTableLoadOptimal (CATDIR, SKY_TABLE, GSCFILE, TRUE, SKY_DEPTH, VERBOSE);
+  SkyTableSetFilenames (sky, CATDIR, "cpt");
+  
+  // determine the populated SkyRegions overlapping the requested area (default depth)
+  skylist = SkyListByPatch (sky, -1, &UserPatch);
+
+  // load data from each region file, only use bright stars
+  for (i = 0; i < skylist[0].Nregions; i++) {
+
+    // set up the basic catalog info
+    catalog.filename  = skylist[0].filename[i];
+    catalog.catformat = dvo_catalog_catformat (CATFORMAT);    // set the default catformat from config data
+    catalog.catmode   = dvo_catalog_catmode (CATMODE);        // set the default catmode from config data
+    catalog.catflags  = LOAD_AVES | LOAD_MEAS | LOAD_MISS | LOAD_SECF;
+    catalog.Nsecfilt  = GetPhotcodeNsecfilt ();
+
+    if (!dvo_catalog_open (&catalog, skylist[0].regions[i], VERBOSE, "w")) {
+      fprintf (stderr, "ERROR: failure reading catalog %s\n", catalog.filename);
+      exit (1);
+    }
+    if (!catalog.Naves_disk) {
+      if (VERBOSE) fprintf (stderr, "no data in %s, skipping\n", catalog.filename);
+      dvo_catalog_unlock (&catalog);
+      dvo_catalog_free (&catalog);
+      continue;
+    }
+
+    if (FlagOutlier) {
+      FlagOutliers(&catalog);
+    }
+
+    // XXX consider what gets reset (only ASTROM flags)
+    if (RESET) {
+      for (j = 0; j < catalog.Naverage; j++) {
+	catalog.average[j].flags = 0;
+	m = catalog.average[j].measureOffset;
+	for (k = 0; k < catalog.average[j].Nmeasure; k++) {
+	  catalog.measure[m+k].dbFlags = 0;
+	}
+      }
+    }
+
+    UpdateObjects (&catalog, 1);
+
+    if (!UPDATE) {
+      dvo_catalog_unlock (&catalog);
+      dvo_catalog_free (&catalog);
+      continue;
+    }
+    
+    save_catalogs (&catalog, 1);
+  }
+  
+  return (TRUE);
+}
Index: /branches/eam_branches/relastro.20100326/src/save_catalogs.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/save_catalogs.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/save_catalogs.c	(revision 27487)
@@ -0,0 +1,15 @@
+# include "relastro.h"
+
+void save_catalogs (Catalog *catalog, int Ncatalog) {
+
+  int i;
+
+  /* load data from each region file */
+  for (i = 0; i < Ncatalog; i++) {
+
+    if (VERBOSE) fprintf (stderr, "saving catalog %s\n", catalog[i].filename);
+    dvo_catalog_save (&catalog[i], VERBOSE); 
+    dvo_catalog_unlock (&catalog[i]);
+    dvo_catalog_free (&catalog[i]);
+  }
+}
Index: /branches/eam_branches/relastro.20100326/src/select_images.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/select_images.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/select_images.c	(revision 27487)
@@ -0,0 +1,311 @@
+# include "relastro.h"
+
+/* this function returns a list of all images which overlap the given SkyList (set of
+   SkyRegions).  All images in the image catalog are tested once, so there is no check that an
+   image already has been included.  LineNum stores the locations in the Image database of the
+   list of images */
+
+typedef struct {
+  double Xc[5];
+  double Yc[5];
+  double Rc;
+  double Dc;
+} SkyRegionCoords;
+
+void dsortindex (double *X, off_t *Y, int N);
+off_t getRegionStartByRA (double R, double *Rref, off_t Nregions);
+
+Image *select_images (SkyList *skylist, Image *timage, off_t Ntimage, off_t **LineNumber, off_t *Nimage) {
+  
+  Image *image;
+  off_t i, j, k, m, nStart, iSky, nimage, NIMAGE;
+  off_t *line_number;
+  int InRange, found;
+  double Ri[5], Di[5], Xi[5], Yi[5], dx, dy;
+  Coords tcoords;
+  SkyRegionCoords *skycoords;
+  
+  double *RmaxSky;
+  off_t *index;
+
+  if (skylist[0].Nregions < 1) {
+    *Nimage = 0;
+    *LineNumber = NULL;
+    if (VERBOSE) fprintf (stderr, "no matching sky regions\n");
+    return NULL;
+  }
+
+  // the comparison is made in the catalog local projection. below we set crval1,2
+  tcoords.crpix1 = tcoords.crpix2 = 0.0;
+  tcoords.cdelt1 = tcoords.cdelt2 = 1.0 / 3600.0;
+  tcoords.pc1_1  = tcoords.pc2_2 = 1.0;
+  tcoords.pc1_2  = tcoords.pc2_1 = 0.0;
+  strcpy (tcoords.ctype, "RA---TAN");
+
+  ALLOCATE (skycoords, SkyRegionCoords, skylist[0].Nregions);
+
+  ALLOCATE (RmaxSky, double, skylist[0].Nregions);
+  ALLOCATE (index, off_t, skylist[0].Nregions);
+
+  /* compare with each region file */
+  for (i = 0; i < skylist[0].Nregions; i++) { 
+
+    /* we make positional comparisons in the projection of catalog */
+    skycoords[i].Rc = 0.5*(skylist[0].regions[i][0].Rmax + skylist[0].regions[i][0].Rmin);
+    skycoords[i].Dc = 0.5*(skylist[0].regions[i][0].Dmax + skylist[0].regions[i][0].Dmin);
+    tcoords.crval1 = skycoords[i].Rc;
+    tcoords.crval2 = skycoords[i].Dc;
+
+    /* define catalog corners */
+    RD_to_XY (&skycoords[i].Xc[0], &skycoords[i].Yc[0], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmin, &tcoords);
+    RD_to_XY (&skycoords[i].Xc[1], &skycoords[i].Yc[1], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmin, &tcoords);
+    RD_to_XY (&skycoords[i].Xc[2], &skycoords[i].Yc[2], skylist[0].regions[i][0].Rmax, skylist[0].regions[i][0].Dmax, &tcoords);
+    RD_to_XY (&skycoords[i].Xc[3], &skycoords[i].Yc[3], skylist[0].regions[i][0].Rmin, skylist[0].regions[i][0].Dmax, &tcoords);
+    skycoords[i].Xc[4] = skycoords[i].Xc[0];    
+    skycoords[i].Yc[4] = skycoords[i].Yc[0];    
+
+    RmaxSky[i] = skylist[0].regions[i][0].Rmax;
+    index[i] = i;
+
+    dx = 0.02*(skycoords[i].Xc[2] - skycoords[i].Xc[0]);
+    dy = 0.02*(skycoords[i].Yc[2] - skycoords[i].Yc[0]);
+    skycoords[i].Xc[0] -= dx; skycoords[i].Yc[0] -= dy;
+    skycoords[i].Xc[1] += dx; skycoords[i].Yc[1] -= dy;
+    skycoords[i].Xc[2] += dx; skycoords[i].Yc[2] += dy;
+    skycoords[i].Xc[3] -= dx; skycoords[i].Yc[3] += dy;
+    skycoords[i].Xc[4] -= dx; skycoords[i].Yc[4] -= dy;
+  }
+
+  dsortindex (RmaxSky, index, skylist[0].Nregions);
+
+  if (VERBOSE) fprintf (stderr, "finding images\n");
+  BuildChipMatch (timage, Ntimage);
+
+  nimage = 0;
+  NIMAGE = 100;
+  ALLOCATE (image, Image, NIMAGE);
+  ALLOCATE (line_number, off_t, NIMAGE);
+  
+  // go through the complete list of images, selecting ones which overlap any region
+  for (i = 0; i < Ntimage; i++) {
+      
+    /* select images by photcode, or equiv photcode, if specified */
+    if (NphotcodesKeep > 0) {
+      found = FALSE;
+      for (k = 0; (k < NphotcodesKeep) && !found; k++) {
+	if (photcodesKeep[k][0].code == timage[i].photcode) found = TRUE;
+	if (photcodesKeep[k][0].code == GetPhotcodeEquivCodebyCode(timage[i].photcode)) found = TRUE;
+      }
+      if (!found) continue;
+    }
+    if (NphotcodesSkip > 0) {
+      found = FALSE;
+      for (k = 0; (k < NphotcodesSkip) && !found; k++) {
+	if (photcodesSkip[k][0].code == timage[i].photcode) found = TRUE;
+	if (photcodesSkip[k][0].code == GetPhotcodeEquivCodebyCode(timage[i].photcode)) found = TRUE;
+      }
+      if (found) continue;
+    }
+
+    /* exclude images by time */
+    if (TimeSelect) {
+      if (timage[i].tzero < TSTART) continue;
+      if (timage[i].tzero > TSTOP) continue;
+    }
+    
+    if (!FindMosaicForImage (timage, Ntimage, i)) {
+      fprintf (stderr, "cannot find mosaic for %lld\n", (long long) i);
+      continue;
+    }
+
+    /* define image corners */
+    Xi[0] = 0;            Yi[0] = 0;
+    Xi[1] = timage[i].NX; Yi[1] = 0;
+    Xi[2] = timage[i].NX; Yi[2] = timage[i].NY;
+    Xi[3] = 0;            Yi[3] = timage[i].NY;
+    Xi[4] = 0;            Yi[4] = 0;
+    found = FALSE;
+
+    /* transform corners to ra,dec */
+    double RminImage = 360.0;
+    for (j = 0; j < 5; j++) {
+      XY_to_RD (&Ri[j], &Di[j], Xi[j], Yi[j], &timage[i].coords);
+      RminImage = MIN(RminImage, Ri[j]);
+    }
+
+    // RA(nStart) is guaranteed to be < RminImage:
+    nStart = getRegionStartByRA (RminImage, RmaxSky, skylist[0].Nregions);
+
+    /* compare with each region file */
+    for (iSky = nStart; (iSky < skylist[0].Nregions) && !found; iSky++) { 
+
+      m = index[iSky];
+
+      /* we make positional comparisons in the projection of catalog */
+      tcoords.crval1 = skycoords[m].Rc;
+      tcoords.crval2 = skycoords[m].Dc;
+
+      /* transform corner coords to X,Y in this catalog system */
+      InRange = TRUE;
+      for (j = 0; (j < 5) && InRange; j++) {
+	InRange = RD_to_XY (&Xi[j], &Yi[j], Ri[j], Di[j], &tcoords);
+      }
+      if (!InRange) continue;
+
+      /* check if image corner inside catalog */
+      for (j = 0; (j < 4) && !found; j++) {
+	found = corner_check (&Xi[j], &Yi[j], &skycoords[m].Xc[0], &skycoords[m].Yc[0]);
+      }
+      /* check if catalog corner inside image */
+      for (j = 0; (j < 4) && !found; j++) {
+	found = corner_check (&skycoords[m].Xc[j], &skycoords[m].Yc[j], &Xi[0], &Yi[0]);
+      }
+      /* check if edges cross */
+      for (j = 0; (j < 4) && !found; j++) {
+	for (k = 0; (k < 4) && !found; k++) {
+	  found = edge_check (&Xi[j], &Yi[j], &skycoords[m].Xc[k], &skycoords[m].Yc[k]);
+	}
+      }
+      if (!found) continue;
+
+      image[nimage] = timage[i]; 
+      /* always allow 'few' images to succeed, if possible */
+      if (image[nimage].flags & ID_IMAGE_FEW) { 
+	image[nimage].flags &= ~(ID_IMAGE_FEW | ID_IMAGE_POOR);
+      }
+      if (RESET) {
+	// XXX do we need / want to do this in relastro?
+	assignMcal (&image[nimage], (double *) NULL, -1);
+	image[nimage].dMcal = NAN;
+	image[nimage].flags &= ~ID_IMAGE_POOR;
+      }
+      line_number[nimage] = i;
+      nimage ++;
+      if (nimage == NIMAGE) {
+	NIMAGE += 100;
+	REALLOCATE (image, Image, NIMAGE);
+	REALLOCATE (line_number, off_t, NIMAGE);
+      }
+    }
+  }
+      
+  if (VERBOSE) fprintf (stderr, "found %lld images\n", (long long) nimage);
+
+  REALLOCATE (image, Image, MAX (nimage, 1));
+  REALLOCATE (line_number, off_t, MAX (nimage, 1));
+  free (skycoords);
+
+  *Nimage  = nimage;
+  *LineNumber = line_number;
+  return (image);
+}
+
+/* check if line between points 0 and 1 of x1
+   crosses line between points 0 and 1 of x2 */
+int edge_check (double *x1, double *y1, double *x2, double *y2) {
+
+  double theta1, theta2;
+  double Theta1, Theta2;
+
+  theta1 = opening_angle (x1[0], y1[0], x2[0], y2[0], x1[1], y1[1]); 
+  theta2 = opening_angle (x1[0], y1[0], x2[0], y2[0], x2[1], y2[1]); 
+
+  if (theta1*theta2 < 0.0) {
+    return (FALSE);
+  }
+
+  if (fabs(theta1) < fabs(theta2)) {
+    return (FALSE);
+  }
+
+  Theta1 = theta1;
+  Theta2 = theta2;
+  theta1 = opening_angle (x2[0], y2[0], x1[1], y1[1], x2[1], y2[1]); 
+  theta2 = opening_angle (x2[0], y2[0], x1[1], y1[1], x1[0], y1[0]); 
+  
+ 
+  if (theta1*theta2 < 0.0) {
+    return (FALSE);
+  }
+
+  if (fabs(theta1) < fabs(theta2)) {
+    return (FALSE);
+  }
+
+  return (TRUE);
+
+}
+
+/* check if point x1,y1 is in box formed by x2[0-4] */
+int corner_check (double *x1, double *y1, double *x2, double *y2) {
+
+  int i;
+  double theta;
+
+  theta = 0;
+
+  for (i = 0; i < 4; i++) {
+    theta += opening_angle (x2[i], y2[i], x1[0], y1[0], x2[i+1], y2[i+1]); 
+  }
+  if (fabs(theta) > 6) {
+    return (TRUE);
+  } else {
+    return (FALSE);
+  }
+}
+
+/* returns the opening angle between the three points (2 is in middle) 
+   in range -pi to pi */
+
+double opening_angle (double x1, double y1, double x2, double y2, double x3, double y3) {
+
+  double dx1, dy1, dx2, dy2, ct, st, theta;
+
+  dx1 = x1 - x2;
+  dy1 = y1 - y2;
+  
+  dx2 = x3 - x2;
+  dy2 = y3 - y2;
+  
+  ct = (dx1*dx2 + dy1*dy2);
+  st = (dx1*dy2 - dx2*dy1);
+
+  theta = atan2 (st, ct);
+
+  return (theta);
+
+}
+
+void dsortindex (double *X, off_t *Y, int N) {
+
+# define SWAPFUNC(A,B){ double tmpf; off_t tmpi; \
+  tmpf = X[A]; X[A] = X[B]; X[B] = tmpf; \
+  tmpi = Y[A]; Y[A] = Y[B]; Y[B] = tmpi; \
+}
+# define COMPARE(A,B)(X[A] < X[B])
+
+  OHANA_SORT (N, COMPARE, SWAPFUNC);
+
+# undef SWAPFUNC
+# undef COMPARE
+
+}
+
+off_t getRegionStartByRA (double R, double *Rref, off_t Nregions) {
+
+  // use bisection to find the overlapping mosaic
+
+  off_t Nlo, Nhi, N;
+
+  // find the last mosaic before start
+  Nlo = 0; Nhi = Nregions;
+  while (Nhi - Nlo > 10) {
+    N = 0.5*(Nlo + Nhi);
+    if (Rref[N] < R) {
+      Nlo = MAX(N, 0);
+    } else {
+      Nhi = MIN(N, Nregions);
+    }
+  }
+  return (Nlo);
+}
Index: /branches/eam_branches/relastro.20100326/src/write_coords.c
===================================================================
--- /branches/eam_branches/relastro.20100326/src/write_coords.c	(revision 27487)
+++ /branches/eam_branches/relastro.20100326/src/write_coords.c	(revision 27487)
@@ -0,0 +1,29 @@
+# include "relastro.h"
+
+# define CD_COORDS 1
+
+void write_coords (Header *header, Coords *coords) {
+
+  gfits_modify (header, "CTYPE1",   "%s",  1, "RA---TAN");
+  gfits_modify (header, "CTYPE2",   "%s",  1, "DEC--TAN");
+
+  gfits_modify (header, "CRVAL1",   "%lf", 1, coords[0].crval1);
+  gfits_modify (header, "CRVAL2",   "%lf", 1, coords[0].crval2);  
+
+  gfits_modify (header, "CRPIX1",   "%lf", 1, coords[0].crpix1);
+  gfits_modify (header, "CRPIX2",   "%lf", 1, coords[0].crpix2);
+
+# if (CD_COORDS)  
+  gfits_modify (header, "CD1_1",    "%le", 1, coords[0].pc1_1 * coords[0].cdelt1);
+  gfits_modify (header, "CD2_1",    "%le", 1, coords[0].pc2_1 * coords[0].cdelt1);
+  gfits_modify (header, "CD1_2",    "%le", 1, coords[0].pc1_2 * coords[0].cdelt2);
+  gfits_modify (header, "CD2_2",    "%le", 1, coords[0].pc2_2 * coords[0].cdelt2);
+# else
+  gfits_modify (header, "CDELT1",   "%le", 1, coords[0].cdelt1); 
+  gfits_modify (header, "CDELT2",   "%le", 1, coords[0].cdelt2);
+  gfits_modify (header, "PC001001", "%le", 1, coords[0].pc1_1);
+  gfits_modify (header, "PC001002", "%le", 1, coords[0].pc1_2);
+  gfits_modify (header, "PC002001", "%le", 1, coords[0].pc2_1);
+  gfits_modify (header, "PC002002", "%le", 1, coords[0].pc2_2);
+# endif
+}