１２．航程線とメルカトル図法

航程線とメルカトル図法

地図に線を引いて行程とすることは何の不思議もありません。ただ、陸地ではまっすぐ進めないことが多いので、引いた直線は方向や直線距離を与えることになります。

海上や航空機も、まっすぐ進めるわけではないでしょうが、海図に線を引くのは航程や位置の把握が目的であって、地図に描いた線が距離や方位を示すのは当然です。

航程線が Rhumb line の訳語として使用が始まったのかどうかは分かりませんが、幕末以降の日本人には海図上の線は Rhumb line でした。

1. 地図について

地図はヒッパルコスの時代から経緯度で方眼に描かれていたのだろうと思います。すでにヒッパルコスの時代（紀元前150年ごろ）には、角度は360°で測られ、月食を利用して経度を測り、地球の大きさや月までの距離を倍率で割り出していたと見られています。

地球を球と考えれば、中心角 1° に相当する地表の弧の長さはどこでも同じです。しかし、1° の長さを決定することは1669年から1670年にかけて行われた「ピカールの三角測量」のようです。このとき、緯度 1° は 112.21km と割り出されました。

ヨーロッパは天測や測量技術を失っていたと見られ、プトレマイオスの「アルマゲスト」が12世紀、「地理学」は15世紀になって再認識されました。ヨーロッパの地図は、「地理学」に収録されていた古代の8000点の経緯度データとプトレマイオス図によって再開されました。ローマの時代には既にガリア、ゲルマニア、ブリタニアは重要な地域になっており、昼の長さから緯度を推定するなどの工夫によって広くデータが集められていたようです。

ヒッパルコスの業績はプトレマイオスによる引用で知られるもののようです。三角関数などの数表も含まれます。

ロドス島のパン岬は、経度58°、緯度35°55（35°+1/2+1/3+1/12）のように記録されました。

プトレマイオス図は扇形で、経度差 1° の緯線の長さが北極に向かって短くなることを考慮したものでした。緯線も弧を描いています。

左図のように経度差180°を考えると、緯線の弧の長さは、緯度 φ に応じて、aπcosφ と変化します。

方位磁石の使用は11世紀以降と見られ、それまでの地図は俯瞰に近いほど良いと考えられたものと推測します。

緯度や経度を測定することに比べると方位磁石は遥かに手軽に常時利用できます。

中心角1°に相当する距離は、まだ、正確には知られていませんでしたが比率としては有効であり、方位磁石の使用は地図を実用的なものに変えたのだと思います。

メルカトル図法は1569年に出版された地図に使用された投影法を指すようです。

縄を張って長さを測ることのできない天体の位置を角度で表すことは天を球と見なすことです。地球もまた球であると考えた人々が遠い位置関係を緯度、経度で表しました。

地上を俯瞰することはほとんど不可能だったと思うのですが、天球はいつも曲面を俯瞰しています。星座は見たままに描かれたのだろうと思います。天の北極を中心に放射状に星座を配置した星図が作られるのも想像できます。また、天球儀や渾天儀が知られるように実際に球に描いていました。

しかし、恒星の間を惑星が動くといったことは、方眼に星を配置して考えたに違いありません。地図も同様で、研究者は方眼にデータを描いたと思います。メルカトル図法の説明には円筒投影が引かれていますが、実際に円筒投影で地図を作る必要があったとは思えません。

2. 方位角

方位磁石に関するヨーロッパの最初の記録は12世紀末で、それ以降不可欠な道具となったようです。方位磁石は天の北極を指す訳ではありませんが、ここでは自転軸の北を指すものとして、球体の地球を考えます。

方位はその位置での球の接平面上に表されます。その地点と両極を通る大円が、その地点を通る経線です。この経線が作る平面は球の中心を含みます。この平面と、接平面の交線が北南を示します。

バンクーバーと香港を結ぶ太平洋航路の運航は1887年に開始され横浜は経由地だったようです。

横浜から見たバンクーバーの方位角は、横浜を通る経線が作る平面と、バンクーバーを含んだ横浜における法平面の、2つの平面の成す角です。

球を考えると、球の中心と横浜を結んだ直線は、横浜の接平面の法線です。この直線を含む平面はすべて横浜における法平面で、バンクーバーを含むものを取ります。

方位角は、44.2°ほどと算出されます。

横浜とバンクーバーを結んだ方位を示す弧は、途中で交差する経線とは異なった角度で交差しています。

横浜では方位角44.2°ですが、進むにしたがって方位角は大きくなっています。

目標に向かって進めば、この方位線上を進むことになりますが、現実には目標が見えません。

3. 航程線

実際に航程線を引いてみます。GeographicLib には RhumbSolve と言うコマンドラインツールがあって、使用例に従って、航程線の通過点を列挙できます。球とみなして扁平率 f=0 で計算しています。

青い線が航程線です。赤は2地点を通る大円です。

航程線は、経線と常に同じ角度で交差する「等角航路」と説明されています。

同じ緯度経度のデータを平面に描いてみます。

メルカトル図は高緯度ほど緯度1°の長さが長く描かれるので、方眼に描いた地図より傾きが小さくなっています。

また、方眼に描いた航程線は弧を描いているように見えます。

左図は、方眼の図の航程線を拡大したもので、黒い線は直線です。

2点間の方位角は44.2°ほどですが、RhumbSolveコマンドは、azi12の値として「079:02:51.2」と表示します。

メルカトル図に描いた航程線の傾きは、10.95245° と計算されます。横浜を通る経線の北からは、79.04755°です。

メルカトル図での航程線は、経線を平行に描くので、どこでも、この角度で交差します。

一方、球体での交差角は一様ではありません。メルカトル図で直線となる航程線は、実際には常に曲がり続けることですが RhumbSolveコマンドの出力は一様ではありません。2点間を60等分し、隣り合う次の点への方位角を調べると、78.34716° から78.61467° の範囲の値でした。

隣り合う2点の方位角は法平面の截り線で測ったものですが、航程線はその間でも曲がり続ける線ですから、区間を短くすれば少し大きくなり、メルカトル図での79.04755°に近づくと言うことだと思います。

600等分して、最初の区間の方位角を計算すると79.00293°となります。

4. 横に縦を合わせる

経線を平行に引いた地図を描くことは、最初の図のように、高緯度ほど横が拡大されることです。地球の赤道半径を a とすれば、赤道の経度差180°の弧の長さは aπ です。緯度φでは aπcosφ と短くなります。どちらも経度差180°なので、平行な経線に従って同じ長さに描くと緯度φでは 1/cosφ に引き延ばされることになります。

方眼にプロットして位置関係を把握していたと考えると、縦も同じだけ引き延ばせば良いと思うはずです。

天体と異なり、誰も俯瞰することが不可能に思われた地上は、見た目より定規を当てて測れる図が重要だったと思います。

球面に正方形を2つ描いてみます。正方形の辺は大円の弧で全て同じ長さです。赤道が通る下の正方形は、経度差30°の経線が左右の辺になっています。向かい合う辺を作っている大円は4組とも30°の角度で交差しています。下の正方形の中心は (1,0,0) にあって、これをY軸の周りに60°回転したものが上の正方形です。「球面の正方形」

説明のために、球面上の正方形の中央に十字を描いておきます。

この2つの球面上の正方形をXY直交座標にプロットしてみます。おそらく、方眼にプロットすることは常に行われることです。

