blob: e5360259878fe586cbf29433dc3b7c33c69a45d4 (
plain) (
tree)
|
|
/* See LICENSE file for copyright and license details. */
#include "common.h"
#ifndef TEST
double
libtellurian_sea_level_radians(double latitude)
{
double a = LIBTELLURIAN_EQUATORIAL_RADIUS;
double b = LIBTELLURIAN_POLAR_RADIUS;
double c = cos(latitude);
double s = sin(latitude);
double x = a * c * a;
double y = b * s * b;
double num = fma(x, x, y * y);
double denom = fma(x, c, y * s);
return sqrt(num / denom);
}
#else
int
main(void)
{
double x = pow(LIBTELLURIAN_EQUATORIAL_RADIUS, 4.0) + pow(LIBTELLURIAN_POLAR_RADIUS, 4.0);
x = sqrt(x / (pow(LIBTELLURIAN_EQUATORIAL_RADIUS, 2.0) + pow(LIBTELLURIAN_POLAR_RADIUS, 2.0)));
ASSERT(libtellurian_sea_level_radians(0) == LIBTELLURIAN_EQUATORIAL_RADIUS);
ASSERT(libtellurian_sea_level_radians(D90) == LIBTELLURIAN_POLAR_RADIUS);
ASSERT(libtellurian_sea_level_radians(D180) == LIBTELLURIAN_EQUATORIAL_RADIUS);
ASSERT(libtellurian_sea_level_radians(-D90) == LIBTELLURIAN_POLAR_RADIUS);
ASSERT(libtellurian_sea_level_radians(D45) < LIBTELLURIAN_EQUATORIAL_RADIUS);
ASSERT(libtellurian_sea_level_radians(D45) > LIBTELLURIAN_POLAR_RADIUS);
ASSERT(libtellurian_sea_level_radians(-D45) == libtellurian_sea_level_radians(D45));
ASSERT(fabs(libtellurian_sea_level_radians(D45) / x - 1.0) < 1.0e-12);
return 0;
}
#endif
|