From 4f9f987099224ca41aab91a582b80fca8eb969f8 Mon Sep 17 00:00:00 2001
From: Lynn <62256001+solis-lumine-vorago@users.noreply.github.com>
Date: Mon, 6 May 2024 23:31:51 +0200
Subject: [PATCH] Ellipse functions (#13025)

# Objective

- Add some useful methods to `Ellipse`

## Solution

- Added `Ellipse::perimeter()` and `::focal_length()`

---------

Co-authored-by: IQuick 143 <IQuick143cz@gmail.com>
---
 crates/bevy_math/src/primitives/dim2.rs | 82 +++++++++++++++++++++++++
 1 file changed, 82 insertions(+)

diff --git a/crates/bevy_math/src/primitives/dim2.rs b/crates/bevy_math/src/primitives/dim2.rs
index 14660edc5b03b..e3c613a9fae5c 100644
--- a/crates/bevy_math/src/primitives/dim2.rs
+++ b/crates/bevy_math/src/primitives/dim2.rs
@@ -118,6 +118,73 @@ impl Ellipse {
         (a * a - b * b).sqrt() / a
     }
 
+    #[inline(always)]
+    /// Get the focal length of the ellipse. This corresponds to the distance between one of the foci and the center of the ellipse.
+    ///
+    /// The focal length of an ellipse is related to its eccentricity by `eccentricity = focal_length / semi_major`
+    pub fn focal_length(&self) -> f32 {
+        let a = self.semi_major();
+        let b = self.semi_minor();
+
+        (a * a - b * b).sqrt()
+    }
+
+    #[inline(always)]
+    /// Get an approximation for the perimeter or circumference of the ellipse.
+    ///
+    /// The approximation is reasonably precise with a relative error less than 0.007%, getting more precise as the eccentricity of the ellipse decreases.
+    pub fn perimeter(&self) -> f32 {
+        let a = self.semi_major();
+        let b = self.semi_minor();
+
+        // In the case that `a == b`, the ellipse is a circle
+        if a / b - 1. < 1e-5 {
+            return PI * (a + b);
+        };
+
+        // In the case that `a` is much larger than `b`, the ellipse is a line
+        if a / b > 1e4 {
+            return 4. * a;
+        };
+
+        // These values are  the result of (0.5 choose n)^2 where n is the index in the array
+        // They could be calculated on the fly but hardcoding them yields more accurate and faster results
+        // because the actual calculation for these values involves factorials and numbers > 10^23
+        const BINOMIAL_COEFFICIENTS: [f32; 21] = [
+            1.,
+            0.25,
+            0.015625,
+            0.00390625,
+            0.0015258789,
+            0.00074768066,
+            0.00042057037,
+            0.00025963783,
+            0.00017140154,
+            0.000119028846,
+            0.00008599834,
+            0.00006414339,
+            0.000049109784,
+            0.000038430585,
+            0.000030636627,
+            0.000024815668,
+            0.000020380836,
+            0.000016942893,
+            0.000014236736,
+            0.000012077564,
+            0.000010333865,
+        ];
+
+        // The algorithm used here is the Gauss-Kummer infinite series expansion of the elliptic integral expression for the perimeter of ellipses
+        // For more information see https://www.wolframalpha.com/input/?i=gauss-kummer+series
+        // We only use the terms up to `i == 20` for this approximation
+        let h = ((a - b) / (a + b)).powi(2);
+
+        PI * (a + b)
+            * (0..=20)
+                .map(|i| BINOMIAL_COEFFICIENTS[i] * h.powi(i as i32))
+                .sum::<f32>()
+    }
+
     /// Returns the length of the semi-major axis. This corresponds to the longest radius of the ellipse.
     #[inline(always)]
     pub fn semi_major(&self) -> f32 {
@@ -861,6 +928,21 @@ mod tests {
         assert_eq!(circle.eccentricity(), 0., "incorrect circle eccentricity");
     }
 
+    #[test]
+    fn ellipse_perimeter() {
+        let circle = Ellipse::new(1., 1.);
+        assert_relative_eq!(circle.perimeter(), 6.2831855);
+
+        let line = Ellipse::new(75_000., 0.5);
+        assert_relative_eq!(line.perimeter(), 300_000.);
+
+        let ellipse = Ellipse::new(0.5, 2.);
+        assert_relative_eq!(ellipse.perimeter(), 8.578423);
+
+        let ellipse = Ellipse::new(5., 3.);
+        assert_relative_eq!(ellipse.perimeter(), 25.526999);
+    }
+
     #[test]
     fn triangle_math() {
         let triangle = Triangle2d::new(