Index: /branches/eam_branches/ipp-20120805/psphot/src/psphot.h
===================================================================
--- /branches/eam_branches/ipp-20120805/psphot/src/psphot.h	(revision 34359)
+++ /branches/eam_branches/ipp-20120805/psphot/src/psphot.h	(revision 34360)
@@ -282,5 +282,5 @@
 bool            psphotVisualShowPSFModel (pmReadout *readout, pmPSF *psf);
 bool            psphotVisualPlotRadialProfile (int myKapa, pmSource *source, psImageMaskType maskVal);
-bool            psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources);
+bool            psphotVisualPlotRadialProfiles (psMetadata *recipe, psArray *sources, pmSourceMode showmode);
 bool            psphotVisualShowFlags (psArray *sources);
 bool            psphotVisualShowSourceSize (pmReadout *readout, psArray *sources);
Index: /branches/eam_branches/ipp-20120805/psphot/src/psphotDeblendSatstars.c
===================================================================
--- /branches/eam_branches/ipp-20120805/psphot/src/psphotDeblendSatstars.c	(revision 34359)
+++ /branches/eam_branches/ipp-20120805/psphot/src/psphotDeblendSatstars.c	(revision 34360)
@@ -78,4 +78,6 @@
     psVector *index = psVectorSortIndex (NULL, SN);
     // this results in an index of increasing SN
+
+    psphotVisualPlotRadialProfiles (recipe, sources, PM_SOURCE_MODE_SATSTAR);
 
     // examine sources in decreasing SN order
@@ -219,2 +221,206 @@
     return true;
 }
