#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include <strings.h>
#include <string.h>
#include <math.h>
#include <assert.h>
#include <pslib.h>

#include "pmHDU.h"
#include "pmFPA.h"
#include "pmFPAfile.h"
#include "pmAstrometryObjects.h"
#include "pmAstrometryVisual.h"
#include "pmFPAExtent.h"

# if (HAVE_KAPA)
# include <kapa.h>
# include "pmKapaPlots.h"
#include "pmVisual.h"

//variables to determine when things are plotted
static bool isVisual             = false;
static bool plotGridMatch        = true;
static bool plotTweak            = true;
static bool plotRawStars         = true;
static bool plotRefStars         = false;
static bool plotLumFunc          = true;
static bool plotRemoveClumps     = true;
static bool plotOneChipFit       = true;
static bool plotFixChips         = true;
static bool plotAstromGuessCheck = true;
static bool plotMosaicMatches    = true;
static bool plotCommonScale      = true;
static bool plotMosaicOneChip    = true;

// variables to store plotting window indices
static int kapa = -1;
static int kapa2 = -1;

//helper prototypes
bool residPlot (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe,
                char *title);


/* Initialization Routines  */

bool pmAstromSetVisual (bool mode) {
    isVisual = mode;
    return true;
}


bool pmAstromVisualClose()
{
    if(kapa != -1)
        KiiClose(kapa);
    if(kapa2 != -1)
        KiiClose(kapa2);
    return true;
}


/* Plotting Routines */
bool pmAstromVisualPlotRawStars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe)
{
    // make sure we want to plot this
    if (!plotRawStars || !isVisual) return true;

    //set up plot region
    if (!pmVisualInitWindow (&kapa, "psastro:plots")){
        isVisual = false;
        return false;
    }
    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;

    //initialize and populate plotting vectors
    bool status = false;
    float iMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MIN");
    float iMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.INST.MAG.MAX");

    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);

    //Chip coordinates
    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);
    pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false);

    //Focal Plane Coordinates
    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);
    pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false);

    // Tangent Plane Coordinates
    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);
    pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, zVec, false);

    //sky coordinates
    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);
    pmVisualTriplePlot (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:
    pmVisualAskUser(&plotRawStars, &isVisual);

    psFree (xVec);
    psFree (yVec);
    psFree (zVec);
    return true;
}


bool pmAstromVisualPlotRefStars (psArray *refstars, psMetadata *recipe)
{
    //make sure we want to plot this
    if (!isVisual || !plotRefStars) return true;

    //set up plotting variables
    if (!pmVisualInitWindow (&kapa, "psastro:plots")) {
        isVisual = false;
        return false;
    }

    Graphdata graphdata;
    KapaInitGraph (&graphdata);
    KapaClearSections (kapa);

    graphdata.color = KapaColorByName ("black");
    graphdata.ptype = 7;
    graphdata.size = 0.5;
    graphdata.style = 2;

    //initialize and populate plot vectors
    bool status = false;
    float rMagMin = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MIN");
    float rMagMax = psMetadataLookupF32 (&status, recipe, "PSASTRO.PLOT.REF.MAG.MAX");

    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);
    pmVisualTriplePlot(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:
    pmVisualAskUser(&plotRefStars, &isVisual);

    psFree (xVec);
    psFree (yVec);
    psFree (zVec);
    return true;
}


bool pmAstromVisualPlotLuminosityFunction (psVector *lnMag,   // Log(n) for each magnitude bin
                                          psVector *Mag,     // magnitude bins
                                          pmLumFunc *lumFunc,// Fit to the reference star luminosity function
                                          pmLumFunc *rawFunc // Fit to the raw star luminoisty function
                                          )
{

    // make sure we want to plot this
    if ( !isVisual || !plotLumFunc ) return true;

    //set up plotting variables
    if ( !pmVisualInitWindow (&kapa, "psastro:plots")){
        isVisual = false;
        return false;
    }

    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
    pmVisualScaleGraphdata(&graphdata, Mag, lnMag, false);

    //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:
        pmVisualAskUser (&plotLumFunc, &isVisual);
    }
    return true;
} // end of pmAstromVisualPlotLuminosityFunction


