Index: trunk/psModules/src/objects/pmPSF.c
===================================================================
--- trunk/psModules/src/objects/pmPSF.c	(revision 13034)
+++ trunk/psModules/src/objects/pmPSF.c	(revision 13064)
@@ -6,6 +6,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.21 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-26 01:20:29 $
+ *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-27 22:14:08 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -236,4 +236,84 @@
 }
 
+// New Concept: the PSF modelling function fits the polarization terms e0, e1, e2:
+
+// convert the parameters used in the fitted source model 
+// to the parameters used in the 2D PSF model
+bool pmPSF_FitToModel (psF32 *fittedPar, float minMinorAxis)
+{
+    psEllipsePol pol;
+
+    pol.e0 = fittedPar[PM_PAR_E0];
+    pol.e1 = fittedPar[PM_PAR_E1];
+    pol.e2 = fittedPar[PM_PAR_E2];
+
+    psEllipseAxes axes = psEllipsePolToAxes (pol, minMinorAxis);
+    psEllipseShape shape = psEllipseAxesToShape (axes);
+
+    fittedPar[PM_PAR_SXX] = shape.sx * M_SQRT2;
+    fittedPar[PM_PAR_SYY] = shape.sy * M_SQRT2;
+    fittedPar[PM_PAR_SXY] = shape.sxy;
+
+    return true;
+}
+
+// convert the PSF parameters used in the 2D PSF model fit into the
+// parameters used in the source model
+psEllipsePol pmPSF_ModelToFit (psF32 *modelPar)
+{
+    psEllipseShape shape;
+
+    shape.sx  = modelPar[PM_PAR_SXX] / M_SQRT2;
+    shape.sy  = modelPar[PM_PAR_SYY] / M_SQRT2;
+    shape.sxy = modelPar[PM_PAR_SXY];
+    
+    psEllipsePol pol = psEllipseShapeToPol (shape);
+    
+    return pol;
+}
+
+// convert the parameters used in the fitted source model to the psEllipseAxes representation
+// (major,minor,theta)
+psEllipseAxes pmPSF_ModelToAxes (psF32 *modelPar, double maxAR)
+{
+    psEllipseShape shape;
+    psEllipseAxes axes;
+
+    shape.sx  = modelPar[PM_PAR_SXX] / M_SQRT2;
+    shape.sy  = modelPar[PM_PAR_SYY] / M_SQRT2;
+    shape.sxy = modelPar[PM_PAR_SXY];
+
+    if ((shape.sx == 0) || (shape.sy == 0)) {
+	axes.major = 0.0;
+	axes.minor = 0.0;
+	axes.theta = 0.0;
+    } else {
+	// XXX this is not really consistent with the model fit range above
+	axes = psEllipseShapeToAxes (shape, maxAR);
+    }
+    
+    return axes;
+}
+
+// convert the psEllipseAxes representation (major,minor,theta) to the parameters used in the
+// fitted source model
+bool pmPSF_AxesToModel (psF32 *modelPar, psEllipseAxes axes)
+{
+    if ((axes.major <= 0) || (axes.minor <= 0)) {
+	modelPar[PM_PAR_SXX] = 0.0;
+	modelPar[PM_PAR_SYY] = 0.0;
+	modelPar[PM_PAR_SXY] = 0.0;
+	return true;
+    }	
+    
+    psEllipseShape shape = psEllipseAxesToShape (axes);
+
+    modelPar[PM_PAR_SXX] = shape.sx * M_SQRT2;
+    modelPar[PM_PAR_SYY] = shape.sy * M_SQRT2;
+    modelPar[PM_PAR_SXY] = shape.sxy;
+
+    return true;
+}
+
 // generate a psf model of the requested type, with fixed shape
 pmPSF *pmPSFBuildSimple (char *typeName, float sxx, float syy, float sxy, ...)
@@ -247,13 +327,21 @@
     pmPSF *psf = pmPSFAlloc (type, true, psfTrend);
 
+    psVector *par = psVectorAlloc (psf->params_NEW->n, PS_TYPE_F32);
+    par->data.F32[PM_PAR_SXX] = sxx;
+    par->data.F32[PM_PAR_SYY] = syy;
+    par->data.F32[PM_PAR_SXY] = sxy;
+
+    psEllipsePol pol = pmPSF_ModelToFit(par->data.F32);
+
     // set the psf shape parameters
-    psPolynomial2D *poly = psf->params_NEW->data[PM_PAR_SXX];
-    poly->coeff[0][0] = M_SQRT2*sxx;
-
-    poly = psf->params_NEW->data[PM_PAR_SYY];
-    poly->coeff[0][0] = M_SQRT2*syy;
-
-    poly = psf->params_NEW->data[PM_PAR_SXY];
-    poly->coeff[0][0] = sxy;
+    psPolynomial2D *poly = NULL;
+    poly = psf->params_NEW->data[PM_PAR_E0];
+    poly->coeff[0][0] = pol.e0;
+
+    poly = psf->params_NEW->data[PM_PAR_E1];
+    poly->coeff[0][0] = pol.e1;
+
+    poly = psf->params_NEW->data[PM_PAR_E2];
+    poly->coeff[0][0] = pol.e2;
 
     for (int i = PM_PAR_SXY + 1; i < psf->params_NEW->n; i++) {
@@ -263,4 +351,5 @@
     va_end(ap);
 
+    psFree (par);
     psFree (psfTrend);
     return psf;
