手書き PostScript で円制限三体問題
円制限三体問題を PostScript で数値計算してみた。
円制限三体問題というのは、重力によって 2 個の質点が円軌道で公転している中に、もう 1 個他の質点に及ぼす重力が無視できるくらい小さい質点を置いた時、どう運動するか、という問題です。ラグランジュポイントはこの問題を解くと出てきます。
以下ソース
%!PS %%BoundingBox: 0 0 150 150 % 円制限三体問題 % 4次のシンプレクティック法 /pi 3.141592653589793238 def /m1 10 def /m2 1 def /a 50 def %軌道長半径 /T 20 def %公転周期 /dt 0.1 def %時間きざみ % G*(m1+m2)=a^3*(2pi/T)^2 /G a dup dup mul mul 2 pi mul T div dup mul mul m1 m2 add div def % 初期座標と初期速度 L4 からちょっとずれたところ /t 0 def /x 0.5 a mul m1 m2 sub m1 m2 add div mul def /y a 3 sqrt mul 0.5 mul 1 add def /z 0 def /vx 2 pi mul T div y mul neg def /vy 2 pi mul T div x mul def /vz 0 def /position [t x y z vx vy vz] def % 大きさ /R1 7 def /R2 5 def /R3 3 def % 色 /color1 {1 0 0 setrgbcolor} def /color2 {0 1 0 setrgbcolor} def /color3 {0 0 1 setrgbcolor} def /distance { % x y z x1 y1 z1 3 -1 roll 6 -1 roll sub dup mul % y z y1 z1 (x-x1)^2 3 -1 roll 5 -1 roll sub dup mul % z z1 (x-x1)^2 (y-y1)^2 4 2 roll sub dup mul add add sqrt } def /m1m2position { 2 dict begin /t exch def /phase 360 t mul T div def a m2 mul m1 m2 add div phase cos mul neg %x1 a m2 mul m1 m2 add div phase sin mul neg %y1 0 %z1 a m1 mul m1 m2 add div phase cos mul %x2 a m1 mul m1 m2 add div phase sin mul %y2 0 %z2 end } def /force { 12 dict begin position 0 get /t exch def position 1 get /x exch def position 2 get /y exch def position 3 get /z exch def t m1m2position /z2 exch def /y2 exch def /x2 exch def /z1 exch def /y1 exch def /x1 exch def /r1cube x y z x1 y1 z1 distance dup dup mul mul def /r2cube x y z x2 y2 z2 distance dup dup mul mul def m1 x x1 sub mul r1cube div m2 x x2 sub mul r2cube div add G neg mul %Fx m1 y y1 sub mul r1cube div m2 y y2 sub mul r2cube div add G neg mul %Fy m1 z z1 sub mul r1cube div m2 z z2 sub mul r2cube div add G neg mul %Fz end } def /movex { 1 dict begin /dt exch def position 1 position 1 get position 4 get dt mul add put position 2 position 2 get position 5 get dt mul add put position 3 position 3 get position 6 get dt mul add put end } def /movev { 4 dict begin /dt exch def force /Fz exch def /Fy exch def /Fx exch def position 4 position 4 get Fx dt mul add put position 5 position 5 get Fy dt mul add put position 6 position 6 get Fz dt mul add put end } def /symplectic2 { 1 dict begin /dt exch def 0.5 dt mul movex 0.5 dt mul movev position 0 position 0 get dt add put 0.5 dt mul movev 0.5 dt mul movex end } def /sm 1 2 2 1 3 div exp sub div def /symplectic4 { 1 dict begin /dt exch def sm dt mul symplectic2 1 2 sm mul sub dt mul symplectic2 sm dt mul symplectic2 end } def /drawplanets { 6 dict begin position 0 get m1m2position /z2 exch def /y2 exch def /x2 exch def /z1 exch def /y1 exch def /x1 exch def color1 x1 y1 R1 0 360 arc fill color2 x2 y2 R2 0 360 arc fill color3 position 1 get position 2 get R3 0 360 arc fill end } def /i 0 def 400 { dt symplectic4 i 2 mod 0 eq { 75 75 translate drawplanets showpage } if /i i 1 add def } repeat