/* 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