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


double
libtellurian_elevated_gravity_radians(double gravity, double latitude, double altitude)
{
	double a = LIBTELLURIAN_EQUATORIAL_RADIUS;
	double b = LIBTELLURIAN_POLAR_RADIUS;
	double omega = LIBTELLURIAN_ANGULAR_VELOCITY;
	double GM = LIBTELLURIAN_GEOCENTRIC_GRAVITATIONAL_CONSTANT;
	double f = 1.0 - b / a;
	double m = (omega * a) * (omega * a) * (b / GM);
	double neg_k1 = fma(-2.0, f + m, -2.0) / a;
	double k2 = 4.0 * f / a;
	double neg_k3 = -3.0 / (a * a);
	double sin2_phi = sin(latitude) * sin(latitude);
	double s = fma(neg_k3, altitude, fma(k2, sin2_phi, neg_k1));
	return fma(s * altitude, gravity, gravity);
}


#else


int
main(void)
{
	ASSERT(libtellurian_elevated_gravity_radians(9.2, 0, 0) == 9.2);
	ASSERT(libtellurian_elevated_gravity_radians(4.2, 0, 0) == 4.2);
	ASSERT(libtellurian_elevated_gravity_radians(8.9, 2.0, 0) == 8.9);
	ASSERT(libtellurian_elevated_gravity_radians(8.9, 2.0, 1.0) < 8.9);
	ASSERT(libtellurian_elevated_gravity_radians(8.9, 2.0, 2.0) < libtellurian_elevated_gravity_radians(8.9, 2.0, 1.0));
	ASSERT(libtellurian_elevated_gravity_radians(8.9, 3.0, 1.0) < libtellurian_elevated_gravity_radians(8.9, 2.0, 1.0));
	return 0;
}


#endif