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