高校レベルの数学で歳差周期2万6000年を導く

# はじめに

下記の記事で一太陽年や一朔望月などの天体物理の基礎用語を調べた。

https://qiita.com/TETSURO1999/items/ff4ea53d3363ae00bdc4?1

このうち地球が太陽の引力（潮汐力）によって地軸の向きが緩やかに変わること，すなわち歳差運動により北極星とされる星（現在はこぐま座のα星）が移り変わっていくことはもともと知っていた。しかし一太陽年は地球が太陽の周りを回って元の位置に戻る手前で定義されていること，すなわち一太陽年後に地球は元の位置には戻らないことを知って驚いた。地球の自転軸と太陽の向きの関係が変わらないだけで，地球の位置は少しずつ変わっていくのだ。この通りなら冬や夏の星座と呼ばれているものも2万6000年周期で入れ替わることになる・・・と思って調べてみたらその通りだった。いわゆる知り得た情報が知識として身に付いていなかったということだろう。

さて，今回はこの歳差周期2万6000年を高校レベルの数学で求めてみたい。

:::note
本記事では慣性モーメントや角運動量など剛体の物理を取り扱っており，内容自体は大学レベルだと思うが，あくまで使用している数学は高校レベルという意味である。
:::

# 潮汐力を受ける部分の質量

地球は完全球体ではなく，自転による遠心力により赤道方向に僅かに膨らんだ回転楕円体となっている。地球は太陽および月などによって潮汐力を受けるが，地球が完全な球体ではないことにより，この潮汐力が地軸の傾きを小さくする方向へのトルクとして作用する。

```plantuml
@startuml
<style>
root {
  margin: 0;
  padding: 0;
}
title {
  fontSize: 16;
  fontStyle: normal;
}
</style>
sprite FIGURE <svg width="500" height="400" viewbox="0 0 0 500 400">
<ellipse cx="250" cy="200" rx="120" ry="100" fill="pink" transform="rotate(20,250,200)" />
<ellipse cx="250" cy="200" rx="120" ry="100" fill="none" stroke="black" transform="rotate(20,250,200)" />
<circle cx="250" cy="200" r="100" fill="white" />
<circle cx="250" cy="200" r="100" fill="none" stroke="black" />
<g transform="rotate(20,250,200)">
<path d="M250 50 v10 h1 v-10" fill="black" transform="rotate(20,250,50)" />
<path d="M250 50 v10 h1 v-10" fill="black" transform="rotate(-20,250,50)" />
<path d="M250 50 v300 h1 v-300" fill="black" />
<path d="M425 200 h-10 v1 h10" fill="black" transform="rotate(20,425,200)" />
<path d="M425 200 h-10 v1 h10" fill="black" transform="rotate(-20,425,200)" />
<path d="M75 200 h350 v1 h-350" fill="black" />
</g>
<path d="M156 166 h25 v1 h-25" fill="black" transform="rotate(-20,156,166)" />
<path d="M344 234 h-28 v1 h28" fill="black" transform="rotate(50,344,234)" />
<circle cx="344" cy="234" r="5" fill="red" />
<circle cx="344" cy="234" r="5" fill="none" stroke="black" />
<circle cx="156" cy="166" r="5" fill="red" />
<circle cx="156" cy="166" r="5" fill="none" stroke="black" />
<g fill="blue">
<path d="M375 220 h15 v3 l15 -5 l-15 -5 v3 h-15" />
<path d="M125 180 h-15 v3 l-15 -5 l15 -5 v3 h15" />
</g>
<g fill="lime">
<path d="M245 55  h-15 v3 l-15 -5 l15 -5 v3 h15" transform="rotate(20,250,200)" />
<path d="M255 345 h15 v3 l15 -5 l-15 -5 v3 h-15" transform="rotate(20,250,200)" />
</g>
<text x="375" y="240" fill="black" font-size="16">潮汐力</text>
<text x="80"  y="170" fill="black" font-size="16">潮汐力</text>
<text x="215" y="55"  fill="black" font-size="16">トルク</text>
<text x="235" y="350" fill="black" font-size="16">トルク</text>
<text x="290" y="210" fill="black" font-size="16">リング</text>
<text x="180" y="160" fill="black" font-size="16">リング</text>
</svg>
title
<$FIGURE>
図１　地軸の傾きと潮汐力
end title
@enduml
```


簡単のため地球を完全球体とみなし，膨らんだ部分の余剰質量が細いリングとして赤道上に集約しているものと考える。リングの線密度 $\mu$ とおくと余剰質量の総和 $\Delta M$ は地球の半径 $r_p$ として

```math
\Delta M = 2 \pi \cdot r_p \cdot \mu \tag{1}
```

となる。一方，地球の極半径 $r_p$，赤道半径 $r_e$ とおくと膨らんだ部分の余剰体積 $\Delta V$ は

