sf/geographiclib/Gnomonic.java

/**
 * Implementation of the net.sf.geographiclib.Gnomonic class
 *
 * Copyright (c) BMW Car IT GmbH (2014-2019) <sebastian.mattheis@bmw-carit.de>
 * and Charles Karney (2020) <charles@karney.com>∑ and licensed and licensed
 * under the MIT/X11 License. For more information, see
 * https://geographiclib.sourceforge.io/
 **********************************************************************/
package net.sf.geographiclib;

/**
 * Gnomonic projection.
 * <p>
 * <i>Note: Gnomonic.java has been ported to Java from its C++ equivalent
 * Gnomonic.cpp, authored by C. F. F. Karney and licensed under MIT/X11
 * license.  The following documentation is mostly the same as for its C++
 * equivalent, but has been adopted to apply to this Java implementation.</i>
 * <p>
 * Gnomonic projection centered at an arbitrary position <i>C</i> on the
 * ellipsoid. This projection is derived in Section 8 of
 * <ul>
 * <li>
 * C. F. F. Karney, <a href="https://doi.org/10.1007/s00190-012-0578-z">
 * Algorithms for geodesics</a>, J. Geodesy <b>87</b>, 43&ndash;55 (2013);
 * DOI: <a href="https://doi.org/10.1007/s00190-012-0578-z">
 * 10.1007/s00190-012-0578-z</a>; addenda:
 * <a href="https://geographiclib.sourceforge.io/geod-addenda.html">
 * geod-addenda.html</a>.
 * </li>
 * </ul>
 * <p>
 * The gnomonic projection of a point <i>P</i> on the ellipsoid is defined as
 * follows: compute the geodesic line from <i>C</i> to <i>P</i>; compute the
 * reduced length <i>m12</i>, geodesic scale <i>M12</i>, and &rho; =
 * <i>m12</i>/<i>M12</i>; finally, this gives the coordinates <i>x</i> and
 * <i>y</i> of <i>P</i> in gnomonic projection with <i>x</i> = &rho; sin
 * <i>azi1</i>; <i>y</i> = &rho; cos <i>azi1</i>, where <i>azi1</i> is the
 * azimuth of the geodesic at <i>C</i>. The method
 * {@link Gnomonic#Forward(double, double, double, double)} performs the
 * forward projection and
 * {@link Gnomonic#Reverse(double, double, double, double)} is the
 * inverse of the projection. The methods also return the azimuth
 * <i>azi</i> of the geodesic at <i>P</i> and reciprocal scale
 * <i>rk</i> in the azimuthal direction. The scale in the radial
 * direction is 1/<i>rk</i><sup>2</sup>.
 * <p>
 * For a sphere, &rho; reduces to <i>a</i> tan(<i>s12</i>/<i>a</i>), where
 * <i>s12</i> is the length of the geodesic from <i>C</i> to <i>P</i>, and the
 * gnomonic projection has the property that all geodesics appear as straight
 * lines. For an ellipsoid, this property holds only for geodesics interesting
 * the centers. However geodesic segments close to the center are approximately
 * straight.
 * <p>
 * Consider a geodesic segment of length <i>l</i>. Let <i>T</i> be the point on
 * the geodesic (extended if necessary) closest to <i>C</i>, the center of the
 * projection, and <i>t</i>, be the distance <i>CT</i>. To lowest order, the
 * maximum deviation (as a true distance) of the corresponding gnomonic line
 * segment (i.e., with the same end points) from the geodesic is<br>
 * <br>
 * (<i>K</i>(<i>T</i>) &minus; <i>K</i>(<i>C</i>))
 * <i>l</i><sup>2</sup> <i>t</i> / 32.
 * <br>
 * <br>
 * where <i>K</i> is the Gaussian curvature.
 * <p>
 * This result applies for any surface. For an ellipsoid of revolution,
 * consider all geodesics whose end points are within a distance <i>r</i> of
 * <i>C</i>.  For a given <i>r</i>, the deviation is maximum when the latitude
 * of <i>C</i> is 45&deg;, when endpoints are a distance <i>r</i> away, and
 * when their azimuths from the center are &plusmn; 45&deg; or &plusmn;
 * 135&deg;. To lowest order in <i>r</i> and the flattening <i>f</i>, the
 * deviation is <i>f</i> (<i>r</i>/2<i>a</i>)<sup>3</sup> <i>r</i>.
 * <p>
 * <b>CAUTION:</b> The definition of this projection for a sphere is standard.
 * However, there is no standard for how it should be extended to an ellipsoid.
 * The choices are:
 * <ul>
 * <li>
 * Declare that the projection is undefined for an ellipsoid.
 * </li>
 * <li>
 * Project to a tangent plane from the center of the ellipsoid. This causes
 * great ellipses to appear as straight lines in the projection; i.e., it
 * generalizes the spherical great circle to a great ellipse. This was proposed
 * by independently by Bowring and Williams in 1997.
 * </li>
 * <li>
 * Project to the conformal sphere with the constant of integration chosen so
 * that the values of the latitude match for the center point and perform a
 * central projection onto the plane tangent to the conformal sphere at the
 * center point. This causes normal sections through the center point to appear
 * as straight lines in the projection; i.e., it generalizes the spherical
 * great circle to a normal section. This was proposed by I. G. Letoval'tsev,
 * Generalization of the gnomonic projection for a spheroid and the principal
 * geodetic problems involved in the alignment of surface routes, Geodesy and
 * Aerophotography (5), 271&ndash;274 (1963).
 * </li>
 * <li>
 * The projection given here. This causes geodesics close to the center point
 * to appear as straight lines in the projection; i.e., it generalizes the
 * spherical great circle to a geodesic.
 * </li>
 * </ul>
 * <p>
 * Example of use:
 *
 * <pre>
 * // Example of using the Gnomonic.java class
 * import net.sf.geographiclib.Geodesic;
 * import net.sf.geographiclib.Gnomonic;
 * import net.sf.geographiclib.GnomonicData;
 * public class ExampleGnomonic {
 *   public static void main(String[] args) {
 *     Geodesic geod = Geodesic.WGS84;
 *     double lat0 = 48 + 50 / 60.0, lon0 = 2 + 20 / 60.0; // Paris
 *     Gnomonic gnom = new Gnomonic(geod);
 *     {
 *       // Sample forward calculation
 *       double lat = 50.9, lon = 1.8; // Calais
 *       GnomonicData proj = gnom.Forward(lat0, lon0, lat, lon);
 *       System.out.println(proj.x + &quot; &quot; + proj.y);
 *     }
 *     {
 *       // Sample reverse calculation
 *       double x = -38e3, y = 230e3;
 *       GnomonicData proj = gnom.Reverse(lat0, lon0, x, y);
 *       System.out.println(proj.lat + &quot; &quot; + proj.lon);
 *     }
 *   }
 * }
 * </pre>
 */

