001: /*
002: * GeoTools - OpenSource mapping toolkit
003: * http://geotools.org
004: * (C) 2003-2006, GeoTools Project Managment Committee (PMC)
005: * (C) 2001, Institut de Recherche pour le Développement
006: *
007: * This library is free software; you can redistribute it and/or
008: * modify it under the terms of the GNU Lesser General Public
009: * License as published by the Free Software Foundation;
010: * version 2.1 of the License.
011: *
012: * This library is distributed in the hope that it will be useful,
013: * but WITHOUT ANY WARRANTY; without even the implied warranty of
014: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
015: * Lesser General Public License for more details.
016: *
017: * This package contains documentation from OpenGIS specifications.
018: * OpenGIS consortium's work is fully acknowledged here.
019: */
020: package org.geotools.referencing.datum;
021:
022: // J2SE dependencies and extensions
023: import java.util.Map;
024: import javax.units.Unit;
025:
026: /**
027: * A ellipsoid which is spherical. This ellipsoid implements a faster
028: * {@link #orthodromicDistance} method.
029: *
030: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/referencing/src/main/java/org/geotools/referencing/datum/Spheroid.java $
031: * @version $Id: Spheroid.java 20874 2006-08-07 10:00:01Z jgarnett $
032: * @author Martin Desruisseaux
033: *
034: * @since 2.0
035: */
036: final class Spheroid extends DefaultEllipsoid {
037: /**
038: * Serial number for interoperability with different versions.
039: */
040: private static final long serialVersionUID = 7867565381280669821L;
041:
042: /**
043: * Constructs a new sphere using the specified radius.
044: *
045: * @param properties Set of properties. Should contains at least <code>"name"</code>.
046: * @param radius The equatorial and polar radius.
047: * @param ivfDefinitive {@code true} if the inverse flattening is definitive.
048: * @param unit The units of the radius value.
049: */
050: protected Spheroid(Map properties, double radius,
051: boolean ivfDefinitive, Unit unit) {
052: super (properties, check("radius", radius), radius,
053: Double.POSITIVE_INFINITY, ivfDefinitive, unit);
054: }
055:
056: /**
057: * Returns the orthodromic distance between two geographic coordinates.
058: * The orthodromic distance is the shortest distance between two points
059: * on a sphere's surface. The orthodromic path is always on a great circle.
060: *
061: * @param x1 Longitude of first point (in decimal degrees).
062: * @param y1 Latitude of first point (in decimal degrees).
063: * @param x2 Longitude of second point (in decimal degrees).
064: * @param y2 Latitude of second point (in decimal degrees).
065: * @return The orthodromic distance (in the units of this ellipsoid's axis).
066: */
067: public double orthodromicDistance(double x1, double y1, double x2,
068: double y2) {
069: /*
070: * The calculation of orthodromic distance on an ellipsoidal surface is complex,
071: * subject to rounding errors and has no solution near the poles. In some situation
072: * we use a calculation based on a spherical shape of the earth. A Fortran program
073: * which calculates orthodromic distances on an ellipsoidal surface can be downloaded
074: * from the NOAA site:
075: *
076: * ftp://ftp.ngs.noaa.gov/pub/pcsoft/for_inv.3d/source/
077: */
078: y1 = Math.toRadians(y1);
079: y2 = Math.toRadians(y2);
080: final double dx = Math.toRadians(Math.abs(x2 - x1) % 360);
081: double rho = Math.sin(y1) * Math.sin(y2) + Math.cos(y1)
082: * Math.cos(y2) * Math.cos(dx);
083: assert Math.abs(rho) < 1.0000001 : rho;
084: if (rho > +1)
085: rho = +1; // Catch rounding error.
086: if (rho < -1)
087: rho = -1; // Catch rounding error.
088: final double distance = Math.acos(rho) * getSemiMajorAxis();
089: /*
090: * Compare the distance with the orthodromic distance using ellipsoidal
091: * computation. This should be close to the same.
092: */
093: try {
094: double delta;
095: assert (delta = Math.abs(super .orthodromicDistance(x1, Math
096: .toDegrees(y1), x2, Math.toDegrees(y2))
097: - distance)) < getSemiMajorAxis() / 1E+9 : delta;
098: } catch (ArithmeticException exception) {
099: // The ellipsoidal model do not converge. Give up the assertion test.
100: // Note: the assertion fails for illegal latitudes (i.e. abs(y1)>90°
101: // or abs(y2)>90°).
102: }
103: return distance;
104: }
105: }
|