```math
\begin{align}
\Delta V &= \frac{4}{3}\pi \cdot r_p \cdot r_e^2 - \frac{4}{3}\pi \cdot r_p^3 \\
&= \frac{4}{3} \pi \cdot r_p \cdot ( r_e^2 - r_p^2)  \tag{2}
\end{align}
```

となる。膨らんだ部分は地球の中でも最外層の軽い部分であるから地殻密度 $\rho$ とおくと線密度 $\mu$ は

```math
\begin{align}
\Delta M &= \rho \cdot \Delta V \\
2 \pi \cdot r_p \cdot \mu &= \frac{4}{3} \pi \cdot \rho \cdot r_p \cdot ( r_e^2 - r_p^2) \\
\mu &= \frac{2}{3} \cdot \rho \cdot (r_e^2 - r_p^2) \tag{3}
\end{align}
```

となる。

# 地球の受けるトルク

簡単のため傾いた地軸をZ軸にとる。地球と太陽の距離，すなわち地球の公転軌道半径 $R_e$ とおき，太陽はXZ平面上の座標 $(R_e \cos \varepsilon, 0, R_e \sin \varepsilon)$ に位置するものとする。ただし，太陽は十分遠方にあるものとし，地球上において受ける潮汐力は平行とする。

```plantuml
@startuml
<style>
root {
  margin: 0;
  padding: 0;
}
title {
  fontSize: 16;
  fontStyle: normal;
}
</style>
sprite FIGURE <svg width="680" height="720" viewbox="0 0 680 720">
<path d="M170 200 v-100 A100 100 0 0 0 136 106" fill="cyan" />
<circle cx="170" cy="200" r="100" fill="none" stroke="black" />
<path d="M170 200 h100 A100 25 0 0 0 253 186" fill="yellow" />
<g transform="rotate(-30,170,200)">
<path d="M320 200 h-10 v1 h10" fill="black" transform="rotate(20,320,200)" />
<path d="M320 200 h-10 v1 h10" fill="black" transform="rotate(-20,320,200)" />
<path d="M20 200 h300 v1 h-300" fill="black" />
</g>
<path d="M170 200 h87 v1 h-87" fill="black" transform="rotate(-10,170,200)"/>
<path d="M170 50 v10 h1 v-10" fill="black" transform="rotate(20,170,50)" />
<path d="M170 50 v10 h1 v-10" fill="black" transform="rotate(-20,170,50)" />
<path d="M170 50 v300 h1 v-300" fill="black" />
<path d="M320 200 h-10 v1 h10" fill="black" transform="rotate(20,320,200)" />
<path d="M320 200 h-10 v1 h10" fill="black" transform="rotate(-20,320,200)" />
<path d="M20 200 h300 v1 h-300" fill="black" />
<ellipse cx="170" cy="200" rx="100" ry="25" fill="none" stroke="red" stroke-width="2" />
<circle cx="253" cy="186" r="3" fill="white" />
<circle cx="253" cy="186" r="3" fill="none" stroke="black" />
<text x="325" y="205" fill="black" font-size="16">x</text>
<text x="305" y="125" fill="black" font-size="16">y</text>
<text x="166" y="46"  fill="black" font-size="16">z</text>

<path d="M500 200 h100 A100 100 0 0 0 571 129" fill="yellow" />
<path d="M500 200 h100 v1 h-100" fill="black" transform="rotate(-45,500,200)" />
<path d="M500 50 v10 h1 v-10" fill="black" transform="rotate(20,500,50)" />
<path d="M500 50 v10 h1 v-10" fill="black" transform="rotate(-20,500,50)" />
<path d="M500 50 v300 h1 v-300" fill="black" />
<path d="M650 200 h-10 v1 h10" fill="black" transform="rotate(20,650,200)" />
<path d="M650 200 h-10 v1 h10" fill="black" transform="rotate(-20,650,200)" />
<path d="M350 200 h300 v1 h-300" fill="black" />
<circle cx="500" cy="200" r="100" fill="none" stroke="red" stroke-width="2" />
<circle cx="571" cy="129" r="3" fill="white" />
<circle cx="571" cy="129" r="3" fill="none" stroke="black" />
<circle cx="500" cy="200" r="6" fill="white" />
<circle cx="500" cy="200" r="2" fill="black" />
<circle cx="500" cy="200" r="6" fill="none" stroke="black" />
<path d="M585 130 h15 v3 l15 -5 l-15 -5 v3 h-15" fill="blue" />
<text x="655" y="205" fill="black" font-size="16">x</text>
<text x="496" y="45"  fill="black" font-size="16">y</text>
<text x="485" y="190" fill="black" font-size="16">z</text>
<text x="515" y="195" fill="black" font-size="16">ϕ</text>
<text x="575" y="117" fill="black" font-size="16">潮汐力</text>

<path d="M170 540 v-100 A100 100 0 0 0 136 446" fill="cyan" />
<path d="M170 440 v200 h1 v-200" fill="black" transform="rotate(-20,170,540)" />
<g transform="rotate(-20,250,540)">
<path d="M175 530 h10 v1 h-10" fill="black" />
<path d="M185 530 v10 h1 v-10" fill="black" />
<path d="M250 540 h-75 v1 h75" fill="black" />
</g>
<circle cx="170" cy="540" r="100" fill="none" stroke="black" />
<path d="M170 390 v10 h1 v-10" fill="black" transform="rotate(20,170,390)" />
<path d="M170 390 v10 h1 v-10" fill="black" transform="rotate(-20,170,390)" />
<path d="M170 390 v300 h1 v-300" fill="black" />
<path d="M320 540 h-10 v1 h10" fill="black" transform="rotate(20,320,540)" />
<path d="M320 540 h-10 v1 h10" fill="black" transform="rotate(-20,320,540)" />
<path d="M20 540 h300 v1 h-300" fill="black" />
<path d="M70 539 h200 v3 h-200" fill="red" />
<circle cx="250" cy="540" r="3" fill="white" />
<circle cx="250" cy="540" r="3" fill="none" stroke="black" />
<circle cx="170" cy="540" r="6" fill="white" />
<path d="M170 534 v12 h1 v-12" fill="black" transform="rotate(45,170,540)" />
<path d="M165 540 h12 v1 h-12" fill="black" transform="rotate(45,170,540)" />
<circle cx="170" cy="540" r="6" fill="none" stroke="black" />
<path d="M260 542 h15 v3 l15 -5 l-15 -5 v3 h-15" fill="blue" transform="rotate(-20,250,540)" />
<path d="M248 530 v-15 h-3 l5 -15 l5 15 h-3 v15" fill="lime" />
<text x="325" y="545" fill="black" font-size="16">x</text>
<text x="152" y="555" fill="black" font-size="16">y</text>
<text x="166" y="385" fill="black" font-size="16">z</text>
<text x="290" y="525" fill="black" font-size="16">潮汐力</text>
<text x="226" y="495"  fill="black" font-size="16">トルク</text>
<text x="159" y="500"  fill="black" font-size="16">ε</text>
</svg>
title
<$FIGURE>
図２　潮汐力によるトルク
end title
@enduml
```
球対称の部分はトルクが相殺されるため無視してよい。このため非対称部分，すなわち赤道上に置いたリングに作用するトルクについてのみ考えればよい。

