８．楕円の弧の長さ

楕円の弧の長さ

ある地点で、ある方位に進むことは、垂直截線で表され、その曲率によって角度と長さの関係が決まります。垂直截線のうち経線と卯酉線は、それぞれ曲率の最小、最大のようです。

地球を真球と見なす場合、垂直截線はすべて大円です。球の表面の2点A、Bを考えると、点Aにおける点Bを通る垂直截線は、点Bにおいても垂直截線です。地心と点A、Bは同じ平面にあって、大円の弧長が距離です。

重力だけが作用した球面上のパチンコ玉は大円上を慣性運動し、大円は平面における直線と考えられます。

準拠楕円体では、表面の2点A、Bの一方の垂直截線は他方の垂直截線ではありません。例外は同一経線上の２点と、赤道上の２点です。例外を除くと、点Bに向かって垂直截線上を点Aから出発した途端に、その地点と点Bを通る異なる新たな垂直截線上を進むことになります。

おそらく、これが測地線です。この計算は、短い区間の垂直截線の弧長の積算値です。

1. 同一経線上の２点の距離

準拠楕円体上の経線は楕円を描きます。同一経線上の２点の距離は、一般的な楕円の弧長の計算です。子午線弧長の差とも考えられます。

2. 赤道上の２点の距離

準拠楕円体は極軸を中心に回転した回転楕円体なので、赤道は真円です。赤道上の２点の距離は円の弧の長さです。しかし、180° 離れた2点の場合、極を通過する方が短いことは明らかです。

3. 子午線弧長

楕円の弧の長さを「子午線弧長」と呼ぶ場合は、準拠楕円体の経線上の弧長のことであり、位置を示すのは地理緯度だと考えられます。また、位置を1つだけ指定して常に緯度ゼロ（赤道）からの長さを計算することが習慣のようです。Wikipedia にあった計算です。

#子午線弧
#地理緯度(ただし、弧度)。赤道面からの長さを計算。
#河瀬和重の式(wikipediaの式から)
kawaseArc2 <- function(fa,a,b)
{
J <- 3; J2 <- J * 2 # 次数
fa2 <- fa * 2 # 緯度の2倍
n <- (a - b) / (a + b) # 第3扁平率
s <- 1:J2
C <- (1/s-4*s)*sin(fa2 * s) # C[]
E <- (3/2)*n/(1:J2) - n # E[]
A <- rep(E[1],J) # A[]
i <- 2
while(i <= J) { A[i] <- E[i] * A[i-1]; i = i + 1; } # 総積
A <- A^2 # 総積を2乗
m <- floor((1:J2)/2)
s <- seq(1,J2,by=2)
m[s] <- (-m[s]) # 奇数番目を負に
sum <- fa
for(j in 1:J)
{
j2 <- j * 2
D <- E[j+m[1:j2]]
s <- seq(1,j2,by=2)
D[s] <- 1 / D[s]
for(k in 2:j2) D[k] <- D[k] * D[k-1] # 総積
B <- sum(C[1:j2] * D)
sum <- sum + A[j] * (fa + B)
}
(a+b)/2 * sum # a/(1+n)==(a+b)/2
}

緯度45°、90° を計算してみます。

r$> a <- 6378137
f <- 1/298.257222101
b <- a*(1 - f)
φ <- 45 /180*pi
kawaseArc2(φ,a,b)
[1] 4984944.377857997
φ <- 90 /180*pi
r$> kawaseArc2(φ,a,b)
[1] 10001965.72923046

測量計算サイトでは、緯度45°に対して 4,984,944.378(m)となります。90°は計算しないことにしているようで、緯度89.99999999° に対して、10,001,965.728(m) となります。

経線の赤道から45° の弧の長さは、45°から90°までの45°分の弧の長さ 5017021.351372463 に対して 32076.9735 m 短く曲率の差を反映しています。

おそらく、子午線弧長の計算は、三角関数の積分を級数展開して置くことで非力な計算機でも高速に計算できるように考えられたものです。また、準拠楕円体の経線の弧長を計算するためのもので、様々な形状の垂直截線の弧長が計算できるかどうかは何も説明されていません。