これは、経線を平行に描くことです。緯線も水平の直線に描かれることになります。

緯度の1°と経度の1°を均等に描けば、上の中央の図になります。

球面では極に集まる経線を平行に描くので緯度が大きくなるほど横の方向に拡大して描かれることになります。

球面の正方形の辺は大円の弧で、弧AB、DCはそれぞれ15°、ー15°の経線上にあります。中央の方眼に描いた図では、経線は垂直な平行線となります。赤道上の線分EFの長さは π/6 です。緯度が60°の緯線上では、同じ経線の幅 m は、

$\frac{\pi{}}{6}cos\frac{\pi{}}{3}=\frac{\pi{}}{12}$

と半分の長さを表すことになります。

緯度60°では横方向が2倍に拡大されますが、点T、U は60°の緯線より南にあって、そこ間隔は π/3 より少し短くなっています。①、②とした球面上の正方形の辺と60°の緯線の交点は、60°の緯線上にありますが、もともと弧TUより長いので π/3 より少し長くなります。

緯度は球と同様に均等なので、線分GHとVWは同じ長さに描かれます。

この図に定規を当てて使う場合は、横方向に測った値に cosφ を乗じることになります。

上の右の図はメルカトル図法の変換式でY軸方向を引き伸ばして描いた図です。

図の横方向は 1/cosφ に伸びますが、縦を同じように引き延ばすには緯度ごとに異なった倍率を乗じなければなりません。

線分ABの長さは、以下のように表現されます。

$\overline{AB}=\int\limits_{{\varphi{}}_{A}}^{{\varphi{}}_{B}}\frac{1}{cos\varphi{}}d\varphi{}$ (4.1)

一方、作図はメルカトル図法の変換式によっています。

$y=log\left({tan\left({\frac{\pi{}}{4}+\frac{\varphi{}}{2}}\right)}\right)$ (4.2)

図の線分ABの長さは、以下のように計算されます。

$\overline{AB}=log\left({tan\left({\frac{\pi{}}{4}+\frac{{\varphi{}}_{B}}{2}}\right)}\right)-log\left({tan\left({\frac{\pi{}}{4}+\frac{{\varphi{}}_{A}}{2}}\right)}\right)$

> φA <- -0.2532616
> φB <- 0.2532616
> integrate(function(φ){1/cos(φ)},φA,φB)$value
[1] 0.5120265
> log(tan(pi/4+φB/2))-log(tan(pi/4+φA/2))
[1] 0.5120265

点Bの緯度φBのメルカトル図のy座標の値は、

$y=\int\limits_{0}^{{\varphi{}}_{B}}\frac{1}{cos\varphi{}}d\varphi{}$ (4.3)

と計算されるはずです。

> integrate(function(φ){1/cos(φ)},0,φB)$value
[1] 0.2560133
> log(tan(pi/4+φB/2))
[1] 0.2560133

どうやら、

$y=\int\limits_{0}^{\varphi{}}\frac{1}{cos\varphi{}}d\varphi{}=log\left({tan\left({\frac{\pi{}}{4}+\frac{\varphi{}}{2}}\right)}\right)$

のようです。

また、式(4.2)から、

${e}^{y}=tan\left({\frac{\pi{}}{4}+\frac{\varphi{}}{2}}\right)$

$\frac{\pi{}}{4}+\frac{\varphi{}}{2}=atan\left({{e}^{y}}\right)$

$\varphi{}=2\cdot{}atan\left({{e}^{y}}\right)-\frac{\pi{}}{2}$ (4.4)

と逆変換できます。

球面の正方形の中に描いた十字の縦棒は、2つとも同じ経線上にあります。球の表面の点G、H、V、W の緯度を、それぞれ φG、φH、φV、φW とします。2つの縦棒（線分GH、VW）の長さは同じですが、メルカトル図にすると、VWが長く描かれます。

2つの正方形の中心は緯度ゼロと60°にあるので2倍ほどになるはずです。

> (φH-φG)==pi/6 # 30°
[1] TRUE
> (φW-φV)==pi/6 # 30°
[1] TRUE
> yG <- log(tan(pi/4+φG/2))
> yH <- log(tan(pi/4+φH/2))
> yV <- log(tan(pi/4+φV/2))
> yW <- log(tan(pi/4+φW/2))
> (yH-yG)/pi*180 # 度
[1] 30.34869
> (yW-yV)/pi*180 # 度
[1] 65.67333
> (yW-yV)/(yH-yG)
[1] 2.16396

球面では弧GH、VWが π/6（30°）になるように正方形を描きました。メルカトル図にすると、線分GHは 30.35°ほど、線分VWは 65.67° ほどに描かれ、線分VWの長さはGHの 2.16 倍でした。

メルカトル図に定規を当てて、VとWのy座標の値 yV、yW を知ったとします。

> c(yV, yW)
[1] 0.8813736 2.0275894
> φV <- 2*atan(exp(yV))-pi/2
> φW <- 2*atan(exp(yW))-pi/2
> c(φV, φW)
[1] 0.7853982 1.3089969
> (φW - φV)
[1] 0.5235988
> (φW - φV)/pi*180
[1] 30

逆変換をすれば、φV、φW が得られ、2点間は π/6（30°）であったことが分かります。

もともと20世紀後半まで数表を使って三角関数などを計算していたわけで、数表を用意すれば計算が困難なものではないと思います。

しかし、あまり便利そうには見えません。

5. 「直線を引いて経線となす角度」を読み取る図

方眼に描くことに比べて、メルカトル図の利点は「航程線が直線となる」ことのようです。

図は、経線を平行に描くので、図に引いた直線が垂直な場合を除いて、どの経線とも同じ角度で交差します。

この角度を測って「航程線の方位」とすると言うことは、図の縦横の比は決まっていなければなりません。

緯度は南緯を負で表して、－90° から 90° までを弧度－π/2 から π/2 で表した φ から以下のように計算します。

$y=log\left({tan\left({\frac{\pi{}}{4}+\frac{\varphi{}}{2}}\right)}\right)$

これをXY直交座標の y の値としてプロットします。

この計算は、 φ=±90° のとき、±90° の正接を計算することになってしまいます。±89° なら値が求められますが、±4.741 (±271.659°) の位置に描かれることになります。

φ=±85.05112877° （±1.4844222297451）のとき、概ね y=±π となるので、これをY軸の範囲としてみます。

緯度 φ に対して、目盛りはー85.05112877° から 85.05112877° を描き、y の値としては－π から πに相当します。

経度は東経と西経に分かれていますが、計算に都合が良い連続した値にします。東経を 0° から 180° とし、西経は360°から引いた値として180°から360° と考えます。西経123° は、237° です。

0°を跨ぐ位置関係は連続しませんが、計算に必要なのは、たいてい経度差なので十分です。

XY直交座標の x の値としては、０から 2π と弧度を使うことにします。

経度も緯度も2πの幅なので、縦横1:1の図に描けば図を測った角度と、計算上の傾きが一致することになります。

左の図は、計算に使用している値と目盛りが一致した図です。

右の図は、常用に近い経度、緯度の目盛りにしたものです。

上の図の右の形式で、縦横共に360mmの地図を作ったとします。

1mm は 1° となるので、水平、垂直方向に地図を測った長さ s mm は、 $\frac{\pi{}\cdot{}s}{180}$ が左の地図の目盛りの値となり、計算ができます。

