From b0b2d52a2f25d623be2498e31df9939286383722 Mon Sep 17 00:00:00 2001
From: Dan Goodliffe <dan@randomdan.homeip.net>
Date: Sat, 27 Feb 2021 20:15:28 +0000
Subject: Calculate the radius to join to point+direction vector pairs

This uses a mental formula that was derived using symbolabs.com, it works but there just has
to be simpler form of it!
---
 test/test-maths.cpp |  6 +++++-
 utility/maths.cpp   | 30 ++++++++++++++++++++++++++++++
 utility/maths.h     |  1 +
 3 files changed, 36 insertions(+), 1 deletion(-)

diff --git a/test/test-maths.cpp b/test/test-maths.cpp
index d5848c6..6151c39 100644
--- a/test/test-maths.cpp
+++ b/test/test-maths.cpp
@@ -2,7 +2,6 @@
 
 #include <boost/test/data/test_case.hpp>
 #include <boost/test/unit_test.hpp>
-#include <iomanip>
 #include <stream_support.hpp>
 
 #include <glm/glm.hpp>
@@ -92,3 +91,8 @@ BOOST_DATA_TEST_CASE(test_find_arc_centre,
 	BOOST_CHECK_CLOSE(exp.y, c.first.y, 1);
 	BOOST_CHECK_EQUAL(lr, c.second);
 }
+
+BOOST_AUTO_TEST_CASE(test_find_arcs_radius)
+{
+	BOOST_CHECK_CLOSE(find_arcs_radius({10.32, 26.71}, {0.4, .92}, {2.92, 22.41}, {-0.89, -0.45}), 2.29, 1);
+}
diff --git a/utility/maths.cpp b/utility/maths.cpp
index 866097a..75346a0 100644
--- a/utility/maths.cpp
+++ b/utility/maths.cpp
@@ -83,3 +83,33 @@ find_arc_centre(glm::vec2 as, glm::vec2 ad, glm::vec2 bs, glm::vec2 bd)
 	}
 	throw std::runtime_error("no intersection");
 }
+
+float
+find_arcs_radius(glm::vec2 start, glm::vec2 ad, glm::vec2 end, glm::vec2 bd)
+{
+	// Short name functions for big forula
+	auto sq = [](auto v) {
+		return v * v;
+	};
+	auto sqrt = [](float v) {
+		return std::sqrt(v);
+	};
+
+	// Calculates path across both arcs along the normals... pythagorean theorem... for some known radius r
+	// (2r)^2 = ((m + (X*r)) - (o + (Z*r)))^2 + ((n + (Y*r)) - (p + (W*r)))^2
+	// According to symbolabs.com equation tool, that solves for r to give:
+	// r=(-2 m X+2 X o+2 m Z-2 o Z-2 n Y+2 Y p+2 n W-2 p W-sqrt((2 m X-2 X o-2 m Z+2 o Z+2 n Y-2 Y p-2 n W+2 p W)^(2)-4
+	// (X^(2)-2 X Z+Z^(2)+Y^(2)-2 Y W+W^(2)-4) (m^(2)-2 m o+o^(2)+n^(2)-2 n p+p^(2))))/(2 (X^(2)-2 X Z+Z^(2)+Y^(2)-2 Y
+	// W+W^(2)-4))
+
+	// These exist cos limitations of online formula rearrangement, and I'm OK with that.
+	const auto &m {start.x}, &n {start.y}, &o {end.x}, &p {end.y};
+	const auto &X {ad.x}, &Y {ad.y}, &Z {bd.x}, &W {bd.y};
+
+	return (2 * m * X - 2 * X * o - 2 * m * Z + 2 * o * Z + 2 * n * Y - 2 * Y * p - 2 * n * W + 2 * p * W
+				   - sqrt(sq(-2 * m * X + 2 * X * o + 2 * m * Z - 2 * o * Z - 2 * n * Y + 2 * Y * p + 2 * n * W
+								  - 2 * p * W)
+						   - (4 * (sq(X) - 2 * X * Z + sq(Z) + sq(Y) - 2 * Y * W + sq(W) - 4)
+								   * (sq(m) - 2 * m * o + sq(o) + sq(n) - 2 * n * p + sq(p)))))
+			/ (2 * (sq(X) - 2 * X * Z + sq(Z) + sq(Y) - 2 * Y * W + sq(W) - 4));
+}
diff --git a/utility/maths.h b/utility/maths.h
index 2ce5d94..b0048c3 100644
--- a/utility/maths.h
+++ b/utility/maths.h
@@ -55,5 +55,6 @@ float normalize(float ang);
 
 std::pair<glm::vec2, bool> find_arc_centre(glm::vec2 start, float entrys, glm::vec2 end, float entrye);
 std::pair<glm::vec2, bool> find_arc_centre(glm::vec2 start, glm::vec2 ad, glm::vec2 end, glm::vec2 bd);
+float find_arcs_radius(glm::vec2 start, glm::vec2 ad, glm::vec2 end, glm::vec2 bd);
 
 #endif
-- 
cgit v1.2.3