1function [x, y, gam, k] = tranmerc_fwd(lat0, lon0, lat, lon, ellipsoid) 2%TRANMERC_FWD Forward transverse Mercator projection 3% 4% [x, y] = TRANMERC_FWD(lat0, lon0, lat, lon) 5% [x, y, gam, k] = TRANMERC_FWD(lat0, lon0, lat, lon, ellipsoid) 6% 7% performs the forward transverse Mercator projection of points (lat,lon) 8% to (x,y) 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. The common 13% case of lat0 = 0 is treated efficiently provided that lat0 is specified 14% as a scalar. projdoc defines the projection and gives the restrictions 15% on the allowed ranges of the arguments. The inverse projection is 16% given by tranmerc_inv. The scale on the central meridian is 1. 17% 18% gam and k give metric properties of the projection at (lat,lon); gam is 19% the meridian convergence at the point and k is the scale. 20% 21% lat0, lon0, lat, lon, gam are in degrees. The projected coordinates x, 22% y are in meters (more precisely the units used for the equatorial 23% radius). k is dimensionless. 24% 25% This implementation of the projection is based on the series method 26% described in 27% 28% C. F. F. Karney, Transverse Mercator with an accuracy of a few 29% nanometers, J. Geodesy 85(8), 475-485 (Aug. 2011); 30% Addenda: https://geographiclib.sourceforge.io/tm-addenda.html 31% 32% This extends the series given by Krueger (1912) to sixth order in the 33% flattening. This is a substantially better series than that used by 34% the MATLAB mapping toolbox. In particular the errors in the projection 35% are less than 5 nanometers within 3900 km of the central meridian (and 36% less than 1 mm within 7600 km of the central meridian). The mapping 37% can be continued accurately over the poles to the opposite meridian. 38% 39% See also PROJDOC, TRANMERC_INV, UTMUPS_FWD, UTMUPS_INV, 40% DEFAULTELLIPSOID, FLAT2ECC. 41 42% Copyright (c) Charles Karney (2012-2018) <charles@karney.com>. 43 44 narginchk(4, 5) 45 if nargin < 5, ellipsoid = defaultellipsoid; end 46 try 47 S = size(lat0 + lon0 + lat + lon); 48 catch 49 error('lat0, lon0, lat, lon have incompatible sizes') 50 end 51 if length(ellipsoid(:)) ~= 2 52 error('ellipsoid must be a vector of size 2') 53 end 54 55 Z = zeros(prod(S),1); 56 maxpow = 6; 57 58 a = ellipsoid(1); 59 e2 = real(ellipsoid(2)^2); 60 f = e2 / (1 + sqrt(1 - e2)); 61 e2m = 1 - e2; 62 cc = sqrt(e2m) * exp(eatanhe(1, e2)); 63 n = f / (2 -f); 64 alp = alpf(n); 65 b1 = (1 - f) * (A1m1f(n) + 1); 66 a1 = b1 * a; 67 68 lat = LatFix(lat(:)) + Z; 69 lon = AngDiff(lon0(:), lon(:)) + Z; 70 71 latsign = 1 - 2 * (lat < 0); 72 lonsign = 1 - 2 * (lon < 0); 73 lon = lon .* lonsign; 74 lat = lat .* latsign; 75 backside = lon > 90; 76 latsign(backside & lat == 0) = -1; 77 lon(backside) = 180 - lon(backside); 78 [sphi, cphi] = sincosdx(lat); 79 [slam, clam] = sincosdx(lon); 80 tau = sphi ./ max(sqrt(realmin), cphi); 81 taup = taupf(tau, e2); 82 xip = atan2(taup, clam); 83 etap = asinh(slam ./ hypot(taup, clam)); 84 gam = atan2dx(slam .* taup, clam .* hypot(1, taup)); 85 k = sqrt(e2m + e2 * cphi.^2) .* hypot(1, tau) ./ hypot(taup, clam); 86 c = ~(lat ~= 90); 87 if any(c) 88 xip(c) = pi/2; 89 etap(c) = 0; 90 gam(c) = lon(c); 91 k(c) = cc; 92 end 93 c0 = cos(2 * xip); ch0 = cosh(2 * etap); 94 s0 = sin(2 * xip); sh0 = sinh(2 * etap); 95 a = 2 * complex(c0 .* ch0, -s0 .* sh0); 96 j = maxpow; 97 y0 = complex(Z); y1 = y0; 98 z0 = y0; z1 = y0; 99 if mod(j, 2) 100 y0 = y0 + alp(j); 101 z0 = z0 + 2*j * alp(j); 102 j = j - 1; 103 end 104 for j = j : -2 : 1 105 y1 = a .* y0 - y1 + alp(j); 106 z1 = a .* z0 - z1 + 2*j * alp(j); 107 y0 = a .* y1 - y0 + alp(j-1); 108 z0 = a .* z1 - z0 + 2*(j-1) * alp(j-1); 109 end 110 a = a/2; 111 z1 = 1 - z1 + z0 .* a; 112 a = complex(s0 .* ch0, c0 .* sh0); 113 y1 = complex(xip, etap) + y0 .* a; 114 gam = gam - atan2dx(imag(z1), real(z1)); 115 k = k .* (b1 * abs(z1)); 116 xi = real(y1); eta = imag(y1); 117 xi(backside) = pi - xi(backside); 118 y = a1 * xi .* latsign; 119 x = a1 * eta .* lonsign; 120 gam(backside) = 180 - gam(backside); 121 gam = AngNormalize(gam .* latsign .* lonsign); 122 123 if isscalar(lat0) && lat0 == 0 124 y0 = 0; 125 else 126 [sbet0, cbet0] = sincosdx(LatFix(lat0(:))); 127 [sbet0, cbet0] = norm2((1-f) * sbet0, cbet0); 128 y0 = a1 * (atan2(sbet0, cbet0) + ... 129 SinCosSeries(true, sbet0, cbet0, C1f(n))); 130 end 131 y = y - y0; 132 133 x = reshape(x, S); y = reshape(y, S); 134 gam = reshape(gam, S); k = reshape(k, S); 135end 136 137function alp = alpf(n) 138 persistent alpcoeff 139 if isempty(alpcoeff) 140 alpcoeff = [ ... 141 31564, -66675, 34440, 47250, -100800, 75600, 151200, ... 142 -1983433, 863232, 748608, -1161216, 524160, 1935360, ... 143 670412, 406647, -533952, 184464, 725760, ... 144 6601661, -7732800, 2230245, 7257600, ... 145 -13675556, 3438171, 7983360, ... 146 212378941, 319334400, ... 147 ]; 148 end 149 maxpow = 6; 150 alp = zeros(1, maxpow); 151 o = 1; 152 d = n; 153 for l = 1 : maxpow 154 m = maxpow - l; 155 alp(l) = d * polyval(alpcoeff(o : o + m), n) / alpcoeff(o + m + 1); 156 o = o + m + 2; 157 d = d * n; 158 end 159end 160