Index: trunk/psphot/src/psphotVisual.c
===================================================================
--- trunk/psphot/src/psphotVisual.c	(revision 31452)
+++ trunk/psphot/src/psphotVisual.c	(revision 34404)
@@ -262,4 +262,7 @@
 }
 
+static psImage *posImage = NULL;
+static psImage *delImage = NULL;
+
 bool psphotVisualShowImage (pmReadout *readout) {
 
@@ -277,4 +280,8 @@
     psphotVisualScaleImage (kapa, readout->variance, readout->mask, "variance", 1.0, 1);
     psphotVisualScaleImage (kapa, readout->image, readout->mask, "image", sqrt(factor), 0);
+
+    if (posImage == NULL) {
+	posImage = psImageCopy (NULL, readout->image, PS_TYPE_F32);
+    }
 
     pmVisualAskUser(NULL);
@@ -1139,5 +1146,5 @@
     // after displaying (as an image) the psf stars, we cycle throught them and display their
     // radial profiles:
-    psphotVisualPlotRadialProfiles (recipe, sources);
+    psphotVisualPlotRadialProfiles (recipe, sources, PM_SOURCE_MODE_PSFSTAR);
 
     return true;
@@ -1273,5 +1280,5 @@
 }
 
-static void plotline (int myKapa, Graphdata *graphdata, float x0, float y0, float x1, float y1) 
+void plotline (int myKapa, Graphdata *graphdata, float x0, float y0, float x1, float y1) 
 {
     float x[2], y[2];
@@ -1285,7 +1292,9 @@
 }
 
-bool psphotVisualPlotRadialProfile (int myKapa, pmSource *source, psImageMaskType maskVal) {
+bool psphotVisualPlotRadialProfile (int myKapa, pmSource *source, psImageMaskType maskVal, pmSourceMode showmode) {
 
     Graphdata graphdata;
+
+    float Rmax = (showmode & PM_SOURCE_MODE_SATSTAR) ? 100.0 : 30.0;
 
     bool subtracted = (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED);
@@ -1302,9 +1311,19 @@
     int ng = 0;
     int nb = 0;
-    float Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS] - source->pixels->col0;
-    float Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS] - source->pixels->row0;
+
+    float Xo = NAN;
+    float Yo = NAN;
+
+    if (source->modelPSF) {
+	Xo = source->modelPSF->params->data.F32[PM_PAR_XPOS] - source->pixels->col0;
+	Yo = source->modelPSF->params->data.F32[PM_PAR_YPOS] - source->pixels->row0;
+    } else {
+	Xo = source->moments->Mx - source->pixels->col0;
+	Yo = source->moments->My - source->pixels->row0;
+    }
+
     for (int iy = 0; iy < source->pixels->numRows; iy++) {
 	for (int ix = 0; ix < source->pixels->numCols; ix++) {
-	    if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix]) {
+	    if ((source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskVal)) {
 		rb->data.F32[nb] = hypot (ix + 0.5 - Xo, iy + 0.5 - Yo) ;
 		// rb->data.F32[nb] = hypot (ix - Xo, iy - Yo) ;
@@ -1322,67 +1341,4 @@
     }
 
-    // generate model profiles (major and minor axis):
-    // create a model with theta = 0.0 so major and minor axes are equiv to x and y:
-    psEllipseShape rawShape, rotShape;
-
-    rawShape.sx  = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2;
-    rawShape.sy  = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2;
-    rawShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY];
-
-    psEllipseAxes axes = psEllipseShapeToAxes (rawShape, 20.0);
-
-    axes.theta = 0.0;
-
-    rotShape = psEllipseAxesToShape (axes);
-
-    psVector *params = psVectorAlloc(source->modelPSF->params->n, PS_TYPE_F32);
-    for (int i = 0; i < source->modelPSF->params->n; i++) {
-	params->data.F32[i] = source->modelPSF->params->data.F32[i];
-    }
-    params->data.F32[PM_PAR_SXX] = rotShape.sx * M_SQRT2;
-    params->data.F32[PM_PAR_SYY] = rotShape.sy * M_SQRT2;
-    params->data.F32[PM_PAR_SXY] = rotShape.sxy;
-    params->data.F32[PM_PAR_XPOS] = 0.0;
-    params->data.F32[PM_PAR_YPOS] = 0.0;
-
-    psVector *rmod = psVectorAlloc(300, PS_TYPE_F32);
-    psVector *fmaj = psVectorAlloc(300, PS_TYPE_F32);
-    psVector *fmin = psVectorAlloc(300, PS_TYPE_F32);
-
-    psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
-
-    float r = 0.0;
-    for (int i = 0; i < rmod->n; i++) {
-	r = i*0.1;
-	rmod->data.F32[i] = r;
-
-	coord->data.F32[1] = r;
-	coord->data.F32[0] = 0.0;
-	fmaj->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord));
-
-	coord->data.F32[0] = r;
-	coord->data.F32[1] = 0.0;
-	fmin->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord));
-    }
-    psFree (coord);
-    psFree (params);
-
-    float FWHM_MAJOR = 2.0*source->modelPSF->modelRadius (source->modelPSF->params, 0.5*source->modelPSF->params->data.F32[PM_PAR_I0]);
-    float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
-    if (FWHM_MAJOR < FWHM_MINOR) PS_SWAP (FWHM_MAJOR, FWHM_MINOR); 
-
-    psEllipseMoments emoments;
-    emoments.x2 = source->moments->Mxx;
-    emoments.xy = source->moments->Mxy;
-    emoments.y2 = source->moments->Myy;
-    axes = psEllipseMomentsToAxes (emoments, 20.0);
-    float MOMENTS_MAJOR = 2.355*axes.major;
-    float MOMENTS_MINOR = 2.355*axes.minor;
-
-    float logHM = log10(0.5*source->modelPSF->params->data.F32[PM_PAR_I0]);
-
-    // reset source Add/Sub state to recorded
-    if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-
     KapaInitGraph (&graphdata);
 