赤道リングの経度 $\phi$ において潮汐力を受ける部分の質量 $dm$ を示す。

```math
dm = \mu \cdot r_p \cdot d\phi \tag{4}
```

太陽から受ける潮汐力 $df_s$ は地球の中心を通り，かつ太陽方向と鉛直な平面までの距離 $r_p \cos \phi \cos \varepsilon$ に比例する。

```math
df_s = \frac{2GM_s}{R_e^3} \cdot dm \cdot r_p \cos \phi \cos \varepsilon \tag{5}
```

太陽から受けるトルク $d\tau_s$ は潮汐力 $df_s$ のうち赤道面と垂直な成分 $df_s \sin \varepsilon$ とY軸までの距離 $r_p\cos\phi$ の積である。

```math
\begin{align}
d\tau_s &= df_s \cdot r_p \cos \phi \sin \varepsilon \\
&= \frac{GM_s}{R_e^3} \cdot dm \cdot r_p^2  \cos^2 \phi \sin 2\varepsilon \\
&= \frac{GM_s}{R_e^3} \cdot \mu \cdot r_p^3  \cos^2 \phi \sin 2\varepsilon \cdot d\phi \tag{6}
\end{align}
```

地球全体が太陽から受けるトルク $\tau_s$ は，これを積分したものなので

```math
\begin{align}
\tau_s &= \int d\tau_s \\
&= \int_{-\pi}^{\pi} \frac{GM_s}{R_e^3} \cdot \mu \cdot r_p^3 \cos^2 \phi \sin 2 \varepsilon \cdot d\phi \\
&= \frac{GM_s}{R_e^3} \cdot \mu \cdot r_p^3 \sin 2 \varepsilon \int_{-\pi}^{\pi} \cos^2 \phi \cdot d\phi \\
&= \frac{GM_s}{R_e^3} \cdot \mu \cdot \pi \cdot r_p^3 \sin 2 \varepsilon \\
&= \frac{2GM_s}{3R_e^3} \cdot \rho \cdot \pi \cdot r_p^3 \cdot (r_e^2 - r_p^2) \sin 2 \varepsilon \tag{7}
\end{align}
```

となる。同様にして地球が月から受けるトルク $\tau_m$ を求めると

```math
\tau_m = \frac{2GM_m}{3R_m^3} \cdot \rho \cdot \pi \cdot r_p^3 \cdot (r_e^2 - r_p^2) \sin 2\varepsilon \tag{8}
```

となる。表１より，それぞれのトルクを求めると 

