0001: /*
0002: * GeoTools - OpenSource mapping toolkit
0003: * http://geotools.org
0004: * (C) 2004-2006, GeoTools Project Managment Committee (PMC)
0005: *
0006: * This library is free software; you can redistribute it and/or
0007: * modify it under the terms of the GNU Lesser General Public
0008: * License as published by the Free Software Foundation;
0009: * version 2.1 of the License.
0010: *
0011: * This library is distributed in the hope that it will be useful,
0012: * but WITHOUT ANY WARRANTY; without even the implied warranty of
0013: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
0014: * Lesser General Public License for more details.
0015: *
0016: * Portions of this file is adapted from Fortran code provided by NOAA.
0017: * Programmed for CDC-6600 by LCDR L.Pfeifer NGS ROCKVILLE MD 18FEB75
0018: * Modified for IBM SYSTEM 360 by John G.Gergen NGS ROCKVILLE MD 7507
0019: * Source: ftp://ftp.ngs.noaa.gov/pub/pcsoft/for_inv.3d/source/
0020: */
0021: package org.geotools.referencing;
0022:
0023: // J2SE dependencies
0024: import java.awt.Shape;
0025: import java.awt.geom.GeneralPath;
0026: import java.awt.geom.Line2D;
0027: import java.awt.geom.Point2D;
0028: import java.text.Format;
0029: import java.util.Iterator;
0030: import java.util.List;
0031: import javax.units.NonSI;
0032:
0033: // OpenGIS dependencies
0034: import org.opengis.referencing.datum.Datum;
0035: import org.opengis.referencing.datum.Ellipsoid;
0036: import org.opengis.referencing.datum.GeodeticDatum;
0037: import org.opengis.referencing.operation.TransformException;
0038: import org.opengis.referencing.cs.CoordinateSystemAxis;
0039: import org.opengis.referencing.cs.CoordinateSystem;
0040: import org.opengis.referencing.cs.AxisDirection;
0041: import org.opengis.referencing.crs.CompoundCRS;
0042: import org.opengis.referencing.crs.GeographicCRS;
0043: import org.opengis.referencing.crs.CoordinateReferenceSystem;
0044: import org.opengis.geometry.coordinate.Position;
0045: import org.opengis.geometry.DirectPosition;
0046:
0047: // Geotools dependencies
0048: import org.geotools.measure.Angle;
0049: import org.geotools.measure.Latitude;
0050: import org.geotools.measure.Longitude;
0051: import org.geotools.measure.CoordinateFormat;
0052: import org.geotools.geometry.DirectPosition2D;
0053: import org.geotools.geometry.GeneralDirectPosition;
0054: import org.geotools.geometry.TransformedDirectPosition;
0055: import org.geotools.referencing.datum.DefaultEllipsoid;
0056: import org.geotools.referencing.datum.DefaultPrimeMeridian;
0057: import org.geotools.referencing.datum.DefaultGeodeticDatum;
0058: import org.geotools.referencing.crs.DefaultGeographicCRS;
0059: import org.geotools.referencing.cs.DefaultEllipsoidalCS;
0060: import org.geotools.resources.CRSUtilities;
0061: import org.geotools.resources.i18n.Errors;
0062: import org.geotools.resources.i18n.ErrorKeys;
0063: import org.geotools.resources.i18n.Vocabulary;
0064: import org.geotools.resources.i18n.VocabularyKeys;
0065: import org.geotools.io.TableWriter;
0066:
0067: /**
0068: * Performs geodetic calculations on an {@linkplain Ellipsoid ellipsoid}. This class encapsulates
0069: * a generic ellipsoid and calculates the following properties:
0070: * <p>
0071: * <ul>
0072: * <li>Distance and azimuth between two points.</li>
0073: * <li>Point located at a given distance and azimuth from an other point.</li>
0074: * </ul>
0075: * <p>
0076: * The calculation use the following informations:
0077: * <p>
0078: * <ul>
0079: * <li>The {@linkplain #setStartingPosition starting position}, which is always considered valid.
0080: * It is initially set at (0,0) and can only be changed to another legitimate value.</li>
0081: * <li><strong>Only one</strong> of the following:
0082: * <ul>
0083: * <li>The {@linkplain #setDestinationPosition destination position}, or</li>
0084: * <li>An {@linkplain #setDirection azimuth and distance}.</li>
0085: * </ul>
0086: * The latest one set overrides the other and determines what will be calculated.</li>
0087: * </ul>
0088: * <p>
0089: * Note: This class is not thread-safe. If geodetic calculations are needed in a multi-threads
0090: * environment, create one distinct instance of {@code GeodeticCalculator} for each thread.
0091: *
0092: * @since 2.1
0093: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/referencing/src/main/java/org/geotools/referencing/GeodeticCalculator.java $
0094: * @version $Id: GeodeticCalculator.java 24925 2007-03-27 20:12:08Z jgarnett $
0095: * @author Daniele Franzoni
0096: * @author Martin Desruisseaux
0097: */
0098: public class GeodeticCalculator {
0099: /**
0100: * Tolerance factors from the strictest (<code>TOLERANCE_0</CODE>)
0101: * to the most relax one (<code>TOLERANCE_3</CODE>).
0102: */
0103: private static final double TOLERANCE_0 = 5.0e-15, // tol0
0104: TOLERANCE_1 = 5.0e-14, // tol1
0105: TOLERANCE_2 = 5.0e-13, // tt
0106: TOLERANCE_3 = 7.0e-3; // tol2
0107:
0108: /**
0109: * Tolerance factor for assertions. It has no impact on computed values.
0110: */
0111: private static final double TOLERANCE_CHECK = 1E-8;
0112:
0113: /**
0114: * The transform from user coordinates to geodetic coordinates used for computation,
0115: * or {@code null} if no transformations are required.
0116: */
0117: private final TransformedDirectPosition userToGeodetic;
0118:
0119: /**
0120: * The coordinate reference system for all methods working on {@link Position} objects.
0121: * If {@code null}, will be created the first time {@link #getCoordinateReferenceSystem}
0122: * is invoked.
0123: */
0124: private CoordinateReferenceSystem coordinateReferenceSystem;
0125:
0126: /**
0127: * The coordinate reference system for all methods working on {@link Point2D} objects.
0128: * If {@code null}, will be created the first time {@link #getGeographicCRS} is invoked.
0129: */
0130: private GeographicCRS geographicCRS;
0131:
0132: /**
0133: * The encapsulated ellipsoid.
0134: */
0135: private final Ellipsoid ellipsoid;
0136:
0137: /*
0138: * The semi major axis of the refereced ellipsoid.
0139: */
0140: private final double semiMajorAxis;
0141:
0142: /*
0143: * The semi minor axis of the refereced ellipsoid.
0144: */
0145: private final double semiMinorAxis;
0146:
0147: /*
0148: * The eccenticity squared of the refereced ellipsoid.
0149: */
0150: private final double eccentricitySquared;
0151:
0152: /*
0153: * The maximum orthodromic distance that could be calculated onto the referenced ellipsoid.
0154: */
0155: private final double maxOrthodromicDistance;
0156:
0157: /**
0158: * GPNARC parameters computed from the ellipsoid.
0159: */
0160: private final double A, B, C, D, E, F;
0161:
0162: /**
0163: * GPNHRI parameters computed from the ellipsoid.
0164: *
0165: * {@code f} if the flattening of the referenced ellipsoid. {@code f2},
0166: * {@code f3} and {@code f4} are <var>f<sup>2</sup></var>,
0167: * <var>f<sup>3</sup></var> and <var>f<sup>4</sup></var> respectively.
0168: */
0169: private final double fo, f, f2, f3, f4;
0170:
0171: /**
0172: * Parameters computed from the ellipsoid.
0173: */
0174: private final double T1, T2, T4, T6;
0175:
0176: /**
0177: * Parameters computed from the ellipsoid.
0178: */
0179: private final double a01, a02, a03, a21, a22, a23, a42, a43, a63;
0180:
0181: /**
0182: * The (<var>latitude</var>, <var>longitude</var>) coordinate of the first point
0183: * <strong>in radians</strong>. This point is set by {@link #setStartingGeographicPoint}.
0184: */
0185: private double lat1, long1;
0186:
0187: /**
0188: * The (<var>latitude</var>, <var>longitude</var>) coordinate of the destination point
0189: * <strong>in radians</strong>. This point is set by {@link #setDestinationGeographicPoint}.
0190: */
0191: private double lat2, long2;
0192:
0193: /**
0194: * The distance and azimuth (in radians) from the starting point
0195: * ({@link #long1}, {@link #lat1}) to the destination point
0196: * ({@link #long2}, {@link #lat2}).
0197: */
0198: private double distance, azimuth;
0199:
0200: /**
0201: * Tell if the destination point is valid.
0202: * {@code false} if {@link #long2} and {@link #lat2} need to be computed.
0203: */
0204: private boolean destinationValid;
0205:
0206: /**
0207: * Tell if the azimuth and the distance are valids.
0208: * {@code false} if {@link #distance} and {@link #azimuth} need to be computed.
0209: */
0210: private boolean directionValid;
0211:
0212: /**
0213: * Constructs a new geodetic calculator associated with the WGS84 ellipsoid.
0214: */
0215: public GeodeticCalculator() {
0216: this (DefaultEllipsoid.WGS84);
0217: }
0218:
0219: /**
0220: * Constructs a new geodetic calculator associated with the specified ellipsoid.
0221: * All calculations done by the new instance are referenced to this ellipsoid.
0222: *
0223: * @param ellipsoid The ellipsoid onto which calculates distances and azimuths.
0224: */
0225: public GeodeticCalculator(final Ellipsoid ellipsoid) {
0226: this (ellipsoid, null);
0227: }
0228:
0229: /**
0230: * Constructs a new geodetic calculator expecting coordinates in the supplied CRS.
0231: * The ellipsoid will be inferred from the CRS.
0232: *
0233: * @param crs The reference system for the {@link Position} objects.
0234: *
0235: * @since 2.2
0236: */
0237: public GeodeticCalculator(final CoordinateReferenceSystem crs) {
0238: this (CRS.getEllipsoid(crs), crs);
0239: }
0240:
0241: /**
0242: * For internal use by public constructors only.
0243: */
0244: private GeodeticCalculator(final Ellipsoid ellipsoid,
0245: final CoordinateReferenceSystem crs) {
0246: if (ellipsoid == null) {
0247: throw new IllegalArgumentException(Errors.format(
0248: ErrorKeys.NULL_ARGUMENT_$1, "ellipsoid"));
0249: }
0250: this .ellipsoid = ellipsoid;
0251: this .semiMajorAxis = ellipsoid.getSemiMajorAxis();
0252: this .semiMinorAxis = ellipsoid.getSemiMinorAxis();
0253: if (crs != null) {
0254: coordinateReferenceSystem = crs;
0255: geographicCRS = getGeographicCRS(crs);
0256: /*
0257: * Note: there is no need to set Hints.LENIENT_DATUM_SHIFT to Boolean.TRUE here since
0258: * the target CRS computed by our internal getGeographicCRS(crs) method should
0259: * returns a CRS using the same datum than the specified CRS. If the factory
0260: * fails with a "Bursa-Wolf parameters required" error message, then we probably
0261: * have a bug somewhere.
0262: */
0263: userToGeodetic = new TransformedDirectPosition(crs,
0264: geographicCRS, null);
0265: } else {
0266: userToGeodetic = null;
0267: }
0268:
0269: /* calculation of GPNHRI parameters */
0270: f = (semiMajorAxis - semiMinorAxis) / semiMajorAxis;
0271: fo = 1.0 - f;
0272: f2 = f * f;
0273: f3 = f * f2;
0274: f4 = f * f3;
0275: eccentricitySquared = f * (2.0 - f);
0276:
0277: /* Calculation of GNPARC parameters */
0278: final double E2 = eccentricitySquared;
0279: final double E4 = E2 * E2;
0280: final double E6 = E4 * E2;
0281: final double E8 = E6 * E2;
0282: final double EX = E8 * E2;
0283:
0284: A = 1.0 + 0.75 * E2 + 0.703125 * E4 + 0.68359375 * E6
0285: + 0.67291259765625 * E8 + 0.6661834716796875 * EX;
0286: B = 0.75 * E2 + 0.9375 * E4 + 1.025390625 * E6 + 1.07666015625
0287: * E8 + 1.1103057861328125 * EX;
0288: C = 0.234375 * E4 + 0.41015625 * E6 + 0.538330078125 * E8
0289: + 0.63446044921875 * EX;
0290: D = 0.068359375 * E6 + 0.15380859375 * E8 + 0.23792266845703125
0291: * EX;
0292: E = 0.01922607421875 * E8 + 0.0528717041015625 * EX;
0293: F = 0.00528717041015625 * EX;
0294:
0295: maxOrthodromicDistance = semiMajorAxis * (1.0 - E2) * Math.PI
0296: * A - 1.0;
0297:
0298: T1 = 1.0;
0299: T2 = -0.25 * f * (1.0 + f + f2);
0300: T4 = 0.1875 * f2 * (1.0 + 2.25 * f);
0301: T6 = 0.1953125 * f3;
0302:
0303: final double a = f3 * (1.0 + 2.25 * f);
0304: a01 = -f2 * (1.0 + f + f2) / 4.0;
0305: a02 = 0.1875 * a;
0306: a03 = -0.1953125 * f4;
0307: a21 = -a01;
0308: a22 = -0.25 * a;
0309: a23 = 0.29296875 * f4;
0310: a42 = 0.03125 * a;
0311: a43 = 0.05859375 * f4;
0312: a63 = 5.0 * f4 / 768.0;
0313: }
0314:
0315: ///////////////////////////////////////////////////////////
0316: //////// ////////
0317: //////// H E L P E R M E T H O D S ////////
0318: //////// ////////
0319: ///////////////////////////////////////////////////////////
0320:
0321: /**
0322: * Returns the first two-dimensional geographic CRS using standard axis, creating one if needed.
0323: */
0324: private static GeographicCRS getGeographicCRS(
0325: final CoordinateReferenceSystem crs) {
0326: if (crs instanceof GeographicCRS) {
0327: final CoordinateSystem cs = crs.getCoordinateSystem();
0328: if (cs.getDimension() == 2
0329: && isStandard(cs.getAxis(0), AxisDirection.EAST)
0330: && isStandard(cs.getAxis(1), AxisDirection.NORTH)) {
0331: return (GeographicCRS) crs;
0332: }
0333: }
0334: final Datum datum = CRSUtilities.getDatum(crs);
0335: if (datum instanceof GeodeticDatum) {
0336: return new DefaultGeographicCRS("Geodetic",
0337: (GeodeticDatum) datum,
0338: DefaultEllipsoidalCS.GEODETIC_2D);
0339: }
0340: if (crs instanceof CompoundCRS) {
0341: final List components = ((CompoundCRS) crs)
0342: .getCoordinateReferenceSystems();
0343: for (final Iterator it = components.iterator(); it
0344: .hasNext();) {
0345: final GeographicCRS candidate = getGeographicCRS((CoordinateReferenceSystem) it
0346: .next());
0347: if (candidate != null) {
0348: return candidate;
0349: }
0350: }
0351: }
0352: throw new IllegalArgumentException(Errors
0353: .format(ErrorKeys.ILLEGAL_COORDINATE_REFERENCE_SYSTEM));
0354: }
0355:
0356: /**
0357: * Returns {@code true} if the specified axis is oriented toward the specified direction and
0358: * uses decimal degrees units.
0359: */
0360: private static boolean isStandard(final CoordinateSystemAxis axis,
0361: final AxisDirection direction) {
0362: return direction.equals(axis.getDirection())
0363: && NonSI.DEGREE_ANGLE.equals(axis.getUnit());
0364: }
0365:
0366: /**
0367: * Returns an angle between -{@linkplain Math#PI PI} and {@linkplain Math#PI PI}
0368: * equivalent to the specified angle in radians.
0369: *
0370: * @param alpha An angle value in radians.
0371: * @return The angle between between -{@linkplain Math#PI PI} and {@linkplain Math#PI PI}.
0372: *
0373: */
0374: private static double castToAngleRange(final double alpha) {
0375: return alpha - (2 * Math.PI)
0376: * Math.floor(alpha / (2 * Math.PI) + 0.5);
0377: }
0378:
0379: /**
0380: * Checks the latidude validity. The argument {@code latidude} should be
0381: * greater or equal than -90 degrees and lower or equals than +90 degrees. As
0382: * a convenience, this method returns the latitude in radians.
0383: *
0384: * @param latitude The latitude value in <strong>decimal degrees</strong>.
0385: * @return The latitude value in <strong>radians</strong>.
0386: * @throws IllegalArgumentException if {@code latitude} is not between -90 and +90 degrees.
0387: */
0388: private static double checkLatitude(final double latitude)
0389: throws IllegalArgumentException {
0390: if (latitude >= Latitude.MIN_VALUE
0391: && latitude <= Latitude.MAX_VALUE) {
0392: return Math.toRadians(latitude);
0393: }
0394: throw new IllegalArgumentException(Errors.format(
0395: ErrorKeys.LATITUDE_OUT_OF_RANGE_$1, new Latitude(
0396: latitude)));
0397: }
0398:
0399: /**
0400: * Checks the longitude validity. The argument {@code longitude} should be
0401: * greater or equal than -180 degrees and lower or equals than +180 degrees. As
0402: * a convenience, this method returns the longitude in radians.
0403: *
0404: * @param longitude The longitude value in <strong>decimal degrees</strong>.
0405: * @return The longitude value in <strong>radians</strong>.
0406: * @throws IllegalArgumentException if {@code longitude} is not between -180 and +180 degrees.
0407: */
0408: private static double checkLongitude(final double longitude)
0409: throws IllegalArgumentException {
0410: if (longitude >= Longitude.MIN_VALUE
0411: && longitude <= Longitude.MAX_VALUE) {
0412: return Math.toRadians(longitude);
0413: }
0414: throw new IllegalArgumentException(Errors.format(
0415: ErrorKeys.LONGITUDE_OUT_OF_RANGE_$1, new Longitude(
0416: longitude)));
0417: }
0418:
0419: /**
0420: * Checks the azimuth validity. The argument {@code azimuth} should be
0421: * greater or equal than -180 degrees and lower or equals than +180 degrees.
0422: * As a convenience, this method returns the azimuth in radians.
0423: *
0424: * @param azimuth The azimuth value in <strong>decimal degrees</strong>.
0425: * @return The azimuth value in <strong>radians</strong>.
0426: * @throws IllegalArgumentException if {@code azimuth} is not between -180 and +180 degrees.
0427: */
0428: private static double checkAzimuth(final double azimuth)
0429: throws IllegalArgumentException {
0430: if (azimuth >= -180.0 && azimuth <= 180.0) {
0431: return Math.toRadians(azimuth);
0432: }
0433: throw new IllegalArgumentException(Errors.format(
0434: ErrorKeys.AZIMUTH_OUT_OF_RANGE_$1, new Longitude(
0435: azimuth)));
0436: }
0437:
0438: /**
0439: * Checks the orthodromic distance validity. Arguments {@code orthodromicDistance}
0440: * should be greater or equal than 0 and lower or equals than the maximum orthodromic distance.
0441: *
0442: * @param distance The orthodromic distance value.
0443: * @throws IllegalArgumentException if {@code orthodromic distance} is not between
0444: * 0 and the maximum orthodromic distance.
0445: */
0446: private void checkOrthodromicDistance(final double distance)
0447: throws IllegalArgumentException {
0448: if (!(distance >= 0.0 && distance <= maxOrthodromicDistance)) {
0449: throw new IllegalArgumentException(Errors.format(
0450: ErrorKeys.DISTANCE_OUT_OF_RANGE_$4, new Double(
0451: distance), new Double(0), new Double(
0452: maxOrthodromicDistance), ellipsoid
0453: .getAxisUnit()));
0454: }
0455: }
0456:
0457: /**
0458: * Checks the number of verteces in a curve. Arguments {@code numberOfPoints}
0459: * should be not negative.
0460: *
0461: * @param numberOfPonits The number of verteces in a curve.
0462: * @throws IllegalArgumentException if {@code numberOfVerteces} is negative.
0463: */
0464: private static void checkNumberOfPoints(final int numberOfPoints)
0465: throws IllegalArgumentException {
0466: if (numberOfPoints < 0) {
0467: throw new IllegalArgumentException(Errors.format(
0468: ErrorKeys.ILLEGAL_ARGUMENT_$2, "numberOfPoints",
0469: new Integer(numberOfPoints)));
0470: }
0471: }
0472:
0473: /**
0474: * Returns a localized "No convergence" error message. The error message
0475: * includes informations about starting and destination points.
0476: */
0477: private String getNoConvergenceErrorMessage() {
0478: final CoordinateFormat cf = new CoordinateFormat();
0479: return Errors.format(ErrorKeys.NO_CONVERGENCE_$2, format(cf,
0480: long1, lat1), format(cf, long2, lat2));
0481: }
0482:
0483: /**
0484: * Format the specified coordinates using the specified formatter, which should be an instance
0485: * of {@link CoordinateFormat}.
0486: */
0487: private static String format(final Format cf,
0488: final double longitude, final double latitude) {
0489: return cf.format(new GeneralDirectPosition(Math
0490: .toDegrees(longitude), Math.toDegrees(latitude)));
0491: }
0492:
0493: ///////////////////////////////////////////////////////////////
0494: //////// ////////
0495: //////// G E O D E T I C M E T H O D S ////////
0496: //////// ////////
0497: ///////////////////////////////////////////////////////////////
0498:
0499: /**
0500: * Returns the coordinate reference system for all methods working on {@link Position} objects.
0501: * This is the CRS specified at {@linkplain #GeodeticCalculator(CoordinateReferenceSystem)
0502: * construction time}.
0503: *
0504: * @since 2.2
0505: */
0506: public CoordinateReferenceSystem getCoordinateReferenceSystem() {
0507: if (coordinateReferenceSystem == null) {
0508: coordinateReferenceSystem = getGeographicCRS();
0509: }
0510: return coordinateReferenceSystem;
0511: }
0512:
0513: /**
0514: * Returns the geographic coordinate reference system for all methods working
0515: * on {@link Point2D} objects. This is inferred from the CRS specified at
0516: * {@linkplain #GeodeticCalculator(CoordinateReferenceSystem) construction time}.
0517: *
0518: * @since 2.3
0519: */
0520: public GeographicCRS getGeographicCRS() {
0521: if (geographicCRS == null) {
0522: final String name = Vocabulary
0523: .format(VocabularyKeys.GEODETIC_2D);
0524: ;
0525: geographicCRS = new DefaultGeographicCRS(name,
0526: new DefaultGeodeticDatum(name, getEllipsoid(),
0527: DefaultPrimeMeridian.GREENWICH),
0528: DefaultEllipsoidalCS.GEODETIC_2D);
0529: }
0530: return geographicCRS;
0531: }
0532:
0533: /**
0534: * Returns the referenced ellipsoid.
0535: */
0536: public Ellipsoid getEllipsoid() {
0537: return ellipsoid;
0538: }
0539:
0540: /**
0541: * Set the starting point in geographic coordinates.
0542: * The {@linkplain #getAzimuth() azimuth},
0543: * the {@linkplain #getOrthodromicDistance() orthodromic distance} and
0544: * the {@linkplain #getDestinationGeographicPoint() destination point}
0545: * are discarted. They will need to be specified again.
0546: *
0547: * @param longitude The longitude in decimal degrees between -180 and +180°
0548: * @param latitude The latitude in decimal degrees between -90 and +90°
0549: * @throws IllegalArgumentException if the longitude or the latitude is out of bounds.
0550: *
0551: * @since 2.3
0552: */
0553: public void setStartingGeographicPoint(double longitude,
0554: double latitude) throws IllegalArgumentException {
0555: // Check first in case an exception is raised
0556: // (in other words, we change all or nothing).
0557: longitude = checkLongitude(longitude);
0558: latitude = checkLatitude(latitude);
0559: // Check passed. Now performs the changes in this object.
0560: long1 = longitude;
0561: lat1 = latitude;
0562: destinationValid = false;
0563: directionValid = false;
0564: }
0565:
0566: /**
0567: * Set the starting point in geographic coordinates. The <var>x</var> and <var>y</var>
0568: * coordinates must be the longitude and latitude in decimal degrees, respectively.
0569: *
0570: * This is a convenience method for
0571: * <code>{@linkplain #setStartingGeographicPoint(double,double)
0572: * setStartingGeographicPoint}(x,y)</code>.
0573: *
0574: * @param point The starting point.
0575: * @throws IllegalArgumentException if the longitude or the latitude is out of bounds.
0576: *
0577: * @since 2.3
0578: */
0579: public void setStartingGeographicPoint(final Point2D point)
0580: throws IllegalArgumentException {
0581: setStartingGeographicPoint(point.getX(), point.getY());
0582: }
0583:
0584: /**
0585: * Set the starting position in user coordinates, which doesn't need to be geographic.
0586: * The coordinate reference system is the one specified to the
0587: * {@linkplain #GeodeticCalculator(CoordinateReferenceSystem) constructor}.
0588: *
0589: * @param position The position in user coordinate reference system.
0590: * @throws TransformException if the position can't be transformed.
0591: *
0592: * @since 2.3
0593: */
0594: public void setStartingPosition(final Position position)
0595: throws TransformException {
0596: DirectPosition p = position.getPosition();
0597: if (userToGeodetic != null) {
0598: userToGeodetic.transform(p);
0599: p = userToGeodetic;
0600: }
0601: setStartingGeographicPoint(p.getOrdinate(0), p.getOrdinate(1));
0602: }
0603:
0604: /**
0605: * Returns the starting point in geographic coordinates. The <var>x</var> and <var>y</var>
0606: * coordinates are the longitude and latitude in decimal degrees, respectively. If the
0607: * starting point has never been set, then the default value is (0,0).
0608: *
0609: * @return The starting point in geographic coordinates.
0610: *
0611: * @since 2.3
0612: */
0613: public Point2D getStartingGeographicPoint() {
0614: return new Point2D.Double(Math.toDegrees(long1), Math
0615: .toDegrees(lat1));
0616: }
0617:
0618: /**
0619: * Returns the starting position in user coordinates, which doesn't need to be geographic.
0620: * The coordinate reference system is the one specified to the
0621: * {@linkplain #GeodeticCalculator(CoordinateReferenceSystem) constructor}.
0622: *
0623: * @throws TransformException if the position can't be transformed to user coordinates.
0624: *
0625: * @since 2.3
0626: */
0627: public DirectPosition getStartingPosition()
0628: throws TransformException {
0629: DirectPosition position = userToGeodetic;
0630: if (position == null) {
0631: position = new DirectPosition2D();
0632: }
0633: position.setOrdinate(0, Math.toDegrees(long1));
0634: position.setOrdinate(1, Math.toDegrees(lat1));
0635: if (userToGeodetic != null) {
0636: position = userToGeodetic.inverseTransform();
0637: }
0638: return position;
0639: }
0640:
0641: /**
0642: * Set the destination point in geographic coordinates. The azimuth and distance values
0643: * will be updated as a side effect of this call. They will be recomputed the next time
0644: * {@link #getAzimuth()} or {@link #getOrthodromicDistance()} are invoked.
0645: *
0646: * @param longitude The longitude in decimal degrees between -180 and +180°
0647: * @param latitude The latgitude in decimal degrees between -90 and +90°
0648: * @throws IllegalArgumentException if the longitude or the latitude is out of bounds.
0649: *
0650: * @since 2.3
0651: */
0652: public void setDestinationGeographicPoint(double longitude,
0653: double latitude) throws IllegalArgumentException {
0654: // Check first in case an exception is raised
0655: // (in other words, we change all or nothing).
0656: longitude = checkLongitude(longitude);
0657: latitude = checkLatitude(latitude);
0658: // Check passed. Now performs the changes in this object.
0659: long2 = longitude;
0660: lat2 = latitude;
0661: destinationValid = true;
0662: directionValid = false;
0663: }
0664:
0665: /**
0666: * Set the destination point in geographic coordinates. The <var>x</var> and <var>y</var>
0667: * coordinates must be the longitude and latitude in decimal degrees, respectively.
0668: *
0669: * This is a convenience method for
0670: * <code>{@linkplain #setDestinationGeographicPoint(double,double)
0671: * setDestinationGeographicPoint}(x,y)</code>.
0672: *
0673: * @param point The destination point.
0674: * @throws IllegalArgumentException if the longitude or the latitude is out of bounds.
0675: *
0676: * @since 2.3
0677: */
0678: public void setDestinationGeographicPoint(final Point2D point)
0679: throws IllegalArgumentException {
0680: setDestinationGeographicPoint(point.getX(), point.getY());
0681: }
0682:
0683: /**
0684: * Set the destination position in user coordinates, which doesn't need to be geographic.
0685: * The coordinate reference system is the one specified to the
0686: * {@linkplain #GeodeticCalculator(CoordinateReferenceSystem) constructor}.
0687: *
0688: * @param position The position in user coordinate reference system.
0689: * @throws TransformException if the position can't be transformed.
0690: *
0691: * @since 2.2
0692: */
0693: public void setDestinationPosition(final Position position)
0694: throws TransformException {
0695: DirectPosition p = position.getPosition();
0696: if (userToGeodetic != null) {
0697: userToGeodetic.transform(p);
0698: p = userToGeodetic;
0699: }
0700: setDestinationGeographicPoint(p.getOrdinate(0), p
0701: .getOrdinate(1));
0702: }
0703:
0704: /**
0705: * Returns the destination point. This method returns the point set by the last
0706: * call to a <code>{@linkplain #setDestinationGeographicPoint(double,double)
0707: * setDestinationGeographicPoint}(...)</code>
0708: * method, <strong>except</strong> if
0709: * <code>{@linkplain #setDirection(double,double) setDirection}(...)</code> has been
0710: * invoked after. In this later case, the destination point will be computed from the
0711: * {@linkplain #getStartingGeographicPoint starting point} to the azimuth and distance
0712: * specified.
0713: *
0714: * @return The destination point. The <var>x</var> and <var>y</var> coordinates
0715: * are the longitude and latitude in decimal degrees, respectively.
0716: * @throws IllegalStateException if the azimuth and the distance have not been set.
0717: *
0718: * @since 2.3
0719: */
0720: public Point2D getDestinationGeographicPoint()
0721: throws IllegalStateException {
0722: if (!destinationValid) {
0723: computeDestinationPoint();
0724: }
0725: return new Point2D.Double(Math.toDegrees(long2), Math
0726: .toDegrees(lat2));
0727: }
0728:
0729: /**
0730: * Returns the destination position in user coordinates, which doesn't need to be geographic.
0731: * The coordinate reference system is the one specified to the
0732: * {@linkplain #GeodeticCalculator(CoordinateReferenceSystem) constructor}.
0733: *
0734: * @throws TransformException if the position can't be transformed to user coordinates.
0735: *
0736: * @since 2.2
0737: */
0738: public DirectPosition getDestinationPosition()
0739: throws TransformException {
0740: if (!destinationValid) {
0741: computeDestinationPoint();
0742: }
0743: DirectPosition position = userToGeodetic;
0744: if (position == null) {
0745: position = new DirectPosition2D();
0746: }
0747: position.setOrdinate(0, Math.toDegrees(long2));
0748: position.setOrdinate(1, Math.toDegrees(lat2));
0749: if (userToGeodetic != null) {
0750: position = userToGeodetic.inverseTransform();
0751: }
0752: return position;
0753: }
0754:
0755: /**
0756: * Set the azimuth and the distance from the {@linkplain #getStartingGeographicPoint
0757: * starting point}. The destination point will be updated as a side effect of this call.
0758: * It will be recomputed the next time {@link #getDestinationGeographicPoint()} is invoked.
0759: *
0760: * @param azimuth The azimuth in decimal degrees from -180° to 180°.
0761: * @param distance The orthodromic distance in the same units as the
0762: * {@linkplain #getEllipsoid ellipsoid} axis.
0763: * @throws IllegalArgumentException if the azimuth or the distance is out of bounds.
0764: *
0765: * @see #getAzimuth
0766: * @see #getOrthodromicDistance
0767: */
0768: public void setDirection(double azimuth, final double distance)
0769: throws IllegalArgumentException {
0770: // Check first in case an exception is raised
0771: // (in other words, we change all or nothing).
0772: azimuth = checkAzimuth(azimuth);
0773: checkOrthodromicDistance(distance);
0774: // Check passed. Now performs the changes in this object.
0775: this .azimuth = azimuth;
0776: this .distance = distance;
0777: destinationValid = false;
0778: directionValid = true;
0779: }
0780:
0781: /**
0782: * Returns the azimuth. This method returns the value set by the last call to
0783: * <code>{@linkplain #setDirection(double,double) setDirection}(azimuth,distance)</code>,
0784: * <strong>except</strong> if <code>{@linkplain #setDestinationGeographicPoint(double,double)
0785: * setDestinationGeographicPoint}(...)</code> has been invoked after. In this later case, the
0786: * azimuth will be computed from the {@linkplain #getStartingGeographicPoint starting point}
0787: * to the destination point.
0788: *
0789: * @return The azimuth, in decimal degrees from -180° to +180°.
0790: * @throws IllegalStateException if the destination point has not been set.
0791: */
0792: public double getAzimuth() throws IllegalStateException {
0793: if (!directionValid) {
0794: computeDirection();
0795: }
0796: return Math.toDegrees(azimuth);
0797: }
0798:
0799: /**
0800: * Returns the orthodromic distance. This method returns the value set by the last call to
0801: * <code>{@linkplain #setDirection(double,double) setDirection}(azimuth,distance)</code>,
0802: * <strong>except</strong> if <code>{@linkplain #setDestinationGeographicPoint(double,double)
0803: * setDestinationGeographicPoint}(...)</code> has been invoked after. In this later case, the
0804: * distance will be computed from the {@linkplain #getStartingGeographicPoint starting point}
0805: * to the destination point.
0806: *
0807: * @return The orthodromic distance, in the same units as the
0808: * {@linkplain #getEllipsoid ellipsoid} axis.
0809: * @throws IllegalStateException if the destination point has not been set.
0810: */
0811: public double getOrthodromicDistance() throws IllegalStateException {
0812: if (!directionValid) {
0813: computeDirection();
0814: assert checkOrthodromicDistance() : this ;
0815: }
0816: return distance;
0817: }
0818:
0819: /**
0820: * Computes the orthodromic distance using the algorithm implemented in the Geotools's
0821: * ellipsoid class (if available), and check if the error is smaller than some tolerance
0822: * error.
0823: */
0824: private boolean checkOrthodromicDistance() {
0825: if (ellipsoid instanceof DefaultEllipsoid) {
0826: double check;
0827: final DefaultEllipsoid ellipsoid = (DefaultEllipsoid) this .ellipsoid;
0828: check = ellipsoid.orthodromicDistance(
0829: Math.toDegrees(long1), Math.toDegrees(lat1), Math
0830: .toDegrees(long2), Math.toDegrees(lat2));
0831: check = Math.abs(distance - check);
0832: return check <= (distance + 1) * TOLERANCE_CHECK;
0833: }
0834: return true;
0835: }
0836:
0837: /**
0838: * Computes the destination point from the {@linkplain #getStartingGeographicPoint starting
0839: * point}, the {@linkplain #getAzimuth azimuth} and the {@linkplain #getOrthodromicDistance
0840: * orthodromic distance}.
0841: *
0842: * @throws IllegalStateException if the azimuth and the distance have not been set.
0843: *
0844: * @see #getDestinationGeographicPoint
0845: */
0846: private void computeDestinationPoint() throws IllegalStateException {
0847: if (!directionValid) {
0848: throw new IllegalStateException(Errors
0849: .format(ErrorKeys.DIRECTION_NOT_SET));
0850: }
0851: // Protect internal variables from changes
0852: final double lat1 = this .lat1;
0853: final double long1 = this .long1;
0854: final double azimuth = this .azimuth;
0855: final double distance = this .distance;
0856: /*
0857: * Solution of the geodetic direct problem after T.Vincenty.
0858: * Modified Rainsford's method with Helmert's elliptical terms.
0859: * Effective in any azimuth and at any distance short of antipodal.
0860: *
0861: * Latitudes and longitudes in radians positive North and East.
0862: * Forward azimuths at both points returned in radians from North.
0863: *
0864: * Programmed for CDC-6600 by LCDR L.Pfeifer NGS ROCKVILLE MD 18FEB75
0865: * Modified for IBM SYSTEM 360 by John G.Gergen NGS ROCKVILLE MD 7507
0866: * Ported from Fortran to Java by Daniele Franzoni.
0867: *
0868: * Source: ftp://ftp.ngs.noaa.gov/pub/pcsoft/for_inv.3d/source/forward.for
0869: * subroutine DIRECT1
0870: */
0871: double TU = fo * Math.sin(lat1) / Math.cos(lat1);
0872: double SF = Math.sin(azimuth);
0873: double CF = Math.cos(azimuth);
0874: double BAZ = (CF != 0) ? Math.atan2(TU, CF) * 2.0 : 0;
0875: double CU = 1 / Math.sqrt(TU * TU + 1.0);
0876: double SU = TU * CU;
0877: double SA = CU * SF;
0878: double C2A = 1.0 - SA * SA;
0879: double X = Math.sqrt((1.0 / fo / fo - 1) * C2A + 1.0) + 1.0;
0880: X = (X - 2.0) / X;
0881: double C = 1.0 - X;
0882: C = (X * X / 4.0 + 1.0) / C;
0883: double D = (0.375 * X * X - 1.0) * X;
0884: TU = distance / fo / semiMajorAxis / C;
0885: double Y = TU;
0886: double SY, CY, CZ, E;
0887: do {
0888: SY = Math.sin(Y);
0889: CY = Math.cos(Y);
0890: CZ = Math.cos(BAZ + Y);
0891: E = CZ * CZ * 2.0 - 1.0;
0892: C = Y;
0893: X = E * CY;
0894: Y = E + E - 1.0;
0895: Y = (((SY * SY * 4.0 - 3.0) * Y * CZ * D / 6.0 + X) * D
0896: / 4.0 - CZ)
0897: * SY * D + TU;
0898: } while (Math.abs(Y - C) > TOLERANCE_1);
0899: BAZ = CU * CY * CF - SU * SY;
0900: C = fo * Math.sqrt(SA * SA + BAZ * BAZ);
0901: D = SU * CY + CU * SY * CF;
0902: lat2 = Math.atan2(D, C);
0903: C = CU * CY - SU * SY * CF;
0904: X = Math.atan2(SY * SF, C);
0905: C = ((-3.0 * C2A + 4.0) * f + 4.0) * C2A * f / 16.0;
0906: D = ((E * CY * C + CZ) * SY * C + Y) * SA;
0907: long2 = long1 + X - (1.0 - C) * D * f;
0908: long2 = castToAngleRange(long2);
0909: destinationValid = true;
0910: }
0911:
0912: /**
0913: * Calculates the meridian arc length between two points in the same meridian
0914: * in the referenced ellipsoid.
0915: *
0916: * @param latitude1 The latitude of the first point (in decimal degrees).
0917: * @param latitude2 The latitude of the second point (in decimal degrees).
0918: * @return Returned the meridian arc length between latitude1 and latitude2
0919: */
0920: public double getMeridianArcLength(final double latitude1,
0921: final double latitude2) {
0922: return getMeridianArcLengthRadians(checkLatitude(latitude1),
0923: checkLatitude(latitude2));
0924: }
0925:
0926: /**
0927: * Calculates the meridian arc length between two points in the same meridian
0928: * in the referenced ellipsoid.
0929: *
0930: * @param P1 The latitude of the first point (in radians).
0931: * @param P2 The latitude of the second point (in radians).
0932: * @return Returned the meridian arc length between P1 and P2
0933: */
0934: private double getMeridianArcLengthRadians(final double P1,
0935: final double P2) {
0936: /*
0937: * Latitudes P1 and P2 in radians positive North and East.
0938: * Forward azimuths at both points returned in radians from North.
0939: *
0940: * Source: ftp://ftp.ngs.noaa.gov/pub/pcsoft/for_inv.3d/source/inverse.for
0941: * subroutine GPNARC
0942: * version 200005.26
0943: * written by Robert (Sid) Safford
0944: *
0945: * Ported from Fortran to Java by Daniele Franzoni.
0946: */
0947: double S1 = Math.abs(P1);
0948: double S2 = Math.abs(P2);
0949: double DA = (P2 - P1);
0950: // Check for a 90 degree lookup
0951: if (S1 > TOLERANCE_0 || S2 <= (Math.PI / 2 - TOLERANCE_0)
0952: || S2 >= (Math.PI / 2 + TOLERANCE_0)) {
0953: final double DB = Math.sin(P2 * 2.0) - Math.sin(P1 * 2.0);
0954: final double DC = Math.sin(P2 * 4.0) - Math.sin(P1 * 4.0);
0955: final double DD = Math.sin(P2 * 6.0) - Math.sin(P1 * 6.0);
0956: final double DE = Math.sin(P2 * 8.0) - Math.sin(P1 * 8.0);
0957: final double DF = Math.sin(P2 * 10.0) - Math.sin(P1 * 10.0);
0958: // Compute the S2 part of the series expansion
0959: S2 = -DB * B / 2.0 + DC * C / 4.0 - DD * D / 6.0 + DE * E
0960: / 8.0 - DF * F / 10.0;
0961: }
0962: // Compute the S1 part of the series expansion
0963: S1 = DA * A;
0964: // Compute the arc length
0965: return Math.abs(semiMajorAxis * (1.0 - eccentricitySquared)
0966: * (S1 + S2));
0967: }
0968:
0969: /**
0970: * Computes the azimuth and orthodromic distance from the
0971: * {@linkplain #getStartingGeographicPoint starting point} and the
0972: * {@linkplain #getDestinationGeographicPoint destination point}.
0973: *
0974: * @throws IllegalStateException if the destination point has not been set.
0975: *
0976: * @see #getAzimuth
0977: * @see #getOrthodromicDistance
0978: */
0979: private void computeDirection() throws IllegalStateException {
0980: if (!destinationValid) {
0981: throw new IllegalStateException(Errors
0982: .format(ErrorKeys.DESTINATION_NOT_SET));
0983: }
0984: // Protect internal variables from change.
0985: final double long1 = this .long1;
0986: final double lat1 = this .lat1;
0987: final double long2 = this .long2;
0988: final double lat2 = this .lat2;
0989: /*
0990: * Solution of the geodetic inverse problem after T.Vincenty.
0991: * Modified Rainsford's method with Helmert's elliptical terms.
0992: * Effective in any azimuth and at any distance short of antipodal.
0993: *
0994: * Latitudes and longitudes in radians positive North and East.
0995: * Forward azimuths at both points returned in radians from North.
0996: *
0997: * Programmed for CDC-6600 by LCDR L.Pfeifer NGS ROCKVILLE MD 18FEB75
0998: * Modified for IBM SYSTEM 360 by John G.Gergen NGS ROCKVILLE MD 7507
0999: * Ported from Fortran to Java by Daniele Franzoni.
1000: *
1001: * Source: ftp://ftp.ngs.noaa.gov/pub/pcsoft/for_inv.3d/source/inverse.for
1002: * subroutine GPNHRI
1003: * version 200208.09
1004: * written by robert (sid) safford
1005: */
1006: final double dlon = castToAngleRange(long2 - long1);
1007: final double ss = Math.abs(dlon);
1008: if (ss < TOLERANCE_1) {
1009: distance = getMeridianArcLengthRadians(lat1, lat2);
1010: azimuth = (lat2 > lat1) ? 0.0 : Math.PI;
1011: directionValid = true;
1012: return;
1013: }
1014: /*
1015: * Computes the limit in longitude (alimit), it is equal
1016: * to twice the distance from the equator to the pole,
1017: * as measured along the equator
1018: */
1019: // tests for antinodal difference
1020: final double ESQP = eccentricitySquared
1021: / (1.0 - eccentricitySquared);
1022: final double alimit = Math.PI * fo;
1023: if (ss >= alimit && lat1 < TOLERANCE_3 && lat1 > -TOLERANCE_3
1024: && lat2 < TOLERANCE_3 && lat2 > -TOLERANCE_3) {
1025: // Computes an approximate AZ
1026: final double CONS = (Math.PI - ss) / (Math.PI * f);
1027: double AZ = Math.asin(CONS);
1028: double AZ_TEMP, S, AO;
1029: int iter = 0;
1030: do {
1031: if (++iter > 8) {
1032: throw new ArithmeticException(
1033: getNoConvergenceErrorMessage());
1034: }
1035: S = Math.cos(AZ);
1036: final double C2 = S * S;
1037: // Compute new AO
1038: AO = T1 + T2 * C2 + T4 * C2 * C2 + T6 * C2 * C2 * C2;
1039: final double CS = CONS / AO;
1040: S = Math.asin(CS);
1041: AZ_TEMP = AZ;
1042: AZ = S;
1043: } while (Math.abs(S - AZ_TEMP) >= TOLERANCE_2);
1044:
1045: final double AZ1 = (dlon < 0.0) ? 2.0 * Math.PI - S : S;
1046: azimuth = castToAngleRange(AZ1);
1047: final double AZ2 = 2.0 * Math.PI - AZ1;
1048: S = Math.cos(AZ1);
1049:
1050: // Equatorial - geodesic(S-s) SMS
1051: final double U2 = ESQP * S * S;
1052: final double U4 = U2 * U2;
1053: final double U6 = U4 * U2;
1054: final double U8 = U6 * U2;
1055: final double BO = 1.0 + 0.25 * U2 + 0.046875 * U4
1056: + 0.01953125 * U6 + -0.01068115234375 * U8;
1057: S = Math.sin(AZ1);
1058: final double SMS = semiMajorAxis * Math.PI
1059: * (1.0 - f * Math.abs(S) * AO - BO * fo);
1060: distance = semiMajorAxis * ss - SMS;
1061: directionValid = true;
1062: return;
1063: }
1064:
1065: // the reduced latitudes
1066: final double u1 = Math.atan(fo * Math.sin(lat1)
1067: / Math.cos(lat1));
1068: final double u2 = Math.atan(fo * Math.sin(lat2)
1069: / Math.cos(lat2));
1070: final double su1 = Math.sin(u1);
1071: final double cu1 = Math.cos(u1);
1072: final double su2 = Math.sin(u2);
1073: final double cu2 = Math.cos(u2);
1074: double xy, w, q2, q4, q6, r2, r3, sig, ssig, slon, clon, sinalf, ab = dlon;
1075: int kcount = 0;
1076: do {
1077: if (++kcount > 8) {
1078: throw new ArithmeticException(
1079: getNoConvergenceErrorMessage());
1080: }
1081: clon = Math.cos(ab);
1082: slon = Math.sin(ab);
1083: final double csig = su1 * su2 + cu1 * cu2 * clon;
1084: ssig = Math.sqrt(slon * cu2 * slon * cu2
1085: + (su2 * cu1 - su1 * cu2 * clon)
1086: * (su2 * cu1 - su1 * cu2 * clon));
1087: sig = Math.atan2(ssig, csig);
1088: sinalf = cu1 * cu2 * slon / ssig;
1089: w = (1.0 - sinalf * sinalf);
1090: final double t4 = w * w;
1091: final double t6 = w * t4;
1092:
1093: // the coefficents of type a
1094: final double ao = f + a01 * w + a02 * t4 + a03 * t6;
1095: final double a2 = a21 * w + a22 * t4 + a23 * t6;
1096: final double a4 = a42 * t4 + a43 * t6;
1097: final double a6 = a63 * t6;
1098:
1099: // the multiple angle functions
1100: double qo = 0.0;
1101: if (w > TOLERANCE_0) {
1102: qo = -2.0 * su1 * su2 / w;
1103: }
1104: q2 = csig + qo;
1105: q4 = 2.0 * q2 * q2 - 1.0;
1106: q6 = q2 * (4.0 * q2 * q2 - 3.0);
1107: r2 = 2.0 * ssig * csig;
1108: r3 = ssig * (3.0 - 4.0 * ssig * ssig);
1109:
1110: // the longitude difference
1111: final double s = sinalf
1112: * (ao * sig + a2 * ssig * q2 + a4 * r2 * q4 + a6
1113: * r3 * q6);
1114: double xz = dlon + s;
1115: xy = Math.abs(xz - ab);
1116: ab = dlon + s;
1117: } while (xy >= TOLERANCE_1);
1118:
1119: final double z = ESQP * w;
1120: final double bo = 1.0
1121: + z
1122: * (1.0 / 4.0 + z
1123: * (-3.0 / 64.0 + z
1124: * (5.0 / 256.0 - z * (175.0 / 16384.0))));
1125: final double b2 = z
1126: * (-1.0 / 4.0 + z
1127: * (1.0 / 16.0 + z
1128: * (-15.0 / 512.0 + z * (35.0 / 2048.0))));
1129: final double b4 = z
1130: * z
1131: * (-1.0 / 128.0 + z
1132: * (3.0 / 512.0 - z * (35.0 / 8192.0)));
1133: final double b6 = z * z * z
1134: * (-1.0 / 1536.0 + z * (5.0 / 6144.0));
1135:
1136: // The distance in ellispoid axis units.
1137: distance = semiMinorAxis
1138: * (bo * sig + b2 * ssig * q2 + b4 * r2 * q4 + b6 * r3
1139: * q6);
1140: double az1 = (dlon < 0) ? Math.PI * (3 / 2) : Math.PI / 2;
1141:
1142: // now compute the az1 & az2 for latitudes not on the equator
1143: if ((Math.abs(su1) >= TOLERANCE_0)
1144: || (Math.abs(su2) >= TOLERANCE_0)) {
1145: final double tana1 = slon * cu2
1146: / (su2 * cu1 - clon * su1 * cu2);
1147: final double sina1 = sinalf / cu1;
1148:
1149: // azimuths from north,longitudes positive east
1150: az1 = Math.atan2(sina1, sina1 / tana1);
1151: }
1152: azimuth = castToAngleRange(az1);
1153: directionValid = true;
1154: return;
1155: }
1156:
1157: /**
1158: * Calculates the geodetic curve between two points in the referenced ellipsoid.
1159: * A curve in the ellipsoid is a path which points contain the longitude and latitude
1160: * of the points in the geodetic curve. The geodetic curve is computed from the
1161: * {@linkplain #getStartingGeographicPoint starting point} to the
1162: * {@linkplain #getDestinationGeographicPoint destination point}.
1163: *
1164: * @param numberOfPoints The number of vertex in the geodetic curve.
1165: * <strong>NOTE:</strong> This argument is only a hint and may be ignored
1166: * in future version (if we compute a real curve rather than a list of line
1167: * segments).
1168: * @return The path that represents the geodetic curve from the
1169: * {@linkplain #getStartingGeographicPoint starting point} to the
1170: * {@linkplain #getDestinationGeographicPoint destination point}.
1171: *
1172: * @todo We should check for cases where the path cross the 90°N, 90°S, 90°E or 90°W boundaries.
1173: */
1174: public Shape getGeodeticCurve(final int numberOfPoints) {
1175: checkNumberOfPoints(numberOfPoints);
1176: if (!directionValid) {
1177: computeDirection();
1178: }
1179: if (!destinationValid) {
1180: computeDestinationPoint();
1181: }
1182: final double long2 = this .long2;
1183: final double lat2 = this .lat2;
1184: final double distance = this .distance;
1185: final double deltaDistance = distance / (numberOfPoints + 1);
1186: final GeneralPath path = new GeneralPath(
1187: GeneralPath.WIND_EVEN_ODD, numberOfPoints + 1);
1188: path.moveTo((float) Math.toDegrees(long1), (float) Math
1189: .toDegrees(lat1));
1190: for (int i = 1; i < numberOfPoints; i++) {
1191: this .distance = i * deltaDistance;
1192: computeDestinationPoint();
1193: path.lineTo((float) Math.toDegrees(this .long2),
1194: (float) Math.toDegrees(this .lat2));
1195: }
1196: this .long2 = long2;
1197: this .lat2 = lat2;
1198: this .distance = distance;
1199: path.lineTo((float) Math.toDegrees(long2), (float) Math
1200: .toDegrees(lat2));
1201: return path;
1202: }
1203:
1204: /**
1205: * Calculates the geodetic curve between two points in the referenced ellipsoid.
1206: * A curve in the ellipsoid is a path which points contain the longitude and latitude
1207: * of the points in the geodetic curve. The geodetic curve is computed from the
1208: * {@linkplain #getStartingGeographicPoint starting point} to the
1209: * {@linkplain #getDestinationGeographicPoint destination point}.
1210: *
1211: * @return The path that represents the geodetic curve from the
1212: * {@linkplain #getStartingGeographicPoint starting point} to the
1213: * {@linkplain #getDestinationGeographicPoint destination point}.
1214: */
1215: public Shape getGeodeticCurve() {
1216: return getGeodeticCurve(10);
1217: }
1218:
1219: /**
1220: * Calculates the loxodromic curve between two points in the referenced ellipsoid.
1221: * The loxodromic curve between two points is a path that links together the two
1222: * points with a constant azimuth. The azimuth from every points of the loxodromic
1223: * curve and the second point is constant.
1224: *
1225: * @return The path that represents the loxodromic curve from the
1226: * {@linkplain #getStartingGeographicPoint starting point} to the
1227: * {@linkplain #getDestinationGeographicPoint destination point}.
1228: */
1229: private Shape getLoxodromicCurve() {
1230: if (true) {
1231: throw new UnsupportedOperationException();
1232: }
1233: /*************************************************************************************
1234: ** THE FOLLOWING IS CHECKED FOR COMPILER ERROR, BUT EXCLUDED FROM THE .class FILE. **
1235: ** THIS CODE IS WRONG: LOXODROMIC CURVES ARE STRAIGHT LINES IN MERCATOR PROJECTION, **
1236: ** NOT IT PLAIN (longitude,latitude) SPACE. FURTHERMORE, THE "OUT OF BOUNDS" CHECK **
1237: ** IS UNFINISHED: WHEN THE PATH CROSS THE 180° LONGITUDE, A +360° ADDITION NEED TO **
1238: ** BE PERFORMED ON ONE OF THE SOURCE OR TARGET POINT BEFORE TO COMPUTE THE LINEAR **
1239: ** INTERPOLATION (OTHERWISE, THE SLOPE VALUE IS WRONG). FORMULAS FOR COMPUTING MID- **
1240: ** POINT ON A LOXODROMIC CURVE ARE AVAILABLE THERE: **
1241: ** **
1242: ** http://mathforum.org/discuss/sci.math/a/t/180912 **
1243: ** **
1244: ** LatM = (Lat0+Lat1)/2 **
1245: ** **
1246: ** (Lon1-Lon0)log(f(LatM)) + Lon0 log(f(Lat1)) - Lon1 log(f(Lat0)) **
1247: ** LonM = --------------------------------------------------------------- **
1248: ** log(f(Lat1)/f(Lat0)) **
1249: ** **
1250: ** where log(f(x)) == log(sec(x)+tan(x)) is the inverse Gudermannian function. **
1251: *************************************************************************************/
1252: if (!directionValid) {
1253: computeDirection();
1254: }
1255: if (!destinationValid) {
1256: computeDestinationPoint();
1257: }
1258: final double x1 = Math.toDegrees(long1);
1259: final double y1 = Math.toDegrees(lat1);
1260: final double x2 = Math.toDegrees(long2);
1261: final double y2 = Math.toDegrees(lat2);
1262: /*
1263: * Check if the azimuth is heading from P1 to P2 (TRUE) or in the opposite direction
1264: * (FALSE). Horizontal (X) and vertical (Y) components are checked separatly. A null
1265: * value means "don't know", because the path is perfectly vertical or horizontal or
1266: * because a coordinate is NaN. If both components are not null (unknow), then they
1267: * must be consistent.
1268: */
1269: final Boolean xDirect = (x2 > x1) ? Boolean
1270: .valueOf(azimuth >= 0) : (x2 < x1) ? Boolean
1271: .valueOf(azimuth <= 0) : null;
1272: final Boolean yDirect = (y2 > y1) ? Boolean
1273: .valueOf(azimuth >= -90 && azimuth <= +90)
1274: : (y2 < y1) ? Boolean.valueOf(azimuth <= -90
1275: || azimuth >= +90) : null;
1276: assert xDirect == null || yDirect == null
1277: || xDirect.equals(yDirect) : this ;
1278: if (!Boolean.FALSE.equals(xDirect)
1279: && !Boolean.FALSE.equals(yDirect)) {
1280: return new Line2D.Double(x1, y1, x2, y2);
1281: }
1282: if (Boolean.FALSE.equals(yDirect)) {
1283: /*
1284: * Crossing North or South pole is more complicated than what we do for now: If we
1285: * follow the 0° longitude toward North, then we have to follow the 180° longitude
1286: * from North to South pole and follow the 0° longitude again toward North up to
1287: * the destination point.
1288: */
1289: throw new UnsupportedOperationException(
1290: "Crossing pole is not yet implemented");
1291: }
1292: /*
1293: * The azimuth is heading in the opposite direction of the path from P1 to P2. Computes
1294: * the intersection points at the 90°N / 90°S boundaries, or the 180°E / 180°W boundaries.
1295: * (xout,yout) is the point where the path goes out (initialized to the corner where the
1296: * azimuth is heading); (xin,yin) is the point where the path come back in the opposite
1297: * hemisphere.
1298: */
1299: double xout = (x2 >= x1) ? -180 : +180;
1300: double yout = (y2 >= y1) ? -90 : +90;
1301: double xin = -xout;
1302: double yin = -yout;
1303: final double dx = x2 - x1;
1304: final double dy = y2 - y1;
1305: if (dx == 0) {
1306: xin = xout = x1; // Vertical line.
1307: } else if (dy == 0) {
1308: yin = yout = y1; // Horizontal line.
1309: } else {
1310: /*
1311: * The path is diagonal (neither horizontal or vertical). The following loop
1312: * is executed exactly twice: the first pass computes the "out" point, and
1313: * the second pass computes the "in" point. Each pass computes actually two
1314: * points: the intersection point against the 180°W or 180°E boundary, and
1315: * the intersection point against the 90°N or 90°S boundary. Usually one of
1316: * those points will be out of range and the other one is selected.
1317: */
1318: boolean in = false;
1319: do {
1320: final double meridX, meridY; // The point where the path cross the +/-180° meridian.
1321: final double zonalX, zonalY; // The point where the path cross the +/- 90° parallel.
1322: meridX = in ? xin : xout;
1323: meridY = dy / dx * (meridX - x1) + y1;
1324: zonalY = in ? yin : yout;
1325: zonalX = dx / dy * (zonalY - y1) + x1;
1326: if (Math.abs(meridY) < Math.abs(zonalX) * 0.5) {
1327: if (in) {
1328: xin = meridX;
1329: yin = meridY;
1330: } else {
1331: xout = meridX;
1332: yout = meridY;
1333: }
1334: } else {
1335: if (in) {
1336: xin = zonalX;
1337: yin = zonalY;
1338: } else {
1339: xout = zonalX;
1340: yout = zonalY;
1341: }
1342: }
1343: } while ((in = !in) == false);
1344: }
1345: final GeneralPath path = new GeneralPath(
1346: GeneralPath.WIND_EVEN_ODD, 4);
1347: path.moveTo((float) x1, (float) y1);
1348: path.lineTo((float) xout, (float) yout);
1349: path.moveTo((float) xin, (float) yin);
1350: path.lineTo((float) x2, (float) y2);
1351: return path;
1352: }
1353:
1354: /**
1355: * Returns a string representation of the current state of this calculator.
1356: */
1357: public String toString() {
1358: final Vocabulary resources = Vocabulary.getResources(null);
1359: final TableWriter buffer = new TableWriter(null, " ");
1360: if (coordinateReferenceSystem != null) {
1361: buffer
1362: .write(resources
1363: .getLabel(VocabularyKeys.COORDINATE_REFERENCE_SYSTEM));
1364: buffer.nextColumn();
1365: buffer.write(coordinateReferenceSystem.getName().getCode());
1366: buffer.nextLine();
1367: }
1368: if (ellipsoid != null) {
1369: buffer.write(resources.getLabel(VocabularyKeys.ELLIPSOID));
1370: buffer.nextColumn();
1371: buffer.write(ellipsoid.getName().getCode());
1372: buffer.nextLine();
1373: }
1374: final CoordinateFormat cf = new CoordinateFormat();
1375: final Format nf = cf.getFormat(0);
1376: if (true) {
1377: buffer.write(resources
1378: .getLabel(VocabularyKeys.SOURCE_POINT));
1379: buffer.nextColumn();
1380: buffer.write(format(cf, long1, lat1));
1381: buffer.nextLine();
1382: }
1383: if (destinationValid) {
1384: buffer.write(resources
1385: .getLabel(VocabularyKeys.TARGET_POINT));
1386: buffer.nextColumn();
1387: buffer.write(format(cf, long2, lat2));
1388: buffer.nextLine();
1389: }
1390: if (directionValid) {
1391: buffer.write(resources.getLabel(VocabularyKeys.AZIMUTH));
1392: buffer.nextColumn();
1393: buffer.write(nf.format(new Angle(Math.toDegrees(azimuth))));
1394: buffer.nextLine();
1395: }
1396: if (directionValid) {
1397: buffer.write(resources
1398: .getLabel(VocabularyKeys.ORTHODROMIC_DISTANCE));
1399: buffer.nextColumn();
1400: buffer.write(nf.format(new Double(distance)));
1401: if (ellipsoid != null) {
1402: buffer.write(' ');
1403: buffer.write(ellipsoid.getAxisUnit().toString());
1404: }
1405: buffer.nextLine();
1406: }
1407: return buffer.toString();
1408: }
1409: }
|