「国立科学博物館で学ぶ物理学」メニュー

ジャイロホイール

回折

HR図

国立科学博物館で学ぶ物理学　＜旧暦／月の運行軌道＞

[概説]

[理論解析]

[Ⅰ] 旧暦

(1)太陰太陽暦

皆さんが慣れ親しんでいる暦とは、太陽の運行を基準とした太陽暦ですが、この他に月の運行を基準とした太陰歴という暦があります。イスラム教の世界では未だに純然たる太陰暦が使用されています。日本でも、例えば１月はどの年であっても冬であるという季節が感じられる暦として、太陽の運行も取り入れた太陰太陽歴が用いられていました。江戸時代以前までは中国で作成された暦を使用していましたが、江戸初期に渋川春海により初めて日本での観測に基ずく大和暦が作成されました。それが貞享暦（左下図に示す展示物は貞享暦に従う享保歴です）というものです。江戸末期には、太陽の動きの取り込み方（平気法と定気法）が貞享暦とは若干異なる暦として天保歴が作成されました。旧暦と言い習わしている暦とはこの天保歴のことを言います。

太陰太陽歴の基本は、

(1)定朔法(朔月) 
    新月（朔）の日を月の初日とする。


(2)連大法(大小暦) 
    ひと月は約29.5日のため、月の
    日数は大（３０日）小（２９日）
    を配置する。必要に応じ大の月
　　(または小の月)を連続配置する。


(3)置閏法(閏月)
    ひと月29.5日の12ケ月では太陽年
    に対し11日不足し、同じ月でも
    季節が変わるので、季節感を保つ
    ために、閏月（約３年 に一度）を
    配置する。

こうした基本原則により暦を作成しますが、太陽と月の運行を正確に予測する必要があります。それには観測が不可欠であり、渾天儀という機材を用いてなされていました。渾天儀の観測原理と天体運行原理を示す模型を右上の図に示しますが、これは教育用として使用されたものです。この模型により、月と太陽が一直線に並ぶ現象、すなわち新月・満月、日食・月食が理解し易くなります。

  (1)中心に地球があります。
