Fourth commit

Signed-off-by: Mattias Andrée <m@maandree.se>

1 files changed, 31 insertions, 5 deletions

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