(original) (raw)

/* 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)$ */ /karney@alum.mit.edu