経線となす角度については分度器で測れる他に、水平、垂直に長さを測って比を求めることもできます。

しかし、これだけでは、「経線とのなす角度を測って方位角と見なして進むと目的地に着く」ということは、まだ理解できません。

6. 方位角と航程線の方位

「方位角」は、球面において、その地点における経線の作る平面と法平面の成す角です。メルカトル図の平行に引いた経線と直線の成す角ではありません。メルカトル図の平行に引いた経線と直線の成す角は、単に bearing で azimuth とは言わないようです。

実際の船の進路は、メルカトル図の平行に引いた経線と直線の成す角が示す「方位」で、「方位角」は使われません。

ある経度での交点を考えると、球面でも、メルカトル図でも、他の経度での交差角と同じだと言うことに変わりはなく、私以外は混乱しないようです。

メルカトル図に引いた直線と、方位角の関係を考えます。

メルカトル図から読み取る方位、航程線の方位を考えてみます。これは、目盛りと値が一致し、縦横が1:1の、上の左の図「(x,y)の値」で考えます。

メルカトル図の方位 B は、

$tanB=\frac{{x}_{2}-{x}_{1}}{{y}_{2}-{y}_{1}}$

メルカトル図の上での2点間の距離は、

$S=\sqrt{{\left({{x}_{2}-{x}_{1}}\right)}^{2}+{\left({{y}_{2}-{y}_{1}}\right)}^{2}}$ (6.1)

メルカトル図は経線を平行に描くので、

${x}_{2}-{x}_{1}={\lambda{}}_{2}-{\lambda{}}_{1}$

式(4.2)から、

${y}_{2}-{y}_{1}=log\left({tan\left({\frac{\pi{}}{4}+\frac{{\varphi{}}_{2}}{2}}\right)}\right)-log\left({tan\left({\frac{\pi{}}{4}+\frac{{\varphi{}}_{1}}{2}}\right)}\right)$

対数を結合して、三角関数を展開すると、

$y_2-y_1=\log\left(\frac{\cos\left(\frac{\mathrm{\varphi2}+\mathrm{\varphi1}}{2}\right)+\sin\left(\frac{\mathrm{\varphi2}-\mathrm{\varphi1}}{2}\right)}{\cos\left(\frac{\mathrm{\varphi2}+\mathrm{\varphi1}}{2}\right)-\sin\left(\frac{\mathrm{\varphi2}-\mathrm{\varphi1}}{2}\right)}\right)$

これを、

$Y\left(\varphi_1,\varphi_2\right)=\log\left(\frac{\cos\left(\frac{\mathrm{\varphi2}+\mathrm{\varphi1}}{2}\right)+\sin\left(\frac{\mathrm{\varphi2}-\mathrm{\varphi1}}{2}\right)}{\cos\left(\frac{\mathrm{\varphi2}+\mathrm{\varphi1}}{2}\right)-\sin\left(\frac{\mathrm{\varphi2}-\mathrm{\varphi1}}{2}\right)}\right)$

と記すと、

$tanB=\frac{{\lambda{}}_{2}-{\lambda{}}_{1}}{Y\left({{\varphi{}}_{1}\ ,\ \ {\varphi{}}_{2}}\right)}$

$S=\sqrt{{\left({{\lambda{}}_{2}-{\lambda{}}_{1}}\right)}^{2}+{Y\left({{\varphi{}}_{1},{\varphi{}}_{2}}\right)}^{2}}$ (6.2)

球面では法面による切り線と、2つの経線が作る球面三角形を考えます。辺はすべて大円の弧です。半径1の球として計算すれば弧度で表した中心角と弧の長さは同じ値です。球面三角法の余弦定理から、

$cos\alpha{}=\frac{{sin\varphi{}}_{2}-{cos(s)sin\varphi{}}_{1}}{sin(s){cos\varphi{}}_{1}}$

$cos(s)=cos{\varphi{}}_{1}cos{\varphi{}}_{2}cos\left({{\lambda{}}_{2}-{\lambda{}}_{1}}\right)+sin{\varphi{}}_{1}sin{\varphi{}}_{2}$

正弦定理から、

$\frac{sin\left({{\lambda{}}_{2}-{\lambda{}}_{1}}\right)}{sin(s)}=\frac{sin\alpha{}}{cos{\varphi{}}_{2}}$

$sin(s)=\frac{cos{\varphi{}}_{2}sin\left({{\lambda{}}_{2}-{\lambda{}}_{1}}\right)}{sin\alpha{}}$ (6.3)

したがって、

$\mathrm{\tan\alpha}=\frac{\cos\mathrm{\varphi2}\cdot\sin\left(\mathrm{\lambda2}-\mathrm{\lambda1}\right)}{\cos\mathrm{\varphi_1}\sin\mathrm{\varphi_2}-\sin\mathrm{\varphi_1}\cos\mathrm{\varphi_2}\cdot\cos\left(\mathrm{\lambda2}-\mathrm{\lambda1}\right)}$

横浜、バンクーバー間を計算してみます。

> YO <- c(35.45033, 139.63422)/180*pi
> VA <- c(49.266667, -123.116667+360)/180*pi
> φ1 <- YO[1]; φ2 <- VA[1]
> λ1 <- YO[2]; λ2 <- VA[2]
> Y <- function(φ1,φ2){
+ cs<-cos((φ2+φ1)/2);sn<-sin((φ2-φ1)/2);log((cs+sn)/(cs-sn))}
> log(tan(pi/4+φ2/2))-log(tan(pi/4+φ1/2))
[1] 0.328463
> tanB <- (λ2-λ1)/Y(φ1,φ2)
> B <- atan(tanB)
> S <- sqrt((λ2-λ1)^2+Y(φ1,φ2)^2)
> c(B, S)/pi*180
[1] 79.04756 99.05335
> tanα <- (sin(λ2-λ1)*cos(φ2))/
+ (cos(φ1)*sin(φ2)-cos(λ2-λ1)*sin(φ1)*cos(φ2))
> α <- atan(tanα)
> sins <- cos(φ2)*sin(λ2-λ1)/sin(α)
> s <- asin(sins)
> α/pi*180
[1] 44.22666
> s*6378137
[1] 7584760

地球を球と見ると横浜から見たバンクーバーの方位角は44.2°ほどで、距離は7,584,760mと計算されます。

メルカトル図に描くと、横浜からバンクーバーまで引いた線分は、横浜を通る経線と北から時計回りに79°ほどの角度です。縦横とも360mmに描けば2点間は99mmの長さになります。

メルカトル図に引いた直線を、球面に描いたものが航程線だと考えてみます。メルカトル図から読み取った B は、メルカトル図の「ベアリング」ですが、球面の航程線の「方位角」でした。

7. 球面に航程線を描く

前述の図「(x,y)の値」の目盛りの値で計算します。メルカトル図の緯度は等間隔ではなく、高緯度ほど引き延ばされています。

メルカトル図を縦横360mmで描くとします。XY直交座標のxとなる経度は0から2πで表します。yとなる緯度はーπからπで表します。

図「(x,y)の値」は、左中央が原点(0,0)で、x軸は赤道です。常用の緯度経度(φ、λ)を以下のように変換して考えます。

λが西経の場合、 $\lambda\leftarrow360-\lambda$
$x=\frac{\lambda{}}{180}\pi{}$
φが南緯の場合、 $\varphi{}\leftarrow{}-\varphi{}$
$\varphi\leftarrow\varphi\frac{\pi}{180}$
$y=log\left({tan\left({\frac{\pi{}}{4}+\frac{\varphi{}}{2}}\right)}\right)$

