From b872e19dec0a7c4f3801ef9e35ee45fae7dc54ab Mon Sep 17 00:00:00 2001
From: Jan Petykiewicz <jan@mpxd.net>
Date: Mon, 17 Jul 2023 21:23:44 -0700
Subject: [PATCH] Improve arclength calculation for elliptical arcs

---
 masque/shapes/arc.py | 17 +++++++++++------
 1 file changed, 11 insertions(+), 6 deletions(-)

diff --git a/masque/shapes/arc.py b/masque/shapes/arc.py
index 17ed3fe..c023e73 100644
--- a/masque/shapes/arc.py
+++ b/masque/shapes/arc.py
@@ -202,13 +202,18 @@ class Arc(Shape):
         # Convert from polar angle to ellipse parameter (for [rx*cos(t), ry*sin(t)] representation)
         a_ranges = self._angles_to_parameters()
 
-        # Approximate perimeter
-        # Ramanujan, S., "Modular Equations and Approximations to ,"
-        #  Quart. J. Pure. Appl. Math., vol. 45 (1913-1914), pp. 350-372
+        # Approximate outer perimeter via numerical integration
         a0, a1 = a_ranges[1]    # use outer arc
-        h = ((r1 - r0) / (r1 + r0)) ** 2
-        ellipse_perimeter = pi * (r1 + r0) * (1 + 3 * h / (10 + math.sqrt(4 - 3 * h)))
-        perimeter = abs(a0 - a1) / (2 * pi) * ellipse_perimeter         # TODO: make this more accurate
+        t = numpy.linspace(a0, a1, 10_000)      # NOTE: could probably use an adaptive number of points
+        r0sin = r0 * numpy.sin(t)
+        r1cos = r1 * numpy.cos(t)
+        perimeter = numpy.trapz(numpy.sqrt(r0sin * r0sin + r1cos * r1cos), dx=t[1] - t[0])
+
+        #from scipy.special import ellipeinc
+        #m = 1 - (r1 / r0) ** 2
+        #t1 = ellipeinc(a1 - pi / 2, m)
+        #t0 = ellipeinc(a0 - pi / 2, m)
+        #perimeter2 = r0 * (t1 - t0)
 
         n = []
         if num_vertices is not None: