blob: c4f80aa4dbe8909e72f6b7b9b2d7fa32a05ec209 (
plain) (
blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
|
/* See LICENSE file for copyright and license details. */
#include "common.h"
#ifndef TEST
double
libtellurian_transverse_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 a * pow(fma(neg_e2, s * s, 1.0), -0.5);
}
#else
static int
approx(double a, double b)
{
return fabs(a - b) <= 1e-8 * (0.5 * (a + b));
}
int
main(void)
{
ASSERT(libtellurian_transverse_radius_radians(0) == LIBTELLURIAN_EQUATORIAL_RADIUS);
ASSERT(approx(libtellurian_transverse_radius_radians(D180), LIBTELLURIAN_EQUATORIAL_RADIUS));
ASSERT(approx(libtellurian_transverse_radius_radians(-D180), LIBTELLURIAN_EQUATORIAL_RADIUS));
ASSERT(approx(libtellurian_transverse_radius_radians(D30), libtellurian_transverse_radius_radians(-D30)));
ASSERT(approx(libtellurian_transverse_radius_radians(D45), libtellurian_transverse_radius_radians(-D45)));
ASSERT(approx(libtellurian_transverse_radius_radians(D60), libtellurian_transverse_radius_radians(-D60)));
ASSERT(approx(libtellurian_transverse_radius_radians(D90), libtellurian_transverse_radius_radians(-D90)));
ASSERT(fabs(libtellurian_transverse_radius_radians(D90) - 6399593.625758674) / 6399593.625758674 < 1.0e-8);
ASSERT(fabs(libtellurian_transverse_radius_radians(D30) - 6383480.917690154) / 6383480.917690154 < 1.0e-8);
return 0;
}
#endif
|