上記5.の式が、φ=±85.05112877° （±1.4844222297451）のとき ±π となるので、x軸、y軸の比を1:1とするために、極に近い部分は除外して考えます。

メルカトル図に始点A、終点Bをマークして直線を引きます。始点Aを $A\left({{A}_{x},{A}_{y}}\right)$ と表し、 ${A}_{x}$ は、図の左端から測った長さ(mm)、 ${A}_{y}$ は、赤道(x軸)から測った長さ(mm)です。南半球は負で表します。

終点Bも同様に $B\left({{B}_{x},{B}_{y}}\right)$ と表します。図を縦横360mmで描いたので、

$A\left({{x}_{A},{y}_{A}}\right)=\left({\frac{\pi{}}{180}{A}_{x},\frac{\pi{}}{180}{A}_{y}}\right)$

$B\left({{x}_{B},{y}_{B}}\right)=\left({\frac{\pi{}}{180}{B}_{x},\frac{\pi{}}{180}{B}_{y}}\right)$

点(x,y) は、以下のように常用の緯度経度(φ、λ)に変換されます。

$x\leftarrow{}x\ mod\ 2\pi{}$
のとき、 $x\leftarrow{}2\pi{}+x$
$x\leq{}\pi{}$ のとき、東経を表し、 $\lambda=\frac{180}{\pi}x$
$x>\pi{}$ のとき、西経を表し、 $\lambda=\frac{180}{\pi}\left(2\pi-x\right)$
$\varphi{}=2\cdot{}atan\left({{e}^{y}}\right)-\frac{\pi{}}{2}$
$\varphi{}\leftarrow{}\frac{180}{\pi{}}\varphi{}$
$\varphi{}\geq{}0$ のとき、北緯を表す
$\varphi{}<0$ のとき、南緯を表す

常用の緯度経度(φ、λ) が既知であるか、図から読み取れれば、XYZ座標に変換して、球面に描くことができます。

8. 航程線の長さ

航程線はXY直交座標で計算したものを3Dに描けば良いことが分かりましたが、航程線の球面での長さを考えてみます。

図の赤い線は2点間の距離を表す大円の弧なので球面三角法で計算できました。

青い航程線は大円の弧ではありません。細かく分けて大円の弧に見立てて、航程線を大円の弧を連ねたものと考えてみます。

式(6.3) から、

$sin(ds)=\frac{cos\varphi{}\cdot{}sin\left({d\lambda{}}\right)}{sinB}$

> YO <- c(35.45033, 139.63422)/180*pi
> VA <- c(49.266667, -123.116667+360)/180*pi
> φ1 <- YO[1]; φ2 <- VA[1]
> λ1 <- YO[2]; λ2 <- VA[2]
> Y <- function(φ1,φ2){
+ cs<-cos((φ2+φ1)/2);sn<-sin((φ2-φ1)/2)
+ log((cs+sn)/(cs-sn))}
> B <- atan2((λ2-λ1),Y(φ1,φ2))
> φ <- seq(φ1,φ2,length.out=600000)
> y <- log(tan(pi/4+φ/2))
> x <- (y-y[1])/tan(pi/2-B)
> dλ <- x[2:length(φ)]-x[1:(length(φ)-1)]
> (s <- sum(asin(cos(φ[2:length(φ)])*sin(dλ)/sin(B))))
[1] 1.269201
> s*6378137
[1] 8095135
> library(geosphere)
> YO <- c(35.45033, 139.63422)
> VA <- c(49.266667, -123.116667)
> distRhumb(c(YO[2],YO[1]),c(VA[2],VA[1]),r=1)
[1] 1.269201
> distRhumb(c(YO[2],YO[1]),c(VA[2],VA[1]),r=6378137)
[1] 8095136

60万分割して、geosphereパッケージの結果に近くなりました。

上の図の下方を ds を斜辺とする直角三角形と見なしてみます。

緯度の変化 dφは、

$d\varphi{}=cosB\cdot{}ds$

と表せます。したがって、

$ds=\frac{d\varphi{}}{cosB}$

表し方は良く分かりませんが、dsの和が長さsであり、dφの和は ${\varphi{}}_{2}-{\varphi{}}_{1}$ に他なりません。

$s=\sum_{\varphi=\varphi_1}^{\varphi_2}ds_{\varphi}=\frac{1}{\cos B}\sum_{\varphi=\varphi_1}^{\varphi_2}d\varphi=\frac{\varphi_2-\varphi_1}{\cos B}$

> YO <- c(35.45033, 139.63422)/180*pi
> VA <- c(49.266667, -123.116667+360)/180*pi
> φ1 <- YO[1]; φ2 <- VA[1]
> λ1 <- YO[2]; λ2 <- VA[2]
> Y <- function(φ1,φ2){
+ cs<-cos((φ2+φ1)/2);sn<-sin((φ2-φ1)/2)
+ log((cs+sn)/(cs-sn))}
> B <- atan2((λ2-λ1),Y(φ1,φ2))
> (s <- (φ2-φ1)/cos(B))
[1] 1.269201
> s*6378137
[1] 8095136

9. メルカトル図の直線が球面で等角ということ

航程線の球面での長さは、

$s=\frac{{\varphi{}}_{2}-{\varphi{}}_{1}}{cosB}$

でした。

これは航程線が方位Bと緯度だけで決まることを示しています。

どの経線とも方位Bで交差するので、赤道との交点でも、交点を通る経線と方位Bで交差します。球ですから、赤道との交点がどんな経度でも航程線の形状には変わりがありません。航程線の形状は方位Bだけで決まります。

また、球面に等角螺旋を1つ引ければ、平面図で直線となるような南北の位置を決めることができたと考えられます。

球面で赤道からの長さを測って平面に移し、直線との差から数表や専用の方眼紙を作ることができたと思います。

「ろくろ」は古くから存在し、筆の直線的な動きで螺旋を描いていました。

10. 常に同じだけ曲がる

Wikipediaのメルカトル図法の説明には「出発地と目的地との間に直線を引いて経線となす角度（「舵角」と言う）」という記述があります。他の要素を無視すれば、舵を切ったまま固定して進めば十分広い水面では円を描くことになりそうです。

これは、「進路に対して等角」や「常に同じ曲率」だと思います。基本的には船速によらず、進んだ距離とその間に変わった進路の角度で決まります。

メルカトル図法の説明にある「経線に対して等角である」との違いを考えてみます。

進む方向に対して進んだ距離に応じた一定の角度で曲がることは円を描くことを示しています。

座標(1,0,0)を出発して、1° 進むごとに1°曲がると、左図の黒い線の円を描きます。

左図は、1°進むごとに、4°、2°、1°、1/2°、1/4° 曲った場合の図です。

この話しでは、地球の自転軸や極、経線といったことは関係がありません。球面も、平面と同じように考えられると言うことです。

メルカトル図に直線を引いて、その経線と成す角を保持するというのは舵を固定することではなく、操舵を伴うことになります。

下段の図の①は、地図から読み取った方位に経度90°分進んだ時の経路と到達点です。まっすぐ進む訳ですが、始点と終点を通る大円のこの上を通ると言うことです。

②は、経度30°ごとに舵を切り直しながら経度90°分進んだ時の経路と到達点です。

