diff --git a/masque/utils/curves.py b/masque/utils/curves.py
index 8b3fcc4..2348678 100644
--- a/masque/utils/curves.py
+++ b/masque/utils/curves.py
@@ -69,14 +69,25 @@ def euler_bend(
     num_points_arc = num_points - 2 * num_points_spiral
 
     def gen_spiral(ll_max: float) -> NDArray[numpy.float64]:
-        xx = []
-        yy = []
-        for ll in numpy.linspace(0, ll_max, num_points_spiral):
-            qq = numpy.linspace(0, ll, 1000)        # integrate to current arclength
-            xx.append(trapezoid( numpy.cos(qq * qq / 2), qq))
-            yy.append(trapezoid(-numpy.sin(qq * qq / 2), qq))
-        xy_part = numpy.stack((xx, yy), axis=1)
-        return xy_part
+        if ll_max == 0:
+            return numpy.zeros((num_points_spiral, 2))
+
+        resolution = 100000
+        qq = numpy.linspace(0, ll_max, resolution)
+        dx = numpy.cos(qq * qq / 2)
+        dy = -numpy.sin(qq * qq / 2)
+
+        dq = ll_max / (resolution - 1)
+        ix = numpy.zeros(resolution)
+        iy = numpy.zeros(resolution)
+        ix[1:] = numpy.cumsum((dx[:-1] + dx[1:]) / 2) * dq
+        iy[1:] = numpy.cumsum((dy[:-1] + dy[1:]) / 2) * dq
+
+        ll_target = numpy.linspace(0, ll_max, num_points_spiral)
+        x_target = numpy.interp(ll_target, qq, ix)
+        y_target = numpy.interp(ll_target, qq, iy)
+
+        return numpy.stack((x_target, y_target), axis=1)
 
     xy_spiral = gen_spiral(ll_max)
     xy_parts = [xy_spiral]