/***********************************/
/***Diagnostic plots for psastro****/
/***********************************/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

# include "psastroInternal.h"
# include <string.h>
# include <stdlib.h>
# include <stdio.h>

# if (HAVE_KAPA)
# include <kapa.h>

# define KAPAX  700
# define KAPAY  700

//variables to determine when things are plotted
static bool isVisual             = false;
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 routine prototypes
bool psastroVisualInitGraph (int kapa, KapaSection *section, Graphdata *graphdata);
bool psastroVisualTriplePlot (int kapa, Graphdata *graphdata, psVector *xVec, psVector *yVec, psVector *zVec, bool increasing);
bool psastroVisualScaleGraphdata(Graphdata *graphdata, psVector *xVec, psVector *yVec, bool clip);
bool psastroVisualTripleOverplot (int kapa, Graphdata *graphdata, psVector *xVec, psVector *yVec, psVector *zVec, bool increasing);
bool psastroVisualCreateScaleVec (psVector *zVec, psVector *zScale, bool increasing);
bool psastroVisualResidPlot (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe, char *title);


/*****************************/
/** 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);
        pmAstromSetVisual(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 **/
/***********************/


/**
 * Plot raw stars as determined from first pass astrometry fit
 * Called within psastroAstromGeuss
 */
bool psastroVisualPlotRawStars (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 (!psastroVisualInitWindow (&kapa, "psastro:plots"))
        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);
    psastroVisualTriplePlot (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);
    psastroVisualTriplePlot (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);
    psastroVisualTriplePlot (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);
    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;
}


/**
 * plot the location of references stars over the entire fpa
 * invoked during psastroChooseRefStars
 */
bool psastroVisualPlotRefStars (psArray *refstars, psMetadata *recipe)
{
    //make sure we want to plot this
    if (!isVisual || !plotRefStars) return true;

    //set up plotting variables
    if (!psastroVisualInitWindow (&kapa, "psastro:plots"))
        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);
    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;
}


/**
 * Plot the two luminosity functions created within 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,   ///< 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 ( !psastroVisualInitWindow (&kapa, "psastro:plots"))
        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
    psastroVisualScaleGraphdata(&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:
        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, ///< 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 ( !psastroVisualInitWindow (&kapa, "psastro:plots")) 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
    psastroVisualAskUser (&plotRemoveClumps);
    psFree (xVec);
    psFree (yVec);

    return true;
}


/**
 * Assess the goodness of fit for a signle chip by
 * plotting the fit residuals
 * invoked during psastroOneChipFit
 */
bool psastroVisualPlotOneChipFit (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;

    //plot the residuals
    if (!psastroVisualResidPlot(rawstars, refstars, match, recipe, "Single Chip Fit Residuals (Chip Coordinates)"))
        return false;

    //ask for user input and finish
    psastroVisualAskUser(&plotOneChipFit);
    return true;
}


/**
 * Plots the chip corners in the FP before and after chips with inconsistent solutions have been fixed.
 * Invoked during psastroFixChips
 */
bool psastroVisualPlotFixChips (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(!psastroVisualInitWindow(&kapa, "psastro:plots")) 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
    psastroVisualScaleGraphdata(&graphdata, xOld, yOld, true);
    psastroVisualInitGraph(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");

    psastroVisualAskUser(&plotFixChips);
    psFree(xNew);
    psFree(yNew);
    psFree(view);
    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, ///< 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 ( !psastroVisualInitWindow (&kapa, "psastro:plots"))
        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);

    psastroVisualScaleGraphdata (&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);

    psastroVisualScaleGraphdata (&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);
    psastroVisualAskUser (&plotAstromGuessCheck);
    return true;
}


/**
 * Plots the pixel scales of the fpa before they are
 * equalized in psastroMosaicCommonScale
 */
bool psastroVisualPlotCommonScale (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 false;

    if (!psastroVisualInitWindow(&kapa, "psastro:plots")) 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++;
    }

    //set up plot window
    KapaInitGraph (&graphdata);
    KapaClearPlots (kapa);
    KapaSetSection (kapa, &section);
    KapaSetFont (kapa, "helvetica", 14);
    psastroVisualTriplePlot (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);

    psastroVisualAskUser (&plotCommonScale);
    return true;
}


/**
 *   plot the residuals between raw stars and ref stars after
 *   fitting in psastroMosaicOneChip
 */
bool psastroVisualPlotMosaicOneChip (psArray *rawstars, psArray *refstars,
                                     psArray *match, psMetadata *recipe)
{

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

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

    //ask for user input and finish
    psastroVisualAskUser(&plotMosaicOneChip);

    return true;
}

