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

#include <pslib.h>
#include "pmTrend2D.h"
#include "pmSourceVisual.h"

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

// functions used to visualize the analysis as it goes
// these are invoked by the -visual options

static bool isVisual = false;
static int kapa1 = -1;
// static int kapa2 = -1;
// static int kapa3 = -1;

bool pmSourcePlotPoints3D (int myKapa, Graphdata *graphdata, psVector *xn, psVector *yn, psVector *zn, float theta, float phi);

bool pmSourceSetVisual (bool mode) {

    isVisual = mode;
    return true;
}

bool pmSourceVisualPSFModelResid (pmTrend2D *trend, psVector *x, psVector *y, psVector *param, psVector *mask) {

    KapaSection section;  // put the positive profile in one and the residuals in another? 

    Graphdata graphdata;

    if (!isVisual) return true;

    if (kapa1 == -1) {
        kapa1 = KapaOpenNamedSocket ("kapa", "pmSource:plots");
	if (kapa1 == -1) {
	    fprintf (stderr, "failure to open kapa; visual mode disabled\n");
	    isVisual = false;
	    return false;
	}
    }  

    KapaClearPlots (kapa1);
    KapaInitGraph (&graphdata);

    float min = +1e32;
    float max = -1e32;
    float Min = +1e32;
    float Max = -1e32;

    psVector *resid = psVectorAlloc (x->n, PS_TYPE_F32);
    psVector *model = psVectorAlloc (x->n, PS_TYPE_F32);

    for (int i = 0; i < x->n; i++) {
	model->data.F32[i] = pmTrend2DEval (trend, x->data.F32[i], y->data.F32[i]);
	resid->data.F32[i] = param->data.F32[i] - model->data.F32[i];
	if (mask->data.U8[i]) continue;
	min = PS_MIN (min, resid->data.F32[i]);
	max = PS_MAX (max, resid->data.F32[i]);
	Min = PS_MIN (min, param->data.F32[i]);
	Max = PS_MAX (max, param->data.F32[i]);
    }

    psVector *xn = psVectorAlloc (x->n, PS_TYPE_F32);
    psVector *yn = psVectorAlloc (x->n, PS_TYPE_F32);
    psVector *zn = psVectorAlloc (x->n, PS_TYPE_F32);
    psVector *Zn = psVectorAlloc (x->n, PS_TYPE_F32);
    psVector *Fn = psVectorAlloc (x->n, PS_TYPE_F32);
    for (int i = 0; i < x->n; i++) {
	xn->data.F32[i] = x->data.F32[i] / 5000.0;
	yn->data.F32[i] = y->data.F32[i] / 5000.0;
	zn->data.F32[i] = (resid->data.F32[i] - min) / (max - min);
	Zn->data.F32[i] = (param->data.F32[i] - Min) / (Max - Min);
	Fn->data.F32[i] = (model->data.F32[i] - Min) / (Max - Min);
    }

    // view 1 on resid
    section.dx = 0.5;
    section.dy = 0.33;
    section.x = 0.0;
    section.y = 0.0;
    section.name = NULL;
    psStringAppend (&section.name, "a1");
    KapaSetSection (kapa1, &section);
    psFree (section.name);
    pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, zn, 30.0*PS_RAD_DEG, -15.0*PS_RAD_DEG);

    // view 2 on resid
    section.dx = 0.5;
    section.dy = 0.33;
    section.x = 0.5;
    section.y = 0.0;
    section.name = NULL;
    psStringAppend (&section.name, "a2");
    KapaSetSection (kapa1, &section);
    psFree (section.name);
    pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, zn, -60.0*PS_RAD_DEG, -15.0*PS_RAD_DEG);

    // view 3 on resid
    section.dx = 0.5;
    section.dy = 0.33;
    section.x = 0.0;
    section.y = 0.33;
    section.name = NULL;
    psStringAppend (&section.name, "a3");
    KapaSetSection (kapa1, &section);
    psFree (section.name);
    pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, Zn, 30.0*PS_RAD_DEG, -15.0*PS_RAD_DEG);

    // view 4 on resid
    section.dx = 0.5;
    section.dy = 0.33;
    section.x = 0.5;
    section.y = 0.33;
    section.name = NULL;
    psStringAppend (&section.name, "a4");
    KapaSetSection (kapa1, &section);
    psFree (section.name);
    pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, Zn, -60.0*PS_RAD_DEG, -15.0*PS_RAD_DEG);

    // view 5 on resid
    section.dx = 0.5;
    section.dy = 0.33;
    section.x = 0.0;
    section.y = 0.66;
    section.name = NULL;
    psStringAppend (&section.name, "a5");
    KapaSetSection (kapa1, &section);
    psFree (section.name);
    pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, Fn, 30.0*PS_RAD_DEG, -15.0*PS_RAD_DEG);

    // view 6 on resid
    section.dx = 0.5;
    section.dy = 0.33;
    section.x = 0.5;
    section.y = 0.66;
    section.name = NULL;
    psStringAppend (&section.name, "a6");
    KapaSetSection (kapa1, &section);
    psFree (section.name);
    pmSourcePlotPoints3D (kapa1, &graphdata, xn, yn, Fn, -60.0*PS_RAD_DEG, -15.0*PS_RAD_DEG);

    psFree (resid);

    psFree (xn);
    psFree (yn);
    psFree (zn);
    psFree (Zn);

    // pause and wait for user input:
    // continue, save (provide name), ??
    char key[10];
    fprintf (stdout, "[c]ontinue? ");
    if (!fgets(key, 8, stdin)) {
	psWarning("Unable to read option");
    }
    return true;
}

