Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
Help us understand the problem. What is going on with this article?

楕円弧長計算

More than 1 year has passed since last update.

こんにちは。
楕円弧長の計算、下記の第二種不完全楕円積分です。与えられた緯度に対し赤道からの子午線弧長です。プログラムでは小さい扁平率を仮定し、倍精度浮動小数点数に現れない$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
kkdd
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away