id	summary	reporter	owner	description	type	status	priority	milestone	component	version	resolution	keywords	cc
815	[PATCH] spherical computation of angles	SaschaR	framm	"The method LatLon#heading(LatLon other) has a bug: It lives
in a plane instead on a sphere.

The methode DrawAction#computeHelperLine() computes the
angle between two way segments by subtraction of two
heading angles. This is incorrect if the segments are
very long.


The following code computes the heading by means of
the spherical sine formular.

Proposal 1: Method LatLon#heading.

{{{
    /** earth radius in metres: 6378135 */
    public final static double EARTH_RADIUS = 6378135;

    /**
     * Returns the heading, in radians, that you have to use to get from
     * this lat/lon to another.
     *
     * @param other the ""destination"" position
     * @return heading
     */
    public double heading(LatLon other) {
        double rv;
        double distance = greatCircleDistance(other) / EARTH_RADIUS;
        double sinDistance = Math.sin(distance);
        if (Math.abs(sinDistance) > 0.0) {
            // spherical sine formula: Our triangle consists of
            // P0 = north pole = (lat=PI/2, lon=0),
            // P1 = this = (lan(), lon()),
            // P2 = other = (other.lan(), other.lon()).
            // We are looking for the angle rv at P1.
            // The distance between north pole and other is PI/2-other.lat().
            // The angle at north pole is delta_lon = other.lon()-lon().
            // sine formular:
            // sin(rv) = sin(delta_lon) * sin(PI/2-other.lat()) / sin(distance)
            //         = sin(delta_lon) * cos(other.lat()) / sin(distance).
            rv = Math.asin(Math.sin(Math.toRadians(other.lon()-lon()))
                    * Math.cos(Math.toRadians(other.lat())) / sinDistance);
            if (lat() > other.lat()) {
                rv = Math.PI - rv;
            }
            if (rv < 0) {
                rv += 2*Math.PI;
            }
        } else {
            rv = Double.NaN;
        }
        return rv;
    }
}}}


Proposal 2: Method LatLon#angle(LatLon p1, LatLon p2).

{{{
    /**
     * Computes the angle between the edges (this, p1) and (this, p2)
     * at this by use of cosine law.
     * @param p1 other point of first edge
     * @param p2 other point of second edge
     * @return angle in radians or NaN
     */
    public double angle(LatLon p1, LatLon p2) {
        double lenEdge1 = greatCircleDistance(p1)/LatLon.EARTH_RADIUS;
        double lenEdge2 = greatCircleDistance(p2)/LatLon.EARTH_RADIUS;
        double lenEdge3 = p1.greatCircleDistance(p2)/LatLon.EARTH_RADIUS;
        double numerator = Math.cos(lenEdge3) - Math.cos(lenEdge1)*Math.cos(lenEdge2);
        double denomin = Math.sin(lenEdge1)*Math.sin(lenEdge2);
        double angle;
        if (denomin != 0) {
            angle = Math.acos(numerator / denomin);
        } else {
            angle = Double.NaN;
        }
        return angle;
    }
}}}
Computation in DrawAction#computeHelperLine():

{{{
        if (previousNode != null) {
            angle = Math.toDegrees(currentBaseNode.coor.angle(
                        previousNode.coor, mouseLatLon));
        }
}}}
instead of
{{{
        if (previousNode != null) {
            angle = hdg - Math.toDegrees(previousNode.coor.heading(currentBaseNode.coor));
            if (angle < 0) angle += 360;
        }
}}}

With best regards

Sascha Rogmann
"	defect	closed	major		Core	latest	wontfix	angle Winkel Winkeldings