Gnomonic.java

Go to the documentation of this file.

package net.sf.geographiclib;

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

  public Gnomonic(Geodesic earth) {
    _earth = earth;
    _a = _earth.MajorRadius();
    _f = _earth.Flattening();
  }

  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 = GeoMath.sincosd(inv.azi1);
      fwd.x = rho * p.first;
      fwd.y = rho * p.second;
    }

    return fwd;
  }

  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;
  }

  public double MajorRadius() { return _a; }

  public double Flattening() { return _f; }
}