楕円弧長計算

こんにちは。
楕円弧長は、下記の第二種不完全楕円積分で表されます。与えられた緯度に対し赤道からの子午線弧長です。

プログラムでは小さい扁平率を仮定し、倍精度浮動小数点数に現れない$n^7$（第三扁平率の冪乗）以上の級数展開項は近似で切り捨てています。このベッセルの級数式は優れた収束ですね。これの級数係数の計算については分数の表示（リスト）。参考：https://en.wikipedia.org/wiki/Meridian_arc#Expansions_in_the_third_flattening_.28n.29

\textrm{An elliptic arc length formula from a modification of Bessel (1825)}\\
\begin{aligned}
m(\varphi) &= \int_0^\varphi M(\varphi) d\varphi \\ 
 &=\frac{a}{1+n}\left( \int_0^\varphi \sqrt{1+2n\cos2\varphi+n^2} \, d\varphi - \frac{2n \sin2\varphi}{\sqrt{1 + 2n \cos2\varphi + n^2}} \right) \\
 &=\frac{a}{1+n}\left( c_0 \, \varphi + \sum_{l=1}^\infty c_l \, l^{-1}  \sin 2 l \varphi - \frac{2n \sin2\varphi}{\sqrt{1 + 2n \cos2\varphi + n^2}} \right) \\
M(\varphi) &= \frac{a(1-e^2)}{(1-e^2\sin^2\varphi)^{3/2}} \\
e^2 &= \frac{4 n}{\left(1+n\right)^2} \\
c_l &= n^l \sum_{j=0}^\infty n^{2j} \binom{1/2}{j} \, \binom{1/2}{j+l}  
\end{aligned}

import math
DEGREE = math.pi / 180
PIHALF = math.pi / 2

# Bessel (1825); elliptic (meridian) arc; incomplete integral
def meridanArc(lat, n):
    n2 = n * n
    coef = [1, 1/4., 1/64., 1/256.]
    ma = lat * poly(n2, coef)
    if lat != PIHALF:  # PIHALF: complete integral
        coef = [[1/2., -1/16., -1/128.], [-1/16., 1/64., 5/2048.], [1/48., -5/768.], [-5/512., 7/2048.], [7/1280.], [-7/2048.],]
        nn, cs, cslat2 = 1, csmake(0), csmake(lat*2)
        for c in coef:
            nn *= n  # n**(i+1)
            cs = csadd(cs, cslat2)  # cos(lat*2*(i+1))
            ma += nn * cs['sin'] * poly(n2, c)
        ma -= 2*n*cslat2['sin']/math.sqrt(1+2*n*cslat2['cos']+n2)
    return ma / (1 + n)

def poly(x, coef):
    sum = 0.0
    for a in coef[::-1]:
        sum = a + sum * x
    return sum

def csmake(arg):
    c, s = math.cos(arg), math.sin(arg)
    return {'cos': c, 'sin': s}

def csadd(cs1, cs2):
    c = cs1['cos']*cs2['cos'] - cs1['sin']*cs2['sin']
    s = cs1['sin']*cs2['cos'] + cs1['cos']*cs2['sin']
    return {'cos': c, 'sin': s}

def main():
    RE = 6378137.0  # GRS80
    f = 1/298.257222101  # GRS80
    n = f / (2 - f)   # the third flattening
    for lat in range(0,100,10):
        ma = meridanArc(lat*DEGREE, n)
        print "%4.1f degree: %12.3f m" % (lat, RE*ma)

出力結果

 0.0 degree:        0.000 m
10.0 degree:  1105854.833 m
20.0 degree:  2212366.254 m
30.0 degree:  3320113.398 m
40.0 degree:  4429529.030 m
50.0 degree:  5540847.042 m
60.0 degree:  6654072.819 m
70.0 degree:  7768980.728 m
80.0 degree:  8885139.872 m
90.0 degree: 10001965.729 m