1function [lat, lon, azi, rk] = gnomonic_inv(lat0, lon0, x, y, ellipsoid)
2%GNOMONIC_INV  Inverse ellipsoidal gnomonic projection
3%
4%   [lat, lon] = GNOMONIC_INV(lat0, lon0, x, y)
5%   [lat, lon, azi, rk] = GNOMONIC_INV(lat0, lon0, x, y, ellipsoid)
6%
7%   performs the inverse ellipsoidal gnomonic projection of points (x,y) to
8%   (lat,lon) using (lat0,lon0) as the center of projection.  These input
9%   arguments can be scalars or arrays of equal size.  The ellipsoid vector
10%   is of the form [a, e], where a is the equatorial radius in meters, e is
11%   the eccentricity.  If ellipsoid is omitted, the WGS84 ellipsoid (more
12%   precisely, the value returned by defaultellipsoid) is used.  projdoc
13%   defines the projection and gives the restrictions on the allowed ranges
14%   of the arguments.  The forward projection is given by gnomonic_fwd.
15%
16%   azi and rk give metric properties of the projection at (lat,lon); azi
17%   is the azimuth of the geodesic from the center of projection and rk is
18%   the reciprocal of the azimuthal scale.  The scale in the radial
19%   direction is 1/rk^2.
20%
21%   In principle, all finite x and y are allowed.  However, it's possible
22%   that the inverse projection fails for very large x and y (when the
23%   geographic position is close to the "horizon").  In that case, NaNs are
24%   returned for the corresponding output variables.
25%
26%   lat0, lon0, lat, lon, azi are in degrees.  The projected coordinates x,
27%   y are in meters (more precisely the units used for the equatorial
28%   radius).  rk is dimensionless.
29%
30%   The ellipsoidal gnomonic projection is an azimuthal projection about a
31%   center point.  All geodesics through the center point are projected
32%   into straight lines with the correct azimuth relative to the center
33%   point.  In addition all geodesics are pass close to the center point
34%   are very nearly straight.  The projection is derived in Section 8 of
35%
36%     C. F. F. Karney, Algorithms for geodesics,
37%     J. Geodesy 87, 43-55 (2013);
38%     https://doi.org/10.1007/s00190-012-0578-z
39%     Addenda: https://geographiclib.sourceforge.io/geod-addenda.html
40%
41%   which also includes methods for solving the "intersection" and
42%   "interception" problems using the gnomonic projection.
43%
44%   See also PROJDOC, GNOMONIC_FWD, GEODRECKON, DEFAULTELLIPSOID, FLAT2ECC.
45
46% Copyright (c) Charles Karney (2012-2015) <charles@karney.com>.
47
48  narginchk(4, 5)
49  if nargin < 5, ellipsoid = defaultellipsoid; end
50  try
51    Z = zeros(size(lat0 + lon0 + x + y));
52  catch
53    error('lat0, lon0, x, y have incompatible sizes')
54  end
55  if length(ellipsoid(:)) ~= 2
56    error('ellipsoid must be a vector of size 2')
57  end
58  a = ellipsoid(1);
59  numit = 10;
60  eps1 = a * 0.01 * sqrt(eps);
61
62  lat0 = lat0 + Z; lon0 = lon0 + Z; x = x + Z; y = y + Z;
63  azi0 = atan2dx(x, y);
64  rho = hypot(x, y);
65  s = a * atan(rho / a);
66  little = rho <= a;
67  rho(~little) = 1 ./ rho(~little);
68  g = Z == 0; trip = ~g;
69  lat = Z; lon = Z; azi = Z; m = Z; M = Z; ds = Z;
70  for k = 1 : numit
71    [lat(g), lon(g), azi(g), ~, m(g), M(g)] = ...
72        geodreckon(lat0(g), lon0(g), s(g), azi0(g), ellipsoid);
73    g = g & ~trip;
74    if ~any(g), break, end
75    c = little & g;
76    ds(c) = (m(c) ./ M(c) - rho(c)) .* M(c).^2;
77    c = ~little & g;
78    ds(c) = (rho(c) - M(c) ./ m(c)) .* m(c).^2;
79    s(g) = s(g) - ds(g);
80    trip(g) = ~(abs(ds(g)) >= eps1);
81  end
82  c = ~trip;
83  if any(c)
84    lat(c) = nan;
85    lon(c) = nan;
86    azi(c) = nan;
87    M(c) = nan;
88  end
89  rk = M;
90end
91