③は、経度0.9°ごとに舵を切り直しながら経度450°分進んだ時の経路と到達点です。

航程線は、この操作の極限だと思います。

経線に対して等角に進むと極の周りを回ることになり、螺旋を描いています。

上段の図の「円」は、真円ではありません。1°進むごとに決まった角度で曲がっていますが、1°の間は、その間の始点と終点を通る大円の弧を描いています。作図は小さな大円の弧を連ねて描かれています。

常時舵を調整しないと通過点や到達点が変わってしまうので目安だったのだろうと思います。

11. 等角螺旋を描く

等角螺旋の線上のすべての点で、その点を通る経線と螺旋が同じ角度で交差している訳ですが、螺旋は大円ではありません。したがって角度が何を意味するのかは明瞭ではありません。

球面に等角螺旋が描かれているとして、その方位を測る術を知りません。

可能なのは方位Bが指定されている場合に、経線の作る面と、法線を軸に角度がBで交差する平面を描くことです。

2つの平面と、接平面の交線が方位を示します。

これは角度がBで交差する大円を描くことです。

短い区間、大円上を進むと考えると、その長さで到達点が変わってしまいます。

func <- function(p,α,dλ,col="black",lwd=1){
P0 <- Plane(c(0,0,0),c(0,0,1),p) # pを通る経線の面
L <- Line2(c(0,0,0),c(0,0,1))
p1 <- Rotate(L,dλ,list(x=p[1],y=p[2],z=p[3]))
P1 <- Plane(c(0,0,0),c(0,0,1),c(p1$x,p1$y,p1$z)) # 終点を通る経線の面
# 北極の座標を、pの法線の周りにα回転する
L <- Line2(c(0,0,0),c(p[1],p[2],p[3]))
p1 <- Rotate(L,-α ,list(x=0,y=0,z=1))
P2 <- Plane(c(0,0,0),p,c(p1$x,p1$y,p1$z))
L <- Intersection2Planes(P2, P1)
p1 <- IntersectionLineEllipse(L,0)
DrawEllipseArc3D2(P2,0,list(x=p[1],y=p[2],z=p[3]),
list(x=p1$x[1],y=p1$y[1],z=p1$z[1]),col=col,lwd=lwd)
invisible(list(x=p1$x[1],y=p1$y[1],z=p1$z[1]))
}

ここで等角螺旋を描いているのは航程線とは無関係に、以下の手順で大円の弧を連ねています。

P0に点pを通る経線が作る平面の式の係数を求める。

経線の作る平面は原点と北極を必ず含み、平面は3点で決まる。

点pをz軸の周りにdλ回転する。この点を通る経線の作る面の式の係数をP1に求める。

P1は点pから出発した次の到達点を通る経線の作る面を示す。

北極の座標を点pの法線の周りに方位α回転してp1とする。

原点、p、p1の3点で決まる平面をP2とする。

この平面は点pを通る経線と点pで方位αで交差する。

平面P1とP2の交線を求める。
交線と球面の交点を求める。
点pから、この交点まで、平面P2上に楕円を描く。
この関数は、弧を描いた終点を返すので、それを始点として再度呼び出せば続きを描く。

12. 航程線に関する計算

入力は弧度とし、南緯を負で表します。西経は2πから減じた値とします。

緯度には、φ=±85.05112877° （±1.4844222297451）の限界がありますが特別な扱いをしません。

以下の地点を計算してみます。

YO <- c(35.45033, 139.63422)
VA <- c(49.266667, -123.116667)
φ1 <- YO[1]/180*pi; φ2 <- VA[1]/180*pi
λ1 <- YO[2]/180*pi; λ2 <- (360+VA[2])/180*pi

航程線を平面で直線に描くことから、XY直交座標で考えます。

(φ , λ) と (x , y) は、

$y=log\left({tan\left({\frac{\pi{}}{4}+\frac{\varphi{}}{2}}\right)}\right)$

$x=\lambda{}$

$\varphi{}=2\cdot{}atan\left({{e}^{y}}\right)-\frac{\pi{}}{2}$

また、φ2-φ1 に相当する y の間隔は、

$Y\left(\varphi_1,\varphi_2\right)=\log\left(\tan\left(\frac{\pi}{4}+\frac{\varphi_2}{2}\right)\right)-\log\left(\tan\left(\frac{\pi}{4}+\frac{\varphi_1}{2}\right)\right)$

$Y\left(\varphi_1,\varphi_2\right)=\log\left(\frac{\tan\left(\frac{\pi}{4}+\frac{\varphi_2}{2}\right)}{\tan\left(\frac{\pi}{4}+\frac{\varphi_2}{1}\right)}\right)=\log\left(\frac{\cos\left(\frac{\mathrm{\varphi_2}+\mathrm{\varphi_1}}{2}\right)+\sin\left(\frac{\mathrm{\varphi_2}-\mathrm{\varphi_1}}{2}\right)}{\cos\left(\frac{\mathrm{\varphi_2}+\mathrm{\varphi_1}}{2}\right)-\sin\left(\frac{\mathrm{\varphi_2}-\mathrm{\varphi_1}}{2}\right)}\right)$

.M <- function(φ){ log(tan(pi/4+φ/2)) } # φ->y
.W <- function(y){ 2*atan(exp(y))-pi/2 } # y->φ
# y2-y1(φ1,φ2から)
.Y <- function(φ1,φ2){cs<-cos((φ2+φ1)/2);sn<-sin((φ2-φ1)/2);log((cs+sn)/(cs-sn))}

12.1. 2点間の航程線の方位

メルカトル図法の方位（ベアリング）は、

$B=\frac{\pi{}}{2}-atan(\frac{{y}_{2}-{y}_{1}}{{x}_{2}-{x}_{1}})=atan(\frac{{x}_{2}-{x}_{1}}{{y}_{2}-{y}_{1}})$

p1(φ1,λ1) から p2(φ2,λ2) へ向かうメルカトル図の方位Bは、以下のように計算します。

.B <- function(p1,p2){
φ1 <- p1[1]; λ1 <- p1[2]; φ2 <- p2[1]; λ2 <- p2[2]
B <- atan2((λ2-λ1),.Y(φ1,φ2))
ifelse(B<0,2*pi+B,B)
}

例：

> YO <- c(35.45033, 139.63422)
> VA <- c(49.266667, -123.116667)
> φ1 <- YO[1]/180*pi; φ2 <- VA[1]/180*pi
> λ1 <- YO[2]/180*pi; λ2 <- (360+VA[2])/180*pi
> .B(c(φ1,λ1),c(φ2,λ2))/pi*180
[1] 79.04756

12.2. 2点間の航程線の長さ

前述のごとく

$s=\frac{{\varphi{}}_{2}-{\varphi{}}_{1}}{cosB}$

RhumbLineLength <- function(p1,p2){
φ1 <- p1[1]; λ1 <- p1[2]; φ2 <- p2[1]; λ2 <- p2[2]
B <- .B(p1,p2)
s <- (φ2-φ1)/cos(B)
list(s=s,B=B)
}

他に、φとyの比を使う計算がありました。

> RhumbLineLength2 <- function(p1,p2){
+ φ1 <- p1[1]; λ1<-p1[2]; φ2 <- p2[1]; λ2 <- p2[2]
+ sqrt(((φ2-φ1)/.Y(φ1,φ2)*(λ2-λ1))^2+(φ2-φ1)^2)
+ }
> RhumbLineLength2(c(φ1,λ1),c(φ2,λ2))
[1] 1.269201

