initial checkin

diff --git a/anneal.c b/anneal.c
new file mode 100644 (file)
index 0000000..e13c0b0
--- /dev/null
+++ b/anneal.c
@@ -0,0 +1,117 @@
+/*
+ * We try to find an optimal triangle grid
+ *
+ * Vertices in strip are numbered as follows:
+ *
+ *     ___ L-2 ___ L-1 ___ 0   ___ 1   ___ 2   ___ 3   ___ 4   __
+ *         W-1     W-1      0       0       0       0       0
+ *        /  \    /  \    /  \    /  \    /  \    /  \    /  \
+ *       /    \  /    \  /    \  /    \  /    \  /    \  /    \
+ *     L-3 ___ L-2 ___ L-1 ___ 0   ___ 1   ___ 2   ___  3  ___ 4
+ *     W-2     W-2     W-2      1       1       1       1       1
+ *       \    /  \    /  \    /  \    /  \    /  \    /  \    /
+ *        \  /    \  /    \  /    \  /    \  /    \  /    \  /
+ *     ___ 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
+ *
+ * Node x,y for
+ *   0 <= x < L     x = distance along    L = length
+ *   0 <= y < W     y = distance across   W = width
+ *
+ * Vertices are in reading order from diagram above ie x varies fastest.
+ *
+ * Where necessary, edges are counted anticlockwise from to-the-right:
+ *
+ *                         \2   /1
+ *                          \  /
+ *                       ___ 0   __
+ *                       3    1   0
+ *                          /  \
+ *                        4/   5\
+ *
+ * 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 W 10
+#define N (L*W)
+
+static double energy_function(const double vertices[N][3]) {
+  /* Energy has the following parts right now:
+   *  - worst curvature radius
+   *  - total angle subtended
+   */
+  /*
+   *
+   *               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'
+   *
+   *
+   *         R'.
+   *                    \.
+   *              \.
+   *                 0
+   *                 |
+   *                 |
+   *                        |
+   *                S'
+   *
+   *  Now radius formula (wikipedia `Circle'), P1=R', P2=0, P3=S':
+   *
+   *          |R'| * |S'| * | R' - S' |
+   *     r =  -------------------------
+   *               2 * | R' x S' |
+   *
+   *  Energy
+   *                            2
+   *     E  = C  min           r
+   *      r    r   j in edges   j
+   *
+   */
+  /*
+   * Angle subtended:
+   *
+   *                -1   | R' . S' |
+   *    gamma =  cos     -----------
+   *                     |R'| * |S'|
+   *
+   *    l = | P - Q |
+   *
+   *  Energy
+   *
+   *    E      = C      sum           gamma  * l
+   *     gamma    gamma   j in edges       j    j
+   *
+   */
+
+static gsl_multimin_fminimizer *minimiser;
+
+
+
+{
+  gsl_multimin_fminimizer_alloc(