Index: trunk/psphot/src/psphotModelTest.c
===================================================================
--- trunk/psphot/src/psphotModelTest.c	(revision 5802)
+++ trunk/psphot/src/psphotModelTest.c	(revision 5828)
@@ -1,13 +1,169 @@
 # include "psphot.h"
+# include "psEllipse.h"
 
-int main (int argc, char **argv) {
+bool psphotModelTest (eamReadout *imdata, psMetadata *config) {
 
-    psphotModelGroupInit ();
+    bool status;
+    int modelType;
+    float obsMag, fitMag, value;
+    char name[64];
+    pmPSF *psf;
 
-    psMetadata *config = modelTestArguments (&argc, argv);
-    eamReadout *imdata = psphotSetup (config);
+    psMetadataItem *item  = NULL;
 
-    modelTestFitSource (imdata, config);
+    // what fitting mode to use?
+    char *psfModeWord = psMetadataLookupPtr (&status, config, "TEST_FIT_MODE");
+    if (!status) {
+	psfModeWord = psStringCopy ("FLT");
+    }
+    bool psfMode = !strcasecmp (psfModeWord, "PSF");
 
-    exit (0);
+    // in psfMode, psf sets the model type
+    if (psfMode) {
+	char *psfFile = psMetadataLookupPtr (&status, config, "PSF_INPUT_FILE");
+	if (!status) psAbort ("psphotModelTest", "PSF_INPUT_FILE not supplied");
+	psf = psphotReadPSF (psfFile);
+	modelType = psf->type;
+    } else {
+	// find the model: supplied by user or first in the PSF_MODEL list
+	char *modelName  = psMetadataLookupPtr (&status, config, "TEST_FIT_MODEL");
+	if (modelName == NULL) {
+	    // get the list pointers for the PSF_MODEL entries
+	    psMetadataItem *mdi = psMetadataLookup (config, "PSF_MODEL");
+	    if (mdi == NULL) psAbort ("psphotChoosePSF", "missing PSF_MODEL selection");
+
+	    // take the first list element
+	    psList *list = (psList *) mdi->data.list;
+	    item = psListGet (list, PS_LIST_HEAD);
+	    modelName = item->data.V;
+	}
+	modelType = pmModelSetType (modelName);
+	if (modelType < 0) psAbort ("fitsource", "unknown model %s", modelName);
+    }
+
+    // find the fitting parameters (try test values first)
+    float INNER = psMetadataLookupF32 (&status, config, "TEST_FIT_INNER_RADIUS");
+    if (!status) {
+	INNER = psMetadataLookupF32 (&status, config, "SKY_INNER_RADIUS");
+    }
+
+    float OUTER = psMetadataLookupF32 (&status, config, "TEST_FIT_OUTER_RADIUS");
+    if (!status) {
+	OUTER = psMetadataLookupF32 (&status, config, "SKY_OUTER_RADIUS");
+    }
+
+    float RADIUS = psMetadataLookupF32 (&status, config, "TEST_FIT_RADIUS");
+    if (!status) {
+	RADIUS = psMetadataLookupF32 (&status, config, "PSF_FIT_RADIUS");
+    }
+
+    float mRADIUS = psMetadataLookupF32 (&status, config, "TEST_MOMENTS_RADIUS");
+    if (!status) {
+	mRADIUS = psMetadataLookupF32 (&status, config, "PSF_MOMENTS_RADIUS");
+    }
+
+    // define the source of interest
+    float xObj     = psMetadataLookupF32 (&status, config, "TEST_FIT_X");
+    float yObj     = psMetadataLookupF32 (&status, config, "TEST_FIT_Y");
+
+    // construct the source structures
+    pmSource *source = pmSourceAlloc();
+    source->peak = pmPeakAlloc (xObj, yObj, 0, 0);
+    psphotDefinePixels (source, imdata, xObj, yObj, OUTER);
+
+    // find the local sky
+    status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER);
+    if (!status) psAbort ("fitsource", "pmSourceLocalSky error");
+
+    // get the source moments
+    status = pmSourceMoments (source, mRADIUS);
+    if (!status) psAbort ("fitsource", "pmSourceLocalSky error");
+    source->peak->counts = source->moments->Peak;
+
+    fprintf (stderr, "sum: %f @ (%f, %f)\n", source->moments->Sum, source->moments->x, source->moments->y);
+
+    EllipseMoments moments;
+    moments.x2 = source->moments->Sx;
+    moments.y2 = source->moments->Sy;
+    moments.xy = source->moments->Sxy;
+    EllipseAxes axes = EllipseMomentsToAxes (moments);
+
+    fprintf (stderr, "axes: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);
+
+    // get the initial model parameter guess
+    pmModel *model = pmSourceModelGuess (source, modelType);
+    // if any parameters are defined, use those values
+    int nParams = pmModelParameterCount (modelType);
+    psF32 *params = model->params->data.F32;
+    for (int i = 0; i < nParams; i++) {
+	if (i == 2) {
+	    params[i] = xObj;
+	    continue;
+	}
+	if (i == 3) {
+	    params[i] = yObj;
+	    continue;
+	}
+	sprintf (name, "TEST_FIT_PAR%d", i);
+	value = psMetadataLookupF32 (&status, config, name);
+	if (status) {
+	    params[i] = value;
+	}
+    }
+
+    float area = params[4]*params[5];
+    fprintf (stderr, "peak: %f @ (%f, %f)\n", source->moments->Sum*area, (double)source->peak->x, (double)source->peak->y);
+
+    if (psfMode) {
+	pmModel *modelPSF = pmModelFromPSF (model, psf);
+	psFree (model);
+	model = modelPSF;
+	params = model->params->data.F32;
+    }
+
+    // list model input shape
+    EllipseShape shape;
+    shape.sx  = 1.4 / model->params->data.F32[4];
+    shape.sy  = 1.4 / model->params->data.F32[5];
+    shape.sxy = model->params->data.F32[6];
+    axes = EllipseShapeToAxes (shape);
+
+    fprintf (stderr, "guess: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);
+
+
+    fprintf (stderr, "input parameters: \n");
+    for (int i = 0; i < nParams; i++) {
+	fprintf (stderr, "%d : %f\n", i, params[i]);
+    }
+
+    // define the pixels used for the fit
+    psImageKeepCircle (source->mask, xObj, yObj, RADIUS, "OR", PSPHOT_MASK_MARKED);
+
+    char *fitset = psMetadataLookupPtr (&status, config, "TEST_FIT_SET");
+    if (status) {
+	status = psphotFitSet (source, model, fitset, psfMode);
+	return status;
+    }
+
+    status = pmSourceFitModel (source, model, psfMode);
+
+    // measure the source mags
+    pmSourcePhotometry (&fitMag, &obsMag, model, source->pixels, source->mask);
+    fprintf (stderr, "ap: %f, fit: %f, apmifit: %f\n", obsMag, fitMag, obsMag - fitMag);
+
+    // write out positive object
+    psphotSaveImage (NULL, source->pixels, "object.fits");
+
+    // subtract object, leave local sky
+    pmSourceSubModel (source->pixels, source->mask, model, false, false);
+    
+    fprintf (stderr, "output parameters: \n");
+    for (int i = 0; i < nParams; i++) {
+	fprintf (stderr, "%d : %f\n", i, params[i]);
+    }
+
+    // write out 
+    psphotSaveImage (NULL, source->pixels, "resid.fits");
+    psphotSaveImage (NULL, source->mask, "mask.fits");
+    return true;
 }