/** * Implementation of the net.sf.geographiclib.Geodesic class * * Copyright (c) Charles Karney (2013-2019) and licensed * under the MIT/X11 License. For more information, see * https://geographiclib.sourceforge.io/ **********************************************************************/ package net.sf.geographiclib; /** * Geodesic calculations. *

* The shortest path between two points on a ellipsoid at (lat1, * lon1) and (lat2, lon2) is called the geodesic. Its * length is s12 and the geodesic from point 1 to point 2 has azimuths * azi1 and azi2 at the two end points. (The azimuth is the * heading measured clockwise from north. azi2 is the "forward" * azimuth, i.e., the heading that takes you beyond point 2 not back to point * 1.) *

* Given lat1, lon1, azi1, and s12, we can * determine lat2, lon2, and azi2. This is the * direct geodesic problem and its solution is given by the function * {@link #Direct Direct}. (If s12 is sufficiently large that the * geodesic wraps more than halfway around the earth, there will be another * geodesic between the points with a smaller s12.) *

* Given lat1, lon1, lat2, and lon2, we can * determine azi1, azi2, and s12. This is the * inverse geodesic problem, whose solution is given by {@link #Inverse * Inverse}. Usually, the solution to the inverse problem is unique. In cases * where there are multiple solutions (all with the same s12, of * course), all the solutions can be easily generated once a particular * solution is provided. *

* The standard way of specifying the direct problem is the specify the * distance s12 to the second point. However it is sometimes useful * instead to specify the arc length a12 (in degrees) on the auxiliary * sphere. This is a mathematical construct used in solving the geodesic * problems. The solution of the direct problem in this form is provided by * {@link #ArcDirect ArcDirect}. An arc length in excess of 180° indicates * that the geodesic is not a shortest path. In addition, the arc length * between an equatorial crossing and the next extremum of latitude for a * geodesic is 90°. *

* This class can also calculate several other quantities related to * geodesics. These are: *

* reduced length. If we fix the first point and increase * azi1 by dazi1 (radians), the second point is displaced * m12 dazi1 in the direction azi2 + 90°. The * quantity m12 is called the "reduced length" and is symmetric under * interchange of the two points. On a curved surface the reduced length * obeys a symmetry relation, m12 + m21 = 0. On a flat * surface, we have m12 = s12. The ratio s12/m12 * gives the azimuthal scale for an azimuthal equidistant projection. *
* geodesic scale. Consider a reference geodesic and a second * geodesic parallel to this one at point 1 and separated by a small distance * dt. The separation of the two geodesics at point 2 is M12 * dt where M12 is called the "geodesic scale". M21 is * defined similarly (with the geodesics being parallel at point 2). On a * flat surface, we have M12 = M21 = 1. The quantity * 1/M12 gives the scale of the Cassini-Soldner projection. *
* area. The area between the geodesic from point 1 to point 2 and * the equation is represented by S12; it is the area, measured * counter-clockwise, of the geodesic quadrilateral with corners * (lat1,lon1), (0,lon1), (0,lon2), and * (lat2,lon2). It can be used to compute the area of any * geodesic polygon. *

* The quantities m12, M12, M21 which all specify the * behavior of nearby geodesics obey addition rules. If points 1, 2, and 3 all * lie on a single geodesic, then the following rules hold: *

* s13 = s12 + s23 *
* a13 = a12 + a23 *
* S13 = S12 + S23 *
* m13 = m12 M23 + m23 M21 *
* M13 = M12 M23 − (1 − M12 * M21) m23 / m12 *
* M31 = M32 M21 − (1 − M23 * M32) m12 / m23 *

* The results of the geodesic calculations are bundled up into a {@link * GeodesicData} object which includes the input parameters and all the * computed results, i.e., lat1, lon1, azi1, lat2, * lon2, azi2, s12, a12, m12, M12, * M21, S12. *

* The functions {@link #Direct(double, double, double, double, int) Direct}, * {@link #ArcDirect(double, double, double, double, int) ArcDirect}, and * {@link #Inverse(double, double, double, double, int) Inverse} include an * optional final argument outmask which allows you specify which * results should be computed and returned. If you omit outmask, then * the "standard" geodesic results are computed (latitudes, longitudes, * azimuths, and distance). outmask is bitor'ed combination of {@link * GeodesicMask} values. For example, if you wish just to compute the distance * between two points you would call, e.g., *

 * {@code
 *  GeodesicData g = Geodesic.WGS84.Inverse(lat1, lon1, lat2, lon2,
 *                      GeodesicMask.DISTANCE); }

* Additional functionality is provided by the {@link GeodesicLine} class, * which allows a sequence of points along a geodesic to be computed. *

* The shortest distance returned by the solution of the inverse problem is * (obviously) uniquely defined. However, in a few special cases there are * multiple azimuths which yield the same shortest distance. Here is a * catalog of those cases: *

* lat1 = −lat2 (with neither point at a pole). If * azi1 = azi2, the geodesic is unique. Otherwise there are * two geodesics and the second one is obtained by setting [azi1, * azi2] → [azi2, azi1], [M12, M21] * → [M21, M12], S12 → −S12. * (This occurs when the longitude difference is near ±180° for * oblate ellipsoids.) *
* lon2 = lon1 ± 180° (with neither point at a * pole). If azi1 = 0° or ±180°, the geodesic is * unique. Otherwise there are two geodesics and the second one is obtained * by setting [ azi1, azi2] → [−azi1, * −azi2], S12 → − S12. (This occurs * when lat2 is near −lat1 for prolate ellipsoids.) *
* Points 1 and 2 at opposite poles. There are infinitely many geodesics * which can be generated by setting [azi1, azi2] → * [azi1, azi2] + [d, −d], for arbitrary * d. (For spheres, this prescription applies when points 1 and 2 are * antipodal.) *
* s12 = 0 (coincident points). There are infinitely many geodesics * which can be generated by setting [azi1, azi2] → * [azi1, azi2] + [d, d], for arbitrary d. *

* The calculations are accurate to better than 15 nm (15 nanometers) for the * WGS84 ellipsoid. See Sec. 9 of * arXiv:1102.1215v1 for * details. The algorithms used by this class are based on series expansions * using the flattening f as a small parameter. These are only accurate * for |f| < 0.02; however reasonably accurate results will be * obtained for |f| < 0.2. Here is a table of the approximate * maximum error (expressed as a distance) for an ellipsoid with the same * equatorial radius as the WGS84 ellipsoid and different values of the * flattening.

 *     |f|      error
 *     0.01     25 nm
 *     0.02     30 nm
 *     0.05     10 um
 *     0.1     1.5 mm
 *     0.2     300 mm

* The algorithms are described in *

C. F. F. Karney, * * Algorithms for geodesics, * J. Geodesy 87, 43–55 (2013) * (addenda). *

* Example of use: *

 * {@code
 * // Solve the direct geodesic problem.
 *
 * // This program reads in lines with lat1, lon1, azi1, s12 and prints
 * // out lines with lat2, lon2, azi2 (for the WGS84 ellipsoid).
 *
 * import java.util.*;
 * import net.sf.geographiclib.*;
 * public class Direct {
 *   public static void main(String[] args) {
 *     try {
 *       Scanner in = new Scanner(System.in);
 *       double lat1, lon1, azi1, s12;
 *       while (true) {
 *         lat1 = in.nextDouble(); lon1 = in.nextDouble();
 *         azi1 = in.nextDouble(); s12 = in.nextDouble();
 *         GeodesicData g = Geodesic.WGS84.Direct(lat1, lon1, azi1, s12);
 *         System.out.println(g.lat2 + " " + g.lon2 + " " + g.azi2);
 *       }
 *     }
 *     catch (Exception e) {}
 *   }
 * }}

**********************************************************************/ public class Geodesic { /** * The order of the expansions used by Geodesic. **********************************************************************/ protected static final int GEOGRAPHICLIB_GEODESIC_ORDER = 6; protected static final int nA1_ = GEOGRAPHICLIB_GEODESIC_ORDER; protected static final int nC1_ = GEOGRAPHICLIB_GEODESIC_ORDER; protected static final int nC1p_ = GEOGRAPHICLIB_GEODESIC_ORDER; protected static final int nA2_ = GEOGRAPHICLIB_GEODESIC_ORDER; protected static final int nC2_ = GEOGRAPHICLIB_GEODESIC_ORDER; protected static final int nA3_ = GEOGRAPHICLIB_GEODESIC_ORDER; protected static final int nA3x_ = nA3_; protected static final int nC3_ = GEOGRAPHICLIB_GEODESIC_ORDER; protected static final int nC3x_ = (nC3_ * (nC3_ - 1)) / 2; protected static final int nC4_ = GEOGRAPHICLIB_GEODESIC_ORDER; protected static final int nC4x_ = (nC4_ * (nC4_ + 1)) / 2; private static final int maxit1_ = 20; private static final int maxit2_ = maxit1_ + GeoMath.digits + 10; // Underflow guard. We require // tiny_ * epsilon() > 0 // tiny_ + epsilon() == epsilon() protected static final double tiny_ = Math.sqrt(Double.MIN_NORMAL); private static final double tol0_ = Math.ulp(1.0); // Increase multiplier in defn of tol1_ from 100 to 200 to fix inverse case // 52.784459512564 0 -52.784459512563990912 179.634407464943777557 // which otherwise failed for Visual Studio 10 (Release and Debug) private static final double tol1_ = 200 * tol0_; private static final double tol2_ = Math.sqrt(tol0_); // Check on bisection interval private static final double tolb_ = tol0_ * tol2_; private static final double xthresh_ = 1000 * tol2_; protected double _a, _f, _f1, _e2, _ep2, _b, _c2; private double _n, _etol2; private double _A3x[], _C3x[], _C4x[]; /** * Constructor for a ellipsoid with *

* @param a equatorial radius (meters). * @param f flattening of ellipsoid. Setting f = 0 gives a sphere. * Negative f gives a prolate ellipsoid. * @exception GeographicErr if a or (1 − f ) a is * not positive. **********************************************************************/ public Geodesic(double a, double f) { _a = a; _f = f; _f1 = 1 - _f; _e2 = _f * (2 - _f); _ep2 = _e2 / GeoMath.sq(_f1); // e2 / (1 - e2) _n = _f / ( 2 - _f); _b = _a * _f1; _c2 = (GeoMath.sq(_a) + GeoMath.sq(_b) * (_e2 == 0 ? 1 : (_e2 > 0 ? GeoMath.atanh(Math.sqrt(_e2)) : Math.atan(Math.sqrt(-_e2))) / Math.sqrt(Math.abs(_e2))))/2; // authalic radius squared // The sig12 threshold for "really short". Using the auxiliary sphere // solution with dnm computed at (bet1 + bet2) / 2, the relative error in // the azimuth consistency check is sig12^2 * abs(f) * min(1, 1-f/2) / 2. // (Error measured for 1/100 < b/a < 100 and abs(f) >= 1/1000. For a // given f and sig12, the max error occurs for lines near the pole. If // the old rule for computing dnm = (dn1 + dn2)/2 is used, then the error // increases by a factor of 2.) Setting this equal to epsilon gives // sig12 = etol2. Here 0.1 is a safety factor (error decreased by 100) // and max(0.001, abs(f)) stops etol2 getting too large in the nearly // spherical case. _etol2 = 0.1 * tol2_ / Math.sqrt( Math.max(0.001, Math.abs(_f)) * Math.min(1.0, 1 - _f/2) / 2 ); if (!(GeoMath.isfinite(_a) && _a > 0)) throw new GeographicErr("Equatorial radius is not positive"); if (!(GeoMath.isfinite(_b) && _b > 0)) throw new GeographicErr("Polar semi-axis is not positive"); _A3x = new double[nA3x_]; _C3x = new double[nC3x_]; _C4x = new double[nC4x_]; A3coeff(); C3coeff(); C4coeff(); } /** * Solve the direct geodesic problem where the length of the geodesic * is specified in terms of distance. *

* @param lat1 latitude of point 1 (degrees). * @param lon1 longitude of point 1 (degrees). * @param azi1 azimuth at point 1 (degrees). * @param s12 distance between point 1 and point 2 (meters); it can be * negative. * @return a {@link GeodesicData} object with the following fields: * lat1, lon1, azi1, lat2, lon2, * azi2, s12, a12. *

* lat1 should be in the range [−90°, 90°]. The values * of lon2 and azi2 returned are in the range [−180°, * 180°]. *

* If either point is at a pole, the azimuth is defined by keeping the * longitude fixed, writing lat = ±(90° − ε), * and taking the limit ε → 0+. An arc length greater that * 180° signifies a geodesic which is not a shortest path. (For a * prolate ellipsoid, an additional condition is necessary for a shortest * path: the longitudinal extent must not exceed of 180°.) **********************************************************************/ public GeodesicData Direct(double lat1, double lon1, double azi1, double s12) { return Direct(lat1, lon1, azi1, false, s12, GeodesicMask.STANDARD); } /** * Solve the direct geodesic problem where the length of the geodesic is * specified in terms of distance and with a subset of the geodesic results * returned. *

* @param lat1 latitude of point 1 (degrees). * @param lon1 longitude of point 1 (degrees). * @param azi1 azimuth at point 1 (degrees). * @param s12 distance between point 1 and point 2 (meters); it can be * negative. * @param outmask a bitor'ed combination of {@link GeodesicMask} values * specifying which results should be returned. * @return a {@link GeodesicData} object with the fields specified by * outmask computed. *

* lat1, lon1, azi1, s12, and a12 are * always included in the returned result. The value of lon2 returned * is in the range [−180°, 180°], unless the outmask * includes the {@link GeodesicMask#LONG_UNROLL} flag. **********************************************************************/ public GeodesicData Direct(double lat1, double lon1, double azi1, double s12, int outmask) { return Direct(lat1, lon1, azi1, false, s12, outmask); } /** * Solve the direct geodesic problem where the length of the geodesic * is specified in terms of arc length. *

* @param lat1 latitude of point 1 (degrees). * @param lon1 longitude of point 1 (degrees). * @param azi1 azimuth at point 1 (degrees). * @param a12 arc length between point 1 and point 2 (degrees); it can * be negative. * @return a {@link GeodesicData} object with the following fields: * lat1, lon1, azi1, lat2, lon2, * azi2, s12, a12. *

* lat1 should be in the range [−90°, 90°]. The values * of lon2 and azi2 returned are in the range [−180°, * 180°]. *

* If either point is at a pole, the azimuth is defined by keeping the * longitude fixed, writing lat = ±(90° − ε), * and taking the limit ε → 0+. An arc length greater that * 180° signifies a geodesic which is not a shortest path. (For a * prolate ellipsoid, an additional condition is necessary for a shortest * path: the longitudinal extent must not exceed of 180°.) **********************************************************************/ public GeodesicData ArcDirect(double lat1, double lon1, double azi1, double a12) { return Direct(lat1, lon1, azi1, true, a12, GeodesicMask.STANDARD); } /** * Solve the direct geodesic problem where the length of the geodesic is * specified in terms of arc length and with a subset of the geodesic results * returned. *

* @param lat1 latitude of point 1 (degrees). * @param lon1 longitude of point 1 (degrees). * @param azi1 azimuth at point 1 (degrees). * @param a12 arc length between point 1 and point 2 (degrees); it can * be negative. * @param outmask a bitor'ed combination of {@link GeodesicMask} values * specifying which results should be returned. * @return a {@link GeodesicData} object with the fields specified by * outmask computed. *

* lat1, lon1, azi1, and a12 are always included * in the returned result. The value of lon2 returned is in the range * [−180°, 180°], unless the outmask includes the {@link * GeodesicMask#LONG_UNROLL} flag. **********************************************************************/ public GeodesicData ArcDirect(double lat1, double lon1, double azi1, double a12, int outmask) { return Direct(lat1, lon1, azi1, true, a12, outmask); } /** * The general direct geodesic problem. {@link #Direct Direct} and * {@link #ArcDirect ArcDirect} are defined in terms of this function. *

* @param lat1 latitude of point 1 (degrees). * @param lon1 longitude of point 1 (degrees). * @param azi1 azimuth at point 1 (degrees). * @param arcmode boolean flag determining the meaning of the * s12_a12. * @param s12_a12 if arcmode is false, this is the distance between * point 1 and point 2 (meters); otherwise it is the arc length between * point 1 and point 2 (degrees); it can be negative. * @param outmask a bitor'ed combination of {@link GeodesicMask} values * specifying which results should be returned. * @return a {@link GeodesicData} object with the fields specified by * outmask computed. *

* The {@link GeodesicMask} values possible for outmask are *

* outmask |= {@link GeodesicMask#LATITUDE} for the latitude * lat2; *
* outmask |= {@link GeodesicMask#LONGITUDE} for the latitude * lon2; *
* outmask |= {@link GeodesicMask#AZIMUTH} for the latitude * azi2; *
* outmask |= {@link GeodesicMask#DISTANCE} for the distance * s12; *
* outmask |= {@link GeodesicMask#REDUCEDLENGTH} for the reduced * length m12; *
* outmask |= {@link GeodesicMask#GEODESICSCALE} for the geodesic * scales M12 and M21; *
* outmask |= {@link GeodesicMask#AREA} for the area S12; *
* outmask |= {@link GeodesicMask#ALL} for all of the above; *
* outmask |= {@link GeodesicMask#LONG_UNROLL}, if set then * lon1 is unchanged and lon2 − lon1 indicates * how many times and in what sense the geodesic encircles the ellipsoid. * Otherwise lon1 and lon2 are both reduced to the range * [−180°, 180°]. *

* The function value a12 is always computed and returned and this * equals s12_a12 is arcmode is true. If outmask * includes {@link GeodesicMask#DISTANCE} and arcmode is false, then * s12 = s12_a12. It is not necessary to include {@link * GeodesicMask#DISTANCE_IN} in outmask; this is automatically * included is arcmode is false. **********************************************************************/ public GeodesicData Direct(double lat1, double lon1, double azi1, boolean arcmode, double s12_a12, int outmask) { // Automatically supply DISTANCE_IN if necessary if (!arcmode) outmask |= GeodesicMask.DISTANCE_IN; return new GeodesicLine(this, lat1, lon1, azi1, outmask) . // Note the dot! Position(arcmode, s12_a12, outmask); } /** * Define a {@link GeodesicLine} in terms of the direct geodesic problem * specified in terms of distance with all capabilities included. *