new energy functions

[moebius2.git] / energy.c

diff --git a/energy.c b/energy.c
index 867330cbf2da74a6293cdf55c58ba58edffc77af..b53199e2555bf8d28be3c258ad5ddcf4a6da71a4 100644 (file)
--- a/energy.c
+++ b/energy.c
@@ -40,7 +40,7 @@ static double compute_energy(const Vertices vertices) {
 
   COST(1e2, edgewise_vertex_displacement_cost(vertices));
 //  COST(1e0, graph_layout_cost(vertices,vertex_areas));
-  COST(1e4, noncircular_rim_cost(vertices));
+//  COST(1e4, noncircular_rim_cost(vertices));
 
   if (printing) printf("| total %# e |", energy);
 
@@ -222,7 +222,7 @@ int main(int argc, const char *const *argv) {
    *
    *  Giving energy contribution:
    *
-   *                                   2
+   *                                   3
    *                            l delta
    *    E             =  F   .  --------
    *     vd, edge PQ      vd       d
@@ -231,7 +231,7 @@ int main(int argc, const char *const *argv) {
    *  (The dimensions of this are those of F_vd.)
    *
    *  We calculate delta as  atan2(|AxB|, A.B)
-   *  where A = RM, B = MS
+   *  where A = PQ, B = QR
    *
    *  In practice to avoid division by zero we'll add epsilon to d and
    *  |AxB| and the huge energy ought then to be sufficient for the
@@ -241,22 +241,25 @@ int main(int argc, const char *const *argv) {
 double edgewise_vertex_displacement_cost(const Vertices vertices) {
   static const double axb_epsilon= 1e-6;
 
-  int pi,e,qi,ri,si, k;
-  double m[D3], a[D3], b[D3], axb[D3];
+  int pi,e,qi,ri, k; //,si
+  double  a[D3], b[D3], axb[D3]; //m[D3],
   double total_cost= 0;
 
-  FOR_EDGE(pi,e,qi) {
-    ri= EDGE_END2(pi,(e+1)%V6); if (ri<0) continue;
-    si= EDGE_END2(pi,(e+5)%V6); if (si<0) continue;
+  FOR_EDGE(qi,e,ri) {
+    pi= EDGE_END2(qi,(e+3)%V6); if (pi<0) continue;
 
-    K m[k]= (vertices[pi][k] + vertices[qi][k]) * 0.5;
-    K a[k]= -vertices[ri][k] + m[k];
-    K b[k]= -m[k] + vertices[si][k];
+//    K m[k]= (vertices[pi][k] + vertices[qi][k]) * 0.5;
+    K a[k]= -vertices[pi][k] + vertices[qi][k];
+    K b[k]= -vertices[qi][k] + vertices[ri][k];
 
     xprod(axb,a,b);
     
     double delta= atan2(magnD(axb) + axb_epsilon, dotprod(a,b));
-    double cost= delta * delta;
+    double cost= pow(delta,3);
+
+    if (!e && !(qi & YMASK))
+      cost *= 1e3;
+
     total_cost += cost;
   }
   return total_cost;