221 lines
7.3 KiB
Plaintext
221 lines
7.3 KiB
Plaintext
/*
|
|
Compute series approximations for Transverse Mercator Projection
|
|
|
|
Copyright (c) Charles Karney (2009-2010) <charles@karney.com> and
|
|
licensed under the MIT/X11 License. For more information, see
|
|
https://geographiclib.sourceforge.io/
|
|
|
|
Reference:
|
|
|
|
Charles F. F. Karney,
|
|
Transverse Mercator with an accuracy of a few nanometers,
|
|
J. Geodesy 85(8), 475-485 (Aug. 2011).
|
|
DOI 10.1007/s00190-011-0445-3
|
|
preprint https://arxiv.org/abs/1002.1417
|
|
resource page https://geographiclib.sourceforge.io/tm.html
|
|
|
|
Compute coefficient for forward and inverse trigonometric series for
|
|
conversion from conformal latitude to rectifying latitude. This prints
|
|
out assignments which with minor editing are suitable for insertion into
|
|
C++ code. (N.B. n^3 in the output means n*n*n; 3/5 means 0.6.)
|
|
|
|
To run, start maxima and enter
|
|
|
|
writefile("tmseries.out")$
|
|
load("tmseries.mac")$
|
|
closefile()$
|
|
|
|
With maxpow = 6, the output (after about 5 secs) is
|
|
|
|
A=a/(n+1)*(1+n^2/4+n^4/64+n^6/256);
|
|
alpha[1]=n/2-2*n^2/3+5*n^3/16+41*n^4/180-127*n^5/288+7891*n^6/37800;
|
|
alpha[2]=13*n^2/48-3*n^3/5+557*n^4/1440+281*n^5/630-1983433*n^6/1935360;
|
|
alpha[3]=61*n^3/240-103*n^4/140+15061*n^5/26880+167603*n^6/181440;
|
|
alpha[4]=49561*n^4/161280-179*n^5/168+6601661*n^6/7257600;
|
|
alpha[5]=34729*n^5/80640-3418889*n^6/1995840;
|
|
alpha[6]=+212378941*n^6/319334400;
|
|
beta[1]=n/2-2*n^2/3+37*n^3/96-n^4/360-81*n^5/512+96199*n^6/604800;
|
|
beta[2]=n^2/48+n^3/15-437*n^4/1440+46*n^5/105-1118711*n^6/3870720;
|
|
beta[3]=17*n^3/480-37*n^4/840-209*n^5/4480+5569*n^6/90720;
|
|
beta[4]=4397*n^4/161280-11*n^5/504-830251*n^6/7257600;
|
|
beta[5]=4583*n^5/161280-108847*n^6/3991680;
|
|
beta[6]=+20648693*n^6/638668800;
|
|
|
|
Notation of output matches that of
|
|
|
|
L. Krueger,
|
|
Konforme Abbildung des Erdellipsoids in der Ebene
|
|
Royal Prussian Geodetic Institute, New Series 52, 172 pp. (1912).
|
|
|
|
with gamma replaced by alpha.
|
|
|
|
Alter maxpow to generate more or less terms for the series
|
|
approximations to the forward and reverse projections. This has been
|
|
tested out to maxpow = 30; but this takes a long time (see below).
|
|
|
|
*/
|
|
|
|
/* Timing:
|
|
maxpow time
|
|
4 2s
|
|
6 5s
|
|
8 11s
|
|
10 24s
|
|
12 52s
|
|
20 813s = 14m
|
|
30 13535s = 226m
|
|
*/
|
|
|
|
maxpow:8$ /* Max power for forward and reverse projections */
|
|
/* Notation
|
|
|
|
e = eccentricity
|
|
e2 = e^2 = f*(2-f)
|
|
n = third flattening = f/(2-f)
|
|
phi = (complex) geodetic latitude
|
|
zetap = Gauss-Schreiber TM = complex conformal latitude
|
|
psi = Mercator = complex isometric latitude
|
|
zeta = scaled UTM projection = complex rectifying latitude
|
|
*/
|
|
|
|
taylordepth:6$
|
|
triginverses:'all$
|
|
/*
|
|
revert
|
|
var2 = expr(var1) = series in eps
|
|
to
|
|
var1 = revertexpr(var2) = series in eps
|
|
|
|
Require that expr(var1) = var1 to order eps^0. This throws in a
|
|
trigreduce to convert to multiple angle trig functions.
|
|
*/
|
|
reverta(expr,var1,var2,eps,pow):=block([tauacc:1,sigacc:0,dsig],
|
|
dsig:ratdisrep(taylor(expr-var1,eps,0,pow)),
|
|
dsig:subst([var1=var2],dsig),
|
|
for n:1 thru pow do (tauacc:trigreduce(ratdisrep(taylor(
|
|
-dsig*tauacc/n,eps,0,pow))),
|
|
sigacc:sigacc+expand(diff(tauacc,var2,n-1))),
|
|
var2+sigacc)$
|
|
|
|
/* Expansion for atan(x+eps) for small eps. Equivalent to
|
|
taylor(atan(x + eps), eps, 0, maxpow) but tidied up a bit.
|
|
*/
|
|
atanexp(x,eps):=''(ratdisrep(taylor(atan(x+eps),eps,0,maxpow)))$
|
|
|
|
/* Convert from n to e^2 */
|
|
e2:4*n/(1+n)^2$
|
|
|
|
/* zetap in terms of phi. The expansion of
|
|
atan(sinh( asinh(tan(phi)) + e * atanh(e * sin(phi)) ))
|
|
*/
|
|
zetap_phi:block([psiv,tanzetap,zetapv,qq,e],
|
|
/* Here qq = atanh(sin(phi)) = asinh(tan(phi)) */
|
|
psiv:qq-e*atanh(e*tanh(qq)),
|
|
psiv:subst([e=sqrt(e2),qq=atanh(sin(phi))],
|
|
ratdisrep(taylor(psiv,e,0,2*maxpow)))
|
|
+asinh(sin(phi)/cos(phi))-atanh(sin(phi)),
|
|
tanzetap:subst([abs(cos(phi))=cos(phi),sqrt(sin(phi)^2+cos(phi)^2)=1],
|
|
ratdisrep(taylor(sinh(psiv),n,0,maxpow)))+tan(phi)-sin(phi)/cos(phi),
|
|
zetapv:atanexp(tan(phi),tanzetap-tan(phi)),
|
|
zetapv:subst([cos(phi)=sqrt(1-sin(phi)^2),
|
|
tan(phi)=sin(phi)/sqrt(1-sin(phi)^2)],
|
|
(zetapv-phi)/cos(phi))*cos(phi)+phi,
|
|
zetapv:ratdisrep(taylor(zetapv,n,0,maxpow)),
|
|
expand(trigreduce(zetapv)))$
|
|
|
|
/* phi in terms of zetap */
|
|
phi_zetap:reverta(zetap_phi,phi,zetap,n,maxpow)$
|
|
|
|
/* Mean radius of meridian */
|
|
a1:expand(integrate(
|
|
ratdisrep(taylor((1+n)*(1-e2)/(1-e2*sin(phi)^2)^(3/2),
|
|
n,0,maxpow)),
|
|
phi, 0, %pi/2)/(%pi/2))/(1+n)$
|
|
|
|
/* zeta in terms of phi. The expansion of
|
|
zeta = pi/(2*a1) * int( (1-e^2)/(1-e^2*sin(phi)^2)^(3/2) )
|
|
*/
|
|
zeta_phi:block([zetav],
|
|
zetav:integrate(trigreduce(ratdisrep(taylor(
|
|
(1-e2)/(1-e2*sin(phi)^2)^(3/2)/a1,
|
|
n,0,maxpow))),phi),
|
|
expand(zetav))$
|
|
|
|
/* phi in terms of zeta */
|
|
phi_zeta:reverta(zeta_phi,phi,zeta,n,maxpow)$
|
|
|
|
/* zeta in terms of zetap */
|
|
/* This is slow. The next version speeds it up a little.
|
|
zeta_zetap:expand(trigreduce(ratdisrep(
|
|
taylor(subst([phi=phi_zetap],zeta_phi),n,0,maxpow))))$
|
|
*/
|
|
zeta_zetap:block([phiv:phi_zetap,zetav:expand(zeta_phi),acc:0],
|
|
for i:0 thru maxpow do (
|
|
phiv:ratdisrep(taylor(phiv,n,0,maxpow-i)),
|
|
acc:acc + expand(n^i * trigreduce(ratdisrep(taylor(
|
|
subst([phi=phiv],coeff(zetav,n,i)),n,0,maxpow-i))))),
|
|
acc)$
|
|
|
|
/* zetap in terms of zeta */
|
|
/* This is slow. The next version speeds it up a little.
|
|
zetap_zeta:expand(trigreduce(ratdisrep(
|
|
taylor(subst([phi=phi_zeta],zetap_phi),n,0,maxpow))))$
|
|
*/
|
|
zetap_zeta:block([phiv:phi_zeta,zetapv:expand(zetap_phi),acc:0],
|
|
for i:0 thru maxpow do (
|
|
phiv:ratdisrep(taylor(phiv,n,0,maxpow-i)),
|
|
acc:acc + expand(n^i * trigreduce(ratdisrep(taylor(
|
|
subst([phi=phiv],coeff(zetapv,n,i)),n,0,maxpow-i))))),
|
|
acc)$
|
|
|
|
printa1():=block([],
|
|
print(concat("A=",string(a/(n+1)),"*(",
|
|
string(taylor(a1*(1+n),n,0,maxpow)),");")),
|
|
0)$
|
|
printtxf():=block([alpha:zeta_zetap,t],
|
|
for i:1 thru maxpow do (t:coeff(alpha,sin(2*i*zetap)),
|
|
print(concat("alpha[",i,"]=",string(taylor(t,n,0,maxpow)),";")),
|
|
alpha:alpha-expand(t*sin(2*i*zetap))),
|
|
/* should return zero */
|
|
alpha:alpha-zetap)$
|
|
printtxr():=block([beta:zetap_zeta,t],
|
|
for i:1 thru maxpow do (t:coeff(beta,sin(2*i*zeta)),
|
|
print(concat("beta[",i,"]=",string(taylor(-t,n,0,maxpow)),";")),
|
|
beta:beta-expand(t*sin(2*i*zeta))),
|
|
/* should return zero */
|
|
beta:beta-zeta)$
|
|
|
|
printseries():=[printa1(),printtxf(),printtxr()]$
|
|
|
|
/* Define array functions which are be read by tm.mac */
|
|
defarrayfuns():=block([aa:a1*(1+n),alpha:zeta_zetap,beta:zetap_zeta,t],
|
|
for i:1 thru maxpow do (
|
|
define(a1_a[i](n),ratdisrep(taylor(aa,n,0,i))/(1+n)),
|
|
t:expand(ratdisrep(taylor(alpha,n,0,i))),
|
|
define(zeta_a[i](zp,n),
|
|
zp+sum(coeff(t,sin(2*k*zetap))*sin(2*k*zp),k,1,i)),
|
|
t:diff(t,zetap),
|
|
define(zeta_d[i](zp,n),
|
|
1+sum(coeff(t,cos(2*k*zetap))*cos(2*k*zp),k,1,i)),
|
|
t:expand(ratdisrep(taylor(beta,n,0,i))),
|
|
define(zetap_a[i](z,n),
|
|
z+sum(coeff(t,sin(2*k*zeta))*sin(2*k*z),k,1,i)),
|
|
t:diff(t,zeta),
|
|
define(zetap_d[i](z,n),
|
|
1+sum(coeff(t,cos(2*k*zeta))*cos(2*k*z),k,1,i))))$
|
|
|
|
printseries()$
|
|
(b1:a1,
|
|
for i:1 thru maxpow do (alp[i]:coeff(zeta_zetap,sin(2*i*zetap)),
|
|
bet[i]:coeff(expand(-zetap_zeta),sin(2*i*zeta))))$
|
|
load("polyprint.mac")$
|
|
printtm():=block([macro:GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER],
|
|
value1('(b1*(1+n)),'n,2,0),
|
|
array1('alp,'n,1,0),
|
|
array1('bet,'n,1,0))$
|
|
printtm()$
|
|
/*
|
|
defarrayfuns()$
|
|
save("tmseries.lsp",maxpow,arrays)$
|
|
*/
|