bool pmAstromVisualPlotRemoveClumps (psArray *input, // Array containing the field stars
                                    psImage *count, // A 2D histogram of the field star distribution
                                    int scale,      // The pixel size of the histogram
                                    float limit     // The minimum numuber of stars in a bin flagged as a clump
                                    )
{

    //make sure we want to plot this
    if (!isVisual || !plotRemoveClumps) return true;

    //set up plot variables
    if ( !pmVisualInitWindow (&kapa, "psastro:plots")) {
        isVisual = false;
        return false;
    }

    KapaSection section;
    Graphdata graphdata;
    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");
    KapaClearPlots(kapa);

    //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");
        }
    }

    //ask for user input and finish
    pmVisualAskUser (&plotRemoveClumps, &isVisual);
    psFree (xVec);
    psFree (yVec);

    return true;
}


bool pmAstromVisualPlotOneChipFit (psArray *rawstars, // stars detected in the image
                                  psArray *refstars, // reference stars over the same region
                                  psArray *match,    // contains which rawstars match to which refstars
                                  psMetadata *recipe // data reduction recipe
                                  )
{

    //make sure we want to plot this
    if (!isVisual || !plotOneChipFit)
        return true;

    if(!pmVisualInitWindow(&kapa, "psastro:plots") || !pmVisualInitWindow(&kapa2, "psastro:plots")) {
        isVisual = false;
        return false;
    }

    //plot the residuals
    if (!residPlot(rawstars, refstars, match, recipe, "Single Chip Fit Residuals (Chip Coordinates)")) {
        plotOneChipFit = false;
        return false;
    }

    //ask for user input and finish
    pmVisualAskUser(&plotOneChipFit, &isVisual);
    return true;
}


bool pmAstromVisualPlotFixChips (pmFPAfile *input, // focal plane array file
                                psVector *xOld, // old X location of chip cornerss
                                psVector *yOld // old Y location of chip corners
                                )
{
    //make sure we want to plot this
    if(!isVisual || !plotFixChips) return true;

    if(!pmVisualInitWindow(&kapa, "psastro:plots")) {
        isVisual = false;
        return false;
    }

    KapaSection section = {"s1", .05, .05, .9, .9};
    Graphdata graphdata;
    KapaInitGraph( &graphdata);
    KapaClearPlots (kapa);
    graphdata.ptype = 2;
    graphdata.style = 2;

    psVector *xNew = psVectorAllocEmpty (xOld->n, PS_TYPE_F32);
    psVector *yNew = psVectorAllocEmpty (yOld->n, PS_TYPE_F32);

    // copy of the code in psastroFixChips that generated xOld, yOld, but for xNew, yNew
    pmFPAview *view = pmFPAviewAlloc (0);

    pmChip *obsChip = NULL;
    while ((obsChip = pmFPAviewNextChip (view, input->fpa, 1)) != NULL) {
        if (!obsChip->process || !obsChip->file_exists || !obsChip->data_exists) { continue; }

        psRegion *region = pmChipPixels(obsChip);
        psPlane ptCP, ptFP;

        ptCP.x = region->x0; ptCP.y = region->y0;
        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
        psVectorAppend (xNew, ptFP.x);
        psVectorAppend (yNew, ptFP.y);

        ptCP.x = region->x0; ptCP.y = region->y1;
        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
        psVectorAppend (xNew, ptFP.x);
        psVectorAppend (yNew, ptFP.y);

        ptCP.x = region->x1; ptCP.y = region->y1;
        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
        psVectorAppend (xNew, ptFP.x);
        psVectorAppend (yNew, ptFP.y);

        ptCP.x = region->x1; ptCP.y = region->y0;
        psPlaneTransformApply (&ptFP, obsChip->toFPA, &ptCP);
        psVectorAppend (xNew, ptFP.x);
        psVectorAppend (yNew, ptFP.y);

        psFree (region);
    }

    //set up graph
    pmVisualScaleGraphdata(&graphdata, xOld, yOld, true);
    pmVisualInitGraph(kapa, &section, &graphdata);
    KapaSendLabel (kapa, "L (FP)", KAPA_LABEL_XM);
    KapaSendLabel (kapa, "M (FP)", KAPA_LABEL_YM);
    KapaSendLabel (kapa, "Chip corners before (black) and after (red) FixChips", KAPA_LABEL_XP);
    KapaPrepPlot (kapa, xOld->n, &graphdata);
    KapaPlotVector (kapa, xOld->n, xOld->data.F32, "x");
    KapaPlotVector (kapa, xOld->n, yOld->data.F32, "y");

    graphdata.ptype = 1;
    graphdata.color = KapaColorByName("red");
    KapaPrepPlot (kapa, xNew->n, &graphdata);
    KapaPlotVector (kapa, xNew->n, xNew->data.F32, "x");
    KapaPlotVector (kapa, xNew->n, yNew->data.F32, "y");

    pmVisualAskUser(&plotFixChips, &isVisual);
    psFree(xNew);
    psFree(yNew);
    psFree(view);
    return true;
}


