X-Git-Url: http://source.jalview.org/gitweb/?a=blobdiff_plain;f=srcjar%2Ffr%2Forsay%2Flri%2Fvarna%2Fmodels%2Fgeom%2FComputeEllipseAxis.java;fp=srcjar%2Ffr%2Forsay%2Flri%2Fvarna%2Fmodels%2Fgeom%2FComputeEllipseAxis.java;h=b8d699b0802ce113602ef3c41c7f6ca1c2e426f3;hb=ec8f3cedf60fb1feed6d34de6b49f6bfa78b9dd8;hp=0000000000000000000000000000000000000000;hpb=056dad85a910551cc95e44d451a61f6b8c4dd35d;p=jalview.git

diff --git a/srcjar/fr/orsay/lri/varna/models/geom/ComputeEllipseAxis.java b/srcjar/fr/orsay/lri/varna/models/geom/ComputeEllipseAxis.java
new file mode 100644
index 0000000..b8d699b
--- /dev/null
+++ b/srcjar/fr/orsay/lri/varna/models/geom/ComputeEllipseAxis.java
@@ -0,0 +1,92 @@
+/**
+ * File written by Raphael Champeimont
+ * UMR 7238 Genomique des Microorganismes
+ */
+package fr.orsay.lri.varna.models.geom;
+
+public class ComputeEllipseAxis {
+	public static boolean debug = false;
+	
+	/**
+	 * Given one axis half-length b and the circumference l,
+	 * find the other axis half-length a.
+	 * Returns 0 if there is no solution.
+	 * Implemented with Newton's method.
+	 */
+	public static double computeEllipseAxis(double b, double l) {
+		if (l/4 <= b || b <= 0 || l <= 0) {
+			// No such ellipse can exist.
+			return 0;
+		} else {
+			int steps = 0;
+			double x_n = 10;
+			double x_n_plus_1, f_x_n, f_x_n_plus_1;
+			f_x_n = f(x_n,b) - l;
+			while (true) {
+				//System.out.println("x_n = " + x_n + "  f(x_n)=" + f_x_n);
+				x_n_plus_1 = x_n - f_x_n/fprime(x_n,b);
+				f_x_n_plus_1 = f(x_n_plus_1,b) - l;
+				if (x_n_plus_1 < 0) {
+					System.out.println("ComputeEllipseAxis: x_n < 0 => returning 0");
+					return 0;
+				}
+				// We want a precision of 0.01 on arc length
+				if (Math.abs(f_x_n_plus_1 - f_x_n) < 0.01) {
+					if (debug) System.out.println("#steps = " + steps);
+					return x_n_plus_1;
+				}
+				x_n = x_n_plus_1;
+				f_x_n = f_x_n_plus_1;
+				steps++;
+			}
+		}
+	}
+	
+	
+	
+	private static double f(double a, double b) {
+		// This is Ramanujan's approximation of an ellipse circumference
+		// Nice because it is fast to compute (no need to compute an integral)
+		// and the derivative is also simple (and fast to compute).
+		return Math.PI*(3*(a+b) - Math.sqrt(10*a*b + 3*(a*a + b*b)));
+	}
+	
+	private static double fprime(double a, double b) {
+		return Math.PI*(3 - (5*b + 3*a)/Math.sqrt(10*a*b + 3*(a*a + b*b)));
+	}
+	
+	
+	/*
+	private static void test(double a, double b, double l) {
+		double a2 = computeEllipseAxis(b, l);
+		System.out.println("true a=" + a + " l=" + l + "   estimated a=" + a2 + " l(from true a)=" + f(a,b));
+	}
+	
+	private static void test(double b, double l) {
+		double a2 = computeEllipseAxis(b, l);
+		System.out.println("true l=" + l + "   estimated a=" + a2);
+	}
+	
+	
+	public static void main(String[] args) {
+		double b = 4;
+		test(100, b, 401.3143);
+		test(7, b, 35.20316);
+		test(4, b, 25.13274);
+		test(1, b, 17.15684);
+		test(0.5, b, 16.37248);
+		test(0.25, b, 16.11448);
+		test(0.1, b, 16.02288);
+		test(0.01, b, 16.00034);
+		test(0, b, 16);
+		test(10000, b, 40000.03);
+	}
+	*/
+	
+	/*
+	public static void main(String[] args) {
+		test(222.89291150684895, 2240);
+	}
+	*/
+	
+}