path: root/libtellurian_coarse_distance_radians.c

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

#define R LIBTELLURIAN_AUTHALIC_RADIUS

#ifndef TEST


const double libtellurian_coarse_distance_radius = R;


double
libtellurian_coarse_distance_radians(double latitude1, double longitude1,
                                     double latitude2, double longitude2)
{
	double h = fma(haversin(longitude2 - longitude1),
	               cos(latitude1) * cos(latitude2),
	               haversin(latitude2 - latitude1));
	return 2.0 * R * asin(sqrt(h));
}


#else


static int
check(double a, double b, double c, double d, double expect)
{
	double got = libtellurian_coarse_distance_radians(a, b, c, d);
	double rgot = libtellurian_coarse_distance_radians(c, d, a, b);
	double err = expect ? fabs(got / (R * expect) - 1.0) : got;
	return err <= 1.0e-8 && got == rgot;
}

int
main(void)
{
	ASSERT(check(0, 0, 0, 0, 0));
	ASSERT(check(0, 0, D90, 0, M_PI / 2.0));
	ASSERT(check(0, 0, 0, D90, M_PI / 2.0));
	ASSERT(check(0, 0, D180, 0, M_PI));
	ASSERT(check(0, 0, 0, D180, M_PI));
	ASSERT(check(0, 0, D45, 0, M_PI / 4.0));
	ASSERT(check(0, 0, 0, D45, M_PI / 4.0));
	ASSERT(check(-D90, 0, D90, 0, M_PI));
	ASSERT(check(0, -D90, 0, D90, M_PI));
	ASSERT(check(0, 0, D60, D90, M_PI / 2.0));
	ASSERT(check(0, 0, D30, D90, M_PI / 2.0));
	ASSERT(check(0, 0, D90, D90, M_PI / 2.0));
	ASSERT(check(0, 0, D45, D45, M_PI / 3.0));
	ASSERT(check(-D45, -D45, D45, D45, M_PI / 1.5));
	ASSERT(check(D45, -D45, D45, D45, M_PI / 3.0));
	ASSERT(check(-D45, D45, D45, D45, M_PI / 2.0));
	ASSERT(libtellurian_coarse_distance_radius == R);
	return 0;
}


#endif