@@ -1392,5 +1348,5 @@
     // examine sources to set data range
     graphdata.xmin =  -0.05;
-    graphdata.xmax = +30.05;
+    graphdata.xmax = Rmax + 0.05;
     graphdata.ymin = -0.05;
     graphdata.ymax = +8.05;
@@ -1418,33 +1374,4 @@
     KapaPlotVector (myKapa, nb, fb->data.F32, "y");
 
-    graphdata.color = KapaColorByName ("blue");
-    graphdata.ptype = 0;
-    graphdata.size = 0.0;
-    graphdata.style = 0;
-    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
-    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
-    KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y");
-    plotline (myKapa, &graphdata, 0.0, logHM, 30.0, logHM);
-    plotline (myKapa, &graphdata, 0.5*FWHM_MINOR, 0.0, 0.5*FWHM_MINOR, 5.0);
-    graphdata.ltype = 1;
-    plotline (myKapa, &graphdata, 0.5*MOMENTS_MINOR, 0.0, 0.5*MOMENTS_MINOR, 5.0);
-    graphdata.ltype = 0;
-	
-    graphdata.color = KapaColorByName ("green");
-    graphdata.ptype = 0;
-    graphdata.size = 0.0;
-    graphdata.style = 0;
-    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
-    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
-    KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y");
-    plotline (myKapa, &graphdata, 0.5*FWHM_MAJOR, 0.0, 0.5*FWHM_MAJOR, 5.0);
-    graphdata.ltype = 1;
-    plotline (myKapa, &graphdata, 0.5*MOMENTS_MAJOR, 0.0, 0.5*MOMENTS_MAJOR, 5.0);
-    graphdata.ltype = 0;
-	
-    for (int i = 0; i < rmod->n; i++) {
-	rmod->data.F32[i] = log10(rmod->data.F32[i]);
-    }
-
     // ** loglog **
     KapaSelectSection (myKapa, "loglog");
@@ -1452,5 +1379,5 @@
     // examine sources to set data range
     graphdata.xmin = -1.51;
-    graphdata.xmax = +1.51;
+    graphdata.xmax = log10(Rmax) + 0.02;
     graphdata.ymin = -0.05;
     graphdata.ymax = +8.05;
@@ -1479,23 +1406,123 @@
     KapaPlotVector (myKapa, nb, fb->data.F32, "y");
 