```math
\left \{ \begin{aligned}
\tau_s &= 1.147 \times 10^{37} \\
\tau_m &= 2.421 \times 10^{38}
\end{aligned} \right. \tag{9}
```

となる。

:::note
潮汐力（トルク）に関しては太陽より月の影響のほうが２倍以上大きい。
:::

<table>
<caption>表１　定数表</caption>
<thead>
<tr><th style="text-align:center">項目</th>
    <th style="text-align:center">記号</th>
    <th style="text-align:center">数値</th>
    <th style="text-align:center">単位</th></tr>
</thead>
<tbody>
<tr><td>地球の赤道半径</td>
    <td style="text-align:center">$r_e$</td>
    <td>$6,378,178$</td>
    <td>[m]</td></tr>
<tr><td>地球の極半径</td>
    <td style="text-align:center">$r_p$</td>
    <td>$6,356,752$</td>
    <td>[m]</td></tr>
<tr><td>太陽質量</td>
    <td style="text-align:center">$M_s$</td>
    <td>$1.989 \times 10^{30}$</td>
    <td>[kg]</td></tr>
<tr><td>月質量</td>
    <td style="text-align:center">$M_m$</td>
    <td>$7.346 \times 10^{22}$</td>
    <td>[kg]</td></tr>
<tr><td>重力定数</td>
    <td style="text-align:center">$G$</td>
    <td>$6.674 \times 10^{-11}$</td>
    <td>[m³･kg⁻¹･s⁻²]</td></tr>
<tr><td>地球公転半径</td>
    <td style="text-align:center">$R_e$</td>
    <td>$1.496 \times 10^{11}$</td>
    <td>[m]</td></tr>
<tr><td>月公転半径</td>
    <td style="text-align:center">$R_m$</td>
    <td>$3.884 \times 10^8$</td>
    <td>[m]</td></tr>
<tr><td>地球の平均地殻密度</td>
    <td style="text-align:center">$\rho$</td>
    <td>$2.7 \times 10^3$</td>
    <td>[kg/m³]</td></tr>
<tr><td>地球の地軸の傾き</td>
    <td style="text-align:center">$\varepsilon$</td>
    <td>$23.44$</td>
    <td>[deg]</td></tr>
</tbody>
</table>

# 年間を通じた平均トルク

先ほど求めたのは（敢えて言明しなかったが）冬至のとき地球が受ける最大トルクである。春分・秋分のときは太陽が $(0,\pm R_e,0)$ の位置にあると考えると地球の受けるトルクはゼロとなる。一方，夏至のときは太陽が $(-R_e\cos(-\varepsilon),0,-R_e\sin(-\varepsilon))$ の位置にあると考えるとトルクは冬至のときと等しく最大値を取る。

```plantuml
@startuml
<style>
root {
  margin: 0;
  padding: 0;
}
title {
  fontSize: 16;
  fontStyle: normal;
}
</style>
sprite FIGURE <svg width="500" height="300" viewbox="0 0 500 300">
<path d="M50 250 c40 0 60 -200 100 -200 c40 0 60 200 100 200 c40 0 60 -200 100 -200 c40 0 60 200 100 200" fill="pink" />
<path d="M20 250 h460 v1 h-460" fill="black" />
<path d="M50 50 v10 h1 v-10" fill="black" transform="rotate(20,50,50)" />
<path d="M50 50 v10 h1 v-10" fill="black" transform="rotate(-20,50,50)" />
<path d="M50 250 v-200 h1 v200" fill="black" />
<path d="M50 148 h400 v2 h-400" fill="red" />
<text x="27"  y="40" fill="black" font-size="16">トルク</text>
<text x="35"  y="270" fill="black" font-size="16">春分</text>
<text x="135" y="270" fill="black" font-size="16">夏至</text>
<text x="235" y="270" fill="black" font-size="16">秋分</text>
<text x="335" y="270" fill="black" font-size="16">冬至</text>
<text x="435" y="270" fill="black" font-size="16">春分</text>
<text x="455" y="155" fill="black" font-size="16">平均</text>
</svg>
title
<$FIGURE>
図３　年間を通じた平均トルク（太陽）
end title
@enduml
```
<path d="M50 250 q22 0 50 -100 q28 -100 50 -100 q22 0 50 100 q 28 100 50 100 q22 0 50 -100 q28 -100 50 -100 q22 0 50 100 q 28 100 50 100" fill="pink" />

最大値と最小値の間を直線で結ぶのかスムーズに正弦波状に繋げるべきかは不問として，年間を通じた平均トルクは夏至・冬至のとき（最大値）のちょうど半分になると考えられる。月も同様に考えると年間を通じた平均トルク $\bar{\tau}$ は

```math
\bar{\tau} = \frac{1}{2} ( \tau_s + \tau_m ) \tag{10}
```

となると考えられる。

# 慣性モーメントを求める

