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_gaussian_radius_radians.c | 45 +++++++++++++++++++++++++++++++--- 1 file changed, 42 insertions(+), 3 deletions(-) (limited to 'libtellurian_gaussian_radius_radians.c') diff --git a/libtellurian_gaussian_radius_radians.c b/libtellurian_gaussian_radius_radians.c index d61319a..02d0586 100644 --- a/libtellurian_gaussian_radius_radians.c +++ b/libtellurian_gaussian_radius_radians.c @@ -1,5 +1,6 @@ /* See LICENSE file for copyright and license details. */ #include "common.h" +#ifndef TEST double @@ -7,10 +8,48 @@ libtellurian_gaussian_radius_radians(double latitude) { double a = LIBTELLURIAN_EQUATORIAL_RADIUS; double b = LIBTELLURIAN_POLAR_RADIUS; - double f = 1.0 - b / a; - double neg_e2 = (f - 2.0) * f; - double neg_e2_plus_1 = fma(f - 2.0, f, 1.0); + double c = b / a; + double neg_e2 = fma(c, c, -1.0); + double neg_e2_plus_1 = c * c; double s = sin(latitude); double denom = fma(neg_e2, s * s, 1.0); return a * sqrt(neg_e2_plus_1) / denom; } + + +#else + + +static int +check(double lat) +{ + double g = libtellurian_gaussian_radius_radians(lat); + double m = libtellurian_meridian_radius_radians(lat); + double n = libtellurian_transverse_radius_radians(lat); + double e = sqrt(m * n); + double err = fabs((g - e) / e); + return err <= 1.0e-8; +} + +int +main(void) +{ + ASSERT(check(0)); + ASSERT(check(D30)); + ASSERT(check(D45)); + ASSERT(check(D60)); + ASSERT(check(D90)); + ASSERT(check(D180)); + ASSERT(check(-D30)); + ASSERT(check(-D45)); + ASSERT(check(-D60)); + ASSERT(check(-D90)); + ASSERT(check(-D180)); + ASSERT(check(1.0)); + ASSERT(check(2.0)); + ASSERT(check(3.0)); + return 0; +} + + +#endif -- cgit v1.2.3-70-g09d2