-    graphdata.color = KapaColorByName ("blue");
-    graphdata.ptype = 0;
-    graphdata.size = 0.0;
-    graphdata.style = 0;
-    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
-    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
-    KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y");
-
-    graphdata.color = KapaColorByName ("green");
-    graphdata.ptype = 0;
-    graphdata.size = 0.0;
-    graphdata.style = 0;
-    KapaPrepPlot   (myKapa, rmod->n, &graphdata);
-    KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
-    KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y");
-
-    psFree (rmod);
-    psFree (fmin);
-    psFree (fmaj);
+    if (source->modelPSF) {
+	// generate model profiles (major and minor axis):
+	// create a model with theta = 0.0 so major and minor axes are equiv to x and y:
+	psEllipseShape rawShape, rotShape;
+
+	rawShape.sx  = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2;
+	rawShape.sy  = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2;
+	rawShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY];
+
+	psEllipseAxes axes = psEllipseShapeToAxes (rawShape, 20.0);
+
+	axes.theta = 0.0;
+
+	rotShape = psEllipseAxesToShape (axes);
+
+	psVector *params = psVectorAlloc(source->modelPSF->params->n, PS_TYPE_F32);
+	for (int i = 0; i < source->modelPSF->params->n; i++) {
+	    params->data.F32[i] = source->modelPSF->params->data.F32[i];
+	}
+	params->data.F32[PM_PAR_SXX] = rotShape.sx * M_SQRT2;
+	params->data.F32[PM_PAR_SYY] = rotShape.sy * M_SQRT2;
+	params->data.F32[PM_PAR_SXY] = rotShape.sxy;
+	params->data.F32[PM_PAR_XPOS] = 0.0;
+	params->data.F32[PM_PAR_YPOS] = 0.0;
+
+	psVector *rmod = psVectorAlloc(Rmax*10, PS_TYPE_F32);
+	psVector *fmaj = psVectorAlloc(Rmax*10, PS_TYPE_F32);
+	psVector *fmin = psVectorAlloc(Rmax*10, PS_TYPE_F32);
+
+	psVector *coord = psVectorAlloc(2, PS_TYPE_F32);
+
+	float r = 0.0;
+	for (int i = 0; i < rmod->n; i++) {
+	    r = i*0.1;
+	    rmod->data.F32[i] = r;
+
+	    coord->data.F32[1] = r;
+	    coord->data.F32[0] = 0.0;
+	    fmaj->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord));
+
+	    coord->data.F32[0] = r;
+	    coord->data.F32[1] = 0.0;
+	    fmin->data.F32[i] = log10(source->modelPSF->modelFunc (NULL, params, coord));
+	}
+	psFree (coord);
+	psFree (params);
+
+	float FWHM_MAJOR = 2.0*source->modelPSF->modelRadius (source->modelPSF->params, 0.5*source->modelPSF->params->data.F32[PM_PAR_I0]);
+	float FWHM_MINOR = FWHM_MAJOR * (axes.minor / axes.major);
+	if (FWHM_MAJOR < FWHM_MINOR) PS_SWAP (FWHM_MAJOR, FWHM_MINOR); 
+
+	psEllipseMoments emoments;
+	emoments.x2 = source->moments->Mxx;
+	emoments.xy = source->moments->Mxy;
+	emoments.y2 = source->moments->Myy;
+	axes = psEllipseMomentsToAxes (emoments, 20.0);
+	float MOMENTS_MAJOR = 2.355*axes.major;
+	float MOMENTS_MINOR = 2.355*axes.minor;
+
+	float logHM = log10(0.5*source->modelPSF->params->data.F32[PM_PAR_I0]);
+
+	// reset source Add/Sub state to recorded
+	if (subtracted) pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+
+	// ** linlog **
+	KapaSelectSection (myKapa, "linlog");
+	KapaGetGraphData (myKapa, &graphdata);
+	graphdata.color = KapaColorByName ("blue");
+	graphdata.ptype = 0;
+	graphdata.size = 0.0;
+	graphdata.style = 0;
+	KapaPrepPlot   (myKapa, rmod->n, &graphdata);
+	KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
+	KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y");
+	plotline (myKapa, &graphdata, graphdata.xmin, logHM, graphdata.xmax, logHM);
+	plotline (myKapa, &graphdata, 0.5*FWHM_MINOR, graphdata.ymin, 0.5*FWHM_MINOR, graphdata.ymax);
+	graphdata.ltype = 1;
+	plotline (myKapa, &graphdata, 0.5*MOMENTS_MINOR, graphdata.ymin, 0.5*MOMENTS_MINOR, graphdata.ymax);
+	graphdata.ltype = 0;
+	
+	graphdata.color = KapaColorByName ("green");
+	graphdata.ptype = 0;
+	graphdata.size = 0.0;
+	graphdata.style = 0;
+	KapaPrepPlot   (myKapa, rmod->n, &graphdata);
+	KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
+	KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y");
+	plotline (myKapa, &graphdata, 0.5*FWHM_MAJOR, graphdata.ymin, 0.5*FWHM_MAJOR, graphdata.ymax);
+	graphdata.ltype = 1;
+	plotline (myKapa, &graphdata, 0.5*MOMENTS_MAJOR, graphdata.ymin, 0.5*MOMENTS_MAJOR, graphdata.ymax);
+	graphdata.ltype = 0;
+	
+	for (int i = 0; i < rmod->n; i++) {
+	    rmod->data.F32[i] = log10(rmod->data.F32[i]);
+	}
+
+	// ** loglog **
+	KapaSelectSection (myKapa, "loglog");
+	KapaGetGraphData (myKapa, &graphdata);
+	graphdata.color = KapaColorByName ("blue");
+	graphdata.ptype = 0;
+	graphdata.size = 0.0;
+	graphdata.style = 0;
+	KapaPrepPlot   (myKapa, rmod->n, &graphdata);
+	KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
+	KapaPlotVector (myKapa, rmod->n, fmin->data.F32, "y");
+
+	graphdata.color = KapaColorByName ("green");
+	graphdata.ptype = 0;
+	graphdata.size = 0.0;
+	graphdata.style = 0;
+	KapaPrepPlot   (myKapa, rmod->n, &graphdata);
+	KapaPlotVector (myKapa, rmod->n, rmod->data.F32, "x");
+	KapaPlotVector (myKapa, rmod->n, fmaj->data.F32, "y");
+
+	psFree (rmod);
+	psFree (fmin);
+	psFree (fmaj);
+    }
 
     psFree (rg);