+
+bool psphotDeblendSatstarsReadoutOld (pmConfig *config, const pmFPAview *view, const char *filerule, int fileIndex) {
+
+    int N;
+    pmSource *source;
+    bool status;
+
+    psTimerStart ("psphot.deblend.sat");
+
+    int Nblend = 0;
+    float SAT_MIN_RADIUS = 5.0;
+
+    // find the currently selected readout
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, fileIndex); // File of interest
+    psAssert (file, "missing file?");
+
+    pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+    psAssert (readout, "missing readout?");
+
+    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+    psAssert (detections, "missing detections?");
+
+    psArray *sources = detections->newSources;
+    psAssert (sources, "missing sources?");
+
+    if (!sources->n) {
+	psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping satstar blend");
+	return true;
+    }
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    pmCell *cell = readout->parent;
+
+    float SATURATION = NAN;
+
+    // do not completely trust the values in the header...
+    float CELL_SATURATION = psMetadataLookupF32 (&status, cell->concepts, "CELL.SATURATION");
+    float MIN_SATURATION = psMetadataLookupF32 (&status, recipe, "DEBLEND_MIN_SATURATION");
+    if (!status || !isfinite(MIN_SATURATION)) {
+	MIN_SATURATION = 40000.0;
+    }
+    if (!isfinite(CELL_SATURATION)) {
+	SATURATION = MIN_SATURATION;
+    } else {
+	SATURATION = PS_MAX(MIN_SATURATION, CELL_SATURATION);
+    }
+    float SAT_TEST_LEVEL = 0.5*SATURATION;
+
+    // we need sources spatially-sorted to find overlaps
+    sources = psArraySort (sources, pmSourceSortByY);
+
+    // source analysis is done in peak order (brightest first)
+    // we use an index for this so the spatial sorting is kept
+    psVector *SN = psVectorAlloc (sources->n, PS_DATA_F32);
+    for (int i = 0; i < SN->n; i++) {
+        source = sources->data[i];
+        SN->data.F32[i] = source->peak->rawFlux;
+    }
+    psVector *index = psVectorSortIndex (NULL, SN);
+    // this results in an index of increasing SN
+
+    // examine sources in decreasing SN order
+    for (int i = sources->n - 1; i >= 0; i--) {
+        N = index->data.U32[i];
+        source = sources->data[N];
+
+        // XXX filter? if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
+        if (source->mode & PM_SOURCE_MODE_BLEND) continue;
+        if (source->peak->rawFlux < SAT_TEST_LEVEL) continue;
+
+	// save these for reference below
+	int xPeak = source->peak->x;
+	int yPeak = source->peak->y;
+
+        // generate a basic contour (set of x,y coordinates at-or-below flux level)
+        psArray *contour = pmSourceContour (source->pixels, xPeak, yPeak, SAT_TEST_LEVEL);
+        if (contour == NULL) continue;
+
+	// contour consists of a set of X,Y coords giving the boundary
+	psVector *xVec = contour->data[0];
+	psVector *yVec = contour->data[1];
+	if (xVec->n < 5) {
+	    psFree(contour);
+	    continue;
+	}
+
+	// find the center of the contour (let's just use mid[x,y])
+	int xMin = xVec->data.F32[0];
+	int xMax = xVec->data.F32[0];
+	int yMin = yVec->data.F32[0];
+	int yMax = yVec->data.F32[0];
+	for (int j = 0; j < xVec->n; j++) {
+	    xMin = PS_MIN (xMin, xVec->data.F32[j]);
+	    xMax = PS_MAX (xMax, xVec->data.F32[j]);
+	    yMin = PS_MIN (yMin, yVec->data.F32[j]);
+	    yMax = PS_MAX (yMax, yVec->data.F32[j]);
+	}	
+	int xCenter = 0.5*(xMin + xMax);
+	int yCenter = 0.5*(yMin + yMax);
+	psFree (contour);
+
+	psAssert (xCenter >= source->pixels->col0, "invalid shift in object center");
+	psAssert (xCenter <  source->pixels->col0 + source->pixels->numCols, "invalid shift in object center");
+	psAssert (yCenter >= source->pixels->row0, "invalid shift in object center");
+	psAssert (yCenter <  source->pixels->row0 + source->pixels->numRows, "invalid shift in object center");
+
+	// reset the peak for this source to the value of the center pixel
+	source->peak->x = xCenter;
+	source->peak->y = yCenter;
+	source->peak->xf = xCenter;
+	source->peak->yf = yCenter;
+	
+	// temporary array for overlapping objects we find
+        psArray *overlap = psArrayAllocEmpty (100);
+
+        // search backwards for overlapping sources
+        for (int j = N - 1; j >= 0; j--) {
+            pmSource *testSource = sources->data[j];
+            if (testSource->peak->x <  source->pixels->col0) continue;
+            if (testSource->peak->x >= source->pixels->col0 + source->pixels->numCols) continue;
+            if (testSource->peak->y <  source->pixels->row0) break;
+            if (testSource->peak->y >= source->pixels->row0 + source->pixels->numRows) {
+                fprintf (stderr, "warning: invalid condition\n");
+                continue;
+            }
+            psArrayAdd (overlap, 100, testSource);
+        }
+
+        // search forwards for overlapping sources
+        for (int j = N + 1; j < sources->n; j++) {
+            pmSource *testSource = sources->data[j];
+            if (testSource->peak->x <  source->pixels->col0) continue;
+            if (testSource->peak->x >= source->pixels->col0 + source->pixels->numCols) continue;
+            if (testSource->peak->y <  source->pixels->row0) {
+                fprintf (stderr, "warning: invalid condition\n");
+                continue;
+            }
+            if (testSource->peak->y >= source->pixels->row0 + source->pixels->numRows) break;
+            psArrayAdd (overlap, 100, testSource);
+        }
+
+        if (overlap->n == 0) {
+            psFree (overlap);
+            continue;
+        }
+
+	// now find the contour which is at 0.5*SAT_TEST_LEVEL (xPeak, yPeak is valid high pixel)
+        contour = pmSourceContour (source->pixels, xPeak, yPeak, 0.5*SAT_TEST_LEVEL);
+        if (contour == NULL) {
+	    psFree (overlap);
+	    continue; 
+	}
+	
+	// any peaks within this contour should be dropped (mark as blend, drop after loop is done)
+	// also drop any peaks which are too close to this peak
+        psVector *xv = contour->data[0];
+        psVector *yv = contour->data[1];
+        for (int k = 0; k < overlap->n; k++) {
+            pmSource *testSource = overlap->data[k];
+	    float radius = hypot((testSource->peak->x - xCenter), (testSource->peak->y - yCenter));
+	    if (radius < SAT_MIN_RADIUS) {
+                testSource->mode |= PM_SOURCE_MODE_BLEND;
+                Nblend ++;
+		continue;
+	    }
+            for (int j = 0; j < xv->n; j+=2) {
+                if (fabs(yv->data.F32[j] - testSource->peak->y) > 0.5) continue;
+                if (xv->data.F32[j+0] > testSource->peak->x) break;
+                if (xv->data.F32[j+1] < testSource->peak->x) break;
+                testSource->mode |= PM_SOURCE_MODE_BLEND;
+                Nblend ++;
+                j = xv->n; // skip rest of contour
+            }
+        }
+        psFree (overlap);
+        psFree (contour);
+    }
+
+    // drop the sources marked as BLEND
+    for (int i = 0; i < sources->n;) {
+	pmSource *source = sources->data[i];
+
+        if (!(source->mode & PM_SOURCE_MODE_BLEND)) {
+	    i++;
+	    continue;
+	}
+
+	// shuffle the remaining sources forward
+	for (int j = i; j < sources->n - 1; j++) {
+	    sources->data[j] = sources->data[j+1];
+	}
+	psFree (source);
+	sources->n --;
+    }
+
+    psFree (SN);
+    psFree (index);
+
+    psLogMsg ("psphot", PS_LOG_INFO, "found %d satstar blend peaks, leaving %ld sources: %f sec\n", Nblend, sources->n, psTimerMark ("psphot.deblend.sat"));
+    return true;
+}
Index: /branches/eam_branches/ipp-20120805/psphot/src/psphotReadout.c
===================================================================
--- /branches/eam_branches/ipp-20120805/psphot/src/psphotReadout.c	(revision 34359)
+++ /branches/eam_branches/ipp-20120805/psphot/src/psphotReadout.c	(revision 34360)
@@ -129,8 +129,8 @@
 
     // find blended neighbors of very saturated stars (detections->newSources)