/** psastroVisualPlotMosaicMatches
 * Plot the matches between raw and reference stars during psastroVisualMosaicSetMatch
 */
bool psastroVisualPlotMosaicMatches( 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 (!psastroVisualResidPlot(rawstars, refstars, match, recipe, title))
        return false;

    //ask for user input
    psastroVisualAskUser (&plotMosaicMatches);
    return true;
}


/*********************/
/** Helper Routines **/
/*********************/


/** psastroVisualInitGraph (kapa, *section, *graphdata)
 * Initializes graph, sets the section, sets the font,
 * sets the limits, draws the box
 */
bool psastroVisualInitGraph (int kapa, KapaSection *section, Graphdata *graphdata)
{
    KapaSetSection (kapa, section);
    KapaSetFont (kapa, "helvetica", 14);
    KapaSetLimits (kapa, graphdata);
    KapaBox (kapa, graphdata);
    return true;
}


/** 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, true);
    //printf("%f %f %f %f \n",graphdata->xmin, graphdata->xmax, graphdata->ymin, graphdata->ymax);
    // 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->ptype = 7;
    graphdata->style = 2;
    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, bool clip) {

    graphdata->xmin = +FLT_MAX;
    graphdata->xmax = -FLT_MAX;
    graphdata->ymin = +FLT_MAX;
    graphdata->ymax = -FLT_MAX;

    //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 = +FLT_MAX, xlo = -FLT_MAX, yhi = +FLT_MAX, ylo = -FLT_MAX;
    if (clip) {
        xhi  = statsX->sampleMedian + 3 *statsX->sampleStdev;
        xlo = statsX->sampleMedian - 3 *statsX->sampleStdev;
        yhi = statsY->sampleMedian + 3 *statsY->sampleStdev;
        ylo = statsY->sampleMedian - 3 *statsY->sampleStdev;
    }

    // abort if there is no good data
    if (!isfinite(xhi) || !isfinite(xlo) || !isfinite(yhi) || !isfinite(ylo)) {
        graphdata->xmin = -1;
        graphdata->ymin  = -1;
        graphdata->xmax = 1;
        graphdata->ymax = 1;
        psFree(statsX);
        psFree(statsY);
        return false;
    }

    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;
        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(xVec->data.F32[i])) continue;
        if (yVec->data.F32[i] > yhi || yVec->data.F32[i] < ylo) continue;
        graphdata->ymin = PS_MIN (graphdata->ymin, yVec->data.F32[i]);
        graphdata->ymax = PS_MAX (graphdata->ymax, yVec->data.F32[i]);
    }

    // add a whitespace border
    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;
}


/** psastroVisualResdiPlot
 * Plots the residuals between matched raw and reference stars. Used to assess the quality of an
 * astrometry solution
 */
bool psastroVisualResidPlot (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe, char *title) {

    //set up the first window
    if (!psastroVisualInitWindow( &kapa, "psastro:plots")) return false;

    //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);
    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, title, 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, 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;
    psastroVisualTripleOverplot (kapa2, &graphdata, xVec, yVec, zVec, false);

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

// END OF PROGRAM
#else

bool psastroSetVisual (bool mode) {return true};
bool psastroVisualClose() {return true};
bool psastroVisualPlotLuminosityFunction (psVector *lnMag, psVector *Mag, pmLumFunc *lumFunc, pmLumFunc *rawFunc) {return true};
bool psastroVisualPlotRawStars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe) {return true};
bool psastroVisualPlotRefStars (psArray *refstars, psMetadata *recipe) {return true};
bool psastroVisualPlotRemoveClumps (psArray *input, psImage *count, int scale, float limit) {return true};
bool psastroVisualPlotFixChips (pmFPAfile *input, psVector *xOld, psVector *yOld) {return true};
bool psastroVisualPlotOneChipFit (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe) {return true};
bool psastroVisualPlotAstromGuessCheck (psVector *cornerPo, psVector *cornerQo, psVector *cornerPn, psVector *cornerQn, psVector *cornerPd, psVector *cornerQd) {return true};
bool psastroVisualPlotMosaicOneChip (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe) {return true};
bool psastroVisualPlotCommonScale (pmFPA *fpa, psVector *oldScale) {return true};
bool psastroVisualPlotMosaicMatches (psArray *rawstars, psArray *refstars, psArray *match, int iteration, psMetadata *recipe) {return true};

# endif

