From 530270b5afd3a88a01ce184717665c247ab3b3b0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= <marek@necada.org>
Date: Thu, 6 Dec 2018 20:42:11 +0000
Subject: [PATCH] lattices2d.c implement finding the "base vector hexagon"
 inscribed sphere radius

Former-commit-id: 075a19b07a8066d31378990470a90c4a7e8689ce
---
 qpms/lattices.h         |  3 +++
 qpms/lattices2d.c       | 41 +++++++++++++++++++++++++++++++++++++++--
 tests/test_lattices2d.c | 18 ++++++++++++++++++
 3 files changed, 60 insertions(+), 2 deletions(-)
 create mode 100644 tests/test_lattices2d.c

diff --git a/qpms/lattices.h b/qpms/lattices.h
index b35ea1e..266ae75 100644
--- a/qpms/lattices.h
+++ b/qpms/lattices.h
@@ -227,6 +227,9 @@ void l2d_reciprocalBasis1(cart2_t b1, cart2_t b2, cart2_t *rb1, cart2_t *rb2);
 void l2d_reciprocalBasis2pi(cart2_t b1, cart2_t b2, cart2_t *rb1, cart2_t *rb2);
 
 
+// returns the radius of inscribed circle of a hexagon (or rectangle/square if applicable) created by the shortest base triple
+double l2d_hexWebInCircleRadius(cart2_t b1, cart2_t b2);
+
 /*
  * THE MORE OR LESS OK PART
  * ========================
diff --git a/qpms/lattices2d.c b/qpms/lattices2d.c
index b1a1886..6f4e1b4 100644
--- a/qpms/lattices2d.c
+++ b/qpms/lattices2d.c
@@ -704,13 +704,36 @@ bool l2d_is_obtuse_r(cart2_t b1, cart2_t b2, double rtol) {
 }
 
 
-# if 0
 /*
  * TODO doc
  * return value is 4 or 6.
  */
-int l2d_shortestBase46(cart2_t i1, cart2_t i2,  cart2_t *o1, cart2_t *o2, cart2_t *o3, cart2_t *o4, cart2_t *o5, cart2_t *o6, double rtol);
+int l2d_shortestBase46(const cart2_t i1, const cart2_t i2,  cart2_t *o1, cart2_t *o2, cart2_t *o3, cart2_t *o4, cart2_t *o5, cart2_t *o6, double rtol){
+  cart2_t b1, b2, b3;
+  l2d_reduceBasis(i1, i2, &b1, &b2);
+  const double b1s = cart2_normsq(b1);
+  const double b2s = cart2_normsq(b2);
+  b3 = cart2_substract(b2, b1);
+  const double b3s = cart2_normsq(b3);
+  const double eps = rtol * (b1s + b2s); // TODO check the same as in l2d_is_obtuse_r
+  if(fabs(b3s-b2s-b1s) < eps) {
+    *o1 = b1; *o2 = b2; *o3 = cart2_scale(-1, b1); *o4 = cart2_scale(-1, b2);
+    return 4;
+  }
+  else {
+    if (b3s-b2s-b1s > eps) { //obtuse 
+      b3 = b2;
+      b2 = cart2_add(b2, b1);
+    }
+    *o1 = b1; *o2 = b2; *o3 = b3;
+    *o4 = cart2_scale(-1, b1);
+    *o5 = cart2_scale(-1, b2);
+    *o6 = cart2_scale(-1, b3);
+    return 6;
+  }
+}
 
+# if 0
 // variant
 int l2d_shortestBase46_arr(cart2_t i1, cart2_t i2,  cart2_t *oarr, double rtol);
 
@@ -741,4 +764,18 @@ int l2d_cellCornersWS_arr(cart2_t i1, cart2_t i2,  cart2_t *oarr, double rtol);
 // Reciprocal bases
 void l2d_reciprocalBasis1(cart2_t b1, cart2_t b2, cart2_t *rb1, cart2_t *rb2);
 void l2d_reciprocalBasis2pi(cart2_t b1, cart2_t b2, cart2_t *rb1, cart2_t *rb2);
+
 #endif
+
+// returns the radius of inscribed circle of a hexagon (or rectangle/square if applicable) created by the shortest base triple
+double l2d_hexWebInCircleRadius(cart2_t i1, cart2_t i2) {
+  cart2_t b1, b2, b3;
+  l2d_shortestBase3(i1, i2, &b1, &b2, &b3);
+  const double r1 = cart2norm(b1), r2 = cart2norm(b2), r3 = cart2norm(b3);
+  const double p = (r1+r2+r3)*0.5;
+  return 2*sqrt(p*(p-r1)*(p-r2)*(p-r3))/r3; // CHECK is r3 guaranteed to be longest?
+}
+
+
+
+
diff --git a/tests/test_lattices2d.c b/tests/test_lattices2d.c
new file mode 100644
index 0000000..dbf60f5
--- /dev/null
+++ b/tests/test_lattices2d.c
@@ -0,0 +1,18 @@
+// c99 -o l2dtest -ggdb -I ../ test_lattices2d.c ../qpms/lattices2d.c -lm
+#include <qpms/lattices.h>
+#include <stdio.h>
+
+int main(int argc, char **argv) {
+  cart2_t b1 = {1.2,0.1};
+  cart2_t b2 = {0.5, 0.5};
+  cart2_t b3;
+
+  printf("original b1 = (%g, %g), b2 = (%g, %g)\n",
+      b1.x, b1.y, b2.x, b2.y);
+  printf("Inscribed circle radius: %g\n",
+      l2d_hexWebInCircleRadius(b1, b2));
+  l2d_shortestBase3(b1, b2, &b1, &b2, &b3);
+  printf("shortestBase3: b1 = (%g, %g), b2 = (%g, %g), b3 = (%g, %g)\n",
+      b1.x, b1.y, b2.x, b2.y, b3.x, b3.y);
+  return 0;
+}