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