From: Ian Jackson <ian@davenant.relativity.greenend.org.uk>
Date: Thu, 8 Nov 2007 01:27:04 +0000 (+0000)
Subject: new energy md^2/l giving sum proportional to Lambda * Gamma
X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ian/git?p=moebius2.git;a=commitdiff_plain;h=4492246818e419bdc258c5842b381b0edbaad5df

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

diff --git a/anneal.c b/anneal.c
index e13c0b0..a3a4258 100644
--- 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(
+