1/*
2
3Arbitrary precision Transverse Mercator Projection
4
5Copyright (c) Charles Karney (2009-2017) <charles@karney.com> and
6licensed under the MIT/X11 License.  For more information, see
7https://geographiclib.sourceforge.io/
8
9Reference:
10
11   Charles F. F. Karney,
12   Transverse Mercator with an accuracy of a few nanometers,
13   J. Geodesy 85(8), 475-485 (Aug. 2011).
14   DOI 10.1007/s00190-011-0445-3
15   preprint https://arxiv.org/abs/1002.1417
16   resource page https://geographiclib.sourceforge.io/tm.html
17
18The parameters for the transformation are set by
19
20setparams(a,f,k0)$
21    sets the equatorial radius, inverse flattening, and central scale
22    factor. The default is
23      setparams(6378137b0, 1/298.257223563b0, 0.9996b0)$
24    appropriate for UTM applications.
25
26tm(lat,lon);
27    takes lat and lon args (in degrees) and returns
28      [x, y, convergence, scale]
29    [x, y] do not include false eastings/northings but do include the
30    scale factor k0.  convergence is in degrees.
31
32ll(x,y);
33    takes x and y args (in meters) and returns
34      [lat, lon, convergence, scale].
35
36Example:
37
38$ maxima
39Maxima 5.15.0 http://maxima.sourceforge.net
40Using Lisp CLISP 2.43 (2007-11-18)
41Distributed under the GNU Public License. See the file COPYING.
42Dedicated to the memory of William Schelter.
43The function bug_report() provides bug reporting information.
44(%i1) load("tm.mac")$
45(%i2) tm(10b0,20b0);
46(%o2) [2.235209504622466691587930831718465965864199221939781808953597771095103\
476690000464b6, 1.17529734503138466792126931904154130080533935727351398258511134\
4868541970512119385b6, 3.6194756227592979778565787394402350354250845160819430786\
49093514889500602612857052b0, 1.062074627142564335518604915718789933200854739344\
508664109599248189291146283796933b0]
51(%i3) ll(%[1],%[2]);
52(%o3) [1.0b1, 2.0b1, 3.6194756227592979778565787394402350354250845160819430786\
53093514889500602612857053b0, 1.062074627142564335518604915718789933200854739344\
548664109599248189291146283796933b0]
55(%i4) float(%o2);
56(%o4) [2235209.504622467, 1175297.345031385, 3.619475622759298,
57                                                             1.062074627142564]
58(%i5) float(%o3);
59(%o5)         [10.0, 20.0, 3.619475622759298, 1.062074627142564]
60
61This implements GeographicLib::TransverseMercatorExact (i.e., Lee, 1976)
62using bfloats.  However fewer changes from Lee 1976 have been made since
63we rely more heavily on the high precision to deal with problem cases.
64
65To change the precision, change fpprec below and reload.
66
67*/
68
69fpprec:80$
70load("ellint.mac")$ /* Load elliptic functions */
71
72tol:0.1b0^fpprec$
73tol1:0.1b0*sqrt(tol)$ /* For Newton's method */
74tol2:sqrt(0.01*tol*tol1)$ /* Also for Newton's method but more conservative */
75ahypover:log(10b0^fpprec)+2$
76
77pi:bfloat(%pi)$
78degree:pi/180$
79ratprint:false$
80debugprint:false$
81setparams(a1,f1,k1):=(a:bfloat(a1),f:bfloat(f1),k0:bfloat(k1),
82  e2:f*(2-f),
83  e:sqrt(e2),
84  kcu:kc(e2),
85  kcv:kc(1-e2),
86  ecu:ec(e2),
87  ecv:ec(1-e2),
88  n:f/(2-f),
89  'done)$
90setparams(6378137b0, 1/298.257223563b0, 0.9996b0)$  /* WGS 84 */
91/* setparams(6378388b0, 1/297b0, 0.9996b0)$ International */
92/* setparams(1/ec(0.01b0), 1/(30*sqrt(11b0)+100), 1b0)$ testing, eps = 0.1*/
93
94/*
95Interpret x_y(y) as x <- y, i.e., "transform quantity y to quantity x"
96
97Let
98
99phi = geodetic latitude
100psi = isometric latitude ( + i * lambda )
101sigma = TM coords
102thom = Thompson coords
103
104*/
105
106/* sqrt(x^2 + y^2) -- Real only */
107hypot(x,y):=sqrt(x^2 + y^2)$
108
109/* log(1 + x) -- Real only */
110log1p(x) := block([y : 1b0+x],
111  if y = 1b0 then x else x*log(y)/(y - 1))$
112
113/* Real only */
114/* Some versions of Maxima have a buggy atanh
115atnh(x) := block([y : abs(x)],
116  y : log1p(2 * y/(1 - y))/2,
117  if x < 0 then -y else y)$ */
118atnh(x) := atanh(x)$
119
120/* exp(x)-1 -- Real only */
121expm1(x) := block([y : exp(bfloat(x)),z],
122  z : y - 1b0,
123  if abs(x) > 1b0 then z else if z = 0b0 then x else x * z/log(y))$
124
125/* Real only */
126/* Some versions of Maxima have a buggy sinh */
127snh(x) := block([u : expm1(x)],
128  (u / (u + 1)) * (u + 2) /2);
129
130/* Real only */
131psi_phi(phi):=block([s:sin(phi)],
132  asinh(s/max(cos(phi),0.1b0*tol)) - e * atnh(e * s))$
133
134/* Real only */
135phi_psi(psi):=block([q:psi,t,dq],
136  for i do (
137    t:tanh(q),
138    dq : -(q - e * atnh(e * t) - psi) * (1 - e2 * t^2) / (1 - e2),
139    q : q + dq,
140    if debugprint then print(float(q), float(dq)),
141    if abs(dq) < tol1 then return(false)),
142  atan(snh(q)))$
143
144psi_thom_comb(w):=block([jacr:sncndn(bfloat(realpart(w)),1-e2),
145  jaci:sncndn(bfloat(imagpart(w)),e2),d,d1,d2],
146  d:(1-e2)*(jaci[2]^2 + e2 * (jacr[1] * jaci[1])^2)^2,
147  d1:sqrt(jacr[2]^2 + (1-e2) * (jacr[1] * jaci[1])^2),
148  d2:sqrt(e2 * jacr[2]^2 + (1-e2) * jaci[2]^2),
149[
150  (if d1 > 0b0 then asinh(jacr[1]*jaci[3]/ d1) else signnum(snu) * ahypover)
151  - (if d2 > 0b0 then e * asinh(e * jacr[1] / d2) else signnum(snu) * ahypover)
152  + %i * (if d1 > 0b0 and d2 > 0b0 then
153    atan2(jacr[3]*jaci[1],jacr[2]*jaci[2])
154    - e * atan2(e*jacr[2]*jaci[1],jacr[3]*jaci[2]) else 0),
155  jacr[2]*jacr[3]*jaci[3]*(jaci[2]^2-e2*(jacr[1]*jaci[1])^2)/d
156  -%i * jacr[1]*jaci[1]*jaci[2]*((jacr[3]*jaci[3])^2+e2*jacr[2]^2)/d]
157)$
158
159psi_thom(w):=block([tt:psi_thom_comb(w)],tt[1])$
160inv_diff_psi_thom(w):=block([tt:psi_thom_comb(w)],tt[2])$
161
162w0a(psi):=block([lam:bfloat(imagpart(psi)),psia:bfloat(realpart(psi))],
163  rectform(kcu/(pi/2)*( atan2(snh(psia),cos(lam))
164      +%i*asinh(sin(lam)/sqrt(cos(lam)^2 + snh(psia)^2)))))$
165
166w0c(psi):=block([m,a,dlam],
167  dlam:bfloat(imagpart(psi))-pi/2*(1-e),
168  psi:bfloat(realpart(psi)),
169  m:sqrt(psi^2+dlam^2)*3/(1-e2)/e,
170  a:if m = 0b0 then 0 else atan2(dlam-psi, psi+dlam) - 0.75b0*pi,
171  m:m^(1/3), a:a/3,
172  m*cos(a)+%i*(m*sin(a)+kcv))$
173
174w0d(psi):=block([psir:-realpart(psi)/e+1b0,lam:(pi/2-imagpart(psi))/e,uu,vv],
175  uu:asinh(sin(lam)/sqrt(cos(lam)^2+snh(psir)^2))*(1+e2/2),
176  vv:atan2(cos(lam), snh(psir)) *(1+e2/2),
177  (-uu+kcu) + %i * (-vv+kcv))$
178
179w0m(psi):=if realpart(psi)<-e/2*pi/2 and
180imagpart(psi)>pi/2*(1-2*e) and
181realpart(psi) < imagpart(psi)-(pi/2*(1-e)) then w0d(psi) else
182if realpart(psi)<e*pi/2 and imagpart(psi)>pi/2*(1-2*e)
183then w0c(psi) else w0a(psi)$
184w0(psi):=w0m(psi)$
185
186thom_psi(psi):=block([w:w0(psi),dw,v,vv],
187if not(abs(psi-pi/2*(1-e)*%i) < e * tol^0.6b0) then
188  for i do (
189    if i > 100 then error("too many iterations"),
190    vv:psi_thom_comb(w),
191    v:vv[1],
192    dw:-rectform((v-psi)*vv[2]),
193    w:w+dw,
194    dw:abs(dw),
195    if debugprint then print(float(w),float(dw)),
196    /* error is approx dw^2/2 */
197    if dw < tol2 then return(false)
198    ),
199  w
200  )$
201
202sigma_thom_comb(w):=block([u:bfloat(realpart(w)),v:bfloat(imagpart(w)),
203  jacr,jaci,phi,iu,iv,den,den1,er,ei,dnr,dni],
204  jacr:sncndn(u,1-e2),jaci:sncndn(v,e2),
205  er:eirx(jacr[1],jacr[2],jacr[3],e2,ecu),
206  ei:eirx(jaci[1],jaci[2],jaci[3],1-e2,ecv),
207  den:e2*jacr[2]^2+(1-e2)*jaci[2]^2,
208  den1:(1-e2)*(jaci[2]^2 + e2 * (jacr[1] * jaci[1])^2)^2,
209  dnr:jacr[3]*jaci[2]*jaci[3],
210  dni:-e2*jacr[1]*jacr[2]*jaci[1],
211[  er - e2*jacr[1]*jacr[2]*jacr[3]/den
212  + %i*(v - ei + (1-e2)*jaci[1]*jaci[2]*jaci[3]/den),
213  (dnr^2-dni^2)/den1 + %i * 2*dnr*dni/den1])$
214
215sigma_thom(w):=block([tt:sigma_thom_comb(w)],tt[1])$
216inv_diff_sigma_thom(w):=block([tt:sigma_thom_comb(w)],tt[2])$
217
218wx0a(sigma):=rectform(sigma*kcu/ecu)$
219wx0b(sigma):=block([m,aa],
220  sigma:rectform(sigma-%i*(kcv-ecv)),
221  m:abs(sigma)*3/(1-e2),
222  aa:atan2(imagpart(sigma),realpart(sigma)),
223  if aa<-pi/2 then aa:aa+2*pi,
224  aa:aa-pi,
225  rectform(m^(1/3)*(cos(aa/3b0)+%i*sin(aa/3b0))+%i*kcv))$
226wx0c(sigma):=rectform(1/(sigma-(ecu+%i*(kcv-ecv))) + kcu+%i*kcv)$
227wx0m(sigma):=block([eta:bfloat(imagpart(sigma)),
228  xi:bfloat(realpart(sigma))],
229  if eta > 1.25b0 * (kcv-ecv) or (xi < -0.25*ecu and xi < eta-(kcv-ecv)) then
230  wx0c(sigma) else
231  if (eta > 0.75b0 * (kcv-ecv) and xi < 0.25b0 * ecu) or
232  eta > kcv-ecv or xi < 0 then wx0b(sigma) else wx0a(sigma))$
233wx0(sigma):=wx0m(sigma)$
234thom_sigma(sigma):=block([w:wx0(sigma),dw,v,vv],
235  for i do (
236    if i > 100 then error("too many iterations"),
237    vv:sigma_thom_comb(w),
238    v:vv[1],
239    dw:-rectform((v-sigma)*vv[2]),
240    w:w+dw,
241    dw:abs(dw),
242    if debugprint then print(float(w),float(dw)),
243    /* error is approx dw^2/2 */
244    if dw < tol2 then return(false)
245    ),
246  w
247  )$
248
249/* Lee/Thompson's method forward */
250
251tm(phi,lam):=block([psi,thom,jacr,jaci,sigma,gam,scale,c],
252  phi:phi*degree,
253  lam:lam*degree,
254  psi:psi_phi(phi),
255  thom:thom_psi(psi+%i*lam),
256  jacr:sncndn(bfloat(realpart(thom)),1-e2),
257  jaci:sncndn(bfloat(imagpart(thom)),e2),
258  sigma:sigma_thom(thom),
259  c:cos(phi),
260  if c > tol1 then (
261    gam:atan2((1-e2)*jacr[1]*jaci[1]*jaci[2],
262      jacr[2]*jacr[3]*jaci[3]),
263    scale:sqrt(1-e2 + e2 * c^2)/c*
264    sqrt(((1-e2)*jaci[1]^2 + (jacr[2]*jaci[3])^2)/
265      (e2*jacr[2]^2 + (1-e2)*jaci[2]^2)))
266  else (gam : lam, scale : 1b0),
267  [imagpart(sigma)*k0*a,realpart(sigma)*k0*a,gam/degree,k0*scale])$
268
269/* Lee/Thompson's method reverse */
270
271ll(x,y):=block([sigma,thom,jacr,jaci,psi,lam,phi,gam,scale,c],
272  sigma:y/(a*k0)+%i*x/(a*k0),
273  thom:thom_sigma(sigma),
274  jacr:sncndn(bfloat(realpart(thom)),1-e2),
275  jaci:sncndn(bfloat(imagpart(thom)),e2),
276  psi:psi_thom(thom),
277  lam:bfloat(imagpart(psi)),
278  psi:bfloat(realpart(psi)),
279  phi:phi_psi(psi),
280  c:cos(phi),
281  if c > tol1 then (
282    gam:atan2((1-e2)*jacr[1]*jaci[1]*jaci[2],
283      jacr[2]*jacr[3]*jaci[3]),
284    scale:sqrt(1-e2 + e2 * c^2)/c*
285    sqrt(((1-e2)*jaci[1]^2 + (jacr[2]*jaci[3])^2)/
286      (e2*jacr[2]^2 + (1-e2)*jaci[2]^2)))
287  else (gam : lam, scale : 1b0),
288  [phi/degree,lam/degree,gam/degree,k0*scale])$
289
290/* Return lat/lon/x/y for a point specified in Thompson coords */
291/* Pick u in [0, kcu] and v in [0, kcv] */
292lltm(u,v):=block([jacr,jaci,psi,lam,phi,c,gam,scale,sigma,x,y],
293  u:bfloat(u), v:bfloat(v),
294  jacr:sncndn(u,1-e2),
295  jaci:sncndn(v,e2),
296  psi:psi_thom(u+%i*v),
297  sigma:sigma_thom(u+%i*v),
298  x:imagpart(sigma)*k0*a,y:realpart(sigma)*k0*a,
299  lam:bfloat(imagpart(psi)),
300  psi:bfloat(realpart(psi)),
301  phi:phi_psi(psi),
302  c:cos(phi),
303  if c > tol1 then (
304    gam:atan2((1-e2)*jacr[1]*jaci[1]*jaci[2],
305      jacr[2]*jacr[3]*jaci[3]),
306    scale:sqrt(1-e2 + e2 * c^2)/c*
307    sqrt(((1-e2)*jaci[1]^2 + (jacr[2]*jaci[3])^2)/
308      (e2*jacr[2]^2 + (1-e2)*jaci[2]^2)))
309  else (gam : lam, scale : 1b0),
310  [phi/degree,lam/degree,x,y,gam/degree,k0*scale])$
311
312/* Gauss-Krueger series to order n^i forward
313
314Uses the array functions
315
316 a1_a[i](n), zeta_a[i](z,n), zeta_d[i](z,n), zetap_a[i](s,n), zetap_d[i](s,n),
317
318defined in tmseries.mac.
319*/
320
321tms(phi,lam,i):=block([psi,xip,etap,z,sigma,sp,gam,k,b1],
322  phi:phi*degree,
323  lam:lam*degree,
324  psi:psi_phi(phi),
325  xip:atan2(snh(psi), cos(lam)),
326  etap:asinh(sin(lam)/hypot(snh(psi),cos(lam))),
327  k:sqrt(1 - e2*sin(phi)^2)/(cos(phi)*hypot(snh(psi),cos(lam))),
328  gam:atan(tan(xip)*tanh(etap)),
329  z:xip+%i*etap,
330  b1:a1_a[i](n),
331  sigma:rectform(b1*zeta_a[i](z,n)),
332  sp:rectform(zeta_d[i](z,n)),
333  gam : gam - atan2(imagpart(sp),realpart(sp)),
334  k : k * b1 * cabs(sp),
335  [imagpart(sigma)*k0*a,realpart(sigma)*k0*a,gam/degree,k*k0])$
336
337/* Gauss-Krueger series to order n^i reverse */
338
339lls(x,y,i):=block([sigma,b1,s,z,zp,xip,etap,s,c,r,gam,k,lam,psi,phi],
340  sigma:y/(a*k0)+%i*x/(a*k0),
341  b1:a1_a[i](n),
342  s:rectform(sigma/b1),
343  z:rectform(zetap_a[i](s,n)),
344  zp:rectform(zetap_d[i](s,n)),
345  gam : atan2(imagpart(zp), realpart(zp)),
346  k : b1 / cabs(zp),
347  xip:realpart(z),
348  etap:imagpart(z),
349  s:snh(etap),
350  c:cos(xip),
351  r:hypot(s, c),
352  lam:atan2(s, c),
353  psi : asinh(sin(xip)/r),
354  phi :phi_psi(psi),
355  k : k * sqrt(1 - e2*sin(phi)^2) * r/cos(phi),
356  gam : gam + atan(tan(xip) * tanh(etap)),
357  [phi/degree,lam/degree,gam/degree,k*k0])$
358
359/* Approx geodesic distance valid for small displacements */
360
361dist(phi0,lam0,phi,lam):=block([dphi,dlam,nn,hlon,hlat],
362  dphi:(phi-phi0)*degree,
363  dlam:(lam-lam0)*degree,
364  phi0:phi0*degree,
365  lam0:lam0*degree,
366  nn : 1/sqrt(1 - e2 * sin(phi0)^2),
367  hlon : cos(phi0) * nn,
368  hlat : (1 - e2) * nn^3,
369  a * hypot(dphi*hlat, dlam*hlon))$
370
371/* Compute truncation errors for all truncation levels */
372
373check(phi,lam):=block([vv,x,y,gam,k,vf,vb,errf,errr,err2,errlist],
374  phi:min(90-0.01b0,phi),
375  lam:min(90-0.01b0,lam),
376  vv:tm(phi,lam),
377  errlist:[],
378  x:vv[1], y:vv[2], gam:vv[3], k:vv[4],
379  for i:1 thru maxpow do (
380    vf:tms(phi,lam,i),
381    errf:hypot(vf[1]-x,vf[2]-y)/k,
382    errfg:abs(vf[3]-gam),
383    errfk:abs((vf[4]-k)/k),
384    vb:lls(x,y,i),
385    errr:dist(phi,lam,vb[1],vb[2]),
386    errrg:abs(vb[3]-gam),
387    errrk:abs((vb[4]-k)/k),
388    errlist:append(errlist,
389      [max(errf, errr), max(errfg, errrg), max(errfk, errrk)])),
390  errlist)$
391
392/* Max of output of check over a set of points */
393
394checka(lst):=block([errlist:[],errx],
395  for i:1 thru 3*maxpow do errlist:cons(0b0,errlist),
396  for vv in lst do (
397    errx:check(vv[1],vv[2]),
398    for i:1 thru 3*maxpow do errlist[i]:max(errlist[i],errx[i])),
399  errlist)$
400