@@ -1508,5 +1535,5 @@
 }
 
-bool psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources) {
+bool psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources, pmSourceMode showmode) {
 
     KapaSection section;  // put the positive profile in one and the residuals in another?
@@ -1549,7 +1576,27 @@
 
 	pmSource *source = sources->data[i];
-	if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
-
-	psphotVisualPlotRadialProfile (myKapa, source, maskVal);
+	if (!(source->mode & showmode)) continue;
+
+	psphotVisualPlotRadialProfile (myKapa, source, maskVal, showmode);
+
+	if (pmVisualTestLevel("psphot.image", 1)) {
+	  int display = psphotKapaChannel (1);
+	  KiiCenter (display, source->peak->xf, source->peak->yf, 1);
+	  KiiOverlay overlay;
+	  overlay.x = source->peak->xf;
+	  overlay.y = source->peak->yf;
+	  overlay.dx = 5;
+	  overlay.dy = 5;
+	  overlay.angle = 0.0;
+	  overlay.type = KiiOverlayTypeByName ("circle");
+	  KiiLoadOverlay (display, &overlay, 1, "red");
+	  overlay.x = source->moments->Mx;
+	  overlay.y = source->moments->My;
+	  overlay.dx = 8;
+	  overlay.dy = 8;
+	  overlay.angle = 0.0;
+	  overlay.type = KiiOverlayTypeByName ("circle");
+	  KiiLoadOverlay (display, &overlay, 1, "blue");
+	}
 
 	// pause and wait for user input:
@@ -2540,9 +2587,16 @@
     }
 
-    if (reshow) {
+    if (false && reshow) {
 	psphotVisualShowMask (myKapa, readout->mask, "mask", 2);
 	psphotVisualScaleImage (myKapa, readout->variance, readout->mask, "variance", 1.0, 1);
     }
-    psphotVisualScaleImage (myKapa, readout->image, readout->mask, "resid", sqrt(factor), 0);
+
+    if (posImage) {
+	delImage = (psImage *) psBinaryOp(delImage, posImage, "-", readout->image);
+	psphotVisualScaleImage (myKapa, posImage, readout->mask, "posimage", sqrt(factor), 0);
+	psphotVisualScaleImage (myKapa, delImage, readout->mask, "delimage", sqrt(factor), 2);
+    }
+
+    psphotVisualScaleImage (myKapa, readout->image, readout->mask, "resid", sqrt(factor), 1);
 
     pmVisualAskUser(NULL);