Index: trunk/psLib/src/math/psEllipse.c
===================================================================
--- trunk/psLib/src/math/psEllipse.c	(revision 7766)
+++ trunk/psLib/src/math/psEllipse.c	(revision 9768)
@@ -20,21 +20,23 @@
 }
 
+// ellipse rotation (major, minor, theta) -> (sx, sy, sxy)
 psEllipseShape psEllipseAxesToShape(psEllipseAxes axes)
 {
     psEllipseShape shape;
 
-    double r1 = 1.0 / PS_SQR(axes.major) + 1.0 / PS_SQR(axes.minor);
-    double r2 = 1.0 / PS_SQR(axes.major) - 1.0 / PS_SQR(axes.minor);
+    double f1 = 1.0 / PS_SQR(axes.major) + 1.0 / PS_SQR(axes.minor);
+    double f2 = 1.0 / PS_SQR(axes.major) - 1.0 / PS_SQR(axes.minor);
 
-    double sxr = r1 + r2*cos(2*axes.theta);
-    double syr = r1 - r2*cos(2*axes.theta);
+    double sxr = 0.5*f1 + 0.5*f2*cos(2*axes.theta);
+    double syr = 0.5*f1 - 0.5*f2*cos(2*axes.theta);
 
     shape.sx = 1.0 / sqrt(sxr);
     shape.sy = 1.0 / sqrt(syr);
-    shape.sxy = r2*sin(2*axes.theta);
+    shape.sxy = 0.5*r2*sin(2*axes.theta);
 
     return shape;
 }
 
+// ellipse derotation (sx, sy, sxy) -> (major, minor, theta)
 psEllipseAxes psEllipseShapeToAxes(psEllipseShape shape)
 {
@@ -43,18 +45,14 @@
     double f1 = 1.0 / PS_SQR(shape.sx) + 1.0 / PS_SQR(shape.sy);
     double f2 = 1.0 / PS_SQR(shape.sx) - 1.0 / PS_SQR(shape.sy);
+    double f3 = PS_SQR(f2) + 4*PS_SQR(shape.sxy);
 
-    // force the axis ratio to be less than 10
-    double r1 = 0.5*0.95*sqrt (PS_SQR(f1) - PS_SQR(f2));
-
-    shape.sxy = PS_MIN(PS_MAX(shape.sxy, -r1), r1);
-
-    axes.theta = atan2 (-2.0*shape.sxy, f2) / 2.0;
-
-    double Ar = 0.25*f1 + 0.25*sqrt(PS_SQR(f2) + 4*PS_SQR(shape.sxy));
-    double Br = 0.25*f1 - 0.25*sqrt(PS_SQR(f2) + 4*PS_SQR(shape.sxy));
-
-    axes.minor = 1.0 / sqrt (Ar);
-    axes.major = 1.0 / sqrt (Br);
+    axes.theta = 0.5 * atan2 (-2.0*shape.sxy, f2) / 2.0;
+    axes.major = sqrt (2.0 / (f1 - f3));
+    axes.minor = sqrt (2.0 / (f1 + f3));
 
     return axes;
 }
+
+// XXX keep this construction?
+// force the axis ratio to be less than 10
+// double r1 = 0.5*0.95*sqrt (PS_SQR(f1) - PS_SQR(f2));