[PATCH] spherical computation of angles

[SaschaR]

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

[framm]

This is technically true, however we do not usually have segments longer than a few km (and even that is a stretch). We "live on a plane" in other parts of JOSM as well, e.g. when doing "align nodes in circle/line/rectangle" etc.; this patch would introduce half spherical and half planar treatment.

If someone thinks otherwise then please comment and re-open the ticket; I'll set it to wontfix now. Again, Sasche is technically correct, I just doubt that this "bug" has any meaningful consequences for the tasks JOSM is used for.