public class Gnomonic {
  private static final double eps_ = 0.01 * Math.sqrt(Math.ulp(1.0));
  private static final int numit_ = 10;
  private Geodesic _earth;
  private double _a, _f;

  /**
   * Constructor for Gnomonic.
   * <p>
   * @param earth the {@link Geodesic} object to use for geodesic
   *   calculations.
   */
  public Gnomonic(Geodesic earth) {
    _earth = earth;
    _a = _earth.EquatorialRadius();
    _f = _earth.Flattening();
  }

  /**
   * Forward projection, from geographic to gnomonic.
   * <p>
   * @param lat0 latitude of center point of projection (degrees).
   * @param lon0 longitude of center point of projection (degrees).
   * @param lat latitude of point (degrees).
   * @param lon longitude of point (degrees).
   * @return {@link GnomonicData} object with the following fields:
   *   <i>lat0</i>, <i>lon0</i>, <i>lat</i>, <i>lon</i>, <i>x</i>, <i>y</i>,
   *   <i>azi</i>, <i>rk</i>.
   * <p>
   * <i>lat0</i> and <i>lat</i> should be in the range [&minus;90&deg;,
   * 90&deg;] and <i>lon0</i> and <i>lon</i> should be in the range
   * [&minus;540&deg;, 540&deg;). The scale of the projection is
   * 1/<i>rk<sup>2</sup></i> in the "radial" direction, <i>azi</i> clockwise
   * from true north, and is 1/<i>rk</i> in the direction perpendicular to
   * this. If the point lies "over the horizon", i.e., if <i>rk</i> &le; 0,
   * then NaNs are returned for <i>x</i> and <i>y</i> (the correct values are
   * returned for <i>azi</i> and <i>rk</i>). A call to Forward followed by a
   * call to Reverse will return the original (<i>lat</i>, <i>lon</i>) (to
   * within roundoff) provided the point in not over the horizon.
   */
  public GnomonicData Forward(double lat0, double lon0, double lat, double lon)
  {
    GeodesicData inv =
      _earth.Inverse(lat0, lon0, lat, lon,
                     GeodesicMask.AZIMUTH | GeodesicMask.GEODESICSCALE |
                     GeodesicMask.REDUCEDLENGTH);
    GnomonicData fwd =
      new GnomonicData(lat0, lon0, lat, lon, Double.NaN, Double.NaN,
                       inv.azi2, inv.M12);

    if (inv.M12 > 0) {
      double rho = inv.m12 / inv.M12;
      Pair p = new Pair();
      GeoMath.sincosd(p, inv.azi1);
      fwd.x = rho * p.first;
      fwd.y = rho * p.second;
    }

    return fwd;
  }