$s=\sqrt{{\left({\frac{{\varphi{}}_{2}-{\varphi{}}_{1}}{{y}_{2}-{y}_{1}}\left({{\lambda{}}_{2}-{\lambda{}}_{1}}\right)}\right)}^{2}+{\left({{\varphi{}}_{2}-{\varphi{}}_{1}}\right)}^{2}}$

12.3. 長さs進んだ時の到達点

球面の航程線を考えて p1(φ1,λ1) から方位Bへ s 進んだとします。このとき、メルカトル図法の図では p1(φ1,λ1) から方位Bへ直線上をS進んだことになるとします。

球面で航程線の長さは、

$s=\frac{{\varphi{}}_{2}-{\varphi{}}_{1}}{cosB}$

なので、

${\varphi{}}_{2}={\varphi{}}_{1}+s\cdot{}cosB$ (12.3.1)

また、メルカトル図では、到達点の (x,y) は、

$\left({{x}_{2}\ ,\ {y}_{2}}\right)=\left({{\lambda{}}_{1}+S\cdot{}sinB\ ,\ M\left({{\varphi{}}_{1}}\right)+S\cdot{}cosB}\right)$

ですが、 ${y}_{2}=M\left({{\varphi{}}_{2}}\right)=M\left({{\varphi{}}_{1}}\right)+S\cdot{}cosB$ ですから、

$S\cdot{}cosB=M\left({{\varphi{}}_{2}}\right)-M\left({{\varphi{}}_{1}}\right)$

$S=\frac{M\left({{\varphi{}}_{2}}\right)-M\left({{\varphi{}}_{1}}\right)}{cosB}=\frac{Y\left({{\varphi{}}_{2}-{\varphi{}}_{1}}\right)}{cosB}$

x2 は、

${x}_{2}={\lambda{}}_{1}+S\cdot{}sinB={\lambda{}}_{1}+tanB\cdot{}Y\left({{\varphi{}}_{1},\ {\varphi{}}_{2}}\right)$

したがって、

${\varphi{}}_{2}={\varphi{}}_{1}+s\cdot{}cosB$

${\lambda{}}_{2}={\lambda{}}_{1}+tanB\cdot{}Y\left({{\varphi{}}_{2}-{\varphi{}}_{1}}\right)$

$\left({{x}_{2}\ ,\ {y}_{2}}\right)=\left({{\lambda{}}_{2},\ M\left({{\varphi{}}_{2}}\right)}\right)$

となります。

RhumbLineDestination <- function(p1,B,s){
φ1 <- p1[1]; λ1 <- p1[2];
φ2 <- φ1 +s*cos(B)
λ2 <- λ1 + tan(B)*.Y(φ1,φ2)
list(φ=φ2, λ=λ2)
}

例：

> YO <- c(35.45033, 139.63422)
> VA <- c(49.266667, -123.116667)
> φ1 <- YO[1]/180*pi; φ2 <- VA[1]/180*pi
> λ1 <- YO[2]/180*pi; λ2 <- (360+VA[2])/180*pi
> g <- RhumbLineLength(c(φ1,λ1),c(φ2,λ2))
> str(u)
List of 2
$ φ: num 0.86
$ λ: num 4.13
> u <- RhumbLineDestination(c(φ1,λ1),g$B,g$s)
> u$φ/pi*180
[1] 49.26667
> u$λ/pi*180 - 360
[1] -123.1167
> ds <- seq(0,g$s,length.out=5)
> v <- RhumbLineDestination(c(φ1,λ1),g$B,ds)
> str(v)
List of 2
$ φ: num [1:5] 0.619 0.679 0.739 0.8 0.86
$ λ: num [1:5] 2.44 2.83 3.24 3.67 4.13

12.4. 航程線の点列

球面に描いたり、方眼やメルカトル図に描くので、緯度、経度での座標列を作ります。

これは、上の例のように、到達点を複数与えれば実現できます。

2点間に描くのであれば、先に航程線の長さを求めます。

12.5. 赤道と航程線の交点

ある点p1(φ1,λ1) を方位Bで通る航程線と、赤道の交点を求めます。

点p1 は、XY直交座標で (λ1, M(φ1))

$\left({{x}_{1}\ ,\ {y}_{1}}\right)=\left({{\lambda{}}_{1}\ ,M\left({\ {\varphi{}}_{1}}\right)}\right)$

です。航程線は、XY平面で直線なので、その傾き a は、

$a=tan(\frac{\pi{}}{2}-B)=\frac{1}{tanB}$

です。切片bは、

$b=M\left({{\varphi{}}_{1}}\right)+\frac{{\lambda{}}_{1}}{tanB}$

したがって、y=0 のとき、

${x}_{0}={\lambda{}}_{1}-tanB\cdot{}M\left({\ {\varphi{}}_{1}}\right)$

RhumbLineEquatorPoint <- function(p1,B){
φ1 <- p1[1]; λ1 <- p1[2];
λ0 <- λ1 - tan(B)*.M(φ1)
list(φ=0, λ=λ0)
}

YO <- c(35.45033, 139.63422)
VA <- c(49.266667, -123.116667)
φ1 <- YO[1]/180*pi; φ2 <- VA[1]/180*pi
λ1 <- YO[2]/180*pi; λ2 <- (360+VA[2])/180*pi
d <- map("world2Hires",plot=FALSE)
x <- d$x /180*pi
y <- .M(d$y/180*pi)
par(mar=c(3,3,1,1),mgp=c(1.5, 0.5, 0))
plot(x,y,t="l",col="blue",ylim=c(-pi,pi))
abline(v=0:6,h=-3:3,col="gray80",lty=2)
points(c(λ1,λ2),.M(c(φ1,φ2)),pch=20,col="red")
# 大円
yo <- φλ2XYZ(φ1,λ1,EARTH$f)
va <- φλ2XYZ(φ2,λ2,EARTH$f)
P <- Plane(c(0,0,0),c(yo$x,yo$y,yo$z),c(va$x,va$y,va$z))
u <- DrawEllipse3D(P,EARTH$f, draw=FALSE)
φλ <- XYZ2φλ(u$x,u$y,u$z,EARTH$f)
p <- cbind(φλ$λ,φλ$φ)
q <- p[order(p[,1]),]
lines(q[,1],.M(q[,2]),col="red",lwd=1)
# 航程線
u <- RhumbLineLength(c(φ1,λ1),c(φ2,λ2))
ds <- seq(0,u$s,length.out=100)
p <- RhumbLineDestination(c(φ1,λ1),u$B,ds)
lines(p$λ,.M(p$φ),col="navy",lwd=5)
# 長い航程線
q <- RhumbLineEquatorPoint(c(φ1,λ1),u$B)
ds <- seq(0,2*pi,length.out=361)
p <- RhumbLineDestination(c(0,q$λ),u$B,ds)
lines(p$λ,.M(p$φ),col="green",lwd=2)
title("長い航程線",cex.main=0.8)

13. 等角螺旋の図

航程線の特徴には、螺旋を描き、極から放射状に開く経線と、どこでも等角で交差することが上げられています。

メルカトル図には縦横比を維持するために緯度の限界がありますが、球面に描くには90°に近い値でも良いはずです。ただし、±90° の正弦は ±∞ ですから計算不能です。始点を緯度89.9°、経度ゼロとして描きます。

