src2/fr/orsay/lri/varna/models/geom/ComputeEllipseAxis.java

   1 /**
   2  * File written by Raphael Champeimont
   3  * UMR 7238 Genomique des Microorganismes
   4  */
   5 package fr.orsay.lri.varna.models.geom;
   6
   7 public class ComputeEllipseAxis {
   8         public static boolean debug = false;
   9
  10         /**
  11          * Given one axis half-length b and the circumference l,
  12          * find the other axis half-length a.
  13          * Returns 0 if there is no solution.
  14          * Implemented with Newton's method.
  15          */
  16         public static double computeEllipseAxis(double b, double l) {
  17                 if (l/4 <= b || b <= 0 || l <= 0) {
  18                         // No such ellipse can exist.
  19                         return 0;
  20                 } else {
  21                         int steps = 0;
  22                         double x_n = 10;
  23                         double x_n_plus_1, f_x_n, f_x_n_plus_1;
  24                         f_x_n = f(x_n,b) - l;
  25                         while (true) {
  26                                 //System.out.println("x_n = " + x_n + "  f(x_n)=" + f_x_n);
  27                                 x_n_plus_1 = x_n - f_x_n/fprime(x_n,b);
  28                                 f_x_n_plus_1 = f(x_n_plus_1,b) - l;
  29                                 if (x_n_plus_1 < 0) {
  30                                         System.out.println("ComputeEllipseAxis: x_n < 0 => returning 0");
  31                                         return 0;
  32                                 }
  33                                 // We want a precision of 0.01 on arc length
  34                                 if (Math.abs(f_x_n_plus_1 - f_x_n) < 0.01) {
  35                                         if (debug) System.out.println("#steps = " + steps);
  36                                         return x_n_plus_1;
  37                                 }
  38                                 x_n = x_n_plus_1;
  39                                 f_x_n = f_x_n_plus_1;
  40                                 steps++;
  41                         }
  42                 }
  43         }
  44
  45
  46
  47         private static double f(double a, double b) {
  48                 // This is Ramanujan's approximation of an ellipse circumference
  49                 // Nice because it is fast to compute (no need to compute an integral)
  50                 // and the derivative is also simple (and fast to compute).
  51                 return Math.PI*(3*(a+b) - Math.sqrt(10*a*b + 3*(a*a + b*b)));
  52         }
  53
  54         private static double fprime(double a, double b) {
  55                 return Math.PI*(3 - (5*b + 3*a)/Math.sqrt(10*a*b + 3*(a*a + b*b)));
  56         }
  57
  58
  59         /*
  60         private static void test(double a, double b, double l) {
  61                 double a2 = computeEllipseAxis(b, l);
  62                 System.out.println("true a=" + a + " l=" + l + "   estimated a=" + a2 + " l(from true a)=" + f(a,b));
  63         }
  64
  65         private static void test(double b, double l) {
  66                 double a2 = computeEllipseAxis(b, l);
  67                 System.out.println("true l=" + l + "   estimated a=" + a2);
  68         }
  69
  70
  71         public static void main(String[] args) {
  72                 double b = 4;
  73                 test(100, b, 401.3143);
  74                 test(7, b, 35.20316);
  75                 test(4, b, 25.13274);
  76                 test(1, b, 17.15684);
  77                 test(0.5, b, 16.37248);
  78                 test(0.25, b, 16.11448);
  79                 test(0.1, b, 16.02288);
  80                 test(0.01, b, 16.00034);
  81                 test(0, b, 16);
  82                 test(10000, b, 40000.03);
  83         }
  84         */
  85
  86         /*
  87         public static void main(String[] args) {
  88                 test(222.89291150684895, 2240);
  89         }
  90         */
  91
  92 }