手書き Postscript その2

二重振り子
ソースは

%!PS
%%BoundingBox: 0 0 220 230
gsave
% 振り子の支点
/cx 110 def /cy 220 def
% 棒の長さ
/l1 100 def /l2 100 def
% 初期角度 radian
/theta1 0.5 def /theta2 0.3 def
% 初期角速度
/dthdt1 0 def /dthdt2 0 def
% 重りの質量
/m1 10 def /m2 10 def
% 重力
/g 30 def
% 時間ステップ
/dt 0.3 def

/rsin {180 mul 3.1415926535897932384626 div sin} def
/rcos {180 mul 3.1415926535897932384626 div cos} def
/ft1 {4 2 roll pop pop pop} def
/ft2 {4 1 roll pop pop pop} def
/fdt1 {
    8 dict begin
    /dthdt2 exch def /dthdt1 exch def
    /theta2 exch def /theta1 exch def
    /diffthsin theta1 theta2 sub rsin def
    /diffthcos theta1 theta2 sub rcos def
    /denom diffthsin dup mul m2 mul m1 add l1 mul def
    /numer m1 m2 add theta1 rsin mul neg
        m2 theta2 rsin mul diffthcos mul add g mul
        l2 dthdt2 dup mul mul  l1 dthdt1 dup mul mul diffthcos mul add
        m2 mul diffthsin mul sub def
    numer denom div
    end
} def
/fdt2 {
    8 dict begin
    /dthdt2 exch def /dthdt1 exch def
    /theta2 exch def /theta1 exch def
    /diffthsin theta1 theta2 sub rsin def
    /diffthcos theta1 theta2 sub rcos def
    /denom diffthsin dup mul m2 mul m1 add l2 mul def
    /numer m1 m2 add g mul theta1 rcos mul
        m1 m2 add l1 mul dthdt1 dup mul mul add
        m2 l2 mul dthdt2 dup mul mul diffthcos mul add diffthsin mul def
    numer denom div
    end
} def
148 {
    /rth1 theta1 theta2 dthdt1 dthdt2 ft1 dt mul def
    /rth2 theta1 theta2 dthdt1 dthdt2 ft2 dt mul def
    /rdth1 theta1 theta2 dthdt1 dthdt2 fdt1 dt mul def
    /rdth2 theta1 theta2 dthdt1 dthdt2 fdt2 dt mul def
    /dth1 theta1 0.5 rth1 mul add theta2 0.5 rth2 mul add dthdt1 0.5 rdth1 mul add dthdt2 0.5 rdth2 mul add ft1 dt mul def
    /dth2 theta1 0.5 rth1 mul add theta2 0.5 rth2 mul add dthdt1 0.5 rdth1 mul add dthdt2 0.5 rdth2 mul add ft2 dt mul def
    /ddth1 theta1 0.5 rth1 mul add theta2 0.5 rth2 mul add dthdt1 0.5 rdth1 mul add dthdt2 0.5 rdth2 mul add fdt1 dt mul def
    /ddth2 theta1 0.5 rth1 mul add theta2 0.5 rth2 mul add dthdt1 0.5 rdth1 mul add dthdt2 0.5 rdth2 mul add fdt2 dt mul def
    /theta1 theta1 dth1 add def
    /theta2 theta2 dth2 add def
    /dthdt1 dthdt1 ddth1 add def
    /dthdt2 dthdt2 ddth2 add def
    0 0 0 setrgbcolor
    cx cy moveto
    l1 theta1 rsin mul l1 theta1 rcos mul neg rlineto
    l2 theta2 rsin mul l2 theta2 rcos mul neg rlineto stroke
    0 1 0 setrgbcolor
    cx l1 theta1 rsin mul add cy l1 theta1 rcos mul sub 10 0 360 arc fill
    0 0 1 setrgbcolor
    cx l1 theta1 rsin mul add l2 theta2 rsin mul add
    cy l1 theta1 rcos mul sub l2 theta2 rcos mul sub 10 0 360 arc fill
    showpage
} repeat
grestore

これを ImageMagick で gif に変換
もっと簡潔に書けるはず。