From: Ian Jackson <ian@davenant.relativity.greenend.org.uk>
Date: Sun, 19 Oct 2008 21:57:14 +0000 (+0100)
Subject: new rim_twist_cost as yet unused
X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ian/git?p=moebius2.git;a=commitdiff_plain;h=820c80f5cc6696dbf4fef24fefe9d908df88ce31

new rim_twist_cost as yet unused
---

diff --git a/energy.c b/energy.c
index 72c6fd3..aefdd8b 100644
--- a/energy.c
+++ b/energy.c
@@ -497,6 +497,46 @@ double noncircular_rim_cost(const Vertices vertices, int section) {
   return cost;
 }
 
+/*---------- rim contact angle rotation ----------*/
+
+/*
+ *
+ *               vq ---- vs
+ *              /       /
+ *           e0/       /e1
+ *           vp 0--- vr
+ *
+ */
+
+double rim_twist_cost(const Vertices vertices, int section) {
+  const double our_epsilon=1e-6;
+  int vpy,vpx,vpi, vqi,vri,vsi, e0,e1, k;
+  double total_cost= 0.0;
+  double pq[D3], rs[D3], pq_x_rs[D3];
+
+  FOR_RIM_VERTEX(vpy,vpx,vpi, OUTER) {
+    FOR_VPEDGE(e0) {
+      if (e0==0 || e0==3) continue;
+      vqi= EDGE_END2(vpi,e0);  if (vqi<0) continue;
+      vri= EDGE_END2(vpi,0);  assert(vri>=0);
+      e1= !vertices_span_join_p(vpi,vri) ? e0 : V6 - e0;
+      vsi= EDGE_END2(vri,e1);  assert(vsi>=0);
+
+      K pq[k]= -vertices[vpi][k] + vertices[vqi][k];
+      K rs[k]= -vertices[vri][k] + vertices[vsi][k];
+
+      xprod(pq_x_rs, pq,rs);
+      double magndiv= edge_lengths[vpi][e0] * edge_lengths[vri][e1];
+
+      double cost= magnD2(pq_x_rs) / (magndiv*magndiv + our_epsilon);
+
+      total_cost += cost;
+    }
+  }
+
+  return total_cost;
+}
+
 /*---------- overly sharp edge cost ----------*/
 
   /*
diff --git a/minimise.h b/minimise.h
index f4eddfb..4999fd0 100644
--- a/minimise.h
+++ b/minimise.h
@@ -27,6 +27,7 @@ double noncircular_rim_cost(const Vertices vertices, int section);
 double edge_length_variation_cost(const Vertices vertices, int section);
 double prop_edge_length_variation_cost(const Vertices vertices, int section);
 double rim_proximity_cost(const Vertices vertices, int section);
+double rim_twist_cost(const Vertices vertices, int section);
 double edge_angle_cost(const Vertices vertices, int section);
 double small_triangles_cost(const Vertices vertices, int section);
 double nonequilateral_triangles_cost(const Vertices vertices, int section);