Fourth commit

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

1 files changed, 45 insertions, 0 deletions

diff --git a/libtellurian_meridian_radius_radians.c b/libtellurian_meridian_radius_radians.c
new file mode 100644
index 0000000..f550289
--- /dev/null
+++ b/libtellurian_meridian_radius_radians.c
@@ -0,0 +1,45 @@
+/* See LICENSE file for copyright and license details. */
+#include "common.h"
+#ifndef TEST
+
+
+double
+libtellurian_meridian_radius_radians(double latitude)
+{
+	double a = LIBTELLURIAN_EQUATORIAL_RADIUS;
+	double b = LIBTELLURIAN_POLAR_RADIUS;
+	double neg_e2 = fma(b / a, b / a, -1.0);
+	double s = sin(latitude);
+	return fma(a, neg_e2, a) / pow(fma(neg_e2, s * s, 1.0), 1.5);
+}
+
+
+#else
+
+
+static int
+approx(double a, double b)
+{
+        return fabs(a - b) <= 1e-8 * (0.5 * (a + b));
+}
+
+int
+main(void)
+{
+	double x = 0.9933056200098025 * LIBTELLURIAN_EQUATORIAL_RADIUS;
+	ASSERT(fabs(libtellurian_meridian_radius_radians(0) - x) / x < 1.0e-8);
+	ASSERT(approx(libtellurian_meridian_radius_radians(D30), libtellurian_meridian_radius_radians(-D30)));
+	ASSERT(approx(libtellurian_meridian_radius_radians(D45), libtellurian_meridian_radius_radians(-D45)));
+	ASSERT(approx(libtellurian_meridian_radius_radians(D60), libtellurian_meridian_radius_radians(-D60)));
+	ASSERT(approx(libtellurian_meridian_radius_radians(D90), libtellurian_meridian_radius_radians(-D90)));
+	ASSERT(approx(libtellurian_meridian_radius_radians(D180), libtellurian_meridian_radius_radians(0)));
+	ASSERT(approx(libtellurian_meridian_radius_radians(-D180), libtellurian_meridian_radius_radians(0)));
+	ASSERT(approx(libtellurian_meridian_radius_radians(D90 + D180), libtellurian_meridian_radius_radians(D90)));
+	ASSERT(approx(libtellurian_meridian_radius_radians(-D90), libtellurian_meridian_radius_radians(D90)));
+	ASSERT(fabs(libtellurian_meridian_radius_radians(D90) - 6399593.625758674) / 6399593.625758674 < 1.0e-8);
+	ASSERT(fabs(libtellurian_meridian_radius_radians(D30) - 6351377.103715289) / 6351377.103715289 < 1.0e-8);
+	return 0;
+}
+
+
+#endif