重力プログラミング入門 第6回 ～理論としての物理・数学をプログラミングして、動かし、世界を楽しむ～ 惑星の位置・軌道を計算する 別名： 「2D天体シミュレータをグルグル回す」 ～単純なケプラー軌道から「惑星軌道の永年変化 (VSOP)」へ～ 2018/07/18 ver 0.5作成 2018/07/20 ver 0.9作成

1 1. そもそもブラックホールだけが関心事だったが… 2. 天体シミュレーション基礎：惑星の軌道計算 3. 復習：ケプラーの楕円軌道 4. ケプラーの楕円軌道の問題 5. 「VSOP82」から「VSOP87」へ 目次

2 １．そもそもブラックホールだけが関心事だったが…

3 １．そもそもブラックホールだけが関心事だったが… 天体シミュレーションしたい のですが、プログラミングの やり方が分からなくて…

4 自分はブラックホールと 重力しか興味が無くて、

5 天体全体の研究となると、 だいぶノリが違うけど…

6 １．そもそもブラックホールだけが関心事だったが… この本で、「SPH法」という、 個々の粒子の動きを演算する 解析手法を知り…

7 １．そもそもブラックホールだけが関心事だったが… 流体力学と重力物理の関係 を説明する論文を読み…

8 １．そもそもブラックホールだけが関心事だったが…

9

10 「ナビエ-ストークス方程式」という粘性の ある流体の運動を表す偏微分方程式 これを解決すると100万ドルもらえる

11 流体力学とか、ブラックホール とも関連性ありそうなので…

12 ま、ついでにやってみるかｗ

13 ２．天体シミュレーション基礎：惑星の軌道計算

14 ２．天体シミュレーション基礎：惑星の軌道計算

15 ２．天体シミュレーション基礎：惑星の軌道計算 「VSOP87」は、世の中に最も普及した、惑星軌道計算するため のアルゴリズムで、メジャーな天体シミュレータである以下2つでも 採用されています ➢ Celestia (セレスティア) http://orbit.medphys.ucl.ac.uk/ ➢ Orbiter (オービター) http://orbit.medphys.ucl.ac.uk/

16 ３．復習：ケプラーの楕円軌道

17 重力プログラミング入門 第1回がちょうどその内容で

18 ほぼ1年前に登壇しました

2017/09/20 プログラマのための数学勉強会@福岡#6 in LINE Fukuoka 重力プログラミング入門 ～理論としての数学をプログラミングして、動かし、世界を楽しむ～ 再掲

20

21 １．「横に投げたボール」と「人工衛星」は同じ ボールを真横に投げると、途中までは水平に飛んで行きますが、 やがて水平方向の勢いを失って、重力に従って地上方向に落ち ていき、最後は地上に着きます GPSを積んだ人工衛星は、このボールの水平飛行が続いたよう なもので、重力に従って地上方向に落ち続けてはいるが、速度 が維持されているため、地上に着かず、周回し続ける訳です この速度のことを、「第1宇宙速度」と言います = 再掲

22 １．「人工衛星」と「公転する惑星」も実は同じ GPSを積んだ人工衛星が、地球を中心に周回するのと、地球 含めた太陽系惑星が、太陽の周りを公転する理屈は、実は全く 同じなのです 地上に着くことも無く、かといって、どこかに飛び去る訳でも無く、 太陽や地球を中心に、楕円を描き続ける特徴も同じです このような公転の軌道は計算可能で、「万有引力の法則」と、 その元である「ケプラーの第3法則」から導き出すことができます = 再掲

23 ２．「万有引力の法則」 あらゆる2つの物体間には、お互いを引き合う力、つまり「引力」 が発生しており、引力は、以下2つの条件が成り立ちます ➢ 各物体の質量の乗算に比例する (重い程強く引き合う) ➢ 物体間の距離の2乗に反比例する (離れると弱まる) 質量当りにかかる引力の強さを「G 」※、物体の質量を「M 」「m」、 物体間の距離を「r 」とすると、物体間の引力は、以下となります ※「万有引力定数」と呼びます 【Column】 「重力」とは、引力に「自転の遠心力」を合わせた力を指します 本スライドでは、「自転の遠心力」の考慮を割愛し、この2つを同一視します F = G Mm r 2 再掲

24 ２．人工衛星の位置を示すxy軸座標 以下は、地球の中心を0地点とした、人工衛星のxy 座標です 人工衛星の軌道は、「運動方程式」により、以下の通りです ma = F 加速度a を、x 軸・y 軸に分けることで､移動速度を出しましょう x α θ y F Fcosθ Fsinθ r Position( x, y ) 人工衛星の座標 0 ※m は人工衛星の質量、a は加速度 再掲

