球内に閉じ込められた 30 個の電荷

立体視できます。

• 面倒だったので、加速度ではなく速度が力に比例するようにしました。電流が電圧に比例している様子をイメェジしてください。たぶん抵抗が大きいのでしょう。
• 球の表面での反発係数は 0 です。面倒だったので。
• 途中から時間ステップが 20 倍になります。一旦ほぼ静止してからまた動き出すところです。

%!PS
%%BoundingBox: 0 0 200 100
/N 30 def
/R 47 def
/surface? false def
/dtsmall 5 def /dtlarge 100 def
/eps 0.01 def
/dotsize 3 def
/thetaL 90 def /phiL -3 def
/thetaR 90 def /phiR 3 def

/charges [ N {3 array} repeat ] def
0 srand

/force {
    4 dict begin
        /charges exch def
        [
            0 1 N 1 sub {
                /i exch def
                /xyz charges i get def
                [0 0 0]
                0 1 N 1 sub {
                    /j exch def
                    i j ne {
                        xyz charges j get arraysub
                        dup dup prod eps add 1 exch div
                        exch arraymul
                        arrayadd
                    } if
                } for
            } for
        ]
    end
} def
/normalize {
    4 dict begin
        /points exch def
        0 1 N 1 sub {
            /i exch def
            /xyz points i get def
            /r xyz dup prod sqrt def
            surface? {
                R r div xyz arraymul
                points exch i exch put
            }{
                r R gt {
                    R r div xyz arraymul
                    points exch i exch put
                }if
            }ifelse
        } for
    end
} def
/prod {
    3 dict begin
        /b exch def
        /a exch def
        0
        0 1 a length 1 sub {
            /i exch def
            a i get b i get mul add
        } for
    end
} def
/arraymul {
    1 dict begin
        exch
        /c exch def
        [exch {c mul} forall]
    end
} def
/arrayadd {
    3 dict begin
        /b exch def
        /a exch def
        [
            0 1 a length 1 sub {
                /i exch def
                a i get b i get add
            } for
        ]
    end
} def
/arraysub {
    3 dict begin
        /b exch def
        /a exch def
        [
            0 1 a length 1 sub {
                /i exch def
                a i get b i get sub
            } for
        ]
    end
} def
/draw {
    4 dict begin
        /phi exch def
        /theta exch def
        /points exch def
        /rotation [
            [theta cos phi cos mul  theta cos phi sin mul  theta sin neg]
            [phi sin neg  phi cos  0]
            [theta sin phi cos mul  theta sin phi sin mul  theta cos]
        ] def
        0 0 R 0 360 arc stroke
        0 1 N 1 sub {
            /i exch def
            /xyz points i get def
            rotation 1 get xyz prod
            rotation 0 get xyz prod neg
            rotation 2 get xyz prod R div 3 add dotsize mul 0.25 mul
            0 360 arc fill
        } for
    end
} def
/chargecenter {
    1 N div
    [0 0 0]
    0 1 N 1 sub {
        charges exch get arrayadd
    } for
    arraymul
} def

% 乱数で初期位置を決める
0 1 N 1 sub {
    /i exch def
    charges i get
    0 1 2 {
        /j exch def
        dup j rand 2 30 exp div 1 sub R mul put
%        dup j rand 2 31 exp div R mul put
    } for
    pop
} for
charges normalize
/center chargecenter def

0 1 300 {
    150 lt {dtsmall} {dtlarge} ifelse /dt exch def
    R dotsize add R dotsize add translate
    charges thetaL phiL draw
    R dotsize add 2 mul 0 translate
    charges thetaR phiR draw showpage

    /oldcenter center def
    /displacement charges force def
    0 1 N 1 sub {
        /i exch def
        charges i get dt displacement i get arraymul arrayadd
        charges exch i exch put
    } for
    charges normalize
    /center chargecenter def
    center oldcenter arraysub dup prod sqrt 1e-8 le {
        exit
    } if
} for

0 1 N 1 sub {
    /i exch def
    charges i get dup prod sqrt 0.999 R mul lt {
        charges i get ==
    } if
} for