Index: /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/relastroVisual.c
===================================================================
--- /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/relastroVisual.c	(revision 22725)
+++ /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/relastroVisual.c	(revision 22725)
@@ -0,0 +1,248 @@
+/** Diagnostic plots for relastro
+ * @author Chris Beaumont, IfA
+ */
+
+#include "kapa.h"
+#include "relastro.h"
+
+#define KAPAX 700
+#define KAPAY 700
+
+static int kapa = -1;
+static int kapa2 = -1;
+
+static int isVisual = 1;
+static int plotRawRef = 1;
+
+/** Spawn a kapa window */
+static int initWindow(int *kapid) {
+  if (*kapid == -1) {
+    *kapid = KapaOpenNamedSocket("kapa", "relastro_plots");
+    if (*kapid == -1) {
+      fprintf(stderr, "Failure to open kapa.\n");
+      isVisual = 0;
+      return 0;
+    }
+    KapaResize (*kapid, KAPAX, KAPAY);
+  }
+  return 1;
+}
+
+/** Ask the user how to proceed */
+static int askUser(int *flag) {
+  char key[10];
+  fprintf (stdout, "[c]ontinue? [s]kip the rest of these plots? [a]bort all visual plots?");
+  if (!fgets(key, 8, stdin)) {
+    fprintf(stderr, "warning: Unable to read option");
+  }
+  if (key[0] == 's') {
+    *flag = 0;
+  }
+  if (key[0] == 'a') {
+    isVisual = 0;
+  }
+  return 1;
+}
+
+/** Size graphdata to encompass all points */
+static void scaleGraphdata(float x[], float y[], Graphdata *graphdata, int n) {
+  float xlo = 9999999, xhi = -9999999, ylo = 9999999, yhi = -9999999;
+  int i;
+  for(i = 0; i < n; i++) {
+    if(x[i] < xlo) xlo = x[i];
+    if(x[i] > xhi) xhi = x[i];
+    if(y[i] < ylo) ylo = y[i];
+    if(y[i] > yhi) yhi = y[i];
+  }
+  graphdata->xmin = xlo;
+  graphdata->ymin = ylo;
+  graphdata->xmax = xhi;
+  graphdata->ymax = yhi;
+}
+
+/**Plot a vector field*/
+static int plotVectorField(float x[], float y[],
+                           float xVec[], float yVec[],
+                           int npts, int kapaID,
+                           Graphdata *graphdata, int overplot) {
+
+    float singleX[2], singleY[2];
+    float maxVecLength, vecScaleFactor;
+    float graphSize;
+    int i;
+
+    if (!initWindow(&kapaID)) return 0;
+
+    if (!overplot) scaleGraphdata(x, y, graphdata, npts);
+
+    //scale the raw xvec, yvecs down to fit into the window
+    maxVecLength = 0;
+    for(i = 0; i < npts; i++) {
+        if (fabs(xVec[i]) > maxVecLength) {
+            maxVecLength = fabs(xVec[i]);
+        }
+        if (fabs(yVec[i]) > maxVecLength) {
+            maxVecLength = fabs(yVec[i]);
+        }
+    }
+    graphSize = graphdata->xmax - graphdata->xmin;
+    if ((graphdata->ymax - graphdata->ymin) > graphSize) {
+        graphSize = graphdata->ymax - graphdata->ymin;
+    }
+    vecScaleFactor = graphSize * 0.1 / maxVecLength;
+
+    //plot the points as circles
+    if (!overplot) {
+        KapaSetFont (kapaID, "helvetica", 14);
+        KapaSetLimits(kapaID, graphdata);
+        KapaBox(kapaID, graphdata);
+    }
+
+    graphdata->ptype = 7;
+    graphdata->style = 2;
+    KapaPrepPlot(kapaID, npts, graphdata);
+    KapaPlotVector(kapaID, npts, x, "x");
+    KapaPlotVector(kapaID, npts, y, "y");
+
+    //plot each vector individually
+    graphdata->ptype = 0;
+    graphdata->style = 0;
+    for(i = 0; i < npts; i++) {
+        singleX[0] = x[i];
+        singleY[0] = y[i];
+        singleX[1] = xVec[i] * vecScaleFactor;
+        singleY[1] = yVec[i] * vecScaleFactor;
+
+        KapaPrepPlot(kapaID, 2, graphdata);
+        KapaPlotVector(kapaID, 2, singleX, "x");
+        KapaPlotVector(kapaID, 2, singleY, "y");
+    }
+    return 1;
+}
+
+
+/** plot raw vs ref (L, M). Only those whose distance is < drMax are used in fit*/
+int relastroVisualPlotRawRef(StarData *raw, StarData *ref, double dRmax, int numObj) {
+
+  if( !isVisual || !plotRawRef) return 1;
+  if( !initWindow(&kapa)) return 0;
+
+  float *rawX, *rawY,  *refX, *refY;
+  float *rawXfit, *rawYfit, *refXfit, *refYfit;
+  float *magRaw, *magRef, *magRawfit, *magReffit;
+  float *xVec, *yVec;
+  int numFit = 0, numNoFit = 0;
+  double dL, dM, dR;
+
+  ALLOCATE(rawX,      float, numObj);
+  ALLOCATE(rawY,      float, numObj);
+  ALLOCATE(refX,      float, numObj);
+  ALLOCATE(refY,      float, numObj);
+  ALLOCATE(rawXfit,   float, numObj);
+  ALLOCATE(rawYfit,   float, numObj);
+  ALLOCATE(refXfit,   float, numObj);
+  ALLOCATE(refYfit,   float, numObj);
+  ALLOCATE(magRaw,    float, numObj);
+  ALLOCATE(magRef,    float, numObj);
+  ALLOCATE(magRawfit, float, numObj);
+  ALLOCATE(magReffit, float, numObj);
+
+  int i;
+  for(i = 0; i < numObj; i++) {
+    if  (raw[i].mask) continue;
+
+    dL = raw[i].L - ref[i].L;
+    dM = raw[i].M - ref[i].M;
+    dR = hypot (dL, dM);
+    if (dR > dRmax) {
+      rawX[numNoFit] = raw[i].X;
+      rawY[numNoFit] = raw[i].Y;
+      refX[numNoFit] = ref[i].X;
+      refY[numNoFit] = ref[i].Y;
+      magRaw[numNoFit] = raw[i].Mag;
+      magRef[numNoFit] = ref[i].Mag;
+      numNoFit++;
+    } else {
+      rawXfit[numFit] = raw[i].X;
+      rawYfit[numFit] = raw[i].Y;
+      refXfit[numFit] = ref[i].X;
+      refYfit[numFit] = ref[i].Y;
+      magRaw[numFit] = raw[i].Mag;
+      magRef[numFit] = ref[i].Mag;
+      numFit++;
+    }
+  }
+
+  ALLOCATE(xVec, float, numFit);
+  ALLOCATE(yVec, float, numFit);
+
+  Graphdata graphdata;
+
+  KapaInitGraph(&graphdata);
+  KapaClearPlots(kapa);
+
+  graphdata.ptype = 7;
+  graphdata.style = 2;
+
+  scaleGraphdata(rawXfit, rawYfit, &graphdata, numFit);
+  //  scaleGraphdata(refXfit, refYfit, &graphdata, numFit);
+
+  //plot the fitted objects as vectors
+  for(i = 0; i < numFit; i++) {
+      xVec[i] = rawXfit[i] - refXfit[i];
+      yVec[i] = rawYfit[i] - refYfit[i];
+  }
+
+  plotVectorField(refX, refY, xVec, yVec, numFit, kapa2,
+                  &graphdata, 0);
+
+  KapaSetFont(kapa, "helvetica", 14);
+  KapaSetLimits(kapa, &graphdata);
+  KapaBox(kapa, &graphdata);
+  KapaSendLabel( kapa, "X", KAPA_LABEL_XM);
+  KapaSendLabel( kapa, "Y", KAPA_LABEL_YM);
+  KapaSendLabel( kapa, "orange, red, green, blue: (raw, ref), (nofit, fit)",
+                 KAPA_LABEL_XP);
+
+  graphdata.color = KapaColorByName("orange");
+  graphdata.size = 1;
+  KapaPrepPlot(kapa, numNoFit, &graphdata);
+  KapaPlotVector(kapa, numNoFit, rawX, "x");
+  KapaPlotVector(kapa, numNoFit, rawY, "y");
+
+  graphdata.color = KapaColorByName("red");
+  graphdata.size = 2;
+  KapaPrepPlot(kapa, numNoFit, &graphdata);
+  KapaPlotVector(kapa, numNoFit, refX, "x");
+  KapaPlotVector(kapa, numNoFit, refY, "y");
+
+  graphdata.color = KapaColorByName("green");
+  graphdata.size = 1;
+  KapaPrepPlot(kapa, numFit,  &graphdata);
+  KapaPlotVector(kapa, numFit, rawXfit, "x");
+  KapaPlotVector(kapa, numFit, rawYfit, "y");
+
+  graphdata.color = KapaColorByName("blue");
+  graphdata.size = 2;
+  KapaPrepPlot(kapa, numFit, &graphdata);
+  KapaPlotVector(kapa, numFit, refXfit, "x");
+  KapaPlotVector(kapa, numFit, refYfit, "y");
+
+
+  askUser(&plotRawRef);
+
+  FREE(rawX);
+  FREE(rawY);
+  FREE(refX);
+  FREE(refY);
+  FREE(rawXfit);
+  FREE(rawYfit);
+  FREE(refXfit);
+  FREE(refYfit);
+  FREE(magRaw);
+  FREE(magRef);
+  FREE(magRawfit);
+  FREE(magReffit);
+
+  return 1;
+}
Index: /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/relastroVisual.h
===================================================================
--- /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/relastroVisual.h	(revision 22725)
+++ /branches/cnb_branches/cnb_branch_20090301/Ohana/src/relastro/src/relastroVisual.h	(revision 22725)
@@ -0,0 +1,6 @@
+#ifndef RELASTRO_VISUAL_H
+#define RELASTROPVISUAL_H
+
+int relastroVisualPlotRawRef(StarData *raw, StarData *ref, double dRmax, int numObj);
+
+#endif