諸行無常は流転しない

諸行無常である、というこのこと自体は、常である。『万物流転の法則』自体は流転しない。そんな流転しない物事をできるだけ見ていたいね、ということで数学や物理について学んだりしようと思っているブログです。

回転楕円体の子午線曲率半径と卯酉(ぼうゆう)線曲率半径の導出の仕方

 

国土地理院の作業規定の準則の計算式集にある

「1.2 楕円体の諸公式」

にある子午線曲率半径 Mおよび、卯酉(ぼうゆう)線曲率半径 Nを導出する記事である。

 

これから、ガウス・クリューゲル等角投影法などを理解しようと試みる人への手助けである。

 

 

子午線および卯酉線

子午線曲率半径および卯酉線曲率半径の意味

微分幾何学の空間曲線論の話になる。

数学では、曲率を以下のように定義している。

 

空間内のある曲線のある点での曲率とは、

その点において、曲線を円で近似したときの半径の逆数である。

 

円で近似したとき、その円の半径を曲率半径といい、

その円を、曲率円という。

 

 

つまり、

曲率=1/(曲率半径)

である。

 

例としては、

直線は、曲率0である。

曲率半径は、無限大。

つまり、直線は、半径無限大の円と考えるわけである。

 

真円は、

曲率円はその円であり、曲率半径はその円の半径である。

 

道路などにあるR=400は、

次のカーブは半径400mの円の円弧ですよ、

という意味である。

 

逆数の1/400が曲率になる。

 

子午線曲率半径と卯酉(ぼうゆう)線曲率半径の導出

それでは、(地理学的、測地学的)緯度 \phiにおける子午線曲率半径 Mを求めていく。

(地理学的)緯度 \phiは、地球楕円体上の点 Pにおける接平面に直行する線(法線)と赤道面とのなす角である。

楕円をパラメトライズするときの、パラメータ \thetaとは異なるので注意。

これが話をややこしくしている。

(なんで緯度と楕円のパラメタをわざわざ別のものにするのだろう?)

 

地球を回転楕円体と見立てて、

下図のように座標を入れて考える。

 

地球回転楕円体

子午線・平行圏・卯酉線(ぼうゆうせん) | 地理空間情報技術ミュージアム Museum of GIS Technolog.

から拝借

 

しかし、子午線曲率半径を考えるには、子午線だけを考えればよく、

結局のところ、楕円だけを考えればいいことになる。

 

というわけで、楕円の曲率を求める。

楕円の曲率を導出する

弧長パラメタとは限らない任意のパラメタに関して、