簡単のため地球内部の密度を一様とすると自転の慣性モーメント $I_e$ は

```math
I_e = \frac{2}{5} M_e \cdot r_p^2 = 9.653 \times 10^{37} \tag{11}
```

となるが，さすがに誤差が大きい。※実際の値に比べて２割ほど大きい。

<table>
<caption>表２　定数表</caption>
<thead>
<tr><th style="text-align:center">項目</th>
    <th style="text-align:center">記号</th>
    <th style="text-align:center">数値</th>
    <th style="text-align:center">単位</th></tr>
</thead>
<tbody>
<tr><td>地球の質量</td>
    <td style="text-align:center">$M_e$</td>
    <td>$5.9724 \times 10^{24}$</td>
    <td>[kg]</td></tr>
</tbody>
</table>

地球の内部構造は詳しく分かっているので精密に求めてみよう。地球の内部密度は以下のように地球の中心に近づくほど高くなり，層の境界で不連続点を持つ。

```plantuml
@startuml
<style>
root {
  margin: 0;
  padding: 0;
}
title {
  fontStyle: normal;
  fontSize: 16;
}
</style>
sprite FIGURE <svg width="690" height="500" viewbox="0 0 690 500">
<path d="M30 450 L30 47 L152 57 L152 450" fill="yellow" />
<path d="M152 450 L152 75 L378 143 L378 450" fill="orange" />
<path d="M378 450 L378 273 L600 309 L600 450" fill="tomato" />
<path d="M600 450 L600 320 L664 339 L664 450" fill="red" />
<path d="M664 450 L664 353 L666 362 L666 450" fill="brown" />
<path d="M666 450 L666 409 L667 409 L667 450" fill="blue" />
<path d="M30 450 h650 v1 h-650" fill="black" />
<path d="M30 450 v-420 h1 v420" fill="black" />
<path d="M130 450 v-5 h1 v5" fill="black" />
<path d="M230 450 v-5 h1 v5" fill="black" />
<path d="M330 450 v-5 h1 v5" fill="black" />
<path d="M430 450 v-5 h1 v5" fill="black" />
<path d="M530 450 v-5 h1 v5" fill="black" />
<path d="M630 450 v-5 h1 v5" fill="black" />
<path d="M30 420 h5 v1 h-5" fill="black" />
<path d="M30 390 h5 v1 h-5" fill="black" />
<path d="M30 360 h5 v1 h-5" fill="black" />
<path d="M30 330 h5 v1 h-5" fill="black" />
<path d="M30 300 h5 v1 h-5" fill="black" />
<path d="M30 270 h5 v1 h-5" fill="black" />
<path d="M30 240 h5 v1 h-5" fill="black" />
<path d="M30 210 h5 v1 h-5" fill="black" />
<path d="M30 180 h5 v1 h-5" fill="black" />
<path d="M30 150 h5 v1 h-5" fill="black" />
<path d="M30 120 h5 v1 h-5" fill="black" />
<path d="M30  90 h5 v1 h-5" fill="black" />
<path d="M30  60 h5 v1 h-5" fill="black" />
<path d="M30  30 h5 v1 h-5" fill="black" />
<text x="26" y="470" fill="black" font-size="16">0</text>
<text x="106" y="470" fill="black" font-size="16">1,000</text>
<text x="206" y="470" fill="black" font-size="16">2,000</text>
<text x="306" y="470" fill="black" font-size="16">3,000</text>
<text x="406" y="470" fill="black" font-size="16">4,000</text>
<text x="506" y="470" fill="black" font-size="16">5,000</text>
<text x="606" y="470" fill="black" font-size="16">6,000</text>
<text x="15" y="456" fill="black" font-size="16">0</text>
<text x="15" y="426" fill="black" font-size="16">1</text>
<text x="15" y="396" fill="black" font-size="16">2</text>
<text x="15" y="366" fill="black" font-size="16">3</text>
<text x="15" y="336" fill="black" font-size="16">4</text>
<text x="15" y="306" fill="black" font-size="16">5</text>
<text x="15" y="276" fill="black" font-size="16">6</text>
<text x="15" y="246" fill="black" font-size="16">7</text>
<text x="15" y="216" fill="black" font-size="16">8</text>
<text x="15" y="186" fill="black" font-size="16">9</text>
<text x="5"  y="156" fill="black" font-size="16">10</text>
<text x="5"  y="126" fill="black" font-size="16">11</text>
<text x="5"  y="96"  fill="black" font-size="16">12</text>
<text x="5"  y="66"  fill="black" font-size="16">13</text>
<text x="5"  y="36"  fill="black" font-size="16">14</text>
<text x="270" y="490" fill="black" font-size="16">地球中心からの距離[km]</text>
<text x="30" y="25" fill="black" font-size="16">密度[g･cm³]</text>
<text x="75" y="260"  fill="black" font-size="16">内核</text>
<text x="245" y="280" fill="black" font-size="16">外核</text>
<text x="440" y="375" fill="black" font-size="16">下部マントル</text>
<text x="590" y="300" fill="black" font-size="16">上部マントル</text>
<path d="M620 305 v30 h1 v-30" fill="black" transform="rotate(20,620,305)" />
</svg>
title
<$FIGURE>
図４　地球内部の密度
end title
@enduml
```