bool pmAstromVisualPlotAstromGuessCheck (psVector *cornerPo, // P coordinates of chip corners before fitting
                                        psVector *cornerQo, // Q coordinates of chip corners before fitting
                                        psVector *cornerPn, // P coordinates of chip corners after fitting
                                        psVector *cornerQn, // Q coordinates of chip corners after fitting
                                        psVector *cornerPd, // P coordinate residuals of fit from old to new coordinates
                                        psVector *cornerQd  // Q coordinate residuals of fit from old to new coordinates
                                        )
{

    //make sure we want to plot this
    if (!isVisual || !plotAstromGuessCheck) return true;

    //set up graph window
    if ( !pmVisualInitWindow (&kapa, "psastro:plots")) {
        isVisual = false;
        return false;
    }

    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);

    pmVisualScaleGraphdata (&graphdata, cornerPo, cornerPo, true);
    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]);
        yResid->data.F32[i] = (cornerQd->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);

    pmVisualScaleGraphdata (&graphdata, xResid, yResid, true);
    KapaSetLimits (kapa, &graphdata);
    KapaBox (kapa, &graphdata);
    KapaSendLabel (kapa, "dP", KAPA_LABEL_XM);
    KapaSendLabel (kapa, "dQ", 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);
    pmVisualAskUser (&plotAstromGuessCheck, &isVisual);
    return true;
}


bool pmAstromVisualPlotCommonScale (pmFPA *fpa,         // the fpa
                                   psVector *oldScale  // the old pixel scale of each chip in the fpa
                                   )
{
    //make sure we want to plot this
    if (!isVisual || !plotCommonScale) return true;

    if (!pmVisualInitWindow(&kapa, "psastro:plots")){
        isVisual = false;
        return false;
    }

    KapaSection section = {"s1", .05, .05, .9, .9};
    Graphdata graphdata;
    psPlane ptCH, ptFP;
    ptCH.x = 0;
    ptCH.y = 0;
    psVector *xVec = psVectorAlloc (oldScale->n, PS_TYPE_F32);
    psVector *yVec = psVectorAlloc (oldScale->n, PS_TYPE_F32);

    int nobj = 0;

    //project each chip corner to the Focal Plane
    for(int i = 0; i < fpa->chips->n; i++) {
        pmChip *chip = fpa->chips->data[i];
        if (!chip->process || !chip->file_exists) { continue; }
        if (!chip->toFPA) { continue; }

        psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
        xVec->data.F32[nobj] = ptFP.x;
        yVec->data.F32[nobj] = ptFP.y;
        nobj++;
        if (nobj == oldScale->n) break;
    }

    //set up plot window
    KapaInitGraph (&graphdata);
    KapaClearPlots (kapa);
    KapaSetSection (kapa, &section);
    KapaSetFont (kapa, "helvetica", 14);
    pmVisualTriplePlot (kapa, &graphdata, xVec, yVec, oldScale, false);
    KapaSendLabel (kapa, "L (FP)", KAPA_LABEL_XM);
    KapaSendLabel (kapa, "M (FP)", KAPA_LABEL_YM);
    KapaSendLabel (kapa, "Old Pixel Scale of FPA Chips (Not to Scale)", KAPA_LABEL_XP);

    pmVisualAskUser (&plotCommonScale, &isVisual);

    psFree(xVec);
    psFree(yVec);

    return true;
}


