🛰️

宇宙からの地球観測１章

2024/04/11に公開

宇宙からの地球観測　という本は素晴らしいです。
章末問題がついているからです。　私は宇宙も数学も一般人レベルなので数式がならんでいるだけの本だと読み進めるのに時間がかかります。
その点章末問題がついていてくれれば、まず問題を読んでわからないところを狙って読むという受験テクニックを応用できます。

さっそく第一章の問題と、それを読んでわからなかったところをまとめていきたいと思います。

問題1.1

地心緯度を用いて、地球上の任意の位置の半径（地球中心からの距離）rが次式で表現されることを証明しなさい
$r={R_b \over \sqrt{1-e^2cos^2\phi}}$
ここにeは離心率であり　 $\sqrt{R^2_a-R^2_b}/R_a$ である。 $R_a$ は赤道半径、 $R_b$ は極半径、 $\phi$ は地心緯度。

いきなり理解につまづく言葉が並んでいます。日本語なのでなんとなく想像はつきますが説明できない。
さっそく言葉の意味を一章本文から探していきます。

言葉の理解

地心緯度
- 対象物から地球中心を向く方向と赤道面のなす角 $\phi$
- 注意：地球は楕円なので地球中心の方向は重力の方向（真下）と一致しない
赤道半径
- 赤道面の半径　約6378km
- 地球は自転しているのでわずかに赤道面が張り出す回転楕円体
極半径
- 北極、南極と地球中心を結ぶ線　約6356km 赤道半径より22km 短い
離心率
- こちらは一般の数学の話なので美しい解説参考: https://manabitimes.jp/math/2352
- 点Fからの距離と直線lからの距離の比が一定 $e$ である点の軌跡は二次曲線になる。特に，
  - $0<e<1$ のとき楕円
  - $e=1$ のとき放物線
  - $1<e$ のとき双曲線
- 楕円の公式: ${x^2 \over a^2} + {y^2 \over b^2} = 1 (a>b>0)$
- 焦点の公式: $(\sqrt{a^2-b^2},0)$
- 離心率の公式: ${\sqrt{a^2-b^2} \over {a}}$

回答1.1

赤道面からの角度がわかれば赤道半径と、極半径と離心率をつかって半径を表せるということを証明する問題ですね。

まず楕円面を表す式は下記になる: https://manabitimes.jp/math/2644

{x^2 \over a^2} + {y^2 \over b^2} + {z^2 \over c^2} = 1

ここで、地球の準拠楕円体のある位置を $(x_s, y_s, z_s)$ として
地球は回転楕円体なので上記の式の $a=b=R_a$ :赤道半径、 $c=R_b$ 極半径になり、
式1.1: ${x_s^2 \over R^2_a} + {y_s^2 \over R^2_a} + {z_s^2 \over R^2_b} = 1$
XY平面に、 $x_s,y_s$ を投影して中心からの距離を $r_s$ とすると三平方定理から
式1.2: $r_s^2=x_s^2+y_s^2$
地球の中心からある地点までの距離rと地心緯度 $\phi$ を使い $r_s,z$ は以下のように表せる
式1.3: $r*\cos{\phi}=r_s$ , $r*\sin{\phi}=z_s$
式1.1 と 1.2 より
式1.4: ${r_s^2 \over R^2_a} + {z_s^2 \over R^2_b} = 1$
式1.3 と 1.4 より
式1.5: ${(r*\cos{\phi})^2 \over R^2_a} + {(r*\sin{\phi})^2 \over R^2_b} = 1$
展開すると $r^2 {\cos{\phi}^2 \over R^2_a} + r^2 {\sin{\phi}^2 \over R^2_b} = 1$
rでまとめる　 $r^2( \cos^2{\phi} R^2_b + {\sin^2{\phi} R^2_a}) = R_a^2 R_b^2$
三角関数の定理: $\cos^2\theta + \sin^2\theta = 1$ より
$r^2( \cos^2\phi R^2_b + (1-\cos^2\phi) R^2_a) = R_a^2 R_b^2$
r以外を右辺に移行し、ルートを取る　 $r = {R_a R_b \over \sqrt{\cos^2\phi R_b^2 + (1-\cos^2\phi) R_a^2}}$
右辺をcosでまとめ、 $r = {R_a R_b \over \sqrt{\cos^2\phi (R_b^2 -R_a^2) + R_a^2 }}$
分母、分子をR_a で割る式1.6:　 $r = {R_b \over \sqrt{{-\cos^2\phi (R_a^2 -R_b^2)\over R_a^2} + 1 }}$

また、離心率の定義より地球の離心率をeとすると
式1.7: $e={\sqrt{R_a^2-R_b^2} \over R_a}$
式1.6, 1.7 より
$r = {R_b \over \sqrt{1-e\cos^2\phi}}$

問題1.2

地心緯度が30度、経度が145度, 高さが0m の場合に以下を計算しなさい
ただし準拠楕円体はGRS80とする。
(1) X-Y-Z 直交座標系での数値
(2) 測地緯度の値

言葉の理解

GRS80
- 測地基準系 1980年代につくられた
  - 赤道半径　6378137 m
  - 極半径 6356752.5 m
  - 扁平率　1/298.257
