1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
|
/*
Compute series approximations for Transverse Mercator Projection
Copyright (c) Charles Karney (2009-2010) <karney@alum.mit.edu> 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)$
*/
|