bool pmAstromVisualPlotMosaicOneChip (psArray *rawstars, psArray *refstars,
                                     psArray *match, psMetadata *recipe)
{

    //make sure we want to plot this
    if (!isVisual || !plotMosaicOneChip) return true;

    if(!pmVisualInitWindow(&kapa, "psastro:plots") || !pmVisualInitWindow(&kapa2, "psastro:plots")) {
        isVisual = false;
        return false;
    }

    //plot the residuals
    if (!residPlot(rawstars, refstars, match, recipe, "Single Chip Fit Residuals - Mosaic Mode")) {
        isVisual = false;
        return false;
    }

    //ask for user input and finish
    pmVisualAskUser(&plotMosaicOneChip, &isVisual);

    return true;
}


bool pmAstromVisualPlotMosaicMatches( psArray *rawstars, psArray *refstars,
                                    psArray *match, int iteration,
                                    psMetadata *recipe)
{
    //make sure we want to plot this
    if (!isVisual || !plotMosaicMatches) return true;

    char title[60];
    sprintf(title, "Matches found during psastroMosaicSetMatch iteration %d", iteration);


    if(!pmVisualInitWindow(&kapa, "psastro:plots") || !pmVisualInitWindow(&kapa2, "psastro:plots")) {
        isVisual = false;
        return false;
    }

    if (!residPlot(rawstars, refstars, match, recipe, title)){
        isVisual = false;
        return false;
    }

    //ask for user input
    pmVisualAskUser (&plotMosaicMatches, &isVisual);
    return true;
}


