1 files changed, 79 insertions, 0 deletions
diff --git a/szilagyi/_nomogram.py b/szilagyi/_nomogram.py
new file mode 100644
index 0000000..4dcd141
--- /dev/null
+++ b/szilagyi/_nomogram.py
@@ -0,0 +1,79 @@
+import math
+from collections import deque
+
+from . import _dataset
+
+
+class NomogramEdgeCase(Exception):
+	def __init__(self, value):
+		self.value = value
+
+
+def look_downwards(data, x, start):
+	for i in range(start, 0, -1):
+		if data[i - 1].x < x:
+			break
+	else:
+		raise IndexError
+	return i - 1
+
+
+def look_upwards(data, x, start):
+	for i in range(start, len(data)):
+		if data[i + 1].x > x:
+			break
+	else:
+		raise IndexError
+	return i
+
+
+def find_segment(data, x):
+	width = data[-1].x - data[0].x
+	relative = x - data[0].x
+	candidate = math.floor(relative / width * len(data))
+	look = look_downwards if data[candidate].x > x else look_upwards  # May raise IndexError
+	candidate = look(data, x, candidate)
+	return candidate, candidate + 1
+
+
+def find_boundary_curves(swis, x, y):
+	segments = deque()
+	for index, data in swis:
+		i, j = find_segment(data, x)
+		if data[i].y > y and data[j].y > y:
+			segments.append((index, data, i, j))
+			break
+		if data[i].y < y and data[j].y < y:
+			if segments:
+				segments.popleft()
+		segments.append((index, data, i, j))
+	if len(segments) == 3:
+		middle = segments[1][1]
+		run = middle[j].x - middle[i].x
+		if run == 0:
+			raise NomogramEdgeCase(segments[1][0])  # Tidy up the dataset nonetheless
+		slope = (middle[j].y - middle[i].y) / run
+		intercept = middle[j].y - slope * middle[j].x
+		value = slope * x + intercept
+		if value == y:
+			raise NomogramEdgeCase(segments[1][0])
+		if value < y:
+			segments.popleft()
+		else:
+			segments.pop()
+	if len(segments) == 1:
+		raise NomogramEdgeCase(-10)
+	return segments
+
+
+def calculate_swi(x, y):
+	try:
+		low, high = find_boundary_curves(_dataset.INDICES, x, y)
+		vec = _dataset.Vector(x, y)
+		dist_to_low = min(abs(vec - p) for p in (low[1][low[2]], low[1][low[2]]))
+		dist_to_high = min(abs(vec - p) for p in (high[1][high[2]], high[1][high[2]]))
+		return dist_to_low / (dist_to_low + dist_to_high) * (high[0] - low[0]) + low[0]
+	except NomogramEdgeCase as edge:
+		return edge.value
+	except IndexError:
+		return -10  # Reduce the amount of cases in which it may occur