Index: /branches/czw_branch/20160809/psModules/src/astrom/pmAstrometryObjects.c
===================================================================
--- /branches/czw_branch/20160809/psModules/src/astrom/pmAstrometryObjects.c	(revision 39684)
+++ /branches/czw_branch/20160809/psModules/src/astrom/pmAstrometryObjects.c	(revision 39685)
@@ -225,5 +225,6 @@
     psArray *ref,
     psArray *match,
-    psStats *stats)
+    psStats *stats,
+    const psMetadata *config)
 {
     PS_ASSERT_PTR_NON_NULL(map, NULL);
@@ -232,4 +233,11 @@
     PS_ASSERT_PTR_NON_NULL(match, NULL);
     PS_ASSERT_PTR_NON_NULL(stats, NULL);
+    PS_ASSERT_PTR_NON_NULL(config, NULL);
+
+    // sigma of gaussian window to down-weight photometric outliers (ignored if NAN or 0.0)
+    bool status;
+    double photomWindowSigma  = psMetadataLookupF32 (&status, config, "PSASTRO.PHOTOM.WINDOW.SIGMA");
+    bool   photomWindowApply  = !(isnan(photomWindowSigma) || (fabs(photomWindowSigma < 0.01)));
+    double photomWindowFactor = photomWindowApply ? -0.5/PS_SQR(photomWindowSigma) : 0.0;
 
     // reassign values for clip fit
@@ -252,5 +260,6 @@
         y->data.F32[i] = refStar->FP->y;
 
-        wt->data.F32[i] = 1.01 - exp(-10.0*PS_SQR(refStar->magCal - rawStar->magCal)); // sigma = 0.22 mag
+	// wt is used as an error (sqrt(variance)) in the fit. the 1.01 prevents the weight from going to 0.0 for perfect matches
+        wt->data.F32[i] = 1.01 - exp(photomWindowFactor*PS_SQR(refStar->magCal - rawStar->magCal)); 
     }
 
@@ -776,6 +785,11 @@
     double gridOffset = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.OFFSET");
 
-    // sampling scale of the grid
-    double gridScale  = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.SCALE");
+     // sampling scale of the grid
+     double gridScale  = psMetadataLookupF32 (&status, config, "PSASTRO.GRID.SCALE");
+
+    // sigma of gaussian window to down-weight photometric outliers (ignored if NAN or 0.0)
+    double photomWindowSigma  = psMetadataLookupF32 (&status, config, "PSASTRO.PHOTOM.WINDOW.SIGMA");
+    bool   photomWindowApply  = !(isnan(photomWindowSigma) || (fabs(photomWindowSigma < 0.01)));
+    double photomWindowFactor = photomWindowApply ? -0.5/PS_SQR(photomWindowSigma) : 0.0;
 
     // set the static scaling factors
@@ -818,5 +832,6 @@
             }
 
-	    int Npts = 10 * exp(-10.0*PS_SQR(ob1->magCal - ob2->magCal));  // sigma = 0.22 mag
+	    // XXX should I make the scale factor in front a recipe value?
+	    int Npts = 10 * exp(photomWindowFactor*PS_SQR(ob1->magCal - ob2->magCal)); 
 
             // accumulate bin stats
@@ -1057,4 +1072,9 @@
     double tweakNsigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.TWEAK.NSIGMA");
 
+    // sigma of gaussian window to down-weight photometric outliers (ignored if NAN or 0.0)
+    double photomWindowSigma  = psMetadataLookupF32 (&status, recipe, "PSASTRO.PHOTOM.WINDOW.SIGMA");
+    bool   photomWindowApply  = !(isnan(photomWindowSigma) || (fabs(photomWindowSigma < 0.01)));
+    double photomWindowFactor = photomWindowApply ? -0.5/PS_SQR(photomWindowSigma) : 0.0;
+
     nBin = 2*tweakRange / tweakScale;
     psVector *xHist = psVectorAlloc (nBin, PS_TYPE_F32);
@@ -1083,5 +1103,6 @@
                 continue;
 
-	    int Npts = 10 * exp(-10.0*PS_SQR(ob1->magCal - ob2->magCal));  // sigma = 0.22 mag
+	    // XXX should I make the scale factor in front a recipe value?
+	    int Npts = 10 * exp(photomWindowFactor*PS_SQR(ob1->magCal - ob2->magCal));  // sigma = 0.22 mag
 
             xHist->data.F32[xBin] += Npts;