// send in normalized points
bool pmSourcePlotPoints3D (int myKapa, Graphdata *graphdata, psVector *xn, psVector *yn, psVector *zn, float theta, float phi) {

    psVector *xv = psVectorAlloc (PS_MAX(6, 2*xn->n), PS_TYPE_F32);
    psVector *yv = psVectorAlloc (PS_MAX(6, 2*xn->n), PS_TYPE_F32);
    psVector *zv = psVectorAlloc (PS_MAX(6, 2*xn->n), PS_TYPE_F32);

    graphdata->xmin = +1e32;
    graphdata->xmax = -1e32;
    graphdata->ymin = +1e32;
    graphdata->ymax = -1e32;

    for (int i = 0; i < xn->n; i++) {
	xv->data.F32[2*i+0] = +xn->data.F32[i]*cos(theta) + yn->data.F32[i]*sin(theta)*cos(phi) + zn->data.F32[i]*sin(theta)*sin(phi);
	yv->data.F32[2*i+0] = -xn->data.F32[i]*sin(theta) + yn->data.F32[i]*cos(theta)*cos(phi) + zn->data.F32[i]*cos(theta)*sin(phi);
	zv->data.F32[2*i+0] = -yn->data.F32[i]*sin(phi)   + zn->data.F32[i]*cos(phi);
	xv->data.F32[2*i+1] = +xn->data.F32[i]*cos(theta) + yn->data.F32[i]*sin(theta)*cos(phi);
	yv->data.F32[2*i+1] = -xn->data.F32[i]*sin(theta) + yn->data.F32[i]*cos(theta)*cos(phi);
	zv->data.F32[2*i+1] = -yn->data.F32[i]*sin(phi);
	graphdata->xmin = PS_MIN(graphdata->xmin, xv->data.F32[2*i+0]);
	graphdata->xmax = PS_MAX(graphdata->xmax, xv->data.F32[2*i+0]);
	graphdata->ymin = PS_MIN(graphdata->ymin, zv->data.F32[2*i+0]);
	graphdata->ymax = PS_MAX(graphdata->ymax, zv->data.F32[2*i+0]);
	graphdata->xmin = PS_MIN(graphdata->xmin, xv->data.F32[2*i+1]);
	graphdata->xmax = PS_MAX(graphdata->xmax, xv->data.F32[2*i+1]);
	graphdata->ymin = PS_MIN(graphdata->ymin, zv->data.F32[2*i+1]);
	graphdata->ymax = PS_MAX(graphdata->ymax, zv->data.F32[2*i+1]);
    }
    xv->n = xn->n;

    // examine sources to set data range
    KapaSetLimits (myKapa, graphdata);

    // KapaSetFont (myKapa, "helvetica", 14);
    // KapaBox (myKapa, graphdata);
    // KapaSendLabel (myKapa, "&ss&h_x| (pixels)", KAPA_LABEL_XM);
    // KapaSendLabel (myKapa, "&ss&h_y| (pixels)", KAPA_LABEL_YM);

    graphdata->color = KapaColorByName ("black");
    graphdata->ptype = 100;
    graphdata->size = 0.5;
    graphdata->style = 2;
    KapaPrepPlot (myKapa, xv->n, graphdata);
    KapaPlotVector (myKapa, xv->n, xv->data.F32, "x");
    KapaPlotVector (myKapa, xv->n, zv->data.F32, "y");

    graphdata->color = KapaColorByName ("blue");
    graphdata->ptype = 0;
    graphdata->size = 1.5;
    graphdata->style = 2;
    KapaPrepPlot (myKapa, xv->n, graphdata);
    KapaPlotVector (myKapa, xv->n, xv->data.F32, "x");
    KapaPlotVector (myKapa, xv->n, zv->data.F32, "y");

    xv->n = 6;

    // set the three axis lines
    xv->data.F32[0] = +0.0*cos(theta) + 0.0*sin(theta)*cos(phi) + 0.0*sin(theta)*sin(phi);
    yv->data.F32[0] = -0.0*sin(theta) + 0.0*cos(theta)*cos(phi) + 0.0*cos(theta)*sin(phi);
    zv->data.F32[0] =                 - 0.0*sin(phi)            + 0.0*cos(phi);
    xv->data.F32[1] = +1.0*cos(theta) + 0.0*sin(theta)*cos(phi);
    yv->data.F32[1] = -1.0*sin(theta) + 0.0*cos(theta)*cos(phi);
    zv->data.F32[1] =                 - 0.0*sin(phi)            + 0.0*cos(phi);

    xv->data.F32[2] = +0.0*cos(theta) + 0.0*sin(theta)*cos(phi) + 0.0*sin(theta)*sin(phi);
    yv->data.F32[2] = -0.0*sin(theta) + 0.0*cos(theta)*cos(phi) + 0.0*cos(theta)*sin(phi);
    zv->data.F32[2] =                 - 0.0*sin(phi)            + 0.0*cos(phi);
    xv->data.F32[3] = +0.0*cos(theta) + 1.0*sin(theta)*cos(phi);
    yv->data.F32[3] = -0.0*sin(theta) + 1.0*cos(theta)*cos(phi);
    zv->data.F32[3] =                 - 1.0*sin(phi)            + 0.0*cos(phi);

    xv->data.F32[4] = +0.0*cos(theta) + 0.0*sin(theta)*cos(phi) + 1.0*sin(theta)*sin(phi);
    yv->data.F32[4] = -0.0*sin(theta) + 0.0*cos(theta)*cos(phi) + 1.0*cos(theta)*sin(phi);
    zv->data.F32[4] =                 - 0.0*sin(phi)            + 1.0*cos(phi);
    xv->data.F32[5] = +0.0*cos(theta) + 0.0*sin(theta)*cos(phi);
    yv->data.F32[5] = -0.0*sin(theta) + 0.0*cos(theta)*cos(phi);
    zv->data.F32[5] =                 - 0.0*sin(phi)            + 0.0*cos(phi);

    graphdata->color = KapaColorByName ("red");
    graphdata->ptype = 100;
    graphdata->size = 0.5;
    graphdata->style = 2;
    KapaPrepPlot (myKapa, xv->n, graphdata);
    KapaPlotVector (myKapa, xv->n, xv->data.F32, "x");
    KapaPlotVector (myKapa, xv->n, zv->data.F32, "y");

    psFree (xv);
    psFree (yv);
    psFree (zv);

    return true;
}

# else
# endif
