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