経線との交差角度の方位Bは、北から時計回りに測るので、90°に近い方が、傾きが小さくなり、長い線を描きます。上の図は、メルカトル図で直線の傾きが -5° となる、B=95° で描いています。

長さは11.4πとして見ました。

> EARTH <- DrawEarth(f=0) # 地球の描画
> d <- map("world",plot=F) # 行政区画(x:経度,y:緯度)
> d <- degXYZ(d$y,d$x,EARTH$b) # XYZ座標系に変換
> lines3d(d[,1],d[,2],d[,3],col="gray40",add=T)
> ds <- seq(0,11.4*pi,length.out=1000)
> B <- 95 /180*pi
> p <- RhumbLineDestination(c(89.9/180*pi,0),B,ds)
> d <- φλ2XYZ(p$φ,p$λ,EARTH$f) # XYZ座標系に変換
> lines3d(d$x,d$y,d$z,col="blue",lwd=2)

北極付近から5°ずつ南極に向かって下って行きます。

線には太さがないと言うことでは極に到達しませんが、船は大きく、短距離で曲がることもできません。

14. 地球楕円体の航程線

現在では逐次自身の周りの地図を作り出すことができます。メルカトルの時代、大航海時代がどのように捉えられているのかを知るには航程線を知る必要があるようです。その意味では地球は球であって大きさも定まっていませんでした。

しかし、現在では、地球を準拠楕円体と見なして数値計算するので、その説明も多くあります。地球の表面は絶えず動いていて、準拠楕円体はそれを捉えるために必要です。地表の一点の位置を正確に示すことは出来ない訳ですが、準拠楕円体の上の計算であれば正解があります。

準拠楕円体は楕円体ですが、球が地軸の周りに膨らんだ回転楕円体です。したがって、極方向から見れば真円です。緯線は真円を描きます。経線は地軸を含む平面による截り線で、すべて同じ楕円を成し、赤道半径が長軸半径、極半径が短軸です。

上の図を「4. 横に縦を合わせる」と同じように考えてみます。

更成緯度βで、楕円体上の点Qの位置を表します。赤は赤道、青は点Qを通る緯線、緑は補助球の点Qに対応する点Pを通る緯線です。赤道の半径を1とすると、緯線の作る円はいずれも cosβ なので、メルカトル図の平行な経線では、球と同様に、1/cosβ の長さに描かれることになります。

経線方向も同じように伸ばすので、球では、式(4.3) のように積分しました。球面では中心角の1°は、どの2点間でも同じ距離ですが、回転楕円体の経線の弧長は緯度によって変化します。

点P、Qの近傍の弧を考えると、それぞれ、図のような小さな直角三角形と見なせると言うのが説明のようです。

直角三角形の斜辺は、角度を小さくしていけば、それぞれの接線に近づきます。

メルカトル図に描く赤道から点Pまでの長さ y は、経線方向の変化と、水平方向の変化の積で表され、a=1として、式(4.3) の計算になります。

$y=\int\limits_{0}^{\beta{}}a\cdot{}\frac{1}{cos\beta{}}d\beta{}=\int\limits_{0}^{\beta{}}\frac{1}{cos\beta{}}d\beta{}$

メルカトル図に描く赤道から点Qまでの長さ y は、以下のように考えられます。

$y=\int\limits_{0}^{\beta{}}\sqrt{{a}^{2}{sin}^{2}\beta{}+{b}^{2}{cos}^{2}\beta{}}\cdot{}\frac{1}{cos\beta{}}d\beta{}$

a=1 と考えますが、b/a かどうかが紛らわしいので、離心率に置き換えます。離心率 e は、

${e}^{2}=\frac{{a}^{2}-{b}^{2}}{{a}^{2}}$

なので、

${b}^{2}=1-{e}^{2}$

${cos}^{2}\beta{}=1-{sin}^{2}\beta{}$

を適用して、

$y=\int\limits_{0}^{\beta{}}\sqrt{1-{e}^{2}{cos}^{2}\beta{}}\cdot{}\frac{1}{cos\beta{}}d\beta{}$ (14.1)

です。

次に β を φ にします。

長半径が1、短半径がbの楕円の点 $\left({{x}_{0}\ ,\ {y}_{0}}\right)$ における接線は、 ${x}_{0}x+\frac{{y}_{0}y}{{b}^{2}}=1$ なので、 $y=\frac{-{b}^{2}{x}_{0}}{{{y}_{0}}_{\ }}x+\frac{{b}^{2}}{{{y}_{0}}_{\ }}$ として、傾きが $\frac{-{b}^{2}{x}_{0}}{{{y}_{0}}_{\ }}$ であることが分かります。したがって、点 $\left({{x}_{0}\ ,\ {y}_{0}}\right)$ における法線の傾きは、その逆数の $\frac{{y}_{0}}{{b}^{2}{x}_{0}}$ です。

法線の傾きは地理的緯度を表し、点Q $\left({{x}_{0}\ ,\ {y}_{0}}\right)=\left({cos\beta{},\ b\cdot{}sin\beta{}}\right)$ における法線の傾きは、

$tan\varphi{}=\frac{{y}_{0}}{{b}^{2}{x}_{0}}=\frac{b\cdot{}sin\beta{}}{{b}^{2}cos\beta{}}=\frac{1}{b}tan\beta{}$ (14.2)

2乗した ${tan}^{2}\varphi{}=\frac{1}{{b}^{2}}{tan}^{2}\beta{}$ に、三角関数の公式から $tan\theta{}=\frac{sin\theta{}}{cos\theta{}}$ 、 ${sin}^{2}\theta{}+{cos}^{2}\theta{}=1$ 、 ${b}^{2}=1-{e}^{2}$ を適用して、

${cos}^{2}\beta{}=\frac{{cos}^{2}\varphi{}}{{1-e}^{2}{sin}^{2}\varphi{}}$ (14.3)

水平方向の拡大は、

$\frac{1}{cos\beta{}}=\frac{\sqrt{{1-e}^{2}{sin}^{2}\varphi{}}}{cos\varphi{}}$ (14.4)

式(14.1) の積分対象の式は、経線方向の変化に水平方向の変化を乗じたものなので、

$\sqrt{1-e^2\cos^2\beta}\cdot\frac{1}{\cos\beta} =\sqrt{1-e^2\frac{\cos^2\varphi}{1-e^2\sin^2\varphi}} \cdot\frac{\sqrt{1-e^2\sin^2\varphi}}{\cos\varphi} =\frac{\sqrt{1-e^2}}{\cos\varphi}$

式(14.2) の短半径bを離心率で表して、

$tan\varphi{}=\frac{1}{\sqrt{1-{e}^{2}}}tan\beta{}$ (14.5)

$tan\varphi{}\ ,\ tan\varphi{}$ を微分して、

$\frac{1}{{cos}^{2}\varphi{}}d\varphi{}=\frac{1}{\sqrt{1-{e}^{2}}}\cdot{}\frac{1}{{cos}^{2}\beta{}}d\beta{}$

$d\beta=\cos^2\beta\cdot\frac{\sqrt{1-e^2}}{\cos^2\varphi}d\varphi=\frac{\cos^2\varphi}{1-e^2\sin^2\varphi}\cdot\frac{\sqrt{1-e^2}}{\cos^2\varphi}=\frac{\sqrt{1-e^2}}{1-e^2\sin^2\varphi}$

