Index: /branches/cnb_branch_20081011/psastro/src/psastroVisual.c
===================================================================
--- /branches/cnb_branch_20081011/psastro/src/psastroVisual.c	(revision 20264)
+++ /branches/cnb_branch_20081011/psastro/src/psastroVisual.c	(revision 20264)
@@ -0,0 +1,1018 @@
+/***Diagnostic plots for psastro***/
+
+# include "psastroInternal.h"
+
+# if (HAVE_KAPA)
+# include <kapa.h>
+
+# define KAPAX  800
+# define KAPAY  800
+
+//variables to determine when things are plotted
+static bool isVisual             = false;
+static bool plotRawStars         = false;
+static bool plotRefStars         = false;
+static bool plotLumFunc          = false;
+static bool plotRemoveClumps     = true;
+static bool plotOneChipFit       = false;
+static bool plotAstromGuessCheck = false;
+
+// variables to store plotting window indices
+static int kapa = -1; //for plots
+static int kapa2 = -1; //also for plots
+
+// helper routine prototypes
+bool psastroVisualTriplePlot (int kapa, Graphdata *graphdata, psVector *xVec, psVector *yVec, psVector *zVec, bool increasing);
+bool psastroVisualScaleGraphdata(Graphdata *graphdata, psVector *xVec, psVector *yVec);
+bool psastroVisualTripleOverplot (int kapa, Graphdata *graphdata, psVector *xVec, psVector *yVec, psVector *zVec, bool increasing);
+bool psastroVisualCreateScaleVec (psVector *zVec, psVector *zScale, bool increasing);
+
+
+/** Initialization Routines **/
+
+
+/**
+ *  start or stop plotting
+ */
+bool psastroSetVisual (bool mode) {
+    isVisual = mode;
+    return true;
+}
+
+
+/**
+ * open, name, and resize a window if necessary
+ */
+bool psastroVisualInitWindow( int *kapid, char *name ) {
+    if (*kapid == -1) {
+        *kapid = KapaOpenNamedSocket("kapa", name);
+        if (*kapid == -1) {
+            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
+            isVisual = false;
+            return false;
+        }
+        KapaResize (*kapid, KAPAX, KAPAY);
+    }
+    return true;
+}
+
+
+/**
+ * ask the user how to proceed
+ */
+bool psastroVisualAskUser( bool *plotflag ) {
+    char key[10];
+    fprintf (stdout, "[c]ontinue? [s]kip the rest of these plots? [a]bort all visual plots?");
+    if (!fgets(key, 8, stdin)) {
+        psWarning("Unable to read option");
+    }
+    if (key[0] == 's') {
+        *plotflag = false;
+    }
+    if (key[0] == 'a') {
+        psastroSetVisual(false);
+    }
+    return true;
+}
+
+
+/**
+ * destroy windows at the end of a run
+ */
+bool psastroVisualClose() {
+    if(kapa != -1)
+        KiiClose(kapa);
+    if (kapa2 != -1)
+        KiiClose(kapa2);
+    return true;
+}
+
+/** Plotting Routines **/
+
+/*  psastroVisualPlotRawStars */
+//  Show raw star astrometry as determined from first pass fit
+//  Invoked from within psastroAstromGuess
+bool psastroVisualPlotRawStars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe)
+{
+
+    if (!plotRawStars || !isVisual) return true;
+    if (!psastroVisualInitWindow (&kapa, "psastro:plots"))
+        return false;
+
+    Graphdata graphdata;
+    KapaSection section;
+
+    bool status = false;
+    float iMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MIN");
+    float iMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MAX");
+
+    KapaInitGraph (&graphdata);
+    KapaClearPlots (kapa);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.ptype = 7;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+
+    section.dx = 0.4;
+    section.dy = 0.4;
+
+    psVector *xVec = psVectorAlloc (rawstars->n, PS_TYPE_F32);
+    psVector *yVec = psVectorAlloc (rawstars->n, PS_TYPE_F32);
+    psVector *zVec = psVectorAlloc (rawstars->n, PS_TYPE_F32);
+
+    section.x = 0.0;
+    section.y = 0.5;
+    section.name = NULL;
+    psStringAppend (&section.name, "a0");
+    KapaSetSection (kapa, &section);
+    psFree (section.name);
+    int n = 0;
+    for (int i = 0; i < rawstars->n; i++) {
+        pmAstromObj *raw = rawstars->data[i];
+        if (!isfinite(raw->Mag)) continue;
+        if (raw->Mag < iMagMin) continue;
+        if (raw->Mag > iMagMax) continue;
+
+        xVec->data.F32[n] = raw->chip->x;
+        yVec->data.F32[n] = raw->chip->y;
+        zVec->data.F32[n] = raw->Mag;
+        n++;
+    }
+    xVec->n = yVec->n = zVec->n = n;
+
+    KapaSendLabel(kapa, "Chip", KAPA_LABEL_XP);
+    KapaSendLabel(kapa, "X", KAPA_LABEL_XM);
+    KapaSendLabel(kapa, "Y", KAPA_LABEL_YM);
+    psastroVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false);
+
+    section.x = 0.5;
+    section.y = 0.5;
+    section.name = NULL;
+    psStringAppend (&section.name, "a1");
+    KapaSetSection (kapa, &section);
+    psFree (section.name);
+
+    n = 0;
+    for (int i = 0; i < rawstars->n; i++) {
+        pmAstromObj *raw = rawstars->data[i];
+        if (!isfinite(raw->Mag)) continue;
+        if (raw->Mag < iMagMin) continue;
+        if (raw->Mag > iMagMax) continue;
+
+        xVec->data.F32[n] = raw->FP->x;
+        yVec->data.F32[n] = raw->FP->y;
+        zVec->data.F32[n] = raw->Mag;
+        n++;
+    }
+    xVec->n = yVec->n = zVec->n = n;
+
+    KapaSendLabel (kapa, "Focal Plane", KAPA_LABEL_XP);
+    KapaSendLabel (kapa, "L", KAPA_LABEL_XM);
+    KapaSendLabel (kapa, "M", KAPA_LABEL_YM);
+    psastroVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false);
+
+    section.x = 0.0;
+    section.y = 0.0;
+    section.name = NULL;
+    psStringAppend (&section.name, "a2");
+    KapaSetSection (kapa, &section);
+    psFree (section.name);
+
+    n = 0;
+    for (int i = 0; i < rawstars->n; i++) {
+        pmAstromObj *raw = rawstars->data[i];
+        if (!isfinite(raw->Mag)) continue;
+        if (raw->Mag < iMagMin) continue;
+        if (raw->Mag > iMagMax) continue;
+
+        xVec->data.F32[n] = raw->TP->x;
+        yVec->data.F32[n] = raw->TP->y;
+        zVec->data.F32[n] = raw->Mag;
+        n++;
+    }
+    xVec->n = yVec->n = zVec->n = n;
+
+    KapaSendLabel (kapa, "Tangential Plane", KAPA_LABEL_XP);
+    KapaSendLabel (kapa, "P", KAPA_LABEL_XM);
+    KapaSendLabel (kapa, "Q", KAPA_LABEL_YM);
+    psastroVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false);
+
+    section.x = 0.5;
+    section.y = 0.0;
+    section.name = NULL;
+    psStringAppend (&section.name, "a3");
+    KapaSetSection (kapa, &section);
+    psFree (section.name);
+
+    n = 0;
+    for (int i = 0; i < rawstars->n; i++) {
+        pmAstromObj *raw = rawstars->data[i];
+        if (!isfinite(raw->Mag)) continue;
+        if (raw->Mag < iMagMin) continue;
+        if (raw->Mag > iMagMax) continue;
+
+        xVec->data.F32[n] = DEG_RAD*raw->sky->r;
+        yVec->data.F32[n] = DEG_RAD*raw->sky->d;
+        zVec->data.F32[n] = raw->Mag;
+        n++;
+    }
+    xVec->n = yVec->n = zVec->n = n;
+
+    KapaSendLabel (kapa, "Sky", KAPA_LABEL_XP);
+    KapaSendLabel(kapa, "RA", KAPA_LABEL_XM);
+    KapaSendLabel(kapa, "Dec", KAPA_LABEL_YM);
+    psastroVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false);
+
+    // flip x (East increase to left)
+    SWAP (graphdata.xmin, graphdata.xmax);
+    KapaSetLimits (kapa, &graphdata);
+
+    // plot label
+    section.x = 0.0;
+    section.y = 0.0;
+    section.dx = .95;
+    section.dy = .95;
+    section.name = NULL;
+    psStringAppend (&section.name, "a5");
+    KapaSetSection (kapa, &section);
+    KapaSendLabel (kapa, "Raw Star Selection - Initial Astrometry", KAPA_LABEL_XP);
+    psFree (section.name);
+
+    // pause and wait for user input:
+    psastroVisualAskUser( &plotRawStars );
+
+    psFree (xVec);
+    psFree (yVec);
+    psFree (zVec);
+    return true;
+}
+
+
+/* psastroVisualPlotRefStars */
+// plot the location of references stars over the entire fpa
+// invoked during psastroChooseRefStars
+bool psastroVisualPlotRefStars (psArray *refstars, psMetadata *recipe)
+{
+    if (!isVisual || !plotRefStars) return true;
+    Graphdata graphdata;
+    if (!psastroVisualInitWindow (&kapa, "psastro:plots"))
+        return false;
+
+    bool status = false;
+    float rMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MIN");
+    float rMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MAX");
+
+    KapaInitGraph (&graphdata);
+    KapaClearSections (kapa);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.ptype = 7;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+
+    psVector *xVec = psVectorAlloc (refstars->n, PS_TYPE_F32);
+    psVector *yVec = psVectorAlloc (refstars->n, PS_TYPE_F32);
+    psVector *zVec = psVectorAlloc (refstars->n, PS_TYPE_F32);
+
+    int n = 0;
+    for (int i = 0; i < refstars->n; i++) {
+        pmAstromObj *ref = refstars->data[i];
+        if (!isfinite(ref->Mag)) continue;
+        if (ref->Mag > rMagMax) continue;
+        if (ref->Mag < rMagMin) continue;
+
+        xVec->data.F32[n] = DEG_RAD*ref->sky->r;
+        yVec->data.F32[n] = DEG_RAD*ref->sky->d;
+        zVec->data.F32[n] = ref->Mag;
+        n++;
+    }
+    xVec->n = yVec->n = zVec->n = n;
+
+    KapaSendLabel (kapa, "RA", KAPA_LABEL_XM);
+    KapaSendLabel (kapa, "Dec", KAPA_LABEL_YM);
+    KapaSendLabel (kapa, "Reference Stars", KAPA_LABEL_XP);
+    psastroVisualTriplePlot(kapa, &graphdata, xVec, yVec, zVec, false);
+
+    // flip x (East increase to left)
+    SWAP (graphdata.xmin, graphdata.xmax);
+    KapaSetLimits (kapa, &graphdata);
+
+    // pause and wait for user input:
+    psastroVisualAskUser(&plotRefStars);
+
+    psFree (xVec);
+    psFree (yVec);
+    psFree (zVec);
+    return true;
+}
+
+
+/* psastroVisualPlotLuminosityFunction */
+// Plot the two luminosity functions created.
+// Called from psastroRefStarSubset.
+// The luminosity functions are used to select a subset of reference stars,
+// so we plot the cutoff that defines this subset
+bool psastroVisualPlotLuminosityFunction (psVector *lnMag, psVector *Mag,
+                                          pmLumFunc *lumFunc, pmLumFunc *rawFunc) {
+
+    if ( !isVisual || !plotLumFunc ) return true;
+    if ( !psastroVisualInitWindow (&kapa, "psastro:plots"))
+        return false;
+
+    //Set up graph window
+    Graphdata graphdata;
+    KapaSection section1 = {"s1", .1, .1, .85, .35};
+    KapaSection section2 =  {"s2", .1, .5, .85, .35};
+    KapaSection section;
+    if(rawFunc == NULL)
+        section = section1;
+    else
+        section = section2;
+    if (rawFunc == NULL)
+        KapaClearPlots(kapa);
+    KapaInitGraph(&graphdata);
+
+    //Determine Plot Limits
+    psastroVisualScaleGraphdata(&graphdata, Mag, lnMag);
+
+    //Make a line for the fit
+    float x[2] = {graphdata.xmin, graphdata.xmax};
+    float y[2] = {lumFunc->offset + x[0] * lumFunc->slope,
+                 lumFunc->offset + x[1] * lumFunc->slope};
+
+    //Plot Data
+    KapaSetSection(kapa, &section);
+    KapaSetLimits(kapa, &graphdata);
+
+    KapaSetFont (kapa, "helvetica", 14);
+    KapaBox (kapa, &graphdata);
+    KapaSendLabel (kapa, "Magnitude", KAPA_LABEL_XM);
+    KapaSendLabel (kapa, "Log(N)", KAPA_LABEL_YM);
+    if (rawFunc == NULL)
+        KapaSendLabel (kapa, "Raw Star Luminosity Function", KAPA_LABEL_XP);
+    else
+        KapaSendLabel (kapa,
+                       "Reference Star Luminosity Function, Shifted Raw Fit, and Cutoff",
+                       KAPA_LABEL_XP);
+    graphdata.color=KapaColorByName("black");
+    graphdata.style = 1;
+    KapaPrepPlot (kapa, lnMag->n, &graphdata);
+    KapaPlotVector(kapa, lnMag->n,   Mag->data.F32, "x");
+    KapaPlotVector(kapa, lnMag->n, lnMag->data.F32, "y");
+
+    //Overplot fit
+    graphdata.style=0;
+    KapaPrepPlot(kapa,2,&graphdata);
+    KapaPlotVector(kapa, 2, x, "x");
+    KapaPlotVector(kapa, 2, y, "y");
+
+    //If rawFunc was supplied, overplot the raw star fit + cutoff
+    if( rawFunc != NULL) {
+        double mRef = 0.5*(lumFunc->mMin + lumFunc->mMax);
+        double logRho = mRef * lumFunc->slope + lumFunc->offset;
+        double mRaw = (logRho - rawFunc->offset) / rawFunc->slope;
+        double deltaM = mRef - mRaw;
+        double mRefMax = rawFunc->mMax + deltaM;
+
+        float xraw[2] = {rawFunc->mMin + deltaM, rawFunc->mMax + deltaM};
+        float yraw[2] = {rawFunc->offset + (rawFunc->slope) * rawFunc->mMin,
+                        rawFunc->offset + (rawFunc->slope) * rawFunc->mMax};
+        float x[2] = {mRefMax, mRefMax};
+        float y[2] = {graphdata.ymin, graphdata.ymax};
+        graphdata.color= KapaColorByName("red");
+        KapaPrepPlot(kapa, 2, &graphdata);
+        KapaPlotVector(kapa, 2, x, "x");
+        KapaPlotVector(kapa, 2, y, "y");
+        KapaPrepPlot (kapa, 2, &graphdata);
+        KapaPlotVector (kapa, 2, xraw, "x");
+        KapaPlotVector (kapa, 2, yraw, "y");
+
+        // pause and wait for user input:
+        psastroVisualAskUser (&plotLumFunc);
+    }
+    return true;
+} // end of psastroVisualPlotLuminosityFunction
+
+
+/**
+ * Plot the stars in a region, and indicate which stars are part of 'clumps'
+ * These stars are flagged during astrometric fitting, since dense regions are
+ * harder to cross-match than sparse ones. Called during psastroRemoveClumps.
+ */
+
+bool psastroVisualPlotRemoveClumps (psArray *input, psImage *count, int scale, float limit) {
+    if (!isVisual || !plotRemoveClumps) return true;
+    if ( !psastroVisualInitWindow (&kapa, "psastro:plots")) return false;
+
+    KapaSection section;
+    Graphdata graphdata;
+
+    //initialize graph information
+    section.x = 0;
+    section.dx = .9;
+    section.y = 0;
+    section.dy = .9;
+    section.name = NULL;
+    psStringAppend( &section.name, "a0");
+    KapaInitGraph(&graphdata);
+    KapaSetSection(kapa, &section);
+    psFree(section.name);
+    graphdata.ptype = 7;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+    graphdata.color = KapaColorByName ("black");
+
+
+    //set up plot vectors
+    float Xmin = +FLT_MAX;
+    float Xmax = -FLT_MAX;
+    float Ymin = +FLT_MAX;
+    float Ymax = -FLT_MAX;
+    psVector *xVec = psVectorAlloc (input->n, PS_TYPE_F32);
+    psVector *yVec = psVectorAlloc (input->n, PS_TYPE_F32);
+
+    //determine boundaries for histogram bin calculation
+    int n = 0;
+    for (int i=0; i< input->n; i++) {
+        pmAstromObj *obj = (pmAstromObj *)input->data[i];
+        if (!isfinite(obj->FP->x)) continue;
+        if (!isfinite(obj->FP->y)) continue;
+        xVec->data.F32[n] = obj->FP->x;
+        yVec->data.F32[n] = obj->FP->y;
+        Xmin = PS_MIN (Xmin, xVec->data.F32[n]);
+        Xmax = PS_MAX (Xmax, xVec->data.F32[n]);
+        Ymin = PS_MIN (Ymin, yVec->data.F32[n]);
+        Ymax = PS_MAX (Ymax, yVec->data.F32[n]);
+        n++;
+    }
+    xVec->n = yVec->n = n;
+
+    //plot stars
+    graphdata.xmax = Xmax;
+    graphdata.xmin = Xmin;
+    graphdata.ymax = Ymax;
+    graphdata.ymin = Ymin;
+    KapaSetLimits (kapa, &graphdata);
+    KapaSetFont (kapa, "helvetica", 14);
+
+    KapaBox (kapa, &graphdata);
+
+    KapaSendLabel (kapa, "L (pixels)", KAPA_LABEL_XM);
+    KapaSendLabel (kapa, "M (pixels)", KAPA_LABEL_YM);
+    KapaSendLabel (kapa, "Regions Flagged as Clumps (Red Boxes)",
+                   KAPA_LABEL_XP);
+
+    KapaPrepPlot (kapa, xVec->n, &graphdata);
+    KapaPlotVector (kapa, xVec->n, xVec->data.F32, "x");
+    KapaPlotVector (kapa, yVec->n, yVec->data.F32, "y");
+
+    graphdata.color = KapaColorByName ("red");
+    graphdata.style = 1;
+
+    //overplot clumpy regions excluded from analysis
+    for(int i = 0; i < count->numCols; i++) {
+        for (int j = 0; j < count->numRows; j++) {
+            if(count->data.U32[j][i] <= limit) continue; //not a clump
+            float Xbot = (i - 5) * scale + Xmin;
+            float Ybot = (j - 5) * scale + Ymin;
+            if(Xbot < graphdata.xmin || Xbot > graphdata.xmax ||
+               Ybot < graphdata.ymin || Ybot > graphdata.ymax) continue;
+            float x[5] = {Xbot, Xbot + scale, Xbot + scale, Xbot, Xbot};
+            float y[5] = {Ybot, Ybot, Ybot + scale, Ybot + scale, Ybot};
+            KapaPrepPlot (kapa, 5, &graphdata);
+            KapaPlotVector (kapa, 5, x, "x");
+            KapaPlotVector (kapa, 5, y, "y");
+        }
+    }
+    psastroVisualAskUser (&plotRemoveClumps);
+    psFree (xVec);
+    psFree (yVec);
+
+    return true;
+}
+
+
+/* psastroVisualPlotOneChipFit */
+// assess the goodness of fit for a signle chip by
+// plotting the fit residuals
+// invoked during psastroOneChipFit
+bool psastroVisualPlotOneChipFit (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe) {
+    if (!isVisual || !plotOneChipFit)
+        return true;
+
+  if (!psastroVisualInitWindow (&kapa, "psastro:plots"))
+        return false;
+
+    bool status = false;
+    float iMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MIN");
+    float iMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MAX");
+    float rMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MIN");
+    float rMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MAX");
+
+    //initialize graph information
+    Graphdata graphdata;
+    KapaSection section;
+
+    KapaInitGraph (&graphdata);
+    KapaClearPlots (kapa);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.ptype = 7;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+
+    section.dx = 0.4;
+    section.dy = 0.4;
+
+    psVector *xVec = psVectorAlloc (match->n, PS_TYPE_F32);
+    psVector *yVec = psVectorAlloc (match->n, PS_TYPE_F32);
+    psVector *zVec = psVectorAlloc (match->n, PS_TYPE_F32);
+
+    // X vs dX
+    section.x = 0.0;
+    section.y = 0.5;
+    section.name = NULL;
+    psStringAppend (&section.name, "a0");
+    KapaSetSection (kapa, &section);
+    psFree (section.name);
+
+    int n = 0;
+    for (int i = 0; i < match->n; i++) {
+        pmAstromMatch *pair = match->data[i];
+        pmAstromObj *raw = rawstars->data[pair->raw];
+        pmAstromObj *ref = refstars->data[pair->ref];
+
+        if (!isfinite(raw->Mag)) continue;
+        if (raw->Mag < iMagMin) continue;
+        if (raw->Mag > iMagMax) continue;
+        if (ref->Mag < rMagMin) continue;
+        if (ref->Mag > rMagMax) continue;
+
+        xVec->data.F32[n] = raw->chip->x;
+        yVec->data.F32[n] = raw->chip->x - ref->chip->x;
+        zVec->data.F32[n] = raw->Mag;
+        n++;
+    }
+    xVec->n = yVec->n = zVec->n = n;
+
+    KapaSendLabel (kapa, "X", KAPA_LABEL_XM);
+    KapaSendLabel (kapa, "dX", KAPA_LABEL_YM);
+    psastroVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false);
+
+    // X vs dY
+    section.x = 0.5;
+    section.y = 0.5;
+    section.name = NULL;
+    psStringAppend (&section.name, "a1");
+    KapaSetSection (kapa, &section);
+    psFree (section.name);
+
+    n = 0;
+    for (int i = 0; i < match->n; i++) {
+        pmAstromMatch *pair = match->data[i];
+        pmAstromObj *raw = rawstars->data[pair->raw];
+        pmAstromObj *ref = refstars->data[pair->ref];
+
+        if (!isfinite(raw->Mag)) continue;
+        if (raw->Mag < iMagMin) continue;
+        if (raw->Mag > iMagMax) continue;
+        if (ref->Mag < rMagMin) continue;
+        if (ref->Mag > rMagMax) continue;
+
+        xVec->data.F32[n] = raw->chip->x;
+        yVec->data.F32[n] = raw->chip->y - ref->chip->y;
+        zVec->data.F32[n] = raw->Mag;
+        n++;
+    }
+    xVec->n = yVec->n = zVec->n = n;
+
+    KapaSendLabel (kapa, "X", KAPA_LABEL_XM);
+    KapaSendLabel (kapa, "dY", KAPA_LABEL_YM);
+    psastroVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false);
+
+    // Y vs dX
+    section.x = 0.0;
+    section.y = 0.0;
+    section.name = NULL;
+    psStringAppend (&section.name, "a2");
+    KapaSetSection (kapa, &section);
+    psFree (section.name);
+
+    n = 0;
+    for (int i = 0; i < match->n; i++) {
+        pmAstromMatch *pair = match->data[i];
+        pmAstromObj *raw = rawstars->data[pair->raw];
+        pmAstromObj *ref = refstars->data[pair->ref];
+
+        if (!isfinite(raw->Mag)) continue;
+        if (raw->Mag < iMagMin) continue;
+        if (raw->Mag > iMagMax) continue;
+        if (ref->Mag < rMagMin) continue;
+        if (ref->Mag > rMagMax) continue;
+
+        xVec->data.F32[n] = raw->chip->y;
+        yVec->data.F32[n] = raw->chip->x - ref->chip->x;
+        zVec->data.F32[n] = raw->Mag;
+        n++;
+    }
+    xVec->n = yVec->n = zVec->n = n;
+
+    KapaSendLabel (kapa, "Y", KAPA_LABEL_XM);
+    KapaSendLabel (kapa, "dX", KAPA_LABEL_YM);
+    psastroVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false);
+
+    // Y vs dY
+    section.x = 0.5;
+    section.y = 0.0;
+    section.name = NULL;
+    psStringAppend (&section.name, "a3");
+    KapaSetSection (kapa, &section);
+    psFree (section.name);
+
+    n = 0;
+    for (int i = 0; i < match->n; i++) {
+        pmAstromMatch *pair = match->data[i];
+        pmAstromObj *raw = rawstars->data[pair->raw];
+        pmAstromObj *ref = refstars->data[pair->ref];
+
+        if (!isfinite(raw->Mag)) continue;
+        if (raw->Mag < iMagMin) continue;
+        if (raw->Mag > iMagMax) continue;
+        if (ref->Mag < rMagMin) continue;
+        if (ref->Mag > rMagMax) continue;
+
+        xVec->data.F32[n] = raw->chip->y;
+        yVec->data.F32[n] = raw->chip->y - ref->chip->y;
+        zVec->data.F32[n] = raw->Mag;
+        n++;
+    }
+    xVec->n = yVec->n = zVec->n = n;
+
+    KapaSendLabel (kapa, "Y", KAPA_LABEL_XM);
+    KapaSendLabel (kapa, "dY", KAPA_LABEL_YM);
+    pmKapaPlotVectorTriple_AutoLimits_OpenGraph (kapa, &graphdata, xVec, yVec, zVec, false);
+
+    section.x = 0.0;
+    section.y = 0.0;
+    section.dx = 0.95;
+    section.dy = 0.95;
+    section.name = NULL;
+    psStringAppend (&section.name, "a5");
+    KapaSetSection (kapa, &section);
+    KapaSendLabel (kapa, "Single Chip Fit Residuals (Chip Coordinates)", KAPA_LABEL_XP);
+    psFree (section.name);
+
+
+    // *** X vs Y plot (different window)
+    if (!psastroVisualInitWindow( &kapa2, "psastro:plots"))
+        return false;
+
+    KapaInitGraph (&graphdata);
+    KapaClearPlots (kapa2);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.ptype = 2;
+    graphdata.style = 2;
+
+    psFree (xVec);
+    psFree (yVec);
+    psFree (zVec);
+
+    xVec = psVectorAlloc (rawstars->n, PS_TYPE_F32);
+    yVec = psVectorAlloc (rawstars->n, PS_TYPE_F32);
+    zVec = psVectorAlloc (rawstars->n, PS_TYPE_F32);
+
+    // X vs Y by mag (raw)
+    n = 0;
+    for (int i = 0; i < rawstars->n; i++) {
+        pmAstromObj *raw = rawstars->data[i];
+        if (!isfinite(raw->Mag)) continue;
+        if (raw->Mag < iMagMin) continue;
+        if (raw->Mag > iMagMax) continue;
+
+        xVec->data.F32[n] = raw->chip->x;
+        yVec->data.F32[n] = raw->chip->y;
+        zVec->data.F32[n] = raw->Mag;
+        n++;
+    }
+    xVec->n = yVec->n = zVec->n = n;
+
+    KapaSendLabel (kapa2, "X", KAPA_LABEL_XM);
+    KapaSendLabel (kapa2, "Y", KAPA_LABEL_YM);
+    KapaSendLabel (kapa2,
+                   "Chip Coordinates. Black = Raw Stars. Red = Ref Stars. Blue = Matched Stars"
+                   , KAPA_LABEL_XP);
+    psastroVisualTriplePlot (kapa2, &graphdata, xVec, yVec, zVec, false);
+
+    // X vs Y by mag (ref)
+    psFree (xVec);
+    psFree (yVec);
+    psFree (zVec);
+
+    xVec = psVectorAlloc (refstars->n, PS_TYPE_F32);
+    yVec = psVectorAlloc (refstars->n, PS_TYPE_F32);
+    zVec = psVectorAlloc (refstars->n, PS_TYPE_F32);
+
+    graphdata.color = KapaColorByName ("red");
+    graphdata.ptype = 7;
+    graphdata.style = 2;
+
+    n = 0;
+    for (int i = 0; i < refstars->n; i++) {
+        pmAstromObj *ref = refstars->data[i];
+        if (!isfinite(ref->Mag)) continue;
+        if (ref->Mag < rMagMin) continue;
+        if (ref->Mag > rMagMax) continue;
+
+        xVec->data.F32[n] = ref->chip->x;
+        yVec->data.F32[n] = ref->chip->y;
+        zVec->data.F32[n] = ref->Mag;
+        n++;
+    }
+    xVec->n = yVec->n = zVec->n = n;
+    psastroVisualTripleOverplot (kapa2, &graphdata, xVec, yVec, zVec, false);
+
+    //rescale the graph to include all points
+    float xmin = graphdata.xmin;
+    float ymin = graphdata.ymin;
+    float xmax = graphdata.xmax;
+    float ymax = graphdata.ymax;
+    psastroVisualScaleGraphdata(&graphdata, xVec, yVec);
+    graphdata.xmin = PS_MIN(xmin, graphdata.xmin);
+    graphdata.ymin = PS_MIN(ymin, graphdata.ymin);
+    graphdata.xmax = PS_MAX(xmax, graphdata.xmax);
+    graphdata.ymax = PS_MAX(ymax, graphdata.ymax);
+    KapaSetLimits (kapa2, &graphdata);
+
+    //overplot matched stars in blue
+    psFree (xVec);
+    psFree (yVec);
+    psFree (zVec);
+
+    xVec = psVectorAlloc (match->n, PS_TYPE_F32);
+    yVec = psVectorAlloc (match->n, PS_TYPE_F32);
+    zVec = psVectorAlloc (match->n, PS_TYPE_F32);
+
+    graphdata.color = KapaColorByName ("blue");
+    n = 0;
+    for (int i = 0; i < match->n; i++) {
+        pmAstromMatch *pair = match->data[i];
+        pmAstromObj *raw = rawstars->data[pair->raw];
+        pmAstromObj *ref = refstars->data[pair->ref];
+        if (raw->Mag < iMagMin) continue;
+        if (raw->Mag > iMagMax) continue;
+        if (ref->Mag < rMagMin) continue;
+        if (ref->Mag > rMagMax) continue;
+
+        xVec->data.F32[n] = raw->chip->x;
+        yVec->data.F32[n] = raw->chip->y;
+        zVec->data.F32[n] = iMagMin;
+        n++;
+    }
+    xVec->n = yVec->n = zVec->n = n;
+    psastroVisualTripleOverplot (kapa2, &graphdata, xVec, yVec, zVec, false);
+
+    psastroVisualAskUser (&plotOneChipFit);
+
+    psFree (xVec);
+    psFree (yVec);
+    psFree (zVec);
+    return true;
+}
+
+
+/**
+ *  Plots the fpa chip corners projected on to the tangential plane before and after
+ *  the astrometry solution has been applied. In psastroAstromGuessCheck, the old corners
+ *  are then fit to the new corners to get a sense at how far off the initial WCS info was
+ *  in offset, rotation, and scale. This procedure also plots the residuals of the fit from
+ *  old to new coordinates
+ */
+bool psastroVisualPlotAstromGuessCheck (psVector *cornerPo, psVector *cornerQo,
+                                        psVector *cornerPn, psVector *cornerQn,
+                                        psVector *cornerPd, psVector *cornerQd) {
+
+    if (!isVisual || !plotAstromGuessCheck) return true;
+    if ( !psastroVisualInitWindow (&kapa, "psastro:plots"))
+        return false;
+
+    //set up graph window
+    Graphdata graphdata;
+    KapaSection section;
+    KapaInitGraph(&graphdata);
+    KapaClearPlots (kapa);
+
+    graphdata.color = KapaColorByName ("black");
+    graphdata.ptype = 2;
+    graphdata.size = 0.5;
+    graphdata.style = 2;
+
+    section.dx = 0.4;
+    section.dy = 0.4;
+
+    //Old Corners
+    section.x = 0.30;
+    section.y = 0.50;
+    section.name = NULL;
+    psStringAppend (&section.name, "a0");
+    KapaSetSection (kapa, &section);
+    psFree(section.name);
+
+    psastroVisualScaleGraphdata (&graphdata, cornerPo, cornerPo);
+    KapaSetLimits (kapa, &graphdata);
+    KapaBox (kapa, &graphdata);
+    KapaSendLabel (kapa, "P (Pixels)", KAPA_LABEL_XM);
+    KapaSendLabel (kapa, "Q (Pixels)", KAPA_LABEL_YM);
+    KapaSendLabel (kapa,
+                   "Fiducial Points in the Tangent Plane. Black: Initial Astrometry. Red: Final Astrometry",
+                   KAPA_LABEL_XP);
+    KapaPrepPlot (kapa, cornerPo->n, &graphdata);
+    KapaPlotVector (kapa, cornerPo->n, cornerPo->data.F32, "x");
+    KapaPlotVector (kapa, cornerQo->n, cornerQo->data.F32, "y");
+
+    // New Corners
+    graphdata.color = KapaColorByName("red");
+    graphdata.ptype = 7;
+    graphdata.size = 1.5;
+    KapaPrepPlot (kapa, cornerPn->n, &graphdata);
+    KapaPlotVector (kapa, cornerPn->n, cornerPn->data.F32, "x");
+    KapaPlotVector (kapa, cornerQn->n, cornerQn->data.F32, "y");
+
+    // Residuals
+    psVector *xResid = psVectorAlloc(cornerPn->n, PS_DATA_F32);
+    psVector *yResid = psVectorAlloc(cornerQn->n, PS_DATA_F32);
+    for(int i=0; i < cornerPn->n; i++) {
+        xResid->data.F32[i] = (cornerPd->data.F32[i] / cornerPn->data.F32[i]);
+        yResid->data.F32[i] = (cornerQd->data.F32[i] / cornerQn->data.F32[i]);
+    }
+
+    graphdata.color = KapaColorByName("black");
+    graphdata.size=0.5;
+    section.x = 0.3;
+    section.y = 0.0;
+    section.name = NULL;
+    psStringAppend (&section.name, "a1");
+    KapaSetSection (kapa, &section);
+    psFree(section.name);
+
+    psastroVisualScaleGraphdata (&graphdata, xResid, yResid);
+    KapaSetLimits (kapa, &graphdata);
+    KapaBox (kapa, &graphdata);
+    KapaSendLabel (kapa, "dP / P", KAPA_LABEL_XM);
+    KapaSendLabel (kapa, "dQ / Q", KAPA_LABEL_YM);
+    KapaSendLabel (kapa,
+                   "Residual of the Fit from the Initial Astrometry to the Final Astrometry",
+                   KAPA_LABEL_XP);
+    KapaPrepPlot (kapa, cornerPd->n, &graphdata);
+    KapaPlotVector (kapa, cornerPd->n, xResid->data.F32, "x");
+    KapaPlotVector (kapa, cornerQd->n, yResid->data.F32, "y");
+
+    psFree(xResid);
+    psFree(yResid);
+    psastroVisualAskUser (&plotAstromGuessCheck);
+    return true;
+}
+
+/** Helper Routines **/
+
+
+/* psastroVisualTriplePlot */
+// plot 2 vectors whose point sizes are scaled by a third vector
+// autoscales the plot
+bool psastroVisualTriplePlot (int kapid, Graphdata *graphdata, psVector *xVec, psVector *yVec, psVector *zVec, bool increasing)
+{
+    psastroVisualScaleGraphdata (graphdata, xVec, yVec);
+
+    // set the scale vector
+    psVector *zScale = psVectorAlloc (zVec->n, PS_DATA_F32);
+    psastroVisualCreateScaleVec (zVec, zScale, increasing);
+
+    KapaSetFont (kapid, "helvetica", 14);
+    KapaSetLimits(kapid, graphdata);
+    KapaBox (kapid, graphdata);
+
+    // the point size will be scaled from the z vector
+    graphdata->size = -1;
+    KapaPrepPlot (kapid, xVec->n, graphdata);
+    KapaPlotVector (kapid, xVec->n, xVec->data.F32, "x");
+    KapaPlotVector (kapid, yVec->n, yVec->data.F32, "y");
+    KapaPlotVector (kapid, zVec->n, zScale->data.F32, "z");
+    psFree (zScale);
+    return true;
+}
+
+
+/* psastroVisualTripleOverplot */
+// same as above, but does not rescale the plot or redraw the bounding box
+bool psastroVisualTripleOverplot (int kapid, Graphdata *graphdata, psVector *xVec, psVector *yVec, psVector *zVec, bool increasing) {
+
+    // set the scale vector
+    psVector *zScale = psVectorAlloc (zVec->n, PS_DATA_F32);
+    psastroVisualCreateScaleVec (zVec, zScale, increasing);
+
+    KapaSetFont (kapid, "helvetica", 14);
+
+    // the point size will be scaled from the z vector
+    graphdata->size = -1;
+    KapaPrepPlot (kapid, xVec->n, graphdata);
+    KapaPlotVector (kapid, xVec->n, xVec->data.F32, "x");
+    KapaPlotVector (kapid, yVec->n, yVec->data.F32, "y");
+    KapaPlotVector (kapid, zVec->n, zScale->data.F32, "z");
+    psFree (zScale);
+    return true;
+}
+
+
+/* psastroVisualCreateScaleVec */
+//create a vector of plot point sizes from a vector
+bool psastroVisualCreateScaleVec (psVector *zVec, psVector *zScale, bool increasing) {
+    // set limits based on data values
+    float zmin = +FLT_MAX;
+    float zmax = -FLT_MAX;
+
+    for (int i = 0; i < zVec->n; i++) {
+        zmin = PS_MIN (zmin, zVec->data.F32[i]);
+        zmax = PS_MAX (zmax, zVec->data.F32[i]);
+    }
+
+    float range = zmax - zmin;
+    if (range == 0.0) {
+        psVectorInit (zScale, 1.0);
+    } else {
+        for (int i = 0; i < zVec->n; i++) {
+            if (increasing) {
+                zScale->data.F32[i] = PS_MIN (1.5, PS_MAX(0.05, 1.5*(zVec->data.F32[i] - zmin)/range));
+            } else {
+                zScale->data.F32[i] = PS_MIN (1.5, PS_MAX(0.05, 1.5*(zmax - zVec->data.F32[i])/range));
+            }
+        }
+    }
+    return true;
+}
+
+
+/** psastroVisualScaleGraphdata
+ * Scale the graphdata structure based on x and y coordinates. Use sigma clipping to
+ * prevent outliers from making te plot region too big.
+ */
+bool psastroVisualScaleGraphdata(Graphdata *graphdata, psVector *xVec, psVector *yVec) {
+    graphdata->xmin = +FLT_MAX;
+    graphdata->xmax = -FLT_MAX;
+    graphdata->ymin = +FLT_MAX;
+    graphdata->ymax = -FLT_MAX;
+    bool baddataX = true;
+    bool baddataY = true;
+
+    //determine standard deviation of xVec and yVec
+    psStats *statsX = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+    psStats *statsY = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+    psVectorStats (statsX, xVec, NULL, NULL, 0);
+    psVectorStats (statsY, yVec, NULL, NULL, 0);
+
+    float xhi  = statsX->sampleMedian + 3 *statsX->sampleStdev;
+    float xlo = statsX->sampleMedian - 3 *statsX->sampleStdev;
+    float yhi = statsY->sampleMedian + 3 *statsY->sampleStdev;
+    float ylo = statsY->sampleMedian - 3 *statsY->sampleStdev;
+
+    for(int i = 0; i < xVec->n; i++) {
+        if (!isfinite(xVec->data.F32[i])) continue;
+        if (xVec->data.F32[i] > xhi || xVec->data.F32[i] < xlo) continue;
+        baddataX = false;
+        graphdata->xmin = PS_MIN (graphdata->xmin, xVec->data.F32[i]);
+        graphdata->xmax = PS_MAX (graphdata->xmax, xVec->data.F32[i]);
+    }
+
+    for (int i = 0; i < yVec->n; i++) {
+        if (!isfinite(yVec->data.F32[i])) continue;
+        if (yVec->data.F32[i] > yhi || yVec->data.F32[i] < ylo) continue;
+        baddataY = false;
+        graphdata->ymin = PS_MIN (graphdata->ymin, yVec->data.F32[i]);
+        graphdata->ymax = PS_MAX (graphdata->ymax, yVec->data.F32[i]);
+    }
+
+    // check to see if there was no good data
+    if ( baddataX || baddataY) {
+        graphdata->xmin = -1;
+        graphdata->ymin  = -1;
+        graphdata->xmax = 1;
+        graphdata->ymax = 1;
+        return false;
+    }
+    float range = graphdata->xmax - graphdata->xmin;
+    if (range == 0) range = 1;
+    graphdata->xmin -= .05 * range;
+    graphdata->xmax += .05 * range;
+
+    range = graphdata->ymax - graphdata->ymin;
+    if (range == 0) range = 1;
+    graphdata->ymin -= .05 * range;
+    graphdata->ymax += .05 * range;
+
+    psFree (statsX);
+    psFree (statsY);
+    return true;
+}
+
+# endif
+