Changeset 2570 in josm for trunk/src/org/openstreetmap/josm


Ignore:
Timestamp:
2009-12-04T16:30:52+01:00 (15 years ago)
Author:
bastiK
Message:

Changed the greatCircleDistance formula to haver-sine.
The old cosine formula does not work for distances smaller 8 centimeter.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/org/openstreetmap/josm/data/coor/LatLon.java

    r2512 r2570  
    66import java.text.DecimalFormat;
    77import java.text.NumberFormat;
     8
     9import static java.lang.Math.*;
    810
    911import org.openstreetmap.josm.Main;
     
    121123    /**
    122124     * Computes the distance between this lat/lon and another point on the earth.
    123      * Uses spherical law of cosines formula, not Haversine.
     125     * Uses Haversine formular.
    124126     * @param other the other point.
    125127     * @return distance in metres.
    126128     */
    127129    public double greatCircleDistance(LatLon other) {
    128         return (Math.acos(
    129                 Math.sin(Math.toRadians(lat())) * Math.sin(Math.toRadians(other.lat())) +
    130                 Math.cos(Math.toRadians(lat()))*Math.cos(Math.toRadians(other.lat())) *
    131                 Math.cos(Math.toRadians(other.lon()-lon()))) * 6378135);
     130        double R = 6378135;
     131        double sinHalfLat = sin(toRadians(other.lat() - this.lat()) / 2);
     132        double sinHalfLon = sin(toRadians(other.lon() - this.lon()) / 2);
     133        double d = 2 * R * asin(
     134                            sqrt(sinHalfLat*sinHalfLat +
     135                            cos(toRadians(this.lat()))*cos(toRadians(other.lat()))*sinHalfLon*sinHalfLon));
     136        // For points opposite to each other on the sphere,
     137        // rounding errors could make the argument of asin greater than 1
     138        // (This should almost never happen.)
     139        if (java.lang.Double.isNaN(d)) {
     140            System.err.println("Error: NaN in greatCircleDistance");
     141            d = PI * R;
     142        }
     143        return d;
    132144    }
    133145
Note: See TracChangeset for help on using the changeset viewer.