:::note
参考文献[11]ではもっと細かく層が分かれているが，簡単のため６層にまとめた。また三次式でモデル化しているところ，境界値の間を単純に線形補間（一次近似）した。
:::

まず半径 $r$，密度 $\rho(r)$ の球殻の慣性モーメント $I(r)$ を求める。

```math
dI = \rho(r) \cdot r^2 \cos \theta \cdot d\theta \cdot d\phi \cdot (r \cos \theta)^2 \tag{12}
```

```math
\begin{align}
I(r) &= \iint dI \\
&= \rho(r) \cdot r^4 \int_{-\pi/2}^{\pi/2} \cos^3 \theta \cdot d\theta \int_{-\pi}^{\pi} d\phi \\
&= \frac{8\pi}{3} \cdot \rho(r) \cdot r^4 \tag{13}
\end{align}
```

地球中心からの距離 $r_0 \le r \le r_1$ において密度関数 $\rho(r)$ の境界値が次のように与えられたとする。

```math
\left\{ \begin{aligned}
\rho(r_0) &= \rho0 \\
\rho(r_1) &= \rho1 \\
\end{aligned} \right. \tag{14}
```

区間 $r_0 \le r \le r_1$ の間，密度は線形補間するものとする。

```math
\begin{align}
\rho(r) &= \rho_0 + \frac{\rho_1 - \rho_0}{r_1 - r_0} \cdot (r - r_0) \\
&= \frac{\rho_0 \cdot r_1 - \rho_1 \cdot r_0}{r_1 - r_0} + \frac{\rho_1 - \rho_0}{r_1 - r_0} \cdot r \tag{15}
\end{align}
```

区間 $r_0 \le r \le r_1$ の慣性モーメント $I[r_0, r_1]$ を求めると

```math
\begin{align}
I[r_0, r_1] &= \int_{r_0}^{r_1} I(r) \cdot dr \\
&= \int_{r_0}^{r_1} \frac{8\pi}{3} \cdot \rho(r) \cdot r^4 \cdot dr \\
&= \frac{8\pi}{3} \int_{r_0}^{r_1} \left( \frac{\rho_0 \cdot r_1 - \rho_1 \cdot r_0}{r_1 - r_0} + \frac{\rho_1 - \rho_0}{r_1 - r_0} \cdot r \right)  \cdot r^4 \cdot dr \\
&= \frac{8\pi}{3} \left[ \frac{\rho_0 \cdot r_1 - \rho_1 \cdot r_0}{r_1 - r_0} \cdot \frac{r^5}{5} + \frac{\rho_1 - \rho_0}{r_1 - r_0} \cdot \frac{r^6}{6} \right]_{r_0}^{r_1} \\
&= \frac{8\pi}{3} \left\{ \frac{\rho_0 \cdot r_1 - \rho_1 \cdot r_0}{r_1 - r_0} \cdot \frac{r_1^5 - r_0^5}{5} + \frac{\rho_1 - \rho_0}{r_1 - r_0} \cdot \frac{r_1^6 - r_0^6}{6} \right\} \tag{16}
\end{align}
```

となる。こうして各層ごとの慣性モーメントを求め，これらの総和により地球全体の慣性モーメントを求めたところ $I_e = 8.0359 \times 10^{37}$ を得られた。

<table>
<caption>表３　地球内部密度と慣性モーメント</caption>
<thead>
<tr><th style="text-align:center">名称</th>
    <th style="text-align:center">中心からの距離<br>[km]</th>
    <th style="text-align:center">密度<br>[g/cm³]</th>
    <th style="text-align:center">慣性モーメント<br>[kg･m²]</th></tr>
</thead>
<tbody>
<tr><td style="text-align:center">内核</td>
    <td style="text-align:right">$0 \rightarrow 1,222$</td>
    <td style="text-align:right">$13.089 \rightarrow 12.764$</td>
    <td style="text-align:right">$5.8402 \times 10^{34}$</td></tr>
<tr><td style="text-align:center">外核</td>
    <td style="text-align:right">$1,222 \rightarrow 3,480$</td>
    <td style="text-align:right">$12.166 \rightarrow\:\: 9.903$</td>
    <td style="text-align:right">$8.9096 \times 10^{36}$</td></tr>
<tr><td style="text-align:center">下部マントル</td>
    <td style="text-align:right">$3,480 \rightarrow 5,701$</td>
    <td style="text-align:right">$5.566 \rightarrow\:\: 4.381$</td>
    <td style="text-align:right">$4.4296 \times 10^{37}$</td></tr>