  /**
   * Reverse projection, from gnomonic to geographic.
   * <p>
   * @param lat0 latitude of center point of projection (degrees).
   * @param lon0 longitude of center point of projection (degrees).
   * @param x easting of point (meters).
   * @param y northing of point (meters).
   * @return {@link GnomonicData} object with the following fields:
   *   <i>lat0</i>, <i>lon0</i>, <i>lat</i>, <i>lon</i>, <i>x</i>, <i>y</i>,
   *   <i>azi</i>, <i>rk</i>.
   * <p>
   * <i>lat0</i> should be in the range [&minus;90&deg;, 90&deg;] and
   * <i>lon0</i> should be in the range [&minus;540&deg;, 540&deg;).
   * <i>lat</i> will be in the range [&minus;90&deg;, 90&deg;] and <i>lon</i>
   * will be in the range [&minus;180&deg;, 180&deg;]. The scale of the
   * projection is 1/<i>rk<sup>2</sup></i> in the "radial" direction,
   * <i>azi</i> clockwise from true north, and is 1/<i>rk</i> in the direction
   * perpendicular to this. Even though all inputs should return a valid
   * <i>lat</i> and <i>lon</i>, it's possible that the procedure fails to
   * converge for very large <i>x</i> or <i>y</i>; in this case NaNs are
   * returned for all the output arguments. A call to Reverse followed by a
   * call to Forward will return the original (<i>x</i>, <i>y</i>) (to
   * roundoff).
   */
  public GnomonicData Reverse(double lat0, double lon0, double x, double y) {
    GnomonicData rev =
      new GnomonicData(lat0, lon0, Double.NaN, Double.NaN, x, y, Double.NaN,
                       Double.NaN);

    double azi0 = GeoMath.atan2d(x, y);
    double rho = Math.hypot(x, y);
    double s = _a * Math.atan(rho / _a);
    boolean little = rho <= _a;

    if (!little)
      rho = 1 / rho;

    GeodesicLine line =
      _earth.Line(lat0, lon0, azi0, GeodesicMask.LATITUDE
                  | GeodesicMask.LONGITUDE | GeodesicMask.AZIMUTH
                  | GeodesicMask.DISTANCE_IN | GeodesicMask.REDUCEDLENGTH
                  | GeodesicMask.GEODESICSCALE);

    int count = numit_, trip = 0;
    GeodesicData pos = null;

    while (count-- > 0) {
      pos =
        line.Position(s, GeodesicMask.LONGITUDE | GeodesicMask.LATITUDE
                      | GeodesicMask.AZIMUTH | GeodesicMask.DISTANCE_IN
                      | GeodesicMask.REDUCEDLENGTH
                      | GeodesicMask.GEODESICSCALE);

      if (trip > 0)
        break;

      double ds =
        little ? ((pos.m12 / pos.M12) - rho) * pos.M12 * pos.M12
        : (rho - (pos.M12 / pos.m12)) * pos.m12 * pos.m12;
      s -= ds;

      if (Math.abs(ds) <= eps_ * _a)
        trip++;
    }

    if (trip == 0)
      return rev;

    rev.lat = pos.lat2;
    rev.lon = pos.lon2;
    rev.azi = pos.azi2;
    rev.rk = pos.M12;

    return rev;
  }

  /**
   * @return <i>a</i> the equatorial radius of the ellipsoid (meters).  This is
   *   the value inherited from the Geodesic object used in the constructor.
   **********************************************************************/
  public double EquatorialRadius() { return _a; }

  /**
   * @return <i>f</i> the  flattening of the ellipsoid.  This is
   *   the value inherited from the Geodesic object used in the constructor.
   **********************************************************************/
  public double Flattening() { return _f; }

  /**
   * @deprecated An old name for {@link #EquatorialRadius()}.
   * @return <i>a</i> the equatorial radius of the ellipsoid (meters).
   **********************************************************************/
  @Deprecated
  public double MajorRadius() { return EquatorialRadius(); }
}