/* See LICENSE file for copyright and license details. */
#include "common.h"
#ifndef TEST


double
libtellurian_normal_gravity_radians(double latitude)
{
	double a = LIBTELLURIAN_EQUATORIAL_RADIUS;
	double b = LIBTELLURIAN_POLAR_RADIUS;
	double neg_e2 = fma(b / a, b / a, -1.0);
	double ag = a * LIBTELLURIAN_NORMAL_EQUATORIAL_GRAVITY;
	double bg = b * LIBTELLURIAN_NORMAL_POLAR_GRAVITY;
	double k = bg / ag - 1.0;
	double sin2_phi = sin(latitude) * sin(latitude);
	double num = fma(k, sin2_phi, 1.0);
	double denom2 = fma(neg_e2, sin2_phi, 1.0);
	return LIBTELLURIAN_NORMAL_EQUATORIAL_GRAVITY * num / sqrt(denom2);
}


#else


static int
approx(double a, double b)
{
        return fabs(a / b - 1.0) <= 1e-12;
}

int
main(void)
{
	ASSERT(approx(libtellurian_normal_gravity_radians(D30), libtellurian_normal_gravity_radians(-D30)));
	ASSERT(approx(libtellurian_normal_gravity_radians(D45), libtellurian_normal_gravity_radians(-D45)));
	ASSERT(approx(libtellurian_normal_gravity_radians(D60), libtellurian_normal_gravity_radians(-D60)));
	ASSERT(approx(libtellurian_normal_gravity_radians(D90), libtellurian_normal_gravity_radians(-D90)));
	ASSERT(approx(libtellurian_normal_gravity_radians(D180), libtellurian_normal_gravity_radians(0)));
	ASSERT(approx(libtellurian_normal_gravity_radians(-D180), libtellurian_normal_gravity_radians(0)));
	ASSERT(approx(libtellurian_normal_gravity_radians(2.0), libtellurian_normal_gravity_radians(-2.0)));
	ASSERT(approx(libtellurian_normal_gravity_radians(0), LIBTELLURIAN_NORMAL_EQUATORIAL_GRAVITY));
	ASSERT(approx(libtellurian_normal_gravity_radians(D90), LIBTELLURIAN_NORMAL_POLAR_GRAVITY));
	ASSERT(libtellurian_normal_gravity_radians(D45) < LIBTELLURIAN_NORMAL_POLAR_GRAVITY);
	ASSERT(libtellurian_normal_gravity_radians(D45) > LIBTELLURIAN_NORMAL_EQUATORIAL_GRAVITY);
	return 0;
}


#endif