$y=\int_0^{\beta}\sqrt{1-e^2\cos^2\beta}\cdot\frac{1}{\cos\beta}d\beta=\int_0^{\varphi}\frac{\sqrt{1-e^2}}{\cos\varphi}\cdot\frac{\sqrt{1-e^2}}{1-e^2\cos^2\varphi}d\varphi$

$y=\int\limits_{0}^{\varphi{}}\frac{1-{e}^{2}}{cos\varphi{}\left({{1-e}^{2}{sin}^{2}\varphi{}}\right)}d\varphi{}$ (14.6)

Wikipediaには次の式がありました。

$y={tanh}^{-1}sin\varphi{}-e\cdot{}{tanh}^{-1}\left({e\cdot{}sin\varphi{}}\right)$

$y=\frac{1}{2}log\left({\frac{1+sin\varphi{}}{1-sin\varphi{}}}\right)-e\cdot{}\frac{1}{2}\cdot{}log\left({\frac{1+e\cdot{}sin\varphi{}}{1-e\cdot{}sin\varphi{}}}\right)$

また、以下の式も使用されるようです。

$y=log\left({tan\left({\frac{\pi{}}{4}+\frac{\varphi{}}{2}}\right)}\right)-e\cdot{}\frac{1}{2}\cdot{}log\left({\frac{1+e\cdot{}sin\varphi{}}{1-e\cdot{}sin\varphi{}}}\right)$ (14.7)

この式の右辺の最初の項は球のときの計算と同じです。球の場合、

$y=log\left({tan\left({\frac{\pi{}}{4}+\frac{\varphi{}}{2}}\right)}\right)={tanh}^{-1}sin\varphi{}=\frac{1}{2}\cdot{}log\left({\frac{1+sin\varphi{}}{1-sin\varphi{}}}\right)$

と計算できるようです。また、式(14.7) の第2項は球の場合との差を表し、差だけを考えることもできるもののようです。

$e\cdot{}{tanh}^{-1}\left({e\cdot{}sin\varphi{}}\right)=e\cdot{}\frac{1}{2}\cdot{}log\left({\frac{1+e\cdot{}sin\varphi{}}{1-e\cdot{}sin\varphi{}}}\right)$

> # 扁平率から、短半径、離心率
> f <- 0.4; b <- 1-f; e2 <- f*(2-f); e <- sqrt(e2)
> c(f,b,e2,e)
[1] 0.40 0.60 0.64 0.80
> # 回転楕円体の同じ場所を示す地理的緯度と更成緯度
> φ <- pi/3; β <- atan((b*tan(φ)))
> c(φ,β) /pi*180
[1] 60.00000 46.10211
> # 等長緯度の計算3式
> y1 <- atanh(sin(φ)) - e*atanh(e*sin(φ))
> y2 <- (1/2)*log((1+sin(φ))/(1-sin(φ))) -
+ e*(1/2)*log((1+e*sin(φ))/(1-e*sin(φ)))
> y3 <- log(tan(pi/4+φ/2)) - e*(1/2)*log((1+e*sin(φ))/(1-e*sin(φ)))
> c(y1,y2,y3)
[1] 0.6342705 0.6342705 0.6342705

式(14.1)、(14.6) の積分をして見ます。

> # 扁平率から、短半径、離心率
> f <- 0.4; b <- 1-f; e2 <- f*(2-f); e <- sqrt(e2)
> c(f,b,e2,e)
[1] 0.40 0.60 0.64 0.80
> # 回転楕円体の同じ場所を示す地理的緯度と更成緯度
> φ <- pi/3; β <- atan((b*tan(φ)))
> c(φ,β) /pi*180
[1] 60.00000 46.10211
> # 積分
> integrate(function(β){sqrt(1-e^2*cos(β)^2)/cos(β)},0,β)$value
[1] 0.6342705
> integrate(function(φ){(1-e^2)/(cos(φ)*(1-e^2*sin(φ)^2))},0,φ)$value
[1] 0.6342705

GRS80準拠楕円体の扁平率 1/298.257222101のとき、y=πとなるφは、およそ 85.08406° になります。直線を引いたり、角度を読み取るためには、図の縦横比を維持する必要があるので、縦横同じ長さに描くなら、この緯度が上端になります。

「12.航程線に関する計算」の、球の場合の関数 .M、.W、.Y、.B を回転楕円体の用に書き換えます。

.e_ <- function(f=1/298.257222101){sqrt(f*(2-f))}
.M_ <- function(φ,f=1/298.257222101){ # φ->y
e<-.e_(f); log(tan(pi/4+φ/2)) -0.5*e*log((1+e*sin(φ))/(1-e*sin(φ)))}
.W_ <- function(y,f=1/298.257222101){ # y->φ
e<-.e_(f);
m1 = exp(-y)
old <- pi/2 - 2.0 * atan(m1)
while(TRUE){
sn <- e*sin(old)
m2 <- ((1 - sn) / (1 + sn))^(e/2)
φ <- pi/2 - 2 * atan(m1 * m2)
if(abs(old - φ) <= 1e-12) break
old = φ
}
φ
}
.Y_ <- function(φ1,φ2,f=1/298.257222101){ # y2-y1(φ1,φ2から)
e <- .e_(f)
u <- rep(NA,length(φ2))
for(i in 1:length(φ2))
u[i] <- integrate(function(φ){(1-e^2)/
(cos(φ)*(1-e^2*sin(φ)^2))},φ1,φ2[i])$value
u
}
.B_ <- function(p1,p2,f=1/298.257222101){
φ1 <- p1[1]; λ1 <- p1[2]; φ2 <- p2[1]; λ2 <- p2[2]
B <- atan2((λ2-λ1),.Y_(φ1,φ2,f=f))
ifelse(B<0,2*pi+B,B)
}

後は球と同様に、回転楕円体にも、メルカトル図や等角螺旋が描けます。

f <- 0.4; b <- 1-f
clear3d(); rgl.light();
vectors3d(c(1.3,0,0),origin=c(-1.3,0,0),col="red"); text3d(1.4,0,0,"X",col="red")
vectors3d(c(0,1.3,0),origin=c(0,-1.3,0),col="green"); text3d(0,1.4,0,"Y",col="green")
vectors3d(c(0,0,1.3),origin=c(0,0,-1.3),col="blue"); text3d(0,0,1.4,"Z",col="blue")
ee <- ellipsoid(c=b,col="white",alpha=0.7) # 回転楕円体
shade3d(ee)
Circle3D(c(0,0,0), c(1,0,0), c(0,1,0),col="red",lwd=1) # 赤道
# 経線
for(θ in seq(0,pi-pi/6,by=pi/6)){
P <- Plane(c(0,0,0), c(cos(θ),sin(θ),0), c(0,0,1))
DrawEllipse3D(P,f=f,col="orange",lwd=1)
}
# 航程線
p1 <- c(-85/180*pi,0); p2 <- c(85/180*pi,360/180*pi)
u <- RhumbLineLength_(p1,p2,f=f)
ds <- seq(0,11.4*pi,length.out=1000)
B <- 95 /180*pi
p <- RhumbLineDestination_(c(89.9/180*pi,0),B,ds,f=f)
d <- φλ2XYZ(p$φ,p$λ,f=f) # XYZ座標系に変換
lines3d(d$x,d$y,d$z,col="blue",lwd=2)