25 ２．引力と人工衛星の位置から移動速度を計算 ma はmdx とmdy という偏微分に分割し、F は直交座標(xy 座標)から極座標(方向F と角度θ で表現される座標)へと変換 します mdx = -F cosθ, mdy = -F sinθ F を前述した引力に展開し、cosθ, sinθ を人工衛星位置x, y と人工衛星までの距離r で表現すると、以下の通りです mdx = -G ・ = -G mdy = -G ・ = -G Mm r 2 x r Mmx r 3 Mm r 2 y r Mmy r 3 ※方向F は、0地点から外側に向かって いればプラス、逆方向ならマイナスと なるため、この場合はマイナスとなる 再掲

26 ２．引力と人工衛星の位置から移動速度を計算 両辺のm を削除したものが、移動速度になります dx = -G , dy = -G 地球が0地点なので、r は、x とy から計算できます r = √ x + y G およびM は、以下の通りです ➢ 万有引力定数 G = 6.67259 x 10 m / s ➢ 地球の質量 M = 5.974 x 10 kg GM = 1.267 × 10 km / s Mx r 3 My r 3 2 2 -11 3 2 24 8 3 2 再掲

27 ３．公転軌道をJavaScriptでプログラミング 前章で紐解いた、人工衛星の移動速度を用いて、人工衛星が 地球周回軌道をグルグルと回るシミュレーションをJavaScirptで 作ってみましょう dx = -G , dy = -G , r = √ x + y この移動速度の数式をコード化した以下は、数式そのままです （ 最近、関数型しか書いて無いので、状態を持てる・更新できる言語に拒否反応が(-_-u… ） … // 公転対象の位置から重力影響を計算し、速度を変更後、公転対象の位置を移動 moveOnGravityEffect: function() { var r = Math.sqrt( Math.pow( this.x, 2 ) + Math.pow( this.y, 2 ) ); this.dx = this.dx - G * M * this.x / Math.pow( r, 3 ); // x方向の移動速度 this.dy = this.dy - G * M * this.y / Math.pow( r, 3 ); // y方向の移動速度 this.x = this.x + this.dx; // x位置を、x方向の移動速度で更新 this.y = this.y + this.dy; // y位置を、y方向の移動速度で更新 }, … Mx r 3 My r 3 2 2 再掲

28 ３．人工衛星が地球周回軌道をグルグル回る JSで計算した人工衛星の移動速度 (とxy 座標) を以下画面 で表示して、地球周回軌道をグルグル回る様を見ましょう 再掲

29

30 ４．ケプラーの楕円軌道の問題

31 ここから先は、ちょっとした 歴史のお勉強みたいな感じ

32 天文学は、長時間の探求で 辿り着いた、ザ・歴史の側面

33 ４．ケプラーの楕円軌道の問題 惑星の公転が、ケプラーの楕円軌道に沿っているならば、楕円 軌道の形や向きは、前述した軌道計算の通り、「永遠に不変」と なるはずです しかし実際は、楕円の形や向きが、時間の経過につれて、徐々 に変化していく傾向が観察されています そして、このズレの原因として挙げられるのが、「摂動 (せつどう)」 と呼ばれる現象です

34 ４．ケプラーの楕円軌道の問題 「摂動 (せつどう)」とは、天文学の用語で、ある公転する天体と、 公転の重力源となる母天体による公転軌道運動系に、第3の 物体 (惑星、彗星、太陽フレア等) による重力影響がおよんだ 場合、軌道が乱され、楕円の形や向きにズレが生じる現象を 指します ちょうど、人工衛星が惑星でスイングバイするのと似たように、第3 の物体による重力影響は、加速・減速を引き起こし、ズレます

35 ４．ケプラーの楕円軌道の問題 こうした「摂動」による、ケプラー楕円軌道からのズレを、計算する ためのモデルが、過去、幾つも作成されてきました 1740年のジャック・カッシーニによる観測表が、この分野における 決定打となっており、本日のお話も、この末裔についてとなります (ただし、実際の計算による証明は、1970年代にコンピュータが 発明されるまで、二次項による近似式までに限定されていた) ちなみに、「摂動」を厳密に計算することは、「多体問題」を解く、 ということであり、重力プログラミング入門「第4回：スペースコロ ニーを建造する」でも、お話しましたが、「ラグランジュ点」のような、 限定された特殊解で無いと解けない…という課題があります

36 １．宇宙空間に物体を配置する 地球と月を2惑星として、質量が無視できる小さな第3の物体を スペースコロニーとし、宇宙空間に配置したのが、ガンダムシリーズ における「サイド」で、5点のラグランジュ点に、サイド1～サイド8の 計8つのサイドを配置していました 月の裏側にあるのが ジオン公国発祥の地 となったサイド3 (L2) 月と同一軌道上の サイド7でガンダム 大地に立つ (L3) 再掲

37 ３．ラグランジュ点 L4／L5 【余談】 今回のような、3つの天体が、重力によって、お互いに影響し合う シチュエーションで、それら軌道を決定する問題は、「三体問題」 もしくは「多体問題」と呼ばれ、解析的には解けない …つまり、 方程式の変形で解くことは不可能… と言われています しかし、3体のうち、1体の質量が、他2体と比べ、無視できるほど 小さい場合は、「制限三体問題」と呼ばれ、特殊解を求めること が可能となります 更に、2体の軌道が、円軌道である場合、「円制限三体問題」 と呼ばれます 今回のラグランジュ点は、この「円制限三体問題」に該当します 再掲

38 ４．ケプラーの楕円軌道の問題 そこで、「摂動」そのものを計算するのでは無く、「摂動」の原因と なる第3の物質のうち、主要因となり得る惑星間の相互作用に 焦点を当て、各惑星の「振幅」と「主周期」を使い、「時間」にて 級数展開することで、「時間経過に伴う惑星の位置 (惑星間 の相互作用込み) を直接計算する」という、解析的アプローチ※ を取るように発展したのが、1781年の頃です ※「テイラー展開」や「マクローリン級数」、「微分」は、この解析的 アプローチの代表で、AI・MLを支える基礎でもあります

39 ５．「VSOP82」から「VSOP87」へ

40 ５．「VSOP82」から「VSOP87」へ その後、コンピュータの発展に伴い、実用に耐え得る理論として 完成したのは、1982年に発表の「VSOP82」というアルゴリズム でした 「VSOP82」は、精度として、1,000年～100万年程の先まで 適用できることが分かっていますが、この適用期間・精度の限界 は、「観測にかけた時間」の長さに依存し、短いと、公転周期に 誤りが出てきます また、「VSOP82」は、公転軌道のみの計算しか出せないため、 惑星の位置特定 (と計算時の級数打ち切り) ができません これを解消したのが、「VSOP87」でした

41 ５．「VSOP82」から「VSOP87」へ 「VSOP87」では、水星、金星、地球については、元期2,000 (≒2,000年1月1日 0世界時)の±4,000年、木星と土星は ±2,000年、天王星と海王星は±6,000年の範囲の誤差が、 1秒以内という、非常に高い精度を実現しました また、惑星の位置を直交座標系 (x, y, z) の形で計算すること も可能となりました そして現在、VSOP87用の観測データは、オープンデータとして 提供されています

42 ５．「VSOP82」から「VSOP87」へ 右3列が、下記のような、横3列の時系列行列の観測データで、 これを取り出し、右の式で計算すると、xyz座標が取得できます function earth( t ) { X0 += 0.99982928844*Math.cos( 1.75348568475 + 6283.07584999140 * t ); X0 += 0.00835257300*Math.cos( 1.71034539450 + 12566.15169998280 * t ); X0 += 0.00561144206*Math.cos( 0.00000000000 + 0.00000000000 * t ); … VSOP87 VERSION A1 EARTH VARIABLE 1 (XYZ) *T**0 843 TERMS HELIOCENTRIC DYNAMICAL ECLIPTIC AND EQUINOX J2000 1310 1 0 0 1 0 0 0 0 0 0 0 0 0 -0.00001522262 0.99982928833 0.99982928844 1.75348568475 6283.07584999140 1310 2 0 0 2 0 0 0 0 0 0 0 0 0 0.00814054703 -0.00187001868 0.00835257300 1.71034539450 12566.15169998280 1310 3 0 0 0 0 0 0 0 0 0 0 0 0 0.00000000000 0.00561144206 0.00561144206 0.00000000000 0.00000000000 …

43 ５．「VSOP82」から「VSOP87」へ あとは、各惑星毎のxyz座標を、特定時刻毎に求めれば、惑星 の公転軌道が計算できます 今回開発したシミュレータでは、太陽系8惑星の特定時刻位置 を変化させ、軌道を見れます (太陽系の拡大・縮小も可能)

44 ５．「VSOP82」から「VSOP87」へ 下記の式で計算した座標は、3次元で取得されます… …ということは… function earth( t ) { … var X = X0 + X1 * t + X2 * t2 + X3 * t3 + X4 * t4 + X5 * t5; var Y = Y0 + Y1 * t + Y2 * t2 + Y3 * t3 + Y4 * t4 + Y5 * t5; var Z = Z0 + Z1 * t + Z2 * t2 + Z3 * t3 + Z4 * t4 + Z5 * t5; return [ X, Y, Z ]; }

45 ５．「VSOP82」から「VSOP87」へ CelestiaやOrbiterのような、3D天体シミュレータの可能性が…

46 今後の重力プログラミングに 乞うご期待！