bool pmAstromVisualPlotGridMatch (const psArray *raw,
                                  const psArray *ref,
                                  psImage *gridNP,
                                  double offsetX,
                                  double offsetY,
                                  double maxOffpix,
                                  double Scale,
                                  double Offset)
{
    //make sure we want to plot this
    if (!isVisual || !plotGridMatch) return true;
    if (!pmVisualInitWindow(&kapa, "psastro:plots")){
        isVisual = false;
        return false;
    }

    KapaSection section = {"s1", 0.05, 0.05, .75, .75};
    KapaSection sectionY = {"s2", 0.8, 0.05, .15, .75};
    KapaSection sectionX = {"s3", .05, .8, .75, .15};

    Graphdata graphdata;
    int nplot = raw->n * ref->n;                      // number of points to plot
    psVector *dXplot = psVectorAlloc (nplot, PS_TYPE_F32); // x data points
    psVector *dYplot = psVectorAlloc (nplot, PS_TYPE_F32); // y data points

    pmAstromObj *ob1; pmAstromObj *ob2;               // shortcuts to the data in raw and ref
    psU32 **NP = gridNP->data.U32;                    // shortcut to the gridNP data
    float vertHistSlice[gridNP->numRows];             // vertical histogram slice through peak
    float horizHistSlice[gridNP->numCols];            // horizontal histogram slice through peak
    float horizontalIndices[gridNP->numCols];         // the horizontal offset corresponding to each bin of horizHistSlice
    float verticalIndices[gridNP->numRows];           // the vertical offset corresponding to each bin of vertHistSlice
    int maxHorizontalSlice = 0;                       // peak value of horizHistSlice
    int maxVerticalSlice = 0;                         // peak value of vertHistSlice
    int peakXbin = (int) (offsetX / Scale + Offset);  // X bin index of peak
    int peakYbin = (int) (offsetY / Scale + Offset);  // Y bin index of peak

    // psVector *vertHistSlice = psVectorAlloc (gridNP->numRows, PS_TYPE_F32);

    // set up plot information
    KapaClearPlots(kapa);
    KapaInitGraph(&graphdata);
    KapaSetSection(kapa, &section);

    graphdata.xmin = -1.0 * maxOffpix;
    graphdata.xmax =  1.0 * maxOffpix;
    graphdata.ymin = -1.0 * maxOffpix;
    graphdata.ymax =  1.0 * maxOffpix;
    KapaSetLimits(kapa, &graphdata);

    KapaSetFont(kapa, "helvetica", 14);
    KapaBox(kapa, &graphdata);
    KapaSendLabel (kapa, "X offset (FP)", KAPA_LABEL_XM);
    KapaSendLabel (kapa, "Y offset (FP)", KAPA_LABEL_YM);
    KapaSendLabel (kapa, "pmAstromGridAngle residuals. Box: Correlation Peak.",
                   KAPA_LABEL_XP);
    graphdata.style = 2;
    graphdata.ptype = 0;
    graphdata.size = 0.4;
    graphdata.color = KapaColorByName ("black");

    // calculate the plot points
    float dX, dY;
    for (int i = 0; i < raw->n; i++) {
        ob1 = (pmAstromObj *)raw->data[i];
        for (int j = 0; j < ref->n; j++) {
            ob2 = (pmAstromObj *)ref->data[j];
            dX = ob1->FP->x - ob2->FP->x;
            dY = ob1->FP->y - ob2->FP->y;
            dXplot->data.F32[(i * ref->n) + j] = dX;
            dYplot->data.F32[(i * ref->n) + j] = dY;
        }
    }

    // calculate the points for the profiles
    for (int i = 0; i < gridNP->numRows; i++) {
        vertHistSlice[i] = NP[i][peakXbin];
        verticalIndices[i] = (i - Offset) * Scale;
        if (vertHistSlice[i] > maxVerticalSlice) {
            maxVerticalSlice = vertHistSlice[i];
        }
    }
    for (int i = 0; i < gridNP->numCols; i++) {
        horizHistSlice[i] = NP[peakYbin][i];
        horizontalIndices[i] = (i - Offset) * Scale;
        if (horizHistSlice[i] > maxHorizontalSlice) {
            maxHorizontalSlice = horizHistSlice[i];
        }
    }

    //Plot the offsets
    KapaPrepPlot(kapa, nplot, &graphdata);
    KapaPlotVector (kapa, nplot, dXplot->data.F32, "x");
    KapaPlotVector (kapa, nplot, dYplot->data.F32, "y");

    //Overplot bounding box, peak of distribution
    float xbound[5] = { -maxOffpix, maxOffpix, maxOffpix, -maxOffpix, -maxOffpix};
    float ybound[5] = { -maxOffpix, -maxOffpix, maxOffpix, maxOffpix, -maxOffpix};
    float xbin[5] = {offsetX - 0.5 * Scale, offsetX + 0.5 * Scale,
                     offsetX + 0.5 * Scale, offsetX - 0.5 * Scale,
                     offsetX - 0.5 * Scale};
    float ybin[5] = {offsetY - 0.5 * Scale, offsetY - 0.5 * Scale,
                     offsetY + 0.5 * Scale, offsetY + 0.5 * Scale,
                     offsetY - 0.5 * Scale};
    graphdata.color = KapaColorByName("red");
    graphdata.style = 0;
    graphdata.size = 1.0;
    KapaPrepPlot(kapa, 5, &graphdata);
    KapaPlotVector (kapa, 5, xbound, "x");
    KapaPlotVector (kapa, 5, ybound, "y");
    KapaPrepPlot(kapa, 5, &graphdata);
    KapaPlotVector (kapa, 5, xbin, "x");
    KapaPlotVector (kapa, 5, ybin, "y");

    //plot X profile
    KapaSetSection(kapa, &sectionX);
    graphdata.color = KapaColorByName("black");
    graphdata.ptype = 1;
    graphdata.style = 1;
    graphdata.ymin = 0;
    graphdata.ymax = maxHorizontalSlice + 0.5;
    KapaSetLimits(kapa, &graphdata);

    KapaBox(kapa, &graphdata);
    KapaPrepPlot(kapa, gridNP->numCols, &graphdata);
    KapaPlotVector (kapa, gridNP->numCols, horizontalIndices, "x");
    KapaPlotVector (kapa, gridNP->numCols, horizHistSlice, "y");
    float xslice[2] = {offsetX - Scale / 2., offsetX - Scale / 2.};
    float yslice[2] = {-5, 100};
    graphdata.color = KapaColorByName("red");
    KapaPrepPlot(kapa, 2, &graphdata);
    KapaPlotVector (kapa, 2, xslice, "x");
    KapaPlotVector (kapa, 2, yslice, "y");

    //plot Y profile
    KapaSetSection(kapa, &sectionY);
    graphdata.color = KapaColorByName("black");
    graphdata.ptype = 1;
    graphdata.style = 1;
    graphdata.ymin = -maxOffpix;
    graphdata.ymax = maxOffpix;
    graphdata.xmin = -1.0 ;
    graphdata.xmax = maxVerticalSlice + 0.5;
    KapaSetLimits(kapa, &graphdata);

    KapaBox(kapa, &graphdata);
    KapaPrepPlot(kapa, gridNP->numRows, &graphdata);
    KapaPlotVector (kapa, gridNP->numRows, vertHistSlice, "x");
    KapaPlotVector (kapa, gridNP->numRows, verticalIndices, "y");
    yslice[0] = yslice[1] = offsetY - Scale / 2.;
    xslice[0] = -5; xslice[1] = 100;
    graphdata.color = KapaColorByName("red");
    KapaPrepPlot(kapa, 2, &graphdata);
    KapaPlotVector (kapa, 2, xslice, "x");
    KapaPlotVector (kapa, 2, yslice, "y");

    pmVisualAskUser(&plotGridMatch, &isVisual);
    psFree(dXplot);
    psFree(dYplot);
    return true;
} // end of pmAstromVisualPlotGridMatch


