new energy md^2/l giving sum proportional to Lambda * Gamma

diff --git a/anneal.c b/anneal.c
index e13c0b0cb8faaa4d03515148a5f32cf035c1ecd0..a3a425800bf26ebb1914c868b1f673a31f718c3d 100644 (file)
--- a/anneal.c
+++ b/anneal.c
@@ -49,64 +49,65 @@
 
 static double energy_function(const double vertices[N][3]) {
   /* Energy has the following parts right now:
-   *  - worst curvature radius
-   *  - total angle subtended
    */
   /*
+   *  Edgewise expected vertex displacement
    *
-   *               Q
-   *              / | `-.
-   *             /  |    `-.
-   *            /   P ------ S
-   *           /.-'
-   *          /'
-   *         R
-   */
-  /*
-   * Curvature radius:
-   *
-   *  First flatten into plane _|_ PQ:
    *
-   *   R' = (R-P) x (Q-P)/|Q-P|
-   *   S' = (S-P) x (Q-P)/|Q-P|
-   *   Q' = (Q-P) x (Q-P)/|Q-P| = 0 = P'
+   *                       Q `-_
+   *              / |    `-_
+   *   R' - _ _ _/_ |       `-.
+   *    .       /   M - - - - - S
+   *    .      /    |      _,-'
+   *    .     /     |  _,-'
+   *    .    /    , P '
+   *    .   /  ,-'
+   *    .  /,-'
+   *    . /'
+   *            R
    *
    *
-   *         R'.
-   *                    \.
-   *              \.
-   *                 0
-   *                 |
-   *                 |
-   *                        |
-   *                S'
    *
-   *  Now radius formula (wikipedia `Circle'), P1=R', P2=0, P3=S':
+   *  Find R', the `expected' location of R, by
+   *  reflecting S in M (the midpoint of QP).
    *
-   *          |R'| * |S'| * | R' - S' |
-   *     r =  -------------------------
-   *               2 * | R' x S' |
+   *  Let  d = |RR'|
+   *       m = |PQ|
+   *       l = |RS|
    *
-   *  Energy
-   *                            2
-   *     E  = C  min           r
-   *      r    r   j in edges   j
+   *  Giving energy contribution:
    *
-   */
-  /*
-   * Angle subtended:
+   *                               2
+   *                            m d 
+   *    E             =  F   .  ----
+   *     vd, edge PQ      vd      l
    *
-   *                -1   | R' . S' |
-   *    gamma =  cos     -----------
-   *                     |R'| * |S'|
+   *  By symmetry, this calculation gives the same answer
+   *  with R and S exchanged.  Looking at the projection in
+   *  the RMS plane:
    *
-   *    l = | P - Q |
    *
-   *  Energy
+   *                          S'
+   *                        ,'
+   *                      ,'
+   *           R'               ,'             d" = |SS'| = |RR'| = d
+   *     `-._         ,'
+   *                 `-._   ,'
+   *              ` M
+   *             ,'  `-._
+   *           ,'        `-._
+   *         ,'              ` S
+   *       ,'
+   *     ,'
+   *    R
    *
-   *    E      = C      sum           gamma  * l
-   *     gamma    gamma   j in edges       j    j
+   *  We choose this value for l (rather than |RM|+|MS|, say, or |RM|)
+   *  because we want this symmetry and because we're happy to punish
+   *  bending more than uneveness in the metric.
    *
+   *  In practice to avoid division by zero we'll add epsilon to l and
+   *  the huge energy ought then to be sufficient for the model to
+   *  avoid being close to R=S.
    */
 
 static gsl_multimin_fminimizer *minimiser;
@@ -115,3 +116,4 @@ static gsl_multimin_fminimizer *minimiser;
 
 {
   gsl_multimin_fminimizer_alloc(
+