ローレンツアトラクタ

以下ソース

%PS!
gsave

% パラメータ
% dx/dt=-px+py
% dy/dt=-xz+rx-y
% dz/dt=xy-bz
/p 10 def /r 28 def /b 8 3 div def
% 時間ステップ
/dt 0.003 def
% 繰り返し回数と何回毎に描くか
/skip 25 def
/times skip 360 mul def
% 初期位置
/initx 10 def /inity 0 def /initz 40 def
% 射影の方向
/theta 60 def /phi 30 def
% 一回描く毎の射影方向の回転角
/dphi 1 def
% 図の位置と倍率
/ctm 70 8 translate 2 2 scale 6 array currentmatrix def
/lwidth 0.5 def

/fx{
    % x y z -> (y-x)*p
    pop exch sub p mul
}bind def
/fy{
    % x y z -> x*(r-z)-y
    r exch sub 3 2 roll mul exch sub
}bind def
/fz{
    % x y z -> x*y-b*z
    b mul 3 1 roll mul exch sub
}bind def
/step{
    % x y z -> x+dx y+dy z+dz
    3 dict begin
        /z exch def /y exch def /x exch def
        x x y z fx dt mul add
        y x y z fy dt mul add
        z x y z fz dt mul add
    end
}bind def
/addstep{
    % arr -> newarr
    dup length 1 add array dup dup 4 2 roll copy length
    [ 3 copy pop 1 sub get aload pop step ] put
}bind def
/addstep1{
    % 別バージョン
    [
        exch aload pop
        dup [ exch aload pop step ]
    ]
}bind def
/projection{
    % x y z -> Px Py
    3 dict begin
        /z exch def /y exch def /x exch def
%        y phi cos mul x theta sin mul z theta cos mul add phi sin mul sub
%        z theta sin mul x theta cos mul sub
        y phi cos mul x phi sin mul sub
        z theta sin mul x phi cos mul y phi sin mul add theta cos mul sub
    end
}bind def

[[initx inity initz]]
0 1 times{
    dup skip mod 0 eq{
        times div 1.0 0.75 sethsbcolor
        ctm setmatrix
        lwidth setlinewidth

        initx inity initz projection moveto
        dup{
            aload pop projection
            lineto
        }forall
        stroke

        dup dup length 1 sub get aload pop projection
        0 0 0 setrgbcolor
        lwidth 2 mul 0 360 arc fill

        showpage
        /phi phi dphi add def
    }{
        pop
    }ifelse
    addstep
}bind for

grestore

PostScript から gif への変換は

convert -delay 2 -page 135x115 lorenz.ps lorenz.gif