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_meridian_radius_radians.c | 45 ++++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) create mode 100644 libtellurian_meridian_radius_radians.c (limited to 'libtellurian_meridian_radius_radians.c') 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 -- cgit v1.2.3-70-g09d2