bool pmAstromVisualPlotTweak (psVector *xHist, // Smoothed Horizontal cut through the histogram
                              psVector *yHist, // Smoothed Vertical cut throug the histogram
                              int xBin,        // X Bin index of the histogram peak
                              int yBin         // Y bin index of the histogram peak
    )
{
    //make sure we want to plot this
    if (!isVisual || !plotTweak) return true;
    if (!pmVisualInitWindow(&kapa, "psastro:plots")) {
        isVisual = false;
        return false;
    }

    Graphdata graphdata;
    KapaSection section1 = {"s1", 0.05, 0.05, 0.90, 0.4};
    KapaSection section2 = {"s2", 0.05, 0.5, 0.90, 0.4};

    psVector *xIndices = psVectorAlloc (xHist->n, PS_TYPE_F32);
    psVector *yIndices = psVectorAlloc (yHist->n, PS_TYPE_F32);

    //populate the Indices vectors
    for(int i = 0; i < xHist->n; i++) {
        xIndices->data.F32[i] = i;
    }
    for(int i = 0; i < yHist->n; i++) {
        yIndices->data.F32[i] = i;
    }

    // set up plot information
    KapaClearPlots(kapa);
    KapaInitGraph(&graphdata);

    // plot the X histogram
    pmVisualScaleGraphdata(&graphdata, xIndices, xHist, false);
    KapaSetSection(kapa, &section1);
    KapaSetLimits (kapa, &graphdata);
    KapaSetFont(kapa, "helvetica", 14);
    KapaBox(kapa, &graphdata);
    KapaSendLabel (kapa, "X offset Bin", KAPA_LABEL_XM);
    KapaSendLabel (kapa, "Number of Sources", KAPA_LABEL_YM);
    KapaSendLabel (kapa, "Horizontal Profile",
                   KAPA_LABEL_XP);
    graphdata.style = 1;
    graphdata.ptype = 0;
    graphdata.size = 0.4;
    graphdata.color = KapaColorByName ("black");

    KapaPrepPlot (kapa, xHist->n, &graphdata);
    KapaPlotVector (kapa, xHist->n, xIndices->data.F32, "x");
    KapaPlotVector (kapa, xHist->n, xHist->data.F32, "y");

    //overplot the peak
    float x[2] = {xBin, xBin};
    float y[2] = {-500, 500};
    graphdata.color = KapaColorByName ("red");
    KapaPrepPlot (kapa, 2, &graphdata);
    KapaPlotVector (kapa, 2, x, "x");
    KapaPlotVector (kapa, 2, y, "y");

    //plot the Y histogram
    pmVisualScaleGraphdata(&graphdata, yIndices, yHist, false);
    KapaSetSection(kapa, &section2);
    KapaSetLimits (kapa, &graphdata);
    KapaSetFont(kapa, "helvetica", 14);
    graphdata.color = KapaColorByName ("black");
    KapaBox(kapa, &graphdata);
    KapaSendLabel (kapa, "Y offset Bin", KAPA_LABEL_XM);
    KapaSendLabel (kapa, "Number of Sources", KAPA_LABEL_YM);
    KapaSendLabel (kapa, "Vertical Profile",
                   KAPA_LABEL_XP);
    graphdata.style = 1;
    graphdata.ptype = 0;
    graphdata.size = 0.4;

    KapaPrepPlot (kapa, yHist->n, &graphdata);
    KapaPlotVector (kapa, yHist->n, yIndices->data.F32, "x");
    KapaPlotVector (kapa, yHist->n, yHist->data.F32, "y");

    //overplot the peak
    x[0] = x[1] = yBin;
    graphdata.color = KapaColorByName ("red");
    KapaPrepPlot (kapa, 2, &graphdata);
    KapaPlotVector (kapa, 2, x, "x");
    KapaPlotVector (kapa, 2, y, "y");

    //plot title
    graphdata.color = KapaColorByName("black");
    KapaSection section3 = {"s3", 0, 0, 1, 1};
    KapaSetSection( kapa, &section3);
    KapaSendLabel (kapa, "Tweaking the Astrometry Grid Solution. Smoothed profiles + peak location",
                   KAPA_LABEL_XP);

    pmVisualAskUser(&plotTweak, &isVisual);

    psFree(xIndices);
    psFree(yIndices);
    return true;
} //end of pmAstromPlotTweak


