* \ / \ / \ / \ / \ / \ / \ /
* ___ L-3 ___ L-2 ___ L-1 ___ 0 ___ 1 ___ 2 ___ 3 __
* W-3 W-3 W-3 2 2 2 2
- * / \ / \ / \ / \ / \ / \ / \
- * / \ / \ / \ / \ / \ / \ / \
- * . . . . . . . . . . . . . . . .
*
- * _ L-4 ___ L-3 ___ L-2 ___ L-1 ___ 0 ___ 1 ___ 2 ___ 3
- * 1 1 1 1 W-2 W-2 W-2 W-2
- * / \ / \ / \ / \ / \ / \ / \ /
- * / \ / \ / \ / \ / \ / \ / \ /
- * L-5 ___ L-4 ___ L-3 ___ L-2 ___ L-1 ___ 0 ___ 1 ___ 2 __
- * 0 0 0 0 0 W-1 W-1 W-1
+ * . . . . . . . . . . . . . . .
+ *
+ * L-4 ___ L-3 ___ L-2 ___ L-1 ___ 0 ___ 1 ___ 2 ___ 3
+ * 1 1 1 1 W-2 W-2 W-2 W-2
+ * \ / \ / \ / \ / \ / \ / \ /
+ * \ / \ / \ / \ / \ / \ / \ /
+ * ___ L-4 ___ L-3 ___ L-2 ___ L-1 ___ 0 ___ 1 ___ 2 __
+ * 0 0 0 0 W-1 W-1 W-1
*
* Node x,y for
* 0 <= x < L x = distance along L = length
*
* Vertices are in reading order from diagram above ie x varies fastest.
*
- * Where necessary, edges are counted anticlockwise from to-the-right:
+ * W must be even. The actual location opposite (0,0) is (L-(W-1)/2,0),
+ * and likewise opposite (0,W-1) is ((W-1)/2,0).
+ *
+ * To label edges, we counte anticlockwise[*] from to-the-right:
*
* \2 /1
* \ /
* / \
* 4/ 5\
*
+ * [*] That is, in the direction represented as anticlockwise for
+ * the vertices (0,*)..(4,*) in the diagram above; and of course
+ * that is clockwise for the vertices (L-5,*)..(L-1,*). The
+ * numbering has an actual discontinuity between (L-1,*) and (0,*).
+ *
* When we iterate over edges, we iterate first over vertices and then
* over edges 0 to 2, disregarding edges 3 to 5.
*/
-#define L 10
+#define LN2W 4
+#define L (1<<LN2W)
#define W 10
#define N (L*W)
+#define YMASK (~0u << LN2W)
+#define XMASK (~YMASK)
-static double energy_function(const double vertices[N][3]) {
/* Energy has the following parts right now:
*/
/*
* Giving energy contribution:
*
* 2
- * m d
+ * m d
* E = F . ----
* vd, edge PQ vd l
*
* avoid being close to R=S.
*/
+static int edge_end2(int v1, int e) {
+ if (!(v1 &
+
+ if (e==1 || e==2)
+ if (v2 < L)
+
+
+#define FOR_VERTEX(v) \
+ for ((v)=0; (v)<N; (v)++)
+
+#define FOR_VPEDGE(v,e) \
+ for ((e)=0; (e)<6; (e)++)
+
+#define EDGE_END2 edge_end2
+
+#define FOR_VEDGE(v1,e,v2) \
+ FOR_VPEDGE((v1),(e))
+ if (((v2)= EDGE_END2((v1),(e))) >= 0) ; else
+
+#define FOR_EDGE(v1,e,v2) \
+ FOR_VERTEX((v1)) \
+ FOR_VEDGE((v1),(e),(v2))
+
+static double energy_function(const double vertices[N][3]) {
+ int vertex;
+
+ FOR_VERTEX {
+
+
static gsl_multimin_fminimizer *minimiser;