function [x, y, azi, rk] = gnomonic_fwd(lat0, lon0, lat, lon, ellipsoid) %GNOMONIC_FWD Forward ellipsoidal gnomonic projection % % [x, y] = GNOMONIC_FWD(lat0, lon0, lat, lon) % [x, y, azi, rk] = GNOMONIC_FWD(lat0, lon0, lat, lon, ellipsoid) % % performs the forward ellipsoidal gnomonic projection of points % (lat,lon) to (x,y) using (lat0,lon0) as the center of projection. % These input arguments can be scalars or arrays of equal size. The % ellipsoid vector is of the form [a, e], where a is the equatorial % radius in meters, e is the eccentricity. If ellipsoid is omitted, the % WGS84 ellipsoid (more precisely, the value returned by % defaultellipsoid) is used. projdoc defines the projection and gives % the restrictions on the allowed ranges of the arguments. The inverse % projection is given by gnomonic_inv. % % azi and rk give metric properties of the projection at (lat,lon); azi % is the azimuth of the geodesic from the center of projection and rk is % the reciprocal of the azimuthal scale. The scale in the radial % direction is 1/rk^2. % % If the point lies "over the horizon", i.e., if rk <= 0, then NaNs are % returned for x and y (the correct values are returned for azi and rk). % % lat0, lon0, lat, lon, azi are in degrees. The projected coordinates x, % y are in meters (more precisely the units used for the equatorial % radius). rk is dimensionless. % % The ellipsoidal gnomonic projection is an azimuthal projection about a % center point. All geodesics through the center point are projected % into straight lines with the correct azimuth relative to the center % point. In addition all geodesics are pass close to the center point % are very nearly straight. The projection is derived in Section 8 of % % C. F. F. Karney, Algorithms for geodesics, % J. Geodesy 87, 43-55 (2013); % https://doi.org/10.1007/s00190-012-0578-z % Addenda: https://geographiclib.sourceforge.io/geod-addenda.html % % which also includes methods for solving the "intersection" and % "interception" problems using the gnomonic projection. % % See also PROJDOC, GNOMONIC_INV, GEODDISTANCE, DEFAULTELLIPSOID, % FLAT2ECC. % Copyright (c) Charles Karney (2012-2015) <charles@karney.com>. narginchk(4, 5) if nargin < 5, ellipsoid = defaultellipsoid; end try [~] = lat0 + lon0 + lat + lon; catch error('lat0, lon0, lat, lon have incompatible sizes') end [~, azi0, azi, ~, m, M] = geoddistance(lat0, lon0, lat, lon, ellipsoid); rho = m ./ M; [x, y] = sincosdx(azi0); x = rho .* x; y = rho .* y; rk = M; x(M <= 0) = nan; y(M <= 0) = nan; end