bool residPlot (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe,
                        char *title) {


    //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;

    //initialize and populate the plotting vectors
    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");

    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);
    pmVisualTriplePlot (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);
    pmVisualTriplePlot (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);
    pmVisualTriplePlot (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, title, KAPA_LABEL_XP);
    psFree (section.name);

    //second window

    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);
    pmVisualTriplePlot (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;
    pmVisualTripleOverplot (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;
    pmVisualScaleGraphdata(&graphdata, xVec, yVec, true);
    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;
    pmVisualTripleOverplot (kapa2, &graphdata, xVec, yVec, zVec, false);

    psFree (xVec);
    psFree (yVec);
    psFree (zVec);
    return true;
}




# else

bool pmAstromSetVisual(bool mode) { return true; }
bool pmAstromVisualClose() { return true; }
bool pmAstromVisualPlotGridMatch (const psArray *raw, const psArray *ref, psImage *gridNP, double offsetX, double offsetY, double maxOffpix, double Scale, double Offset) { return true; }
bool pmAstromVisualPlotTweak (psVector *xHist, psVector *yHist, int xBin, int yBin) {return true;}
bool pmAstromVisualPlotLuminosityFunction (psVector *lnMag, psVector *Mag, pmLumFunc *lumFunc, pmLumFunc *rawFunc) {return true;}
bool pmAstromVisualPlotRawStars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe) {return true;}
bool pmAstromVisualPlotRefStars (psArray *refstars, psMetadata *recipe) {return true;}
bool pmAstromVisualPlotRemoveClumps (psArray *input, psImage *count, int scale, float limit) {return true;}
bool pmAstromVisualPlotFixChips (pmFPAfile *input, psVector *xOld, psVector *yOld) {return true;}
bool pmAstromVisualPlotOneChipFit (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe) {return true;}
bool pmAstromVisualPlotAstromGuessCheck (psVector *cornerPo, psVector *cornerQo, psVector *cornerPn, psVector *cornerQn, psVector *cornerPd, psVector *cornerQd) {return true;}
bool pmAstromVisualPlotMosaicOneChip (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe) {return true;}
bool pmAstromVisualPlotCommonScale (pmFPA *fpa, psVector *oldScale) {return true;}
bool pmAstromVisualPlotMosaicMatches (psArray *rawstars, psArray *refstars, psArray *match, int iteration, psMetadata *recipe) {return true;}

# endif
