From 09e6857ae73249bc7433f2971dcf291c70e4c766 Mon Sep 17 00:00:00 2001 From: Mattias Andrée Date: Sun, 20 Oct 2024 17:28:46 +0200 Subject: Fourth commit MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Signed-off-by: Mattias Andrée --- libtellurian_elevated_gravity_radians.c | 36 ++++++++++++++++++++++++++++----- 1 file changed, 31 insertions(+), 5 deletions(-) (limited to 'libtellurian_elevated_gravity_radians.c') diff --git a/libtellurian_elevated_gravity_radians.c b/libtellurian_elevated_gravity_radians.c index ffd1a5b..1c97f1a 100644 --- a/libtellurian_elevated_gravity_radians.c +++ b/libtellurian_elevated_gravity_radians.c @@ -1,14 +1,40 @@ /* See LICENSE file for copyright and license details. */ #include "common.h" +#ifndef TEST double libtellurian_elevated_gravity_radians(double gravity, double latitude, double altitude) { - double k1 = -3.1570429877205807e-07; - double k2 = 2.1026896504579084e-09; - double k3 = -7.374516772941995e-14; + double a = LIBTELLURIAN_EQUATORIAL_RADIUS; + double b = LIBTELLURIAN_POLAR_RADIUS; + double omega = LIBTELLURIAN_ANGULAR_VELOCITY; + double GM = LIBTELLURIAN_GEOCENTRIC_GRAVITATIONAL_CONSTANT; + double f = 1.0 - b / a; + double m = (omega * a) * (omega * a) * (b / GM); + double neg_k1 = fma(-2.0, f + m, -2.0) / a; + double k2 = 4.0 * f / a; + double neg_k3 = -3.0 / (a * a); double sin2_phi = sin(latitude) * sin(latitude); - double f = fma(k3, altitude, fma(k2, sin2_phi, k1)); - return fma(f * altitude, gravity, gravity); + double s = fma(neg_k3, altitude, fma(k2, sin2_phi, neg_k1)); + return fma(s * altitude, gravity, gravity); } + + +#else + + +int +main(void) +{ + ASSERT(libtellurian_elevated_gravity_radians(9.2, 0, 0) == 9.2); + ASSERT(libtellurian_elevated_gravity_radians(4.2, 0, 0) == 4.2); + ASSERT(libtellurian_elevated_gravity_radians(8.9, 2.0, 0) == 8.9); + ASSERT(libtellurian_elevated_gravity_radians(8.9, 2.0, 1.0) < 8.9); + ASSERT(libtellurian_elevated_gravity_radians(8.9, 2.0, 2.0) < libtellurian_elevated_gravity_radians(8.9, 2.0, 1.0)); + ASSERT(libtellurian_elevated_gravity_radians(8.9, 3.0, 1.0) < libtellurian_elevated_gravity_radians(8.9, 2.0, 1.0)); + return 0; +} + + +#endif -- cgit v1.2.3-70-g09d2