src2/fr/orsay/lri/varna/models/geom/ComputeArcCenter.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 ComputeArcCenter {
   8
   9         /**
  10          * Given an arc length (l) and segment length (delta) of the arc,
  11          * find where to put the center, returned as a position of the perpendicular
  12          * bisector of the segment. The positive side is the one where the arc is drawn.
  13          * It works using Newton's method.
  14          */
  15         public static double computeArcCenter(double delta, double l) {
  16                 double x_n = 0;
  17                 double x_n_plus_1, f_x_n, f_x_n_plus_1;
  18                 f_x_n = f(x_n,delta);
  19                 while (true) {
  20                         x_n_plus_1 = x_n - (f_x_n - l)/fprime(x_n,delta);
  21                         f_x_n_plus_1 = f(x_n_plus_1,delta);
  22                         // We want a precision of 0.1 on arc length
  23                         if (x_n_plus_1 == Double.NEGATIVE_INFINITY || Math.abs(f_x_n_plus_1 - f_x_n) < 0.1) {
  24                                 //System.out.println("computeArcCenter: steps = " + steps + "    result = " + x_n_plus_1);
  25                                 return x_n_plus_1;
  26                         }
  27                         x_n = x_n_plus_1;
  28                         f_x_n = f_x_n_plus_1;
  29                 }
  30         }
  31
  32         private static double f(double c, double delta) {
  33                 if (c < 0) {
  34                         return 2*Math.atan(delta/(-2*c)) * Math.sqrt(delta*delta/4 + c*c);
  35                 } else if (c != 0) { // c > 0
  36                         return (2*Math.PI - 2*Math.atan(delta/(2*c))) * Math.sqrt(delta*delta/4 + c*c);
  37                 } else { // c == 0
  38                         return Math.PI * Math.sqrt(delta*delta/4 + c*c);
  39                 }
  40         }
  41
  42         /**
  43          * d/dc f(c,delta)
  44          */
  45         private static double fprime(double c, double delta) {
  46                 if (c < 0) {
  47                         return delta/(c*c + delta/4)*Math.sqrt(delta*delta/4 + c*c) + 2*Math.atan(delta/(-2*c))*c/Math.sqrt(delta*delta/4 + c*c);
  48                 } else if (c != 0) { // c > 0
  49                         return delta/(c*c + delta/4)*Math.sqrt(delta*delta/4 + c*c) + (2*Math.PI - 2*Math.atan(delta/(-2*c)))*c/Math.sqrt(delta*delta/4 + c*c);
  50                 } else { // c == 0
  51                         return 2;
  52                 }
  53         }
  54 }