+ *
+ * Then the (1/delta)th power of the cost is
+ * proportional to |D|, and
+ * inversely proportional to |SM|
+ * except that |D| is measured in a wierd way which counts
+ * distance in the same direction as SM 1/lambda times as much
+ * ie the equipotential surfaces are ellipsoids around R',
+ * lengthened by lambda in the direction of RM.
+ *
+ * So
+ * delta
+ * [ -1 ]
+ * cost = [ lambda . ( D . SM/|SM| ) + | D x SM/|SM| | ]
+ * R,S,M [ ------------------------------------------- ]
+ * [ |SM| ]
+ *
+ */
+
+static double vertex_one_displ_cost(const double r[D3], const double s[D3],
+ const double midpoint[D3],
+ double delta, double inv_lambda) {
+ const double smlen2_epsilon= 1e-12;
+ double sm[D3], d[D3], ddot, dcross[D3];
+ int k;
+
+ K sm[k]= -s[k] + midpoint[k];
+ K d[k]= midpoint[k] + sm[k] - r[k];
+ ddot= dotprod(d,sm);
+ xprod(dcross, d,sm);
+ double smlen2= magnD2(sm);
+ double cost_basis= inv_lambda * ddot + magnD(dcross);
+ double cost= pow(cost_basis / (smlen2 + smlen2_epsilon), delta);
+
+ return cost;
+}
+
+/*---------- displacement of vertices opposite at a vertex ----------*/
+
+ /*
+ * At vertex Q considering edge direction e to R
+ * and corresponding opposite edge to S.
+ *
+ * This is vertex displacement as above with M=Q