Index: /branches/czw_branch/20160809/psModules/src/astrom/pmAstrometryObjects.h
===================================================================
--- /branches/czw_branch/20160809/psModules/src/astrom/pmAstrometryObjects.h	(revision 39684)
+++ /branches/czw_branch/20160809/psModules/src/astrom/pmAstrometryObjects.h	(revision 39685)
@@ -348,5 +348,6 @@
     psArray *ref,
     psArray *match,
-    psStats *stats
+    psStats *stats,
+    const psMetadata *config
 );
 
Index: /branches/czw_branch/20160809/psModules/src/astrom/pmAstrometryWCS.c
===================================================================
--- /branches/czw_branch/20160809/psModules/src/astrom/pmAstrometryWCS.c	(revision 39684)
+++ /branches/czw_branch/20160809/psModules/src/astrom/pmAstrometryWCS.c	(revision 39685)
@@ -277,4 +277,29 @@
     wcs->crpix2 = psMetadataLookupF64 (&status, header, "CRPIX2");
     wcs->toSky = psProjectionAlloc (wcs->crval1*PM_RAD_DEG, wcs->crval2*PM_RAD_DEG, PM_RAD_DEG, PM_RAD_DEG, type);
+
+    // XXX if type == ZPN, look for PV2_%d elements:
+    if (type == PS_PROJ_ZPN) {
+      psVector *maxRadial = psVectorAlloc (21, PS_TYPE_F64);
+      for (int i = 0; i <= 20; i++) {
+	char name[64];
+	snprintf (name, 64, "PV2_%d", i);
+
+	maxRadial->data.F64[i] = 0.0;
+	double value = psMetadataLookupF64 (&status, header, name);
+
+	if (status) {
+	  maxRadial->data.F64[i] = value;
+	  maxRadial->n = i;
+	}
+
+	// PV2_1 is implicit if not present
+	if ((i == 1) && !status) {
+	  maxRadial->data.F64[i] = 1.0;
+	  continue;
+	}
+      }
+      maxRadial->n ++;
+      wcs->toSky->radial = maxRadial;
+    }
 
     // These aren't needed but having them empty is disconcerting
@@ -375,4 +400,13 @@
     psMetadataAddF64 (header, PS_LIST_TAIL, "CRPIX1", PS_META_REPLACE, "", wcs->crpix1);
     psMetadataAddF64 (header, PS_LIST_TAIL, "CRPIX2", PS_META_REPLACE, "", wcs->crpix2);
+
+    if (wcs->toSky->type == PS_PROJ_ZPN) {
+      psAssert (wcs->toSky->radial, "missing radial vector");
+      for (int i = 0; i < wcs->toSky->radial->n; i++) {
+	if (wcs->toSky->radial->data.F64[i] == 0.0) continue;
+	snprintf (name, 16, "PV2_%d", i);
+	psMetadataAddF64 (header, PS_LIST_TAIL, name, PS_META_REPLACE, "", wcs->toSky->radial->data.F64[i]);
+      }
+    }
 
     // XXX make it optional to write out CDi_j terms, or other versions
@@ -485,4 +519,5 @@
     // projection from TPA (linear microns) to SKY (radians)
     psProjection *toSky = psProjectionAlloc (wcs->toSky->R, wcs->toSky->D, PM_RAD_DEG*pdelt1, PM_RAD_DEG*pdelt2, wcs->toSky->type);
+    toSky->radial = psMemIncrRefCounter (wcs->toSky->radial);
 
     if (fpa->toSky == NULL) {
@@ -633,4 +668,5 @@
     // cdelt1,2 has units of degrees/micron
     fpa->toSky = psProjectionAlloc (wcs->toSky->R, wcs->toSky->D, wcs->cdelt1*PM_RAD_DEG, wcs->cdelt2*PM_RAD_DEG, wcs->toSky->type);
+    fpa->toSky->radial = psMemIncrRefCounter (wcs->toSky->radial);
 
     // create transformation with 0,0 reference pixel
@@ -701,4 +737,6 @@
     // convert projection from FPA to SKY into wcs projection (degrees to radians)
     wcs->toSky = psProjectionAlloc (fpa->toSky->R, fpa->toSky->D, PM_RAD_DEG, PM_RAD_DEG, fpa->toSky->type);
+    wcs->toSky->radial = psMemIncrRefCounter (fpa->toSky->radial);
+
     wcs->crval1 = fpa->toSky->R*PS_DEG_RAD;
     wcs->crval2 = fpa->toSky->D*PS_DEG_RAD;