-    // if (!psphotDeblendSatstars (config, view, filerule)) {
-    // psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
-    // return psphotReadoutCleanup (config, view, filerule);
-    // }
+    if (!psphotDeblendSatstars (config, view, filerule)) {
+	psError (PSPHOT_ERR_UNKNOWN, false, "failed on satstar deblend analysis");
+	return psphotReadoutCleanup (config, view, filerule);
+    }
 
     // mark blended peaks PS_SOURCE_BLEND (detections->newSources)
Index: /branches/eam_branches/ipp-20120805/psphot/src/psphotVisual.c
===================================================================
--- /branches/eam_branches/ipp-20120805/psphot/src/psphotVisual.c	(revision 34359)
+++ /branches/eam_branches/ipp-20120805/psphot/src/psphotVisual.c	(revision 34360)
@@ -1139,5 +1139,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;
@@ -1302,6 +1302,16 @@
     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;
+	Yo = source->moments->My;
+    }
+
     for (int iy = 0; iy < source->pixels->numRows; iy++) {
 	for (int ix = 0; ix < source->pixels->numCols; ix++) {
@@ -1322,67 +1332,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);
 
@@ -1417,33 +1364,4 @@
     KapaPlotVector (myKapa, nb, rb->data.F32, "x");
     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 **
@@ -1479,23 +1397,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(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);
+
+	// ** 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 +1526,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,5 +1567,5 @@
 
 	pmSource *source = sources->data[i];
-	if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
+	if (!(source->mode & showmode)) continue;
 
 	psphotVisualPlotRadialProfile (myKapa, source, maskVal);