<tr><td style="text-align:center">上部マントル</td>
    <td style="text-align:right">$5,701 \rightarrow 6,347$</td>
    <td style="text-align:right">$3.992 \rightarrow\:\: 3.381$</td>
    <td style="text-align:right">$2.6248 \times 10^{36}$</td></tr>
<tr><td style="text-align:center">地殻</td>
    <td style="text-align:right">$6,347 \rightarrow 6,368$</td>
    <td style="text-align:right">$2.900 \rightarrow\:\: 2.600$</td>
    <td style="text-align:right">$8.0520 \times 10^{35}$</td></tr>
<tr><td style="text-align:center">海洋</td>
    <td style="text-align:right">$6,368 \rightarrow 6,371$</td>
    <td style="text-align:right">$1.020$</td>
    <td style="text-align:right">$4.2195 \times 10^{34}$</td></tr>
<tr><td style="text-align:center">合計</td>
    <td></td>
    <td></td>
    <td style="text-align:right">$8.0359 \times 10^{37}$</td></tr>
<tr><td colspan=4></td></tr>
</tbody>
</table>

同様にして地球の質量も計算できる。区間 $r_0 \le r \le r_1$ の質量 $M[r_0, r_1]$ は

```math
M[r_0, r_1] = 4\pi \left\{ \frac{\rho_0 \cdot r_1 - \rho_1 \cdot r_0}{r_1 - r_0} \cdot \frac{r_1^3 - r_0^3}{3} + \frac{\rho_1 - \rho_0}{r_1 - r_0} \cdot \frac{r_1^4 - r_0^4}{4} \right\} \tag{17}
```

であるから，各層ごとの質量を求め，これらの総和により地球全体の質量を求めたところ $M_e = 5.9455 \times 10^{24} $ [kg] を得られた。これは表２に示す真値の 99.55% に当たる。

# 歳差周期を求める

地球の自転周期 $T_r = 2\pi / \omega_r$，自転角速度 $\omega_r$，慣性モーメント $I_e$ とおくと地球の自転による角運動量ベクトルの大きさ $L = I_e \cdot \omega_r$ となり，歳差運動による角速度 $\omega_p$ とおくと角運動量ベクトルの時間変化率が地球に加わるトルクと等しいので

```math
\omega_p \cdot L \sin \varepsilon = \bar{\tau} \tag{18}
```

となる。

```plantuml
@startuml
<style>
root {
  margin: 0;
  padding: 0;
}
title {
  fontStyle: normal;
  fontSize: 16;
}
</style>
sprite FIGURE <svg width="500" height="450" viewbox="0 0 500 450">
<path d="M250 100 l-48.2 17.5 A100 20 0 0 0 226.8 119.5" fill="yellow" />
<path d="M250 100 v30.2 h1 v-30.2" fill="silver" transform="rotate(50,250,100)" />
<path d="M250 350 v-50 A 50 50 0 0 1 268.6 303.6" fill="cyan" />
<path d="M250 350 v-205 h1 v205" fill="black" />
<path d="M250 50 v75 h1 v-75" fill="black" />
<path d="M250 50 v10 h1 v-10" fill="black" transform="rotate(20,250,50)" />
<path d="M250 50 v10 h1 v-10" fill="black" transform="rotate(-20,250,50)" />
<path d="M50 350 h400 v1 h-400" fill="black" />
<path d="M450 350 h-10 v1 h10" fill="black" transform="rotate(20,450,350)" />
<path d="M450 350 h-10 v1 h10" fill="black" transform="rotate(-20,450,350)" />
<g transform="rotate(-30,250,350)">
<path d="M50 350 h400 v1 h-400" fill="black" />
<path d="M450 350 h-10 v1 h10" fill="black" transform="rotate(20,450,350)" />
<path d="M450 350 h-10 v1 h10" fill="black" transform="rotate(-20,450,350)" />
</g>
<ellipse cx="250" cy="100" rx="100" ry="20" fill="none" stroke="silver" />
<path d="M150 100 h200 v1 h-200" fill="silver" />
<path d="M200 100 h100 v1 h-100" fill="silver" transform="rotate(-20,250,100)" />
<path d="M250 350 v-269.3 h1 v269.3" fill="silver" transform="rotate(21.8,250,350)" />
<path d="M250 350 v-269.3 h1 v269.3" fill="silver" transform="rotate(-21.8,250,350)" />
<g transform="rotate(-11.7,250,350)">
<path d="M249 350 v-222.3 h-2 l3 -15 l3 15 h-2 v222.3" fill="red" />
</g>
<g transform="rotate(-5.7,250,350)">
<path d="M249 350 v-216.7 h-2 l3 -15 l3 15 h-2 v216.7" fill="red" />
</g>
<path d="M202.8 116.5 h12 v-2 l13 3 l-13 3 v-2 h-12" fill="blue" transform="rotate(4.6,201.8,117.5)" />
<text x="246" y="370" fill="black" font-size="16">0</text>
<text x="455" y="355" fill="black" font-size="16">x</text>
<text x="430" y="250" fill="black" font-size="16">y</text>
<text x="246" y="45"  fill="black" font-size="16">z</text>
<text x="258" y="297" fill="black" font-size="16">ε</text>
<text x="180" y="127" fill="red" font-size="16">→</text>
<text x="183" y="140" fill="red" font-size="16">L</text>
<text x="232" y="127" fill="red" font-size="16">→</text>
<text x="265" y="127" fill="red" font-size="16">→</text>
<text x="235" y="140" fill="red" font-size="16">L+ΔL</text>
<text x="192" y="100" fill="blue" font-size="16">→</text>
<text x="185" y="113" fill="blue" font-size="16">ΔL</text>
</svg>
title
<$FIGURE>
<latex>|\Delta \vec{L}| = \omega_p\, L \sin\varepsilon\Delta t = \bar{τ}\,\Delta t</latex>
図５　角運動量ベクトルの時間変化
end title
@enduml
```

