Tolerance check is itself inaccurate

Description

The orthodromicDistance in org.geotools.referencing.operation.projection.MapProjection is not sufficiently accurate for its purpose.

for the points (-33.56099261594231, 1.480512392340082) and (-33.56099261594231, 1.4805123923400947) the calculated and true distances are:

orthodromicDistance: 0.09504164755344391
True distance: -1.3941023023489318E-9

This is a result of numerical error (loss of precision). The following sample code illustrates the issue.

private static final double EPSILON = 1E-6;
private double semiMajor = 6378137.0;

private double orthodromicDistance(final Point2D source, final Point2D target) {
final double y1 = toRadians(source.getY());
final double y2 = toRadians(target.getY());
final double dx = toRadians(abs(target.getX() - source.getX()) % 360);
double rho = sin(y1)*sin(y2) + cos(y1)*cos(y2)*cos(dx);
if (rho > +1)
{assert rho <= +(1+EPSILON) : rho; rho = +1;}
if (rho < -1)
{assert rho >= -(1+EPSILON) : rho; rho = -1;}
return acos(rho) * semiMajor;
}

public @Test void test()
{
Point2D.Double source = new Point2D.Double(-33.56099261594231, 1.480512392340082);
Point2D.Double target = new Point2D.Double(-33.56099261594231, 1.4805123923400947);

double distance = orthodromicDistance(source, target);
System.out.println("orthodromicDistance: "+distance);
if (source.getX() == target.getX())
System.out.println("True distance: "+semiMajor*(toRadians(source.getY())-toRadians(target.getY())));
}

Environment

None

Assignee

Unassigned

Reporter

codehaus

Triage

None

Components

Fix versions

Priority

Medium
Configure