1 files changed, 51 insertions, 0 deletions
diff --git a/libtellurian_effective_gravity_radians.c b/libtellurian_effective_gravity_radians.c
new file mode 100644
index 0000000..54e2b8e
--- /dev/null
+++ b/libtellurian_effective_gravity_radians.c
@@ -0,0 +1,51 @@
+/* See LICENSE file for copyright and license details. */
+#include "common.h"
+
+#define K (LIBTELLURIAN_EQUATORIAL_RADIUS * LIBTELLURIAN_ANGULAR_VELOCITY * LIBTELLURIAN_ANGULAR_VELOCITY)
+
+#ifndef TEST
+
+
+double
+libtellurian_effective_gravity_radians(double gravity, double latitude)
+{
+	return fma(-K / gravity, cos(latitude) * cos(latitude), gravity);
+}
+
+
+#else
+
+
+static int
+approx(double a, double b)
+{
+        return fabs(a - b) <= 1e-8 * (0.5 * (a + b));
+}
+
+int
+main(void)
+{
+	ASSERT(libtellurian_effective_gravity_radians(1, D90) == 1);
+	ASSERT(libtellurian_effective_gravity_radians(2, D90) == 2);
+	ASSERT(libtellurian_effective_gravity_radians(4, D90) == 4);
+	ASSERT(libtellurian_effective_gravity_radians(1, -D90) == 1);
+	ASSERT(libtellurian_effective_gravity_radians(2, -D90) == 2);
+	ASSERT(libtellurian_effective_gravity_radians(4, -D90) == 4);
+	ASSERT(libtellurian_effective_gravity_radians(1, 0) == fma(1.00, -K, 1.0));
+	ASSERT(libtellurian_effective_gravity_radians(2, 0) == fma(0.50, -K, 2.0));
+	ASSERT(libtellurian_effective_gravity_radians(4, 0) == fma(0.25, -K, 4.0));
+	ASSERT(libtellurian_effective_gravity_radians(1, D45) == fma(0.50, -K, 1.0));
+	ASSERT(libtellurian_effective_gravity_radians(2, D45) == fma(0.25, -K, 2.0));
+	ASSERT(libtellurian_effective_gravity_radians(1, D60) == fma(0.250, -K, 1.0));
+	ASSERT(libtellurian_effective_gravity_radians(2, D60) == fma(0.125, -K, 2.0));
+	ASSERT(libtellurian_effective_gravity_radians(4, D30) < 4);
+	ASSERT(libtellurian_effective_gravity_radians(4, D30) > libtellurian_effective_gravity_radians(4, 0));
+	ASSERT(approx(libtellurian_effective_gravity_radians(LIBTELLURIAN_NORMAL_POLAR_GRAVITY, D90),
+	              LIBTELLURIAN_POLAR_GRAVITY));
+	ASSERT(approx(libtellurian_effective_gravity_radians(LIBTELLURIAN_NORMAL_EQUATORIAL_GRAVITY, 0),
+	              LIBTELLURIAN_EQUATORIAL_GRAVITY));
+	return 0;
+}
+
+
+#endif