:::note
あくまで角運動量ベクトルの大きさ $|\vec{L}| = L$ は不変であり，作用するトルクによってベクトルの向きが変わるだけと理解されたい。
:::

歳差運動の角速度 $\omega_p$ は

```math
\omega_p = \frac{\bar{\tau}}{L \sin \varepsilon} = \frac{1}{2} \cdot \frac{\tau_s + \tau_m}{\omega_r \cdot I_e \sin \varepsilon} \tag{19}
```

より $\omega_p = 7.653 \times 10^{-12}$ となり，歳差運動の周期 $T_p = 2\pi / \omega_r = 8.210 \times 10^{11}$ が得られる。

```math
\begin{aligned}
\frac{8.210 \times 10^{11}\textsf{［秒］}}{24 \times 60^2 \textsf{［秒/日］}} &= 9.503 \times 10^6 \textsf{［日］} \\
\frac{9.503 \times 10^6 \textsf{［日］}}{365.25 \textsf{［日/年］}} &= 26,017 \textsf{［年］} 
\end{aligned} \tag{20}
```

<table>
<caption>表４　定数表</caption>
<thead>
<tr><th style="text-align:center">項目</th>
    <th style="text-align:center">記号</th>
    <th style="text-align:center">数値</th>
    <th style="text-align:center">単位</th>
    <th style="text-align:center">備考</th></tr>
</thead>
<tbody>
<tr><td>地球自転周期</td>
    <td style="text-align:center">$T_r$</td>
    <td>$86,164$</td>
    <td>[s]</td>
    <td>23時間56分04秒</td></tr>
</tbody>
</table>

# まとめ

以上より歳差周期 2万6017年を得ることができた。真値2万5772年と比べると素晴らしい精度（+1.0%）で得られている。

このうち赤道リングモデルについては，地球が受けるトルクを（奇跡的に）ほぼ正しく計算することができたようだ。

地球の慣性モーメントについては地球の中心からの距離と密度のみで計算していることから球対称であり，冒頭で述べたように地球は赤道方向に少し膨らんだ回転楕円体であることを考慮していない。これに加えて粗い線形近似モデルにも関わらず，こちらもまた驚異的な精度を得られている。

# 参考文献

[1] [国立天文台 ＞暦計算室＞暦Wiki＞惑星　太陽系天体の半径](https://eco.mtk.nao.ac.jp/koyomi/wiki/CFC7C0B12FC8BEB7C2.html)
[2] [地球内部物理学 - wikipedia](https://ja.wikipedia.org/wiki/地球内部物理学)
[3] [地球の地殻 - wikipedia](https://ja.wikipedia.org/wiki/地球の地殻)
[4] [地球の構造 - wikipedia](https://ja.wikipedia.org/wiki/地球の構造)
[5] [重力定数 - wikipedia](https://ja.wikipedia.org/wiki/万有引力定数)
[6] [地球 - wikipedia](https://ja.wikipedia.org/wiki/地球)
[7] [太陽 - wikipedia](https://ja.wikipedia.org/wiki/太陽)
[8] [月 - wikipedia](https://ja.wikipedia.org/wiki/月)
[9] [歳差 - wikipedia](https://ja.wikipedia.org/wiki/歳差)
[10] [地圏－層構造 - 名古屋市科学館 <img width=16 src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3490655/66e66c9d-a6f3-f215-951a-9d79cc94a9d6.png">](http://www.ncsm.city.nagoya.jp/exhibit_files/pdf/S608.pdf)
[11] [Preliminary reference Earth model, Physics of the Earth and Planetary Interiors, 25 (1981) <img width=16 src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3490655/66e66c9d-a6f3-f215-951a-9d79cc94a9d6.png">](https://www.cfa.harvard.edu/~lzeng/papers/PREM.pdf)