diff --git a/examples/orbit_viewer.rs b/examples/orbit_viewer.rs
index 62f4535..6005651 100644
--- a/examples/orbit_viewer.rs
+++ b/examples/orbit_viewer.rs
@@ -5,6 +5,7 @@ use bevy::pbr::NotShadowCaster;
 use bevy::prelude::*;
 use bevy_egui::egui::{ComboBox, DragValue, Ui};
 use bevy_egui::{egui, EguiContexts, EguiPlugin};
+use keplerian_elements::astro::standard_gravitational_parameter;
 use keplerian_elements::constants::AU;
 use keplerian_elements::utils::{yup2zup, zup2yup};
 use keplerian_elements::{KeplerianElements, StateVectors};
@@ -99,7 +100,11 @@ fn ui(
                             f32::MIN,
                             f32::MAX,
                         );
-                        value_slider(ui, "Eccentricity", &mut orbit.eccentricity);
+                        value_slider(
+                            ui,
+                            "Eccentricity",
+                            &mut orbit.eccentricity,
+                        );
                         value_slider(ui, "Inclination", &mut orbit.inclination);
                         value_slider(
                             ui,
@@ -111,7 +116,11 @@ fn ui(
                             "Argument of periapsis",
                             &mut orbit.argument_of_periapsis,
                         );
-                        value_slider(ui, "Mean anomaly", &mut orbit.mean_anomaly_at_epoch);
+                        value_slider(
+                            ui,
+                            "Mean anomaly",
+                            &mut orbit.mean_anomaly_at_epoch,
+                        );
                         value_slider(ui, "Epoch", &mut orbit.epoch);
 
                         let v = orbit.true_anomaly_at_epoch(
@@ -123,9 +132,12 @@ fn ui(
                         ui.label(format!("True anomaly: {v}",));
 
                         let e = orbit.eccentricity;
-                        let h = orbit.specific_angular_momentum(state.star_mass);
+                        let h =
+                            orbit.specific_angular_momentum(state.star_mass);
                         let μ =
-                            KeplerianElements::standard_gravitational_parameter(state.star_mass);
+                            standard_gravitational_parameter(
+                                state.star_mass,
+                            );
 
                         let r = (h.powi(2) / μ) / (1.0 + e * v.cos());
                         ui.label(format!("r: {r}"));
@@ -139,11 +151,12 @@ fn ui(
                         ui.label(format!("vp: {vp}",));
                         ui.label(format!("vq: {vq}",));
 
-                        planet.state_vectors = planet.orbit.state_vectors_at_epoch(
-                            state.star_mass,
-                            state.epoch,
-                            state.tolerance,
-                        );
+                        planet.state_vectors =
+                            planet.orbit.state_vectors_at_epoch(
+                                state.star_mass,
+                                state.epoch,
+                                state.tolerance,
+                            );
                     });
 
                     // --- State Vectors ---
@@ -151,7 +164,8 @@ fn ui(
                         let sv = &mut planet.state_vectors;
                         ui.label("Position");
 
-                        let mut p = zup2yup(sv.position * state.distance_scaling);
+                        let mut p =
+                            zup2yup(sv.position * state.distance_scaling);
 
                         value_slider(ui, "X", &mut p.x);
                         value_slider(ui, "Y", &mut p.y);
@@ -161,24 +175,36 @@ fn ui(
 
                         ui.label("Velocity");
 
-                        let mut v =
-                            zup2yup(sv.velocity * state.distance_scaling * state.velocity_scaling);
+                        let mut v = zup2yup(
+                            sv.velocity
+                                * state.distance_scaling
+                                * state.velocity_scaling,
+                        );
 
                         value_slider(ui, "Vx", &mut v.x);
                         value_slider(ui, "Vy", &mut v.y);
                         value_slider(ui, "Vz", &mut v.z);
 
-                        sv.velocity =
-                            yup2zup(v / (state.distance_scaling * state.velocity_scaling));
+                        sv.velocity = yup2zup(
+                            v / (state.distance_scaling
+                                * state.velocity_scaling),
+                        );
 
-                        planet.orbit = sv.to_elements(state.star_mass, state.epoch);
+                        planet.orbit =
+                            sv.to_elements(state.star_mass, state.epoch);
                     });
                 });
             }
         });
 
         ui.collapsing("State", |ui| {
-            value_slider_min_max(ui, "Tolerance", &mut state.tolerance, f32::EPSILON, 100.0);
+            value_slider_min_max(
+                ui,
+                "Tolerance",
+                &mut state.tolerance,
+                f32::EPSILON,
+                100.0,
+            );
             value_slider(ui, "Mass", &mut state.star_mass);
             value_slider(ui, "Epoch", &mut state.epoch);
             value_slider(ui, "Epoch scale", &mut state.epoch_scale);
@@ -197,7 +223,11 @@ fn ui(
                     "Show position & velocity",
                 );
 
-                value_slider_u32(ui, "Orbit subdivisions", &mut state.orbit_subdivisions);
+                value_slider_u32(
+                    ui,
+                    "Orbit subdivisions",
+                    &mut state.orbit_subdivisions,
+                );
             }
 
             ui.checkbox(&mut state.draw_soi, "Draw SOI");
@@ -224,17 +254,33 @@ fn ui(
                 ui.label("Mouse rotate sensitivity");
                 ui.horizontal(|ui| {
                     ui.label("x");
-                    ui.add(DragValue::new(&mut camera.mouse_rotate_sensitivity.x).speed(0.01));
+                    ui.add(
+                        DragValue::new(&mut camera.mouse_rotate_sensitivity.x)
+                            .speed(0.01),
+                    );
                     ui.label("y");
-                    ui.add(DragValue::new(&mut camera.mouse_rotate_sensitivity.y).speed(0.01));
+                    ui.add(
+                        DragValue::new(&mut camera.mouse_rotate_sensitivity.y)
+                            .speed(0.01),
+                    );
                 });
 
                 ui.label("Mouse translate sensitivity");
                 ui.horizontal(|ui| {
                     ui.label("x");
-                    ui.add(DragValue::new(&mut camera.mouse_translate_sensitivity.x).speed(0.01));
+                    ui.add(
+                        DragValue::new(
+                            &mut camera.mouse_translate_sensitivity.x,
+                        )
+                        .speed(0.01),
+                    );
                     ui.label("y");
-                    ui.add(DragValue::new(&mut camera.mouse_translate_sensitivity.y).speed(0.01));
+                    ui.add(
+                        DragValue::new(
+                            &mut camera.mouse_translate_sensitivity.y,
+                        )
+                        .speed(0.01),
+                    );
                 });
             });
         }