測地緯度
- 対象物から重力方向(真下)に引いた線分と赤道面のなす角 $\phi$
楕円体表面の位置
- 公式が存在します：楕円面の媒介変数表示
- ${x^2 \over a^2} + {y^2 \over b^2} + {z^2 \over c^2} = 1$ は媒介変数 $\theta , \phi$ を用いて次のように表せる
  - $x=a\cos\theta\cos\phi$
  - $y=b\cos\theta\sin\phi$
  - $z=c\sin\theta$

回答1.2

問題1.1 から地心緯度30での地心半径:R を求める
上記の公式からx,y,z を求める

import math

def to_radians(degrees):
    return degrees * math.pi / 180

def calcR(Ra,Rb, lat):
   e = math.sqrt(Ra**2 - Rb**2) / Ra
   R = Rb / math.sqrt(1-e**2*(math.cos(lat)**2))
   print(f"e:{e:.6f} R:{R:.2f}")
   return R

def calcXYZ(R, lat, lng) :
  x=R * math.cos(lat)* math.cos(lng)
  y=R * math.cos(lat)* math.sin(lng)
  z=R * math.sin(lat)
  print(f"x:{x:.2f} y:{y:.2f} z:{z:.2f}")
  return x, y, z

Ra = 6378137.0       # 赤道半径 (メートル)
Rb = 6356752.5  # 極半径 (メートル)
geocentric_lat = to_radians(30) # 地心緯度
lng = to_radians(145) # 経度

R = calcR(Ra, Rb,geocentric_lat)
x,y,z = calcXYZ(R, geocentric_lat,lng)
--
e:0.08 R:6372770.65
x:-4520884.79 y:3165557.61 z:3186385.32

本文式1.6:下記の式から直行座標x,y,zから測地緯度 $\psi$ 地心緯度 $\phi$ をもとめる。
$\psi=\tan^{-1}{{R_a^2 \over R_b^2} {z_s \over r_s}}$
$\phi=\tan^{-1}{z_s \over r_s}$
（ここで $r_s$ とは、XY平面上の距離の距離 $r_s^2=x_s^2+y_s^2$ ）

def to_degree(rad):
    return rad * 180 / math.pi 

def calcGeocentricLat(zs,rs):
    phi = math.atan((zs / rs))
    print(f"GeocentricLat:{to_degree(phi):.3f}");
    return phi


def calcGeodeticLat(Ra,Rb,zs,rs):
    psi = math.atan((Ra**2 / Rb**2) * (zs / rs))
    print(f"GeodeticLat:{to_degree(psi):.3f}");
    return psi

rs = math.sqrt(x**2+y**2)
geodetic_lat = calcGeodeticLat(Ra,Rb,z,rs)
---
GeodeticLat:30.167

問題1.3

測地緯度が30度、経度が145度, 高さが0m の場合に以下を計算しなさい
ただし準拠楕円体はGRS80とする。
(1) X-Y-Z 直交座標系での数値
(2) 地心緯度の値

本文式1.6 の後半部分
$\psi=\tan^{-1}{{R_a^2 \over R_b^2} {\tan\phi}}$ から
測地緯度 $\psi$ から地心緯度 $\phi$ をもとめるには、
まず、 $\tan^{-1}$ はアークタンジェント（逆正接）を意味するので、両辺のタンジェントを取ることで式を逆にする

\tan \psi = \frac{R_a^2}{R_b^2} \tan \phi

この式から $\phi$ を単独で求めるためには、両辺を $ \frac{R_a^2}{R_b^2} $ で割ります：

\frac{\tan \psi}{\frac{R_a^2}{R_b^2}} = \tan \phi

ここで、 $\tan \phi$ を求めるには、両辺のアークタンジェント（逆正接）を取ります：

\phi = \tan^{-1} \left( \frac{\tan \psi}{\frac{R_a^2}{R_b^2}} \right)

地心緯度がわかればあとは1.2 と同じ計算でx,y,z 座標がでます。

def calcGeocentricFromGeodeticLat(Ra,Rb,lat):
  phi = math.atan(math.tan(lat)*Rb**2/(Ra**2))
  print(f"GeocentricLat:{to_degree(phi):.3f}")
  return phi

geodetic_lat = to_radians(30)
geocentric_lat = calcGeocentricFromGeodeticLat(Ra,Rb,geodetic_lat)
R = calcR(Ra, Rb,geocentric_lat)
x,y,z = calcXYZ(R, geocentric_lat,lng)
--
GeocentricLat:29.834
e:0.081819 R:6372824.47
x:-4528482.69 y:3170877.72 z:3170373.90

以上です。

途中プログラムを書いても回答と全然ちがったり、radianに変換するのをわすれたり、地心半径Rと、赤道面へ投射したRsを勘違いしていたりと、
ただしく問題をとけるプログラムを書くのは一番勉強になりました。
（ただ1.3の回答が、本書の回答と微妙に1meter ほど違う理由がわかりません、、精度の問題なのか、、　ロケットをちゃんと着陸させるのって本当に大変なんだということが図らずも実感できました。）

問題1.1

言葉の理解

回答1.1

問題1.2

言葉の理解

回答1.2

問題1.3

Discussion