手書き 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