　(2)水平面（水平のリング）は観測地点での地平面です。(江戸は北緯36度ですので自転軸は36度傾いています）
　(3)地球を貫く黒い棒は黄道面（地球軌道の面、あるいは地球に対する太陽の軌道面）の法線です
  (4)赤道（地球の赤道）黄道（太陽の軌道）白道（月の軌道）は右上から左下にかけてのリングで表され、
　　 その操作は左上から右下にかけてのリングで行われます。
　　　　　赤道：外から２番目のリング・・・外から２番目の操作リングをねじると太陽・月が動きます（日周運動）                                                                                      
　　　　　黄道：外から３番目のリング・・・外から３番目の操作リングをねじると太陽が動きます（年周運動）                                                                                     
　　　　　白道：外から４番目のリング・・・外から４番目の操作リングをねじると月の軌道面が動きます
     黄道と赤道の交差角は約23.4度、黄道と白道との交差角は約5度となっています。

観測といっても基準がなければ観測結果を表せません。その基準となるものが星（恒星）で、そのために右に示す天球儀の上での星の位置を正確に求めたわけです（星座そのものにあまり興味は向いていなかったと言えます）。図中、赤い線の大円が赤道で、黒っぽい線の大円が黄道です。白道は絶えず移動するので図中に固定して描くことができません。

(2)太陽・月運行全体感と諸データ

太陰太陽暦は、右図に示すような天球上での太陽・月の運行を元にしています。太陽は赤道に対し約23.4度傾斜した大円（黄道）上を１年かけて動きます。(展示物の黄道と位置がちがうようにみえますが、夏か冬か（または昼か夜）で太陽の位置が違うためにそのように見えます。）

月は、黄道に対し約５度傾斜した大円（白道）上を旧暦１ケ月かけて動きますが、太陽の引力のため白道面は同じ所にはなく、白道面の法線が黄道面法線周りを図示する方向に回転するように変化していきます。そのため、月が黄道面と交差する点である交点は、月の進行方向と逆方向に移動します。

この状況を、赤道座標および黄道座標で表したものが右下の図です。

(天体座標詳細は軌道要素・天球座標を参照してください。）

現代の観測による、太陽・月の運行データは以下の通りです。

   (1)太陽　　　①恒星年T₀　　　：　365.25日　　　②太陽年T₀' 　　：　365.24日
　　　　　　　　③地軸歳差周期T_e：　  約25800年　　　

　 (2)月　　　　①恒星月T₁　　　：　 27.32日　　　②朔望月T₂　　　：　 29.53日
　　　　　　　　③交点月T₃　　　：　 27.21日　　　④近点月T₄　　　：　 27.55日
                ⑤交点逆行周期T_g：   18.61年　　  ⑥近点順行周期T_k：　 8.85年

よく見かける関係式を列記しておきます。(これらは次項(3)で説明するような円軌道の仮定による近似式です。）

   (0)太陽年　: 1/T₀'=1/T₀+1/T_e
　 
   (1)朔望月　: 1/T₂=1/T₁-1/T₀
   (2)交点月　: 1/T₃=1/T₁+1/T_g
   (3)近点月　: 1/T₄=1/T₁-1/T_k
   (4)食年　  : 1/T_s=1/T₃-1/T₂ = 346.62日
   (5)ﾒﾄﾝ周期 : T_m=最小公倍数(T₂,T₀）=6939.6日（ T₂/T₀=235朔望月/19年)　　・・・・(19年に7回の閏月を置く周期と同じ）
   (6)ｻﾛｽ周期 : T_t=最小公倍数(T₂,T_s）=6585.32日（ T₂/T_s=223朔望月/18.03年)・・・・(約18年毎に同じ日食が発生する周期）

(3)朔望月日数のの初歩的説明

朔望月のもっとも初歩的な説明は、左図のように月が地球の周りを一周してもその間に太陽が少し動くので、太陽と月が一直線に並ぶには、さらに月はもう少し進まなければならないため、朔望月（すなわち新月から新月まで、もしくは満月から満月までの日数）は、恒星月（同じ星の方向に見えるように周回する日数）より長くなります。数値的に言うと、太陽も月も同一平面上を円軌道で運動しているとすると、太陽・月の角速度をΩ（=360/365.25)・ω(=360/27.32)とし、また一直線になるのは一度一直線になった時からの時間を t として、(ωt－Ωt)=360°の時であるから、よく見かける関係式として、
　　　　　　　　　 $t=\frac{1}{\frac{1}{27.32}-\frac{1}{365.25}}=29.53\hspace{80mm}(\frac{1}{T_2}=\frac{1}{T_1}-\frac{1}{T_0})$
すなわち、約29.5日毎に新月となるので、大の月(30日）と小の月（29日）を大－小－大－小－大－小－と交互に配置すればよいことになります。しかし現実に大－小－大－大－小－大－小－などと「連大」などが現れることを、説明できません。

皆さんはすでに、太陽も月も地球に対して楕円軌道上を運動していると習われているでしょうから、やはり同一平面上を独立の楕円軌道で運動しているとして考えてみましょう。計算式は少々ややこしいので以下に述べる理論解析／時間計算を参照していただくとして、結果を示すと、１年間の変動（左下図参照：横軸／１年間の周回数、縦軸／朔望月、注釈／２つの楕円の長軸のなす角度θ）は、

　　　　　　　　　　　　　　　　　　t = 29.3～29.7日（θ= 0度), 29.2～29.8日（θ= 180度)

２つの楕円軌道を考えると、朔望月が一定の29.5日ではなく変動すること、すなわち「連大」などが現れることが理解できます。
２つの楕円が独立固定ではなく角度θが次第に変わっていくことで、現実の朔望月の変化（すなわち大－小の並びの変化）が起こることも理解できます。

しかし、ある年は連大（大－大）、別の年では３連大（大－大－大）、またまた別の年では、４連大（大－大－大－大）が現れるということを説明するには、このθの時間変化を具体的かつ詳細に求めなければなりません。



（なお、２つの軌道が交差している（傾斜角約５度を有する）影響は小さいので、双方の議論では無視しています。）

(4)二十四節気と閏月

二十四節気とは、１２の中気と１２の節気を言い、中気の雨水を含む月を旧暦の１月とします。（右図参照）

季節との整合性を要請した太陽・太陰暦は、太陽の位置を定めるに、旧暦と称される天保歴は太陽の春分からの角度を24等分した角度（15度）毎の位置を基準とした定気法を、それ以前の歴法は太陽が冬至から１年の24等分した時間（約15.2日)毎にいる位置を基準とした平気法を用いていました。平気法では、中気の間隔(約30.4日)と朔望月(約29.5日)とほぼ同じ間隔であるので各月の中に中気がほぼ確実にあり、約３年に一度発生する中気のない月を閏月としたのです。

定気法では太陽の速度が楕円軌道のため位置により異なりますので、図示するように節気の角度の間隔は違ってきます。このため、ひどい場合は1ケ月に中気が２つも入るというようなことも有りうるのです。

＜参考＞　図からも想像できるように、春分から秋分までの時間は、秋分から春分までの時間より７日ほど長い。
　　　　　　　（この計算も理論解析／時間計算を参照すれば計算できます）

(5)数値計算による太陽・月運行データ

太陽と地球、または地球と月のように２つの天体のみの場合は２体問題と言い、解析的に運行状態を求めることができますが、太陽・地球・月のように３つ以上の場合は、多体問題と言い、これを解析的に求めることはできません。それでも、月の運動に対する太陽の引力の影響は小さいとして近似的に解析解を得る方法（これを摂動法と言います）がありますが、きわめて複雑です。そこで、数値計算で太陽・地球・月の運動を求めてみました。その一つの結果が以下のデータです。

計算方法詳細は、以下に述べる理論解析／３体問題数値計算を参照してください。実に簡単なプログラムで計算できます。結果は、かなりの精度で月の運行状態を計算できていると思います。是非皆さんも計算してみてください。

太陽・月軌道の黄道面投影図　　　＜初期値：皆既日食状態＞＜表示：1日毎１年間＞ (太陽との距離は100分の1に縮小して描画）	朔望月変化(18年）：約29.3～29.8日　　　 (赤線：実績例・４年分）	近点月変化(18年）：約25.0～28.5日　　　 (赤線：実績例・４年分）

計算データ　vs 観測データ　　計算結果はかなりの精度で実際に一致	近点黄径変化(18年）　近点順行周期　： 8.83年　(=360/0.1116/365.24) （永年摂動　：　θ= 0.1116 d + θ₀ ) (下図２段目は、永年摂動を引いた短周期摂動）	交点黄径変化(18年）　交点逆行周期　： 18.56年　(=360/.0531/365.24) （永年摂動　：　θ= －0.0531 d + θ₀ ) (下図２段目は、永年摂動を引いた短周期摂動）
<計算> <理科年表> 太陽　　①恒星年　　365.26年　　 365.25年　　 ②離心率　　0.165 0.167 月 ①平均距離　38.40万Km 38.44万Km ②離心率　　0.059 0.055 ③傾斜角　　5.00度　　　 5.149度　　　 ①恒星月　　27.32日　　 27.32日　　　　　②朔望月　　29.53日　　 29.53日　　　　　③近点月　　27.56日　　 27.55日　　　　　④交点月　　29.21日　　　29.22日　　　　　　⑤近点順行　 8.83年　　 8.85年　　　　　⑥交点逆行　18.56年　　 18.61年

(5)数値計算でわかる月運行図

太陽に対する地球・月の運行をイメージ的に表すと左下図のように、月は、その軌道面が黄道面に対し歳差運動を行いながら、且つ楕円軌道そのものも回転（長軸回転、近点位置回転）しながら、地球を周回することが理解できると思います。ただし、詳細にみると回転も常に一定方向でなく逆行することもあり、歳差も純粋な歳差ではなく章動という歳差円に沿って波打つ現象も有る（これらを含めて短周期摂動といいます）のです。(右下図の白道法線軌跡の数値計算結果を参照してください。この歳差現象はコマの歳差現象と同じような原理で起こります。）
天球儀で示すと、地球と太陽あるいは月との距離は問題とならないので、それらの軌道は単なる大円として描かれます。

[Ⅱ] 理論解析

(1) 基本物理法則

使用する物理法則は、以下のニュートンの第２法則と、万有引力の経験式のみです。
　　　　　 $\bf F\rm=m\bf a\rm}$
　　　　　 $\bf F\rm=G\frac{m_1m_2}{r^2}\cdot\frac{\bf r\rm}{r}$

(a)ケプラーの法則/２体問題

①基本方程式

右図のように、位置r₁,r₂にある質量ｍ₁,ｍ₂の２つの天体の運動を考えると、万有引力定数をＧとして、
　　　　　 $m_1\frac{d^2\bf r_1\rm}{dt^2}=-Gm_1m_2\frac{\bf r_1\rm - \bf r_2\rm}{r^3}\hspace{80mm}\mbox{where}\hspace{40mm}r=|\bf r\rm_1-\bf r\rm_2|$
　　　　　 $m_2\frac{d^2\bf r_2\rm}{dt^2}=-Gm_1m_2\frac{\bf r_2\rm - \bf r_1\rm}{r^3}$
重心の位置ベクトルをr_G、天体１に対する２の位置ベクトルをrとすると、
　　　　　 $\bf r\rm_G=\frac{m_1\bf r\rm_1+m_2\bf r\rm_2}{(m_1+m_2)}$
　　　　　 $\bf r\rm=\bf r_2\rm - \bf r_1\rm \hspace{40mm}\bf r_1\rm=\bf r\rm_G-\frac{m_2}{m_1+m_2}\bf r\rm \hspace{40mm}\bf r_2\rm=\bf r\rm_G+ \frac{m_1}{m_1+m_2}\bf r\rm$
であるから、冒頭の２式の和と差をとることで、
　　　　　 $\frac{d^2\bf r\rm_G}{dt^2}=\frac{1}{m_1+m_2}(m_1\frac{d^2\bf r_1\rm}{dt^2}+m_2\frac{d^2\bf r_2\rm}{dt^2})=0$
　　　　　 $\frac{d^2\bf r\rm}{dt^2}=-\frac{G(m_1+m_2)}{r^3}\bf r\rm\hspace{100mm}$ ：　基本方程式

すなわち、２体の重心は動かず、２体の相対運動は、上式の第２式で書き表すことができる。

②基本方程式の解

この微分方程式を解くに当たり、M=m₁+m₂、r方向の単位ベクトルをpとしてp=r/rとおくと、この式は次のように書ける。
　　　　　 $\ddot{\bf r\rm}=-\frac{GM}{r^2}\frac{\bf r\rm}{r}=-\frac{GM}{r^2}\bf p\rm$
rとこの式の外積をとると、r×r=0であるから、
　　　　　 $\bf r\rm\times\ddot{\bf r\rm}=-\frac{GM}{r^3}\bf r\rm\times\bf r\rm=0$
積分すると、積分定ベクトルをh=hk(kは単位ベクトル）として、＜参考： d/dt{r×(dr/dt)}=(dr/dt)×(dr/dt)+r×(d²r/dt²)=r×(d²r/dt²)＞
　　　　　 $\bf r\rm\times\dot{\bf r\rm}=\bf h\rm=\mbox{(const.)}\hspace{100mm}$ ：　角運動量定義 $\hspace{40mm}(h=rv\sin\phi))$
すなわち、hは右下図でいえばxy面に垂直なベクトルであって、
　　　　　 $\dot{\bf h\rm}=0\hspace{100mm}$ ：　角運動量保存の法則
　　　　　 $\bf r\rm\cdot{\bf h\rm}=0\hspace{100mm}$ ：　位置ベクトルと角運動量ベクトルの直交性
　　　　　 $\dot{\bf r\rm}\cdot{\bf h\rm}=0\hspace{100mm}$ ：　速度ベクトルと角運動量ベクトルの直交性
保存される角運動量hは、rとv(ie v=dr/dt)の作る平面に垂直であることを示している。
微小時間dｔにrの掃く面積速度dSは、dS = (1/2)r・(rdθ)= (1/2)r・(v dt sin(φ))=(1/2)｜r×v｜dt=(1/2)｜h｜dtであるから、
　　　　　 $h=2\cdot(dS/dt)\hspace{100mm}$ ：　角運動量は面積速度の２倍
すなわち、２体問題の相対運動は、角運動量ベクトルに垂直な平面内を運動し、その面積速度は一定であると言える。これがケプラーの第２法則である。

同様に微分方程式とhとの外積をとると、
　　　　　 $\ddot{\bf r\rm}\times\bf h\rm=-\frac{GM}{r^3}\bf r\rm\times\bf h\rm$
ところが、r=rp、p×p=0を考慮すると、
　　　　　 $\bf h\rm=\bf r\rm\times\dot{\bf r\rm}=(r\bf p\rm)\times(\dot{r}\bf p\rm+r\dot{\bf p\rm})=r\bf p\rm\times r\dot{\bf p\rm}=r^2\bf p\rm\times\dot{\bf p\rm}$
従い、公式A×(B×C)=B(A・C)-C(B・A)、p²=1よりp・(dp/dt)=0、p・p=1を考慮すると、
　　　　　 $\bf r\rm\times\bf h\rm=(r\bf p\rm)\times(r^2\bf p\rm\times\dot{\bf p\rm})=-r^3\bf p\rm\times(\bf p\rm\times\dot{\bf p\rm})=-r^3\dot{\bf p\rm}$
であるから、次のように変形できる。＜参考： (d²r/dt²)×h=d/dt{(dr/dt)×h})＞
　　　　　 $\frac{d}{dt}(\dot{\bf r\rm}\times\bf h\rm)=GM\dot{\bf p\rm}$
積分すると、積分定ベクトルをQ=eq(qは単位ベクトル）として、
　　　　　 $\dot{\bf r\rm}\times\bf h\rm=GM(\bf p\rm+\bf Q\rm)=GM(\bf p\rm+e\bf q\rm)\hspace{100mm}$ ・・・・・・　(*)
ここで、もう一度、rとこの式との内積をとると、
　　　　　 $\bf r\rm\cdot(\dot{\bf r\rm}\times\bf h\rm)=GM\bf r\rm\cdot(\bf p\rm+e\bf q\rm)$
またまた、公式A・(B×C)=B・(C×A)=C・(A×B)、r×h=0、r×(dr/dt)=hを考慮すると、
　　　　　 $\bf r\rm\cdot(\dot{\bf r\rm}\times\bf h\rm)=\bf h\rm\cdot(\bf r\rm\times\dot{\bf r\rm})=\bf h\rm\cdot\bf h\rm=h^2$
rとQとのなす角をθとすると、p・q=cosθであるから、式(*)の右辺は、
　　　　　 $GM\bf r\rm\cdot(\bf p\rm+e\bf q\rm)=GM(r\bf p\rm)\cdot(\bf p\rm+e\bf q\rm)=rGM(1+e\cos\theta)$
結論として、微分方程式の解は、
　　　　　 $r=\frac{h^2}{GM}\cdot\frac{1}{1+e\cos\theta}=\frac{L}{1+e\cos\theta}\hspace{100mm}( L=\frac{h^2}{GM} )\hspace{100mm}$ ：　楕円軌道

すなわち、２体問題の相対運動は、一方の天体を焦点とする楕円軌道で表される。これがケプラーの第１法則である。

楕円の長半径 a は、θ=0のときのｒとθ=πのときのｒとの和の半分であるから、
　　　　　 $a=\{\frac{L}{1+e}+\frac{L}{1-e}\}/2=\frac{L}{1-e^2}\hspace{100mm}$ ：　長半径
短半径 b は、dy*/dx*=0となるθを求めて、
　　　　　 $b=a\sqr{1-e^2}\hspace{100mm}\hspace{100mm}$ ：　短半径

近点(r_min)は、θ=0のときであるから、qは原点から近点へのベクトルである。上の式(*)より、
　　　　　 $e\bf q\rm=\frac{1}{GM}(\dot{\bf r\rm}\times\bf h\rm)-\bf p\rm$
従い、ｘ軸からの近点への角ξを具体的に求めるには、x・ｙ・ｚ軸方向単位ベクトルをi,j,k、またr(x,y)、v(v_x,v_y)として、
　　　　　 $\dot{\bf r\rm}\times\bf h\rm=(v_x\bf i\rm+v_y\bf j\rm)\times(h\bf k\rm)=hv_y\bf i\rm-hv_x\bf j\rm$
　　　　　 $\bf p\rm}=(x/r)\bf i\rm+(y/r)\bf j\rm$
であるから、q=Ai+Bjと整理してtanξ＝B/Aを計算すると、
　　　　　 $\tan\xi=\frac{-(\frac{hv_x}{GM}+\frac{y}{r})}{\frac{hv_y}{GM}-\frac{x}{r}}\hspace{100mm}$ ：　近点のｘ軸からの角

離心率 e は、eq=(1/GM)(dr/dt)×h-pであったから、
　　　　　 $e^2=(e\bf q\rm)^2=(1/GM)^2(\dot{\bf r}\times\bf h\rm)^2-2(1/GM)(\dot{\bf r}\times\bf h\rm)\cdot{\bf p\rm}+(\bf p\rm)^2$
ところが、公式A・(B×C)=B・(C×A)=C・(A×B)、A×(B×C)=B(A・C)-C(B・A)および直交性(dr/dt)・h=0より、
　　　　　 $(\dot{\bf r\rm}\times\bf h\rm)\cdot{\bf p\rm}=(1/r)(\dot{\bf r\rm}\times\bf h\rm)\cdot{\bf r\rm}=(1/r)\bf h\rm\cdot(\bf r\rm\times\dot{\bf r\rm})=(1/r)\bf h\rm\cdot\bf h\rm=h^2/r$
　　　　　 $(\dot{\bf r\rm}\times\bf h\rm)^2=(\dot{\bf r\rm}\times\bf h\rm)\cdot(\dot{\bf r\rm}\times\bf h\rm)=\dot{\bf r\rm}\cdot\{\bf h\rm\times(\dot{\bf r\rm}\times\bf h\rm)\}$
　　　　 $\hspace{70mm}=\dot{\bf r\rm}\{ (\dot{\bf r\rm)\cdot(\bf h\rm\cdot\bf h\rm)-\bf h\rm\cdot(\dot{\bf r\rm}\cdot\bf h\rm) \}$
　　　　 $\hspace{70mm}=(\dot{\bf r\rm})^2(\bf h\rm)^2=v^2h^2$
であるから、この２式とp²=1を代入して、
　　　　　 $e=\sqr{(\frac{h}{GM})^2(v^2-2\frac{GM}{r})+1}\hspace{100mm}$ ：　離心率

初期条件(r₀、v₀、φ₀)による軌道の形は、以下のようになる。＜近点角の方向ξは、前出の式で、r,vをr₀、v₀、φ₀に置き換えれば得られる。＞
　　　　　 $k=v_0^2-2GM/r_0$
　　　　　 $\hspace{100mm}k<0\hspace{40mm}\mbox{ie}\hspace{40mm}e<1\hspace{100mm}$ ：　楕円軌道
　　　　　 $\hspace{100mm}k=0\hspace{30mm}\mbox{ie}\hspace{40mm}e=1\hspace{100mm}$ ：　放物線軌道
　　　　　 $\hspace{100mm}k>0\hspace{40mm}\mbox{ie}\hspace{40mm}e>1\hspace{100mm}$ ：　双曲線軌道

以上に述べたように、初期条件r₀、v₀が決まれば、離心率 e・長半径 a・近点離角（近日点引数）ω（前述の図でｘ軸上に昇交点がある場合のξ）の軌道要素が決まる。
多体問題では、その瞬間瞬間のr、vを初期値と考えて得られる仮想的軌道の軌道要素を接触軌道要素といい、理科年表等では、こうした値をベースとしている。

(b)軌道上の位置と時間

①周期

２体問題の軌道関係式は、以下の通りである。
　　　　　 $\bf r\rm=r\cos\theta\bf i\rm+r\sin\theta\bf j\rm$
　　　　　 $r=\frac{L}{1+e\cos\theta}\hspace{80mm}\mbox{where}\hspace{40mm}L=a(1-e^2)$
従い、速度ｖは、
　　　　　 $\bf v\rm=\frac{d\bf r\rm}{dt}=\frac{L}{(1+e\cos\theta)^2}\{-\sin\theta\bf i\rm+(e+\cos\theta)\bf j\rm \}\frac{d\theta}{dt}$
角運動量h=r×vであるから、
　　　　　 $\bf h\rm=\frac{L^2}{(1+e\cos\theta)^2}\frac{d\theta}{dt}\bf k\rm$
h=hk、L=h²/GMであるから、
　　　　　 $\frac{d\theta}{dt}=\frac{h}{L^2}(1+e\cos\theta)^2=(\frac{GM}{L^3})^{1/2}(1+e\cos\theta)^2=K(1+e\cos\theta)^2$
積分して、周期Ｔは、
　　　　　 $T=(1/K)\int_0^{2\pi}\frac{1}{(1+e\cos\theta)^2}d\theta$
積分は少々面倒であるが、x=tan(θ/2)、s={(1+e)/(1-e)}^1/2とおいて変形・整理すると、
　　　　　 $\int\frac{1}{(1+e\cos\theta)^2}d\theta=\frac{2}{(1-e)^3}\{-2e\int\frac{1}{(s^2+x^2)^2}dx+(1-e)\int\frac{1}{s^2+x^2}}dx \}$
これに、以下の積分公式
　　　　　 $\int_0^\infty\frac{1}{(s^2+x^2)^{n+1}}dx=\frac{(2n-1)!!}{(2n)!!}\cdot\frac{\pi}{2s^{2n+1}}$
を援用して、θ=0～2πで積分すると、
　　　　　 $\int_0^{2\pi}\frac{1}{(1+e\cos\theta)^2}d\theta=2\int_0^{\pi}\frac{1}{(1+e\cos\theta)^2}d\theta=\frac{2\pi}{(1-e^2)^{3/2}$
であるから、
　　　　　 $T=(\frac{L^3}{GM})^{1/2}\cdot\frac{2\pi}{(1-e^2)^{3/2}}$
L=a(1-e²）を考慮すると、結論として、
　　　　　 $T=2\pi(\frac{a^3}{GM})^{1/2}$
　　　　　 $\frac{T^2}{a^3}=\frac{(2\pi)^2}{GM}\hspace{100mm}$ ：　周期 T² ∝　長半径 a³
これが、ケプラーの第３法則である。

②ケプラー方程式／位置と時間

周期を求めた式で、被積分関数をθでテーラー展開して積分を行えば、ｔ= f(θ）として時間が求まるが、このｆ自身が級数となっており、数値計算でも、逆関数 θ = g(t) を得ることは容易ではない。

そこで、右図の関係にある角 u を考えると、θとの角関係は以下の如くになる。
＜楕円は円をｙ軸方向に(1-e²)^1/2倍したものであることに注意＞

　　　　　 $x=a\cos u\hspace{100mm}\hspace{100mm}y=a\sqr{1-e^2}\sin u$
　　　　　 $x^*=r\cos\theta=a(\cos u-e)\hspace{100mm}y^*=r\sin\theta=a\sqr{1-e^2}\sin u$

この式より、
　　　　　 $r=\sqr{x^{*2}+y^{*2}}=a(1-e\cos u)$
　　　　　 $\tan\theta=\frac{y^*}{x^*}=\frac{\sqr{1-e^2}\sin u}{\cos u-e}\hspace{40mm}\mbox{or}\hspace{40mm}\tan\frac{\theta}{2}=\sqr{\frac{1+e}{1-e}}\tan \frac{u}{2}}\hspace{100mm}$ : θ－u 関係

ここに導入した角 u を用いると、面積速度は一定(k)であるので S=ktと書け、１周期Ｔでは楕円面積 S=πab = πa²{(1-e²)}^1/2となるから、
　　　　　 $k=\frac{\pi a^2\sqr{1-e^2}}{T}=\{\frac{a^2\sqr{1-e^2}}{2}\}\cdot\{\frac{2\pi}{T}\}\equiv\{\frac{a^2\sqr{1-e^2}}{2}\}n\hspace{100mm}$ : n=平均運動
θ－uの関係式を微分して整理すると、r(dθ/dt)=a(1-e²)^1/2(du/dt)が得られる。
これと、r=a(1-e cos(u))の２式を、面積速度k=dS/dt=(1/2)r²(dθ/dt)=(1/2)r{(r(dθ/dt)}に代入して整理すると、
　　　　　 $(1-e\cos u)\frac{du}{dt}=n$
積分して、次のケプラー方程式を得る。（初期値は t=t₀で、近点位置 ie u=0である）
　　　　　 $u-e\sin u=n(t-t_0)\hspace{80mm}\mbox{where}\hspace{40mm}n=\frac{2\pi}{T}=\sqr{\frac{GM}{a^3}}\hspace{100mm}$ ：　ケプラー方程式

ここに現れる l=u-e sin(u)= n(t-t₀)を平均近点離角といい、uを離心近点離角・θを真近点離角という。

③時間計算／数値計算法

角θを与えて時間 t を得るには、θ－u 関係式から u を求め、ケプラー方程式に代入すればよい。

時間 t 与えて角θを得るには、先ず、ケプラー方程式をニュートン法(以下の２式）を用いて u_iの収斂計算をしなければならない。
　　　　　 $f=u-e\sin u-n(t-t_0)$
　　　　　 $u_{i+1}=u_i-\frac{f(u_i)}{df(u_i)/dt}$
ここで得られた u をθ－u 関係式に代入して、θを得る。

＜参考＞楕円軌道を動く太陽・月の朔望月の計算法

　　２つの楕円軌道の長軸のなす角は、ケプラー方程式中の平均近点離角 nt₀で与えることができる。　　　　　　　　　すなわち太陽の近点をｘ軸方向に固定し、月の近点を角αだけ回転させると、αをuに変換して nt₀=u-e sin(u) (n:月の平均運動）として与えられる。太陽については　t₀=0である。　　　　　　　　　この条件の下、時々刻々上記の方法で２つのθを求め、その差Δθ＝α＋2mπ(m:整数）(ie 地球－月－太陽一直線）となる時刻を記録すれば、朔望月の変化が計算できる。

(c)太陽・地球・月の運動（３体問題）の数値計算

太陽（質量 m₁,位置 r₁)、地球（質量 m₂,位置 r₂)、月（質量 m₃,位置 r₃）とすると、各天体の運動方程式は、
　　　　　＜太陽＞　　 $m_1\frac{d^2\bf r_1\rm}{dt^2}=-Gm_1m_2\frac{\bf r_1\rm - \bf r_2\rm}{{r_{12}}^3}-Gm_1m_3\frac{\bf r_1\rm - \bf r_3\rm}{{r_{13}}^3}\equiv\bf F\rm_{1}$
　　　　　＜地球＞　　 $m_2\frac{d^2\bf r_2\rm}{dt^2}=-Gm_2m_1\frac{\bf r_2\rm - \bf r_1\rm}{{r_{21}}^3}-Gm_2_m_3\frac{\bf r_2\rm - \bf r_3\rm}{{r_{23}}^3}\equiv\bf F\rm_{2}$
　　　　　＜月　＞　　 $m_3\frac{d^2\bf r_3\rm}{dt^2}=-Gm_3m_1\frac{\bf r_3\rm - \bf r_1\rm}{{r_{31}}^3}-Gm_3m_2\frac{\bf r_3\rm - \bf r_2\rm}{{r_{32}}^3}\equiv\bf F\rm_{3}$
これらの方程式を、差分化して解く。（一般的にはルンゲ・クッタ法がよく用いられるが、簡単で、精度も左程違わないのでこれを採用。精度比較については数値計算法を参照）

一般的には、
　　　　　 $m\frac{d^2x\rm}{dt^2}=F_x$
差分化は、時間刻みΔtに対する、xの過去・現在・未来の値を、x₀・x₁・x₂として、未来値（次のステップの値）は、以下の式で与えられる。
　　　　　 $x_2=2x_1-x_0+(F_x/m)(\Delta t)^2\hspace{80mm}v_1=(x_2-x_0)/(2\Delta t)$
同様に y座標、z座標についても解くと、r(t)(x(t),y(t),z(t))が得られる。このようにして、３つの天体について位置を求めれば３体問題としての運動が得られる。

この結果から各種情報を取り出せば、運動が理解できることになる。





＜参考＞[Ⅰ]旧暦で示した計算例条件と解析法

簡単のため、初期状態は皆既日食の状態とした。地球は原点にあり、太陽・月は(x,y)面内ｘ軸上にあって（朔月）、太陽の速度は(x,y)面内、月はｘ軸上で
昇交点、すなわち速度は(x,y)面と傾斜角を有するものとした。初期値は理科年表に示されている平均値を基本としつつ観測値と整合性のとれる条件として、


　　　　(1)太陽　：　初速=地球公転速度×1.0008、　　　　　　　　　　　　　　　　 　　位置ベクトル・速度ベクトル角=(90-0.95)度        
　　　　(2)月　　：　初速=月公転速度×1.003、  　　 初期距離=月公転半径×1.003、     位置ベクトル・速度ベクトル角=(90-2.7)度        

太陽の軌道面すなわち黄道の変化は小さいので、太陽の初期値である位置ベクトルr,速度ベクトルvの決定する(x,y)平面が不動すなわち黄道は(x,y)面
であり不動のものとして考える。更に、地球の位置は動くが地球を中心として考えることとし、太陽・月の位置は地球の位置を差し引いた値に変換して考える。
この仮定のもとに月に関する以下の情報を取り出す。



　　　太陽（X,Y,Z) 　月（x,y,z)　地球（0,0,0)

        (1)月軌道　　　：　  (x,y)座標そのもの           (月の軌道の黄道への投影）
        (2)恒星月　　　：　　ｘ>0 & y=0　　　　　        (初期状態をｘ軸方向とした為）　　
        (2)朔望月　　　：　　tan^-1(y/x)=tan^-1(Y/X)        (黄経同一）
        (3)交点月　　　：　  z=0　　　　　　　　　　　   (黄道=(x,y)面）
        (4)近点月　　　：　　r=r_min for １周、　r²=x²+y²+z² 
      
        (5)軌道要素（離心率 e・長半径 a・傾斜角ε）
　
　　　　　　　　　　近似的には、隣り合う遠点月と近点月の位置から、次の計算で求めることができる。
　　　　　　　　　　　　　　 e=1-2r_min/(r_max+r_min)
                             a=(r_maxn+r_min)/2
                             ε={tan^-1(y/x)遠点月-tan^-1(y/x)近点月}/2　or ε=tan^-1(v_z/(v_x^2+v_y^2)^(1/2))

                    理科年表データと比較するには、すでに述べた接触軌道要素としての計算が必要である。問題は傾斜角であるが、
　　　　　　　　　　黄道を(x,y)平面と考えたのでその法線n=(0,0,1)、月の接触軌道の法線は現在値h=r×vであるから、
　　　　　　　　　　　　　　ε=π/2-β　　　但し　　cosβ＝n・h/h
                    として得ることができる。

(d)多体問題の摂動論

例えば、地球に対する月の運動を考える場合、月・地球間の引力のみでは、２体問題として解析解が得られるが、月・太陽間の引力も無視し得ない。この微小な影響による２体問題の解からの差を計算する手法を摂動論と言うが、簡単には論じられないので割愛するが、その結果は、平均軌道軌道要素（全体摂動から短周期項を除いたもの）のみについて言うと、日数ｄの多項式（どこまでの微少量を考えるかで１次式にも２次・３次・４次式にもなります）として表される。以下の例は、理科年表等でよく見かける式ですが、このような摂動論がベースになっています。

        月 　	: 　近点離角　　ω= 83.3532 + 0.11140352d－0.0103d²　

　　　　　　　　    昇交点経度　Ω= 125.0446－0.05295376d + 0.0021d²

この式の１次の係数は[Ⅰ]節に図示した数値計算による近点月・交点月の永年摂動の直線係数そのものとなっています。摂動項は、太陽の影響だけではなく、地球の扁平楕円体の影響(これをJ₂項という）・惑星の影響もあり、正確にはこれらも取り込まなければなりません。ニューカムは、４００の摂動項を用いて解析したが、最新のものは４７４４の摂動項を取り込んでいるそうである。

[トップへ戻る]