空間曲線 \gamma = c(t)の曲率 \kappa(t)は以下のように書ける。

 \displaystyle \kappa (t) = \frac{\det ( c'(t), c''(t) ) } { c(t)^{3/2}}

この公式を利用する。

楕円 Eを、

 E(\theta) = ( a \cos \theta , b \sin \theta ), a \gt b \gt 0 , 0 \le \theta \lt 2 \pi

とパラメトライズする。

1階と2階の微分を計算して、

 E'(\theta) = ( - a \sin \theta , b \cos \theta )

 E''(\theta) = ( a \cos \theta , b \sin \theta )

空間曲線の曲率の公式に代入して、

 \displaystyle \kappa (\theta) = \frac{ ab } { ( a^2 \sin ^2  \theta + b^2 \cos ^2 \theta )^{3/2}}

となる。

ここで、

離心率 \displaystyle e = \sqrt{\frac { a^2 -b^2}{a^2}}と、

 \cos ^2 \theta + \sin ^2 \theta = 1

を用いて、ルートの中を整理すると、

 \displaystyle \kappa (\theta) = \frac{ b }{a^2 (1 - e^2 \cos ^2 \theta)^{3/2}}

となる。

ここで、 \sinを消去して \cosで整理したが、これは後で楕円のパラメタ \thetaから緯度 \phiに戻すときに、

作業既定の準則の付録の計算式集の、

 W = \sqrt { 1 - e^2 \sin ^2 \phi}

とつじつまが合うので、どちらで整理してもよい。

 

これで楕円の曲率(したがってその逆数の曲率半径)が求められた。

 

楕円のパラメタと緯度の換算式

この楕円のパラメタ \thetaに関する曲率の式を、

緯度 \phiに関する式に換算するために、 \theta \phiの関係式を求める。( \cos \thetaさえわかればいい)

 

緯度 \phiと回転楕円体の接平面が直行するので、

接ベクトル E'(\theta) = (-a \sin \theta , b \cos \theta)と、

緯度の方向ベクトル (\cos \phi ,\sin \phi)が直行することを使って、

 

緯度とパラメタ

内積 -a \sin \theta \cos \phi + b \cos \theta \sin \phi = 0

となる。

 a \sin \theta \cos \phi = b \cos \theta \sin \phi

だから2乗して、 \sin \theta \cos \theta に直して、

 a^2 (1 - \cos ^2 \theta ) \cos^2 \phi = b^2 \cos ^2 \theta \sin ^2 \phi

より、

 \displaystyle \frac { a^2  \cos ^2  \theta  }{ b^2 \sin ^2 \phi } ( 1- \cos ^2 \theta) = \cos ^2 \theta

だから、

 \displaystyle  \cos ^2 \theta = \frac { \frac {a^2 \cos ^2 \phi }{b^2 \sin ^2 \phi}}{1 + \frac{a^2 \cos ^2 \phi }{b^2 \sin^2 \phi }} = \frac{a^2 \cos ^2 \phi }{a^2 \cos ^2 \phi + b^2 \sin ^2 \phi } = \frac {\cos ^2 \phi}{1- e^2 \sin ^2 \phi}

となる。

 

結局、

 \displaystyle  \cos ^2 \theta = \frac {\cos ^2 \phi}{1- e^2 \sin ^2 \phi}

を得る。

 

(厳密なことを書くと、 \sin \phi = 0のとき、 \theta = \phiとなるので、 \theta \ne \phiとしておく。ゼロ割防止のため)

 

子午線曲率半径 Mを求めよう

以上により得られた楕円の曲率

 \displaystyle \kappa (\theta) = \frac{ b }{a^2 (1 - e^2 \cos ^2 \theta)^{3/2}}

とパラメタと緯度の換算式

 \displaystyle  \cos ^2 \theta = \frac {\cos ^2 \phi}{1- e^2 \sin ^2 \phi}

より、緯度 \phiにおける子午線曲率は、

 \displaystyle \kappa (\phi) = \frac{b}{a^2 (1- e^2 \cdot \frac{\cos ^2 \phi}{1- e^2 \sin ^2 \phi})^{3/2}} = \frac { b (1 - e^2 \sin ^2 \phi)^{3/2} }{a^2 (1-e^2)^{3/2} }=\frac{(1-e^2 \sin ^2\phi)^{3/2}}{a(1-e^2)}

となる。

ここで、

離心率 \displaystyle e = \sqrt{\frac { a^2 -b^2}{a^2}}より、

 \displaystyle   1-e^2 = \frac{b^2}{a^2}を使って、

 b を消去した。

 

卯酉線曲率半径 Nを求める

卯酉(ぼうゆう)線が何かということと、卯酉円曲率半径が以下の図の赤線部になることは前提として、この赤線部の長さを求めていく。

卯酉線曲率半径

経度 \phiと楕円のパラメタ \thetaの換算式が求まった今、それほどややこしいことはない。

流れとしては、赤線を表示するベクトル方程式を求めて、

 y軸との交点を求めて、線長を計算するだけである。

赤線のベクトル方程式は、

楕円上の点 E(\theta)=(a \cos \theta , b \sin \theta)を通り、

傾きは、接ベクトル E'(\theta)=(-a \sin \theta , b \cos \theta)を反時計回りに90度回転させたものになるので、

 (X, Y) = t \cdot ( -a \cos \theta, -b \sin \theta) + (a \cos \theta , b \sin \theta)

となる。( tは任意の実数)

 X  = 0を代入して、

 t= \dfrac{a }{b}を得る。

よって、赤線部の線長 Nは、

 \begin{align} \displaystyle N  = \int_{0}^{\dfrac{a }{b}}  \sqrt{ (\frac{dX}{dt})^2+(\frac{dY}{dt})^2 }dt \\ = \int_{0}^{\dfrac{a }{b}} \sqrt{b^2 \cos ^2 \theta + a^2 \sin ^2 \theta}dt  \\  =  \frac{a}{b} \sqrt{b^2 \cos ^2 \theta + a^2 \sin^2 \theta} \\  = \frac{a^2}{b} \sqrt{1- e^2 \cos ^2 \theta} \\  =  \displaystyle \frac{a^2}{b} \sqrt{1- \frac{ e^2 \cdot\cos^2 \phi}{1 - e^2 \sin ^2 \phi}} \\  =   \frac{a^2}{b} \sqrt{\frac{1 - e^2}{1 - e^2 \sin ^2 \phi}}\\ =  \displaystyle \frac{a}{\sqrt{1 - e^2 \sin ^2 \phi}} \\ \end{align} 

 

(数式の=をそろえられなくてすいません。)

 

まとめ

子午線曲率が、

 \displaystyle \kappa (\phi) = \frac{(1-e^2 \sin ^2\phi)^{3/2}}{a(1-e^2)}

だから、

子午線曲率半径 Mはその逆数

 M = \dfrac{a(1-e^2)}{(1-e^2 \sin ^2\phi)^{3/2}}

卯酉線曲率半径 Nは、

 N =  \displaystyle \frac{a}{\sqrt{1 - e^2 \sin ^2 \phi}}

となる。