@@ -273,7 +319,10 @@ fn ui(
 
                 for (_, name) in &planets {
                     if ui
-                        .selectable_label(current == name.to_string(), &name.to_string())
+                        .selectable_label(
+                            current == name.to_string(),
+                            &name.to_string(),
+                        )
                         .clicked()
                     {
                         state.focus_mode = FocusMode::Planet(name.to_string());
@@ -287,7 +336,13 @@ fn value_slider(ui: &mut Ui, name: &str, value: &mut f32) {
     value_slider_min_max(ui, name, value, f32::MIN, f32::MAX)
 }
 
-fn value_slider_min_max(ui: &mut Ui, name: &str, value: &mut f32, min: f32, max: f32) {
+fn value_slider_min_max(
+    ui: &mut Ui,
+    name: &str,
+    value: &mut f32,
+    min: f32,
+    max: f32,
+) {
     value_slider_min_max_with_speed(ui, name, value, min, max, 0.01);
 }
 
@@ -531,7 +586,8 @@ fn spawn_solar_system(
                 eccentricity: 0.09339410,
                 semi_major_axis: 1.52371034 * AU,
                 inclination: 0.03232349774693498376462675303241,
-                right_ascension_of_the_ascending_node: 0.86760317116638123268876668101569,
+                right_ascension_of_the_ascending_node:
+                    0.86760317116638123268876668101569,
                 argument_of_periapsis: 5.8657025501025428421251399347365,
                 mean_anomaly_at_epoch: 6.2034237603634456152598740984391,
                 epoch: 0.0,
@@ -552,7 +608,8 @@ fn spawn_solar_system(
                 eccentricity: 0.04854,
                 semi_major_axis: 5.2025 * AU,
                 inclination: 0.02267182698340634120423874308267,
-                right_ascension_of_the_ascending_node: 1.7503907068251131326967694717172,
+                right_ascension_of_the_ascending_node:
+                    1.7503907068251131326967694717172,
                 argument_of_periapsis: 0.24905848425959083062701067266333,
                 mean_anomaly_at_epoch: 0.59917153220965334375790304082214,
                 epoch: 0.0,
@@ -573,7 +630,8 @@ fn spawn_solar_system(
                 eccentricity: 0.05551,
                 semi_major_axis: 9.5415 * AU,
                 inclination: 0.04352851154473857964847684776611,
-                right_ascension_of_the_ascending_node: 1.9833921619663561312160821893105,
+                right_ascension_of_the_ascending_node:
+                    1.9833921619663561312160821893105,
                 argument_of_periapsis: 1.6207127434019344451313392476185,
                 mean_anomaly_at_epoch: 0.8740608893987602521233843368591,
                 epoch: 0.0,
@@ -594,7 +652,8 @@ fn spawn_solar_system(
                 eccentricity: 0.04686,
                 semi_major_axis: 19.188 * AU,
                 inclination: 0.01349139511791616762962012964042,
-                right_ascension_of_the_ascending_node: 1.2908455147750061550927616923742,
+                right_ascension_of_the_ascending_node:
+                    1.2908455147750061550927616923742,
                 argument_of_periapsis: 3.0094712292138224894895199921049,
                 mean_anomaly_at_epoch: 5.4838245097661835306942363945912,
                 epoch: 0.0,
@@ -615,7 +674,8 @@ fn spawn_solar_system(
                 eccentricity: 0.00895,
                 semi_major_axis: 30.070 * AU,
                 inclination: 0.03089232776029963351154932660225,
-                right_ascension_of_the_ascending_node: 2.3001694212033269494277320637911,
+                right_ascension_of_the_ascending_node:
+                    2.3001694212033269494277320637911,
                 argument_of_periapsis: 0.81471969483095304650797885073048,
                 mean_anomaly_at_epoch: 5.3096406504171494389172520558961,
                 epoch: 0.0,
@@ -632,12 +692,16 @@ fn update_epoch(time: Res<Time>, mut state: ResMut<State>) {
     }
 }
 
-fn update_planets(mut query: Query<(&mut Transform, &mut Planet)>, state: Res<State>) {
+fn update_planets(
+    mut query: Query<(&mut Transform, &mut Planet)>,
+    state: Res<State>,
+) {
     for (mut transform, mut planet) in query.iter_mut() {
-        planet.state_vectors =
-            planet
-                .orbit
-                .state_vectors_at_epoch(state.star_mass, state.epoch, state.tolerance);
+        planet.state_vectors = planet.orbit.state_vectors_at_epoch(
+            state.star_mass,
+            state.epoch,
+            state.tolerance,
+        );
 
         let position = zup2yup(planet.state_vectors.position);
 
@@ -646,9 +710,13 @@ fn update_planets(mut query: Query<(&mut Transform, &mut Planet)>, state: Res<St
     }
 }
 
-fn update_star(mut query: Query<&mut Transform, With<Star>>, state: Res<State>) {
+fn update_star(
+    mut query: Query<&mut Transform, With<Star>>,
+    state: Res<State>,
+) {
     for mut transform in query.iter_mut() {
-        transform.scale = Vec3::ONE * mass2radius(state.as_ref(), state.star_mass);
+        transform.scale =
+            Vec3::ONE * mass2radius(state.as_ref(), state.star_mass);
     }
 }
 
@@ -692,7 +760,8 @@ fn draw_orbits(
         let color = materials.get(mat).unwrap().base_color;
 
         let first_position =
-            zup2yup(orbit.position_at_true_anomaly(state.star_mass, 0.0)) * state.distance_scaling;
+            zup2yup(orbit.position_at_true_anomaly(state.star_mass, 0.0))
+                * state.distance_scaling;
         let mut prev_position = first_position.clone();
 
         let step = (2.0 * PI) / state.orbit_subdivisions as f32;
@@ -714,8 +783,12 @@ fn draw_orbits(
         let mut debug_arrows = DebugArrows::new(&mut lines, camera_position);
 
         if state.show_position_and_velocity {
-            let StateVectors { position, velocity } =
-                orbit.state_vectors_at_epoch(state.star_mass, state.epoch, state.tolerance);
+            let StateVectors { position, velocity } = orbit
+                .state_vectors_at_epoch(
+                    state.star_mass,
+                    state.epoch,
+                    state.tolerance,
+                );
 
             let position = zup2yup(position);
             let velocity = zup2yup(velocity);
@@ -730,12 +803,14 @@ fn draw_orbits(
         if state.show_nodes {
             debug_arrows.draw_arrow(
                 Vec3::ZERO,
-                zup2yup(orbit.ascending_node(state.star_mass)) * state.distance_scaling,
+                zup2yup(orbit.ascending_node(state.star_mass))
+                    * state.distance_scaling,
                 Color::YELLOW_GREEN,
             );
             debug_arrows.draw_arrow(
                 Vec3::ZERO,
-                zup2yup(orbit.descending_node(state.star_mass)) * state.distance_scaling,
+                zup2yup(orbit.descending_node(state.star_mass))
+                    * state.distance_scaling,
                 Color::YELLOW,
             );
         }
@@ -743,12 +818,14 @@ fn draw_orbits(
         if state.show_peri_and_apo_apsis {
             debug_arrows.draw_arrow(
                 Vec3::ZERO,
-                zup2yup(orbit.periapsis(state.star_mass)) * state.distance_scaling,
+                zup2yup(orbit.periapsis(state.star_mass))
+                    * state.distance_scaling,
                 Color::WHITE,
             );
             debug_arrows.draw_arrow(
                 Vec3::ZERO,
-                zup2yup(orbit.apoapsis(state.star_mass)) * state.distance_scaling,
+                zup2yup(orbit.apoapsis(state.star_mass))
+                    * state.distance_scaling,
                 Color::WHITE,
             );
         }
@@ -772,10 +849,12 @@ fn draw_soi(
     for planet in planets.iter() {
         let r = planet.state_vectors.position.length();
 
-        let soi = keplerian_elements::astro::soi(r, planet.mass, state.star_mass)
-            * state.distance_scaling;
+        let soi =
+            keplerian_elements::astro::soi(r, planet.mass, state.star_mass)
+                * state.distance_scaling;
 
-        let pos = zup2yup(planet.state_vectors.position) * state.distance_scaling;
+        let pos =
+            zup2yup(planet.state_vectors.position) * state.distance_scaling;
 
         let to_camera = (camera_position - pos).normalize();
         let planet_camera_radial = to_camera.cross(pos).normalize();
diff --git a/rustfmt.toml b/rustfmt.toml
new file mode 100644
index 0000000..6156297
--- /dev/null
+++ b/rustfmt.toml
@@ -0,0 +1,3 @@
+max_width = 80
+imports_granularity = "Module"
+group_imports = "StdExternalCrate"
diff --git a/src/astro.rs b/src/astro.rs
index c3e4478..f7199cd 100644
--- a/src/astro.rs
+++ b/src/astro.rs
@@ -1,10 +1,15 @@
 use crate::constants::G;
 use crate::Num;
 
-pub fn soi(r: Num, m1: Num, m2: Num) -> Num {
-    r * (m1 / m2).powf(2.0 / 5.0)
-}
+pub mod elliptic;
+pub mod hyperbolic;
 
+/// https://en.wikipedia.org/wiki/Standard_gravitational_parameter
+#[inline]
 pub fn standard_gravitational_parameter(mass: Num) -> Num {
     G * mass
 }
+
+pub fn soi(r: Num, m1: Num, m2: Num) -> Num {
+    r * (m1 / m2).powf(2.0 / 5.0)
+}
diff --git a/src/astro/elliptic.rs b/src/astro/elliptic.rs
new file mode 100644
index 0000000..74f60d4
--- /dev/null
+++ b/src/astro/elliptic.rs
@@ -0,0 +1,42 @@
+use crate::math::newton_approx;
+use crate::Num;
+
+use super::standard_gravitational_parameter;
+
+/// Eccentric Anomaly (E) is given by the equation:
+/// M = E - e * sin(E)
+/// where
+/// M is the mean anomaly
+/// e is the eccentricity
+///
+/// https://orbital-mechanics.space/time-since-periapsis-and-keplers-equation/elliptical-orbits.html#equation-eq-keplers-equation-ellipse
+pub fn estimate_anomaly(
+    // Mean anomaly
+    M: Num,
+    // Eccentricity
+    e: Num,
+    tolerance: Num,
+) -> Num {
+    newton_approx(
+        // f(E) = E - e*sin(E) - M
+        |E| E - (e * E.sin()) - M,
+        // f'(E) = 1 - e*cos(E)
+        |E| 1.0 - (e * E.cos()),
+        M,
+        tolerance,
+    )
+}
+
+/// Mean motion
+/// https://en.wikipedia.org/wiki/Mean_anomaly
+pub fn mean_motion(h: Num, e: Num, mass: Num) -> Num {
+    let μ = standard_gravitational_parameter(mass);
+
+    (μ.powi(2) / h.powi(3)) * (1.0 - e.powi(2)).powi(3).sqrt()
+}
+
+pub fn true_anomaly(E: Num, e: Num) -> Num {
+    // Circular (practically unattainable), elliptic or parabolic (practically unattainable)
+    // https://orbital-mechanics.space/time-since-periapsis-and-keplers-equation/elliptical-orbits.html#equation-eq-eccentric-anomaly-true-anomaly-ellipse
+    2.0 * ((E / 2.0).tan() / ((1.0 - e) / (1.0 + e)).sqrt()).atan()
+}
diff --git a/src/astro/hyperbolic.rs b/src/astro/hyperbolic.rs
new file mode 100644
index 0000000..3cff3a2
--- /dev/null
+++ b/src/astro/hyperbolic.rs
@@ -0,0 +1,34 @@
+use super::standard_gravitational_parameter;
+use crate::math::newton_approx;
+use crate::Num;
+
+/// Hyperbolic Anomaly (F) is given by the equation:
+/// M = e * sinh(F) - F
+/// where
+/// M is the hyperbolic mean anomaly
+/// e is the eccentricity
+///
+/// https://orbital-mechanics.space/time-since-periapsis-and-keplers-equation/hyperbolic-trajectories.html#equation-eq-hyperbolic-keplers-equation
+pub fn estimate_anomaly(M: Num, e: Num, tolerance: Num) -> Num {
+    newton_approx(
+        // f(F) = e * sinh(F) - F - M
+        |F| (e * F.sinh()) - F - M,
+        // f'(F) = e * cosh(F) - 1
+        |F| e * F.cosh() - 1.0,
+        M,
+        tolerance,
+    )
+}
+
+/// Hyperbolic mean motion
+/// SRC: https://orbital-mechanics.space/time-since-periapsis-and-keplers-equation/hyperbolic-trajectories.html#equation-eq-hyperbolic-mean-anomaly
+pub fn mean_motion(h: Num, e: Num, mass: Num) -> Num {
+    let μ = standard_gravitational_parameter(mass);
+
+    (μ.powi(2) / h.powi(3)) * (e.powi(2) - 1.0).powi(3).sqrt()
+}
+
+pub fn true_anomaly(F: Num, e: Num) -> Num {
+    // https://orbital-mechanics.space/time-since-periapsis-and-keplers-equation/hyperbolic-trajectories.html#equation-eq-eccentric-anomaly-true-anomaly-hyperbola
+    2.0 * ((F / 2.0).tanh() / ((e - 1.0) / (e + 1.0)).sqrt()).atan()
+}
diff --git a/src/elements.rs b/src/elements.rs
index 3966b81..b1bb8a8 100644
--- a/src/elements.rs
+++ b/src/elements.rs
@@ -1,5 +1,5 @@
-use crate::math::newton_approx;
-use crate::{vec3, Mat3, Num, StateVectors, Vec3, G, PI, TWO_PI};
+use crate::astro::{self, standard_gravitational_parameter};
+use crate::{vec3, Mat3, Num, StateVectors, Vec3, PI, TWO_PI};
 
 #[derive(Debug, Clone, Copy)]
 #[cfg_attr(feature = "serde", derive(serde::Serialize, serde::Deserialize))]
@@ -22,13 +22,19 @@ impl KeplerianElements {
         diff += (self.right_ascension_of_the_ascending_node
             - other.right_ascension_of_the_ascending_node)
             .abs();
-        diff += (self.argument_of_periapsis - other.argument_of_periapsis).abs();
-        diff += (self.mean_anomaly_at_epoch - other.mean_anomaly_at_epoch).abs();
+        diff +=
+            (self.argument_of_periapsis - other.argument_of_periapsis).abs();
+        diff +=
+            (self.mean_anomaly_at_epoch - other.mean_anomaly_at_epoch).abs();
 
         diff
     }
 
-    pub fn from_state_vectors(state_vectors: &StateVectors, mass: Num, time: Num) -> Self {
+    pub fn from_state_vectors(
+        state_vectors: &StateVectors,
+        mass: Num,
+        time: Num,
+    ) -> Self {
         state_vectors.to_elements(mass, time)
     }
 
@@ -52,18 +58,13 @@ impl KeplerianElements {
         self.perifocal_to_equatorial(Vec3::Z)
     }
 
-    /// https://en.wikipedia.org/wiki/Standard_gravitational_parameter
-    pub fn standard_gravitational_parameter(mass: Num) -> Num {
-        G * mass
-    }
-
     /// https://en.wikipedia.org/wiki/Orbital_period
     pub fn period(&self, mass: Num) -> Num {
         Self::period_static(self.semi_major_axis, mass)
     }
 
     pub fn period_static(a: Num, mass: Num) -> Num {
-        TWO_PI * (a.powi(3) / Self::standard_gravitational_parameter(mass)).sqrt()
+        TWO_PI * (a.powi(3) / standard_gravitational_parameter(mass)).sqrt()
     }
 
     /// https://en.wikipedia.org/wiki/Mean_anomaly
@@ -73,15 +74,8 @@ impl KeplerianElements {
 
         let epoch_diff = epoch - self.epoch;
 
-        self.mean_anomaly_at_epoch + Self::mean_motion(h, e, mass) * epoch_diff
-    }
-
-    /// Mean motion
-    /// https://en.wikipedia.org/wiki/Mean_anomaly
-    pub fn mean_motion(h: Num, e: Num, mass: Num) -> Num {
-        let μ = Self::standard_gravitational_parameter(mass);
-
-        (μ.powi(2) / h.powi(3)) * (1.0 - e.powi(2)).powi(3).sqrt()
+        self.mean_anomaly_at_epoch
+            + astro::elliptic::mean_motion(h, e, mass) * epoch_diff
     }
 
     /// Hyperbolic mean anomaly
@@ -92,15 +86,8 @@ impl KeplerianElements {
 
         let epoch_diff = epoch - self.epoch;
 
-        self.mean_anomaly_at_epoch + Self::hyperbolic_mean_motion(h, e, mass) * epoch_diff
-    }
-
-    /// Hyperbolic mean motion
-    /// SRC: https://orbital-mechanics.space/time-since-periapsis-and-keplers-equation/hyperbolic-trajectories.html#equation-eq-hyperbolic-mean-anomaly
-    pub fn hyperbolic_mean_motion(h: Num, e: Num, mass: Num) -> Num {
-        let μ = Self::standard_gravitational_parameter(mass);
-
-        (μ.powi(2) / h.powi(3)) * (e.powi(2) - 1.0).powi(3).sqrt()
+        self.mean_anomaly_at_epoch
+            + astro::hyperbolic::mean_motion(h, e, mass) * epoch_diff
     }
 
     /// Eccentric Anomaly (E) is given by the equation:
@@ -110,18 +97,16 @@ impl KeplerianElements {
     /// e is the eccentricity
     ///
     /// https://orbital-mechanics.space/time-since-periapsis-and-keplers-equation/elliptical-orbits.html#equation-eq-keplers-equation-ellipse
-    pub fn estimate_eccentric_anomaly(&self, mass: Num, epoch: Num, tolerance: Num) -> Num {
+    pub fn estimate_eccentric_anomaly(
+        &self,
+        mass: Num,
+        epoch: Num,
+        tolerance: Num,
+    ) -> Num {
         let M = self.mean_anomaly(mass, epoch);
         let e = self.eccentricity;
 
-        newton_approx(
-            // f(E) = E - e*sin(E) - M
-            |E| E - (e * E.sin()) - M,
-            // f'(E) = 1 - e*cos(E)
-            |E| 1.0 - (e * E.cos()),
-            M,
-            tolerance,
-        )
+        astro::elliptic::estimate_anomaly(M, e, tolerance)
     }
 
     /// Hyperbolic Anomaly (F) is given by the equation:
@@ -131,21 +116,24 @@ impl KeplerianElements {
     /// e is the eccentricity
     ///
     /// https://orbital-mechanics.space/time-since-periapsis-and-keplers-equation/hyperbolic-trajectories.html#equation-eq-hyperbolic-keplers-equation
-    pub fn estimate_hyperbolic_anomaly(&self, mass: Num, epoch: Num, tolerance: Num) -> Num {
+    pub fn estimate_hyperbolic_anomaly(
+        &self,
+        mass: Num,
+        epoch: Num,
+        tolerance: Num,
+    ) -> Num {
         let M = self.hyperbolic_mean_anomaly(mass, epoch);
         let e = self.eccentricity;
 
-        newton_approx(
-            // f(F) = e * sinh(F) - F - M
-            |F| (e * F.sinh()) - F - M,
-            // f'(F) = e * cosh(F) - 1
-            |F| e * F.cosh() - 1.0,
-            M,
-            tolerance,
-        )
+        astro::hyperbolic::estimate_anomaly(M, e, tolerance)
     }
 
-    pub fn state_vectors_at_epoch(&self, mass: Num, epoch: Num, tolerance: Num) -> StateVectors {
+    pub fn state_vectors_at_epoch(
+        &self,
+        mass: Num,
+        epoch: Num,
+        tolerance: Num,
+    ) -> StateVectors {
         // Lowercase nu
         let v = self.true_anomaly_at_epoch(mass, epoch, tolerance);
 
@@ -159,7 +147,7 @@ impl KeplerianElements {
     pub fn position_at_true_anomaly(&self, mass: Num, v: Num) -> Vec3 {
         let e = self.eccentricity;
         let h = self.specific_angular_momentum(mass);
-        let μ = Self::standard_gravitational_parameter(mass);
+        let μ = standard_gravitational_parameter(mass);
 
         let r = (h.powi(2) / μ) / (1.0 + e * v.cos());
 
@@ -176,7 +164,7 @@ impl KeplerianElements {
     pub fn velocity_at_true_anomaly(&self, mass: Num, v: Num) -> Vec3 {
         let e = self.eccentricity;
         let h = self.specific_angular_momentum(mass);
-        let μ = Self::standard_gravitational_parameter(mass);
+        let μ = standard_gravitational_parameter(mass);
 
         let vp = -(μ / h) * v.sin();
         let vq = (μ / h) * (e + v.cos());
@@ -200,7 +188,7 @@ impl KeplerianElements {
     }
 
     pub fn specific_angular_momentum(&self, mass: Num) -> Num {
-        let μ = Self::standard_gravitational_parameter(mass);
+        let μ = standard_gravitational_parameter(mass);
         let a = self.semi_major_axis;
         let e = self.eccentricity;
 
@@ -214,20 +202,20 @@ impl KeplerianElements {
     }
 
     /// Calculates true anomaly
-    pub fn true_anomaly_at_epoch(&self, mass: Num, epoch: Num, tolerance: Num) -> Num {
+    pub fn true_anomaly_at_epoch(
+        &self,
+        mass: Num,
+        epoch: Num,
+        tolerance: Num,
+    ) -> Num {
         let e = self.eccentricity;
 
         if self.is_hyperbolic() {
             let F = self.estimate_hyperbolic_anomaly(mass, epoch, tolerance);
-
-            // https://orbital-mechanics.space/time-since-periapsis-and-keplers-equation/hyperbolic-trajectories.html#equation-eq-eccentric-anomaly-true-anomaly-hyperbola
-            2.0 * ((F / 2.0).tanh() / ((e - 1.0) / (e + 1.0)).sqrt()).atan()
+            astro::hyperbolic::true_anomaly(F, e)
         } else {
             let E = self.estimate_eccentric_anomaly(mass, epoch, tolerance);
-
-            // Circular (practically unattainable), elliptic or parabolic (practically unattainable)
-            // https://orbital-mechanics.space/time-since-periapsis-and-keplers-equation/elliptical-orbits.html#equation-eq-eccentric-anomaly-true-anomaly-ellipse
-            2.0 * ((E / 2.0).tan() / ((1.0 - e) / (1.0 + e)).sqrt()).atan()
+            astro::elliptic::true_anomaly(E, e)
         }
     }
 
diff --git a/src/lib.rs b/src/lib.rs
index a9f62e4..bb63546 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -18,7 +18,7 @@ pub mod math;
 pub mod state_vectors;
 pub mod utils;
 
-use constants::{G, PI, TWO_PI};
+use constants::{PI, TWO_PI};
 
 pub use self::elements::KeplerianElements;
 pub use self::state_vectors::StateVectors;
@@ -34,12 +34,17 @@ mod tests {
     const MAX_ABS_DIFF: Num = 0.0001;
     const TOLERANCE: Num = 0.0001;
 
-    fn test_back_and_forth_conversion(original: KeplerianElements, mass: Num, epoch: Num) {
+    fn test_back_and_forth_conversion(
+        original: KeplerianElements,
+        mass: Num,
+        epoch: Num,
+    ) {
         let sv = original.state_vectors_at_epoch(mass, epoch, TOLERANCE);
 
         let elements = KeplerianElements::from_state_vectors(&sv, mass, epoch);
 
-        let sv_converted = elements.state_vectors_at_epoch(mass, epoch, TOLERANCE);
+        let sv_converted =
+            elements.state_vectors_at_epoch(mass, epoch, TOLERANCE);
 
         let pos_diff = sv.position.distance(sv_converted.position);
         assert!(
diff --git a/src/state_vectors.rs b/src/state_vectors.rs
index 7185886..1f6b81a 100644
--- a/src/state_vectors.rs
+++ b/src/state_vectors.rs
@@ -14,7 +14,8 @@ impl StateVectors {
     }
 
     pub fn abs_diff(&self, other: &Self) -> Num {
-        self.position.distance(other.position) + self.velocity.distance(other.velocity)
+        self.position.distance(other.position)
+            + self.velocity.distance(other.velocity)
     }
 
     pub fn to_elements(&self, mass: Num, time: Num) -> KeplerianElements {
@@ -69,8 +70,20 @@ impl StateVectors {
         // Argument of periapsis
         let mut ω = (ev / e).dot(nv).acos();
 
-        if ev.z < 0.0 {
-            ω = TWO_PI - ω;
+        // An edge case for a zero inclination orbit
+        // If an orbit has zero inclination,
+        // the z component of the eccentricity vector
+        // is zero.
+        //
+        // But we can still do a quadrant check using the y component
+        if i.abs() < Num::EPSILON {
+            if ev.y < 0.0 {
+                ω = TWO_PI - ω;
+            }
+        } else {
+            if ev.z < 0.0 {
+                ω = TWO_PI - ω;
+            }
         }
 
         if e == 0.0 {
@@ -92,23 +105,6 @@ impl StateVectors {
             v = TWO_PI - v;
         }
 
-        // Hyperbolic mean anomaly calculation
-        fn calculate_hyperbolic_mean_anomaly(e: Num, v: Num) -> Num {
-            let term1 = (e * (e.powi(2) - 1.0).sqrt() * v.sin()) / (1.0 + e * v.cos());
-            let term2_numerator = (e + 1.0).sqrt() + (e - 1.0).sqrt() * (v / 2.0).tan();
-            let term2_denominator = (e + 1.0).sqrt() - (e - 1.0).sqrt() * (v / 2.0).tan();
-
-            term1 - (term2_numerator / term2_denominator).ln()
-        }
-
-        // Elliptical mean anomaly calculation
-        fn calculate_elliptical_mean_anomaly(e: Num, v: Num) -> Num {
-            let term1 = 2.0 * (((1.0 - e) / (1.0 + e)).sqrt() * (v / 2.0).tan()).atan();
-            let term2 = e * ((1.0 - e.powi(2)).sqrt() * v.sin() / (1.0 + e * v.cos()));
-
-            term1 - term2
-        }
-
         // Mean anomaly calculation
         let M = if is_hyperbolic {
             calculate_hyperbolic_mean_anomaly(e, v)
@@ -127,3 +123,21 @@ impl StateVectors {
         }
     }
 }
+
+// Hyperbolic mean anomaly calculation
+fn calculate_hyperbolic_mean_anomaly(e: Num, v: Num) -> Num {
+    let term1 = (e * (e.powi(2) - 1.0).sqrt() * v.sin()) / (1.0 + e * v.cos());
+    let term2_numerator = (e + 1.0).sqrt() + (e - 1.0).sqrt() * (v / 2.0).tan();
+    let term2_denominator =
+        (e + 1.0).sqrt() - (e - 1.0).sqrt() * (v / 2.0).tan();
+
+    term1 - (term2_numerator / term2_denominator).ln()
+}
+
+// Elliptical mean anomaly calculation
+fn calculate_elliptical_mean_anomaly(e: Num, v: Num) -> Num {
+    let term1 = 2.0 * (((1.0 - e) / (1.0 + e)).sqrt() * (v / 2.0).tan()).atan();
+    let term2 = e * ((1.0 - e.powi(2)).sqrt() * v.sin() / (1.0 + e * v.cos()));
+
+    term1 - term2
+}