4. 単純に

ある区間の楕円の弧の長さ S を、左図のように小区間の和で求めてみます。小区間の弧長 ds は、

$ds=\sqrt{{dx}^{2}+{dy}^{2}}$

また、原点が中心で、傾いていない、x軸半径がa、y軸半径がcの楕円の式は、

$\frac{{x}^{2}}{{a}^{2}}+\frac{{y}^{2}}{{c}^{2}}=1$

です。

また、更成緯度 θ で媒介変数表示すれば、

$x=a\cdot{}cos\theta{}$

$y=c\cdot{}sin\theta{}$

ある区間を分割して dθ ごとに ds を求めます。

三角関数は象限を跨ぐと私には分からなくなるので、 $0\leq{}\theta{}\leq{}\frac{\pi{}}{2}$ を考えます。

EllipseArcLength <- function(φ1,φ2, a=6378137, f=1/298.257222101){
b <- a*(1-f)
β1 <- atan((1-f)*tan(φ1)); β2 <- atan((1-f)*tan(φ2));
β <- seq(β1,β2, length.out=80000)
x <- a*cos(β)
y <- b*sin(β)
x <- x[2:length(x)] - x[2:length(x)-1]
y <- y[2:length(y)] - y[2:length(y)-1]
s <- sqrt(x^2+y^2)
sum(s)
}

関数の引数は地理的緯度 φ で、関数内は更成緯度 β を使用します。

EllipseArcLength(0, pi/4)
[1] 4984944.377837976
EllipseArcLength(0, pi/2)
[1] 10001965.72906979

8万分割していますが、mmまで、前述の計算結果と同じです。

5. 楕円積分

原点が中心で、傾いていない、x軸半径が1、y軸半径がcの楕円の式は、

${x}^{2}+\frac{{y}^{2}}{{c}^{2}}=1$

あるいは媒介変数 θ を使用して、

$x=cos\theta{}$

$y=c\cdot{}sin\theta{}$

また、左図のように楕円の周の長さ S を考えると、

$S=\int\limits_{\ }^{\ }ds=\int\limits_{\ }^{\ }\sqrt{{dx}^{2}+{dy}^{2}}$

のように表せます。

楕円の式の媒介変数表現を微分すると、

$xd=-sin\theta{}$

$yd=c\cdot{}cos\theta{}$

したがって、

$\sqrt{{dx}^{2}+{dy}^{2}}=\sqrt{{sin}^{2}\theta{}+{c}^{2}\cdot{}{cos}^{2}\theta{}}=\sqrt{{1-cos}^{2}\theta{}+{c}^{2}{\cdot{}cos}^{2}\theta{}}=\sqrt{1-\left({{1-c}^{2}}\right){cos}^{2}\theta{}}$

また、x軸半径が1、y軸半径がcの楕円の楕円の離心率 k は、

${k}^{2}=\frac{1-{c}^{2}}{1}=1-{c}^{2}$

したがって、以下のような積分で、離心率 k の0°からφまでの弧長が計算されるものと思います。

$E\left({\varphi{},k}\right)=\int\limits_{0}^{\pi{}/2}\sqrt{1-{k}^{2}{cos}^{2}\theta{}}d\theta{}$

r$> func <- function(β){
f <- 1/298.257222101
k2 <- f*(2-f)
sqrt(1 - k2*cos(β)^2)
}
a <- 6378137
f <- 1/298.257222101
integrate(func,0,atan((1-f)*tan(pi/4)))$value * a
[1] 4984944.377857996
integrate(func,0,pi/2)$value * a
[1] 10001965.72923047

しかし、Wikipediaの「楕円積分」にあるルジャンドルの標準形の第2種楕円積分は、

$E\left({\varphi{},k}\right)=\int\limits_{0}^{\pi{}/2}\sqrt{1-{k}^{2}{sin}^{2}\theta{}}d\theta{}$

です。

更成緯度 β を、極軸から測る、極角だと考えてみます。

楕円は媒介変数 β を使用して、

$x=cos\left({\frac{\pi{}}{2}-\beta{}}\right)=sin\beta{}$

$y=c\cdot{}sin\left({\frac{\pi{}}{2}-\beta{}}\right)=c\cdot{}cos\beta{}$

と表され、

$\sqrt{{dx}^{2}+{dy}^{2}}=\sqrt{{cos}^{2}\theta{}+{c}^{2}\cdot{}{sin}^{2}\theta{}}=\sqrt{1-\left({{1-c}^{2}}\right){sin}^{2}\theta{}}$

おそらく、これがルジャンドルの標準形の第2積分です。

楕円積分関数は C++ の標準ライブラリにもあります。ellint_2()が第二種不完全楕円積分、comp_ellint_2() が第二種完全楕円積分と名前が付いています。それぞれ、以下のような説明が付いています。

$\\E(k, \phi) = \int_0^\phi \mathrm d\theta ~ \sqrt{1 - k^2 \sin^2 \theta} \quad \text{for } |k| \le 1 \\E(k) = E(k, \pi/2) = \int_0^{\pi/2} \mathrm d\theta ~ \sqrt{1 - k^2 \sin^2 \theta} \quad \text{for } |k| \le 1 \\E(k, \pi/2) = E(k)$

完全積分は、積分範囲を指定せず、常に0°から90°までの弧の長さを返します。

プログラミング言語で使用される楕円積分が極角だと考えると上手く行くことをR言語で試しておきます。完全積分の結果から、45°までの不完全積分の結果を引くと、赤道からの45°の弧長に一致します。

> a <- 6378137
> f <- 1/298.257222101
> b <- a*(1-f)
> k <- sqrt(f*(2-f))
> ellint_Ecomp(k) * a
[1] 10001965.72923046
> θ <- pi/2-atan((1-f)*tan(pi/4))
> ellint_E(θ,k) * a
[1] 5017021.351372467
> (ellint_Ecomp(k) - ellint_E(θ,k)) * a
[1] 4984944.377857994

6. 確かめる

離心率 k と、その2乗のどちらを引数にするか混乱したので確認しておきます。

長半径 a=1、扁平率 f=0.3 とします。離心率 k と、その2乗の k2 は以下のようになります。

> f <- 0.3; a <- 1; b <- 1-f
> k2 <- f*(2-f); k <- sqrt(k2)
> c(k2, k)
[1] 0.5100000000 0.7141428429

次のような関数を作って、力ずくで弧長を計算します。

EllipseArcLength0 <- function(a,b, β1, β2, N=1000){
t <- seq(β1, β2,length.out=N)
x <- a*cos(t); y <- b*sin(t)
x <- x[2:N]-x[1:(N-1)]; y <- y[2:N]-y[1:(N-1)];
sum(sqrt(x^2+y^2))
}

細かく区切って弦長を計算し総和を求めて、正しい値の目安にします。

> c(EllipseArcLength0(1,1-f, pi/4, pi/2, N=10000),
+ EllipseArcLength0(1,1-f, 0, pi/4, N=10000))
[1] 0.7474074299 0.5981848151

gslの楕円積分関数は以下のようになります。

> u <- ellint_Ecomp(k)
> v <- ellint_E(pi/4,k)
> c(u,v,u-v)
[1] 1.3455922454 0.7474074300 0.5981848154

完全積分 ellint_Ecomp() は、0°から90°の周長の 1/4 を返します。 ellint_E()に指定する１つの角度は極を0°としていると考えると話しが合います。π/4 を指定して得られる値 v と 1/4周の長さの差は、v より短い値です。

結果の3つの値は普通の緯度の表示で、0°から90°の長さ、45°から90°の長さ、0°から45°の長さと解釈しました。

gslの楕円積分関数は、離心率 k を引数にしています。

普通の積分関数を使うと、象限を跨いでも計算可能だと思います。積分対象の関数に離心率の2乗を渡すとすれば、以下のようになります。

> func <- function(θ, k2){ sqrt(1 - k2*cos(θ)^2) }
> c(integrate(func,pi/4,pi/2,k2)$value,
+ integrate(func,0,pi/4,k2)$value)
[1] 0.7474074300 0.5981848154