cube l

diff --git a/anneal.c b/anneal.c
index ab521f4700da30b1b9fa0d5989ffa6d6d89665f5..4ddddb8934901dce708ad7270094365eaae2f169 100644 (file)
--- a/anneal.c
+++ b/anneal.c
@@ -97,7 +97,8 @@
    *                               2
    *                            b d
    *    E             =  F   .  ----
-   *     vd, edge PQ      vd      l
+   *     vd, edge PQ      vd      3
+   *                             l
    *
    *  By symmetry, this calculation gives the same answer
    *  with R and S exchanged.  Looking at the projection in
@@ -122,10 +123,9 @@
    *  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.
-   */
+   *  In practice to avoid division by zero we'll add epsilon to l^3
+   *  and the huge energy ought then to be sufficient for the model to
+   *  avoid being close to R=S.  */
 
 #define V6 6
 
@@ -204,6 +204,8 @@ static double hypotD(const double p[], const double q[]) {
 # define ffsqu(factor,term) ((factor)*(factor)+(term))
 #endif
 
+static const l3_epsison= 1e-6;
+
 static double energy_function(const double vertices[N][D3]) {
   int pi,e,qi,ri,si, k;
   double m[D3], mprime[D3], b, d2, l, sigma_bd2_l;
@@ -219,8 +221,9 @@ static double energy_function(const double vertices[N][D3]) {
     b= hypotD(vertices[pi], vertices[qi]);
     d2= 0; K d2= ffsqa(m[k] - mprime[k], d2);
     l= hypotD(vertices[ri][k] - vertices[si][k]);
+    l3 = l*l*l + l3_epsilon;
 
-    sigma_bd2_l += b * d2 / l;
+    sigma_bd2_l += b * d2 / l3;
   }