Tcl/Tkで物理シミュレーション

空気抵抗を受ける飛翔体

原点から打ち上げられた飛翔体が、速度に比例する抵抗を受けながら運動する様子をシミュレートするプログラムを書いてみる。

• 数値積分の時間刻みをdt、抵抗の比例係数をalphaとし、入力できるようにしています。原点から60°の角度に初速10m/sで打ち上げられるものとしています。
• moveをチェックすると地面に衝突するまでの軌跡を計算・描画します。1 stepを押すと、1ステップだけ数値積分をして動かします。
• clearで描画した軌跡を全て消去します（なので、これを押さずに、dt等を変更して軌跡の変化を比較したりすることができます）。
• resetで質点は原点に戻ります。
• 初速や打上げ角度などを変更する場合にはコードを直接いじる必要があります。

grvty.tcl

#!/usr/bin/wish
set lx 500
set ly 500
set ifmv 0
#set ifib 1
canvas .c0 -width $lx -height [expr $ly+10]
#menu .m -type menubar
#. configure -menu .m
#.m add cascade -label "File" -under 0 -menu .m.m1
set scl 200
set g 9.8
set alp 0.2
set m 0.1
set dt 0.01
set rad 5
set vxini [expr 10/2]
set vyini [expr 10/2*sqrt(3)]
set xini 0
set yini 0
set  x0 $xini; set  y0 $yini
set vx0 $vxini; set vy0 $vyini
frame .f0
frame .f1
frame .f2
#label .lb2 -text "Initial x"
#entry .e2 -width 10 -textvariable xini
#label .lb3 -text "y"
#entry .e3 -width 10 -textvariable yini
#label .lb0 -text "Vx"
#entry .e0 -width 10 -textvariable vxini
label .lb1 -text "dt"
entry .e1 -width 10 -textvariable dt
label .lb2 -text "alpha"
entry .e2 -width 10 -textvariable alp
bind .e1 <Key> {
    if {%k==36} then {
    reset
    }
}
set  x0f [format "%5.3f"  $x0]
set  y0f [format "%5.3f"  $y0]
set vx0f [format "%5.3f" $vx0]
set vy0f [format "%5.3f" $vy0]
label .lb00 -text "(x,y)"
label .lb01 -width 10 -textvariable x0f
label .lb02 -width 10 -textvariable y0f
label .lb03 -text "(Vx,Vy)"
label .lb04 -width 10 -textvariable vx0f
label .lb05 -width 10 -textvariable vy0f
pack .lb00 .lb01 .lb02 .lb03 .lb04 .lb05 -in .f2 -side left
#pack .lb2 .e2 .lb3 .e3 .lb0 .e0 .lb1 .e1 -in .f1 -side left
pack .lb1 .e1 .lb2 .e2 -in .f1 -side left
checkbutton .cb0 -text "Move" -variable ifmv
#checkbutton .cb1 -text "in Box" -variable ifib
button .b1 -text "Reset" -command "reset"
button .b2 -text "1 step" -command "onestep"
button .b3 -text "Clear" -command {
    .c0 delete all
    set x1 [expr $x0-$rad]
    set x2 [expr $x0+$rad]
    set y1 [expr $ly-$y0-$rad]
    set y2 [expr $ly-$y0+$rad]
    set r(0) [.c0 create oval $x1 $y1 $x2 $y2 -fill brown]
    reset
}
button .b0 -text "Quit" -command "exit"
pack .cb0 .b1 .b2 .b3 .b0 -in .f0 -side left
pack .c0 .f2 .f1 .f0
set xp0 [expr $x0-$vx0*$dt];
set yp0 [expr $y0-$vy0*$dt];
set xn0 $x0;set yn0 $y0
set x1 [expr $x0-$rad]
set x2 [expr $x0+$rad]
set y1 [expr $ly-$y0-$rad]
set y2 [expr $ly-$y0+$rad]
set r(0) [.c0 create oval $x1 $y1 $x2 $y2 -fill brown]

.c0 bind $r(0) <B1-Motion> {
    set ifmv 0
    set x0 [expr %x]
    set y0 [expr $ly-%y]
    set xp0 [expr $x0-$vxini*$dt]; set xn0 $x0;
    set yp0 [expr $y0-$vyini*$dt]; set yn0 $y0;
    set x1 [expr $x0-$rad]
    set x2 [expr $x0+$rad]
    set y1 [expr $ly-$y0-$rad]
    set y2 [expr $ly-$y0+$rad]
    .c0 coords $r(0) $x1 $y1 $x2 $y2
}

proc setcoord {a b} {
    global vxini vyini dt x0 y0 xp0 yp0 xn0 yn0 vxini vyini
    global x0f y0f vx0f vy0f r scl rad ly vx0 vy0
    set x0 $a
    set y0 $b
    set xp0 [expr $x0-$vxini*$dt]; set xn0 $x0;
    set yp0 [expr $y0-$vyini*$dt]; set yn0 $y0;
    set x1 [expr $x0*$scl-$rad]
    set x2 [expr $x0*$scl+$rad]
    set y1 [expr $ly-$y0*$scl-$rad]
    set y2 [expr $ly-$y0*$scl+$rad]
    set  x0f [format "%5.3f"  $x0]
    set  y0f [format "%5.3f"  $y0]
    set vx0f [format "%5.3f" $vx0]
    set vy0f [format "%5.3f" $vy0]
    .c0 coords $r(0) $x1 $y1 $x2 $y2
    update
}

proc reset {} {
    global xini yini
    setcoord $xini $yini
}

proc onestep {} {
    global alp dt m x0 y0 xn0 yn0 xp0 yp0 g
    global vx0 vy0 scl rad ly x0f y0f vx0f vy0f vx0 vy0 r
    set xx [expr $alp*$dt/$m/2]
    set xn0 [expr ($x0*2-(1.0-$xx)*$xp0)/(1+$xx)]
    set yn0 [expr ($y0*2-(1.0-$xx)*$yp0-$g*$dt*$dt)/(1.+$xx)]
    set vx0 [expr ($xn0-$x0)/$dt]
    set vy0 [expr ($yn0-$y0)/$dt]
    set xp0 $x0; set x0 $xn0
    set yp0 $y0; set y0 $yn0
    set x1 [expr $x0*$scl-$rad]
    set x2 [expr $x0*$scl+$rad]
    set y1 [expr $ly-$y0*$scl-$rad]
    set y2 [expr $ly-$y0*$scl+$rad]
    set  x0f [format "%5.3f"  $x0]
    set  y0f [format "%5.3f"  $y0]
    set vx0f [format "%5.3f" $vx0]
    set vy0f [format "%5.3f" $vy0]
    .c0 coords $r(0) $x1 $y1 $x2 $y2
    .c0 create oval $x1 $y1 $x2 $y2
    update
}

while {1} {
#    global ifmv ifib
    if {$ifmv == 1} then {
    #set xn0 [expr $x0*2-$xp0]
    #set yn0 [expr $y0*2-$yp0-$g*$dt*$dt]
    set xx [expr $alp*$dt/$m/2]
    set xn0 [expr ($x0*2-(1.0-$xx)*$xp0)/(1+$xx)]
    set yn0 [expr ($y0*2-(1.0-$xx)*$yp0-$g*$dt*$dt)/(1.+$xx)]

    if {$yn0<0} then {
        set ifmv 0
    }
    set vx0 [expr ($xn0-$x0)/$dt]
    set vy0 [expr ($yn0-$y0)/$dt]
    set xp0 $x0; set x0 $xn0
    set yp0 $y0; set y0 $yn0
    set x1 [expr $x0*$scl-$rad]
    set x2 [expr $x0*$scl+$rad]
    set y1 [expr $ly-$y0*$scl-$rad]
    set y2 [expr $ly-$y0*$scl+$rad]
    set  x0f [format "%5.3f"  $x0]
    set  y0f [format "%5.3f"  $y0]
    set vx0f [format "%5.3f" $vx0]
    set vy0f [format "%5.3f" $vy0]
    .c0 coords $r(0) $x1 $y1 $x2 $y2
    .c0 create oval $x1 $y1 $x2 $y2
    update
    } else { update }
}

3原子の分子動力学（ＭＤ）計算

もっとも単純なＭＤ計算の例として、互いにモースポテンシャルで相互作用する3原子が運動する様子をシミュレーションするプログラム。

• 初期条件が入力できるようになっています（resetでそこから再スタート）。
• moveにチェックを入れるとＭＤがスタートします。3原子の座標と速度が表示されるようになっています。
• 3つの原子はマウスでドラッグして動かすことができます（なので、tcl/tkにおいてキャンバス上に描いた物体をマウスで任意に動かすための例としても参考になると思います）。

md_3atoms.tcl

#!/usr/bin/wish

set lx 600
set ly 600
set ox 150
set oy 150
set ifmv 0
#set ifib 1
canvas .c0 -width $lx -height $ly
#menu .m -type menubar
#. configure -menu .m
#.m add cascade -label "File" -under 0 -menu .m.m1
set natom 3
set dt 1.0e-16
set scl 3.0e11
set ang 1.0e10
set vxini(0)  0e3
set vyini(0)  0e3
set vxini(1)  0e3
set vyini(1)  0e3
set vxini(2)  0e3
set vyini(2)  0e3
set xini(0)  0.0e-10
set yini(0)  0.0e-10
set xini(1)  5.0e-10
set yini(1)  0.0e-10
set xini(2)  5.0e-10
set yini(2)  5.0e-10
set wm [expr 1.6726e-27*1.0]
for {set i 0} {$i < $natom} {incr i} {
    set  x($i)  $xini($i)
    set  y($i)  $yini($i)
    set vx($i) $vxini($i)
    set vy($i) $vyini($i)
    set fx($i) 0
    set fy($i) 0
}
frame .f0
frame .f1
frame .f2
frame .f3
label .lb000 -text "Initial(1) x"
entry .e000 -width 10 -textvariable xini(0)
label .lb001 -text "y"
entry .e001 -width 10 -textvariable yini(0)
label .lb002 -text "Vx"
entry .e002 -width 10 -textvariable vxini(0)
label .lb003 -text "Vy"
entry .e003 -width 10 -textvariable vyini(0)
pack .lb000 .e000 .lb001 .e001 .lb002 .e002 .lb003 .e003 -in .f1 -side left
label .lb100 -text "Initial(2) x"
entry .e100 -width 10 -textvariable xini(1)
label .lb101 -text "y"
entry .e101 -width 10 -textvariable yini(1)
label .lb102 -text "Vx"
entry .e102 -width 10 -textvariable vxini(1)
label .lb103 -text "Vy"
entry .e103 -width 10 -textvariable vyini(1)
pack .lb100 .e100 .lb101 .e101 .lb102 .e102 .lb103 .e103 -in .f2 -side left
label .lb200 -text "Initial(3) x"
entry .e200 -width 10 -textvariable xini(2)
label .lb201 -text "y"
entry .e201 -width 10 -textvariable yini(2)
label .lb202 -text "Vx"
entry .e202 -width 10 -textvariable vxini(2)
label .lb203 -text "Vy"
entry .e203 -width 10 -textvariable vyini(2)
pack .lb200 .e200 .lb201 .e201 .lb202 .e202 .lb203 .e203 -in .f3 -side left

for {set i 0} {$i < $natom} {incr i} {
    set  xf(i) [format "%5.2f" [expr  $x($i)*$ang]]
    set  yf(i) [format "%5.2f" [expr  $y($i)*$ang]]
    set vxf(i) [format "%5.2f" [expr $vx($i)]]
    set vyf(i) [format "%5.2f" [expr $vy($i)]]
}
frame .f00
frame .f01
frame .f02
label .lb00 -text "1:(x,y)\[A\]"
label .lb01 -width 8 -textvariable xf(0)
label .lb02 -width 8 -textvariable yf(0)
label .lb03 -text "(Vx,Vy)\[m/s\]"
label .lb04 -width 10 -textvariable vxf(0)
label .lb05 -width 10 -textvariable vyf(0)
pack .lb00 .lb01 .lb02 .lb03 .lb04 .lb05 -in .f00 -side left
label .lb10 -text "2:(x,y)\[A\]"
label .lb11 -width 8 -textvariable xf(1)
label .lb12 -width 8 -textvariable yf(1)
label .lb13 -text "(Vx,Vy)\[m/s\]"
label .lb14 -width 10 -textvariable vxf(1)
label .lb15 -width 10 -textvariable vyf(1)
pack .lb10 .lb11 .lb12 .lb13 .lb14 .lb15 -in .f01 -side left
label .lb20 -text "3:(x,y)\[A\]"
label .lb21 -width 8 -textvariable xf(2)
label .lb22 -width 8 -textvariable yf(2)
label .lb23 -text "(Vx,Vy)\[m/s\]"
label .lb24 -width 10 -textvariable vxf(2)
label .lb25 -width 10 -textvariable vyf(2)
pack .lb20 .lb21 .lb22 .lb23 .lb24 .lb25 -in .f02 -side left

checkbutton .cb0 -text "Move" -variable ifmv
#checkbutton .cb1 -text "in Box" -variable ifib
button .b1 -text "Reset" -command "reset"
button .b0 -text "Quit" -command "exit"
pack .cb0 .b1 .b0 -in .f0 -side left
pack .c0 .f00 .f01 .f02 .f1 .f2 .f3 .f0
set g 9.8
set rad 25
for {set i 0} {$i < $natom} {incr i} {
    set x1 [expr $x($i)*$scl-$rad+$ox]
    set x2 [expr $x($i)*$scl+$rad+$ox]
    set y1 [expr $ly-$y($i)*$scl-$rad-$oy]
    set y2 [expr $ly-$y($i)*$scl+$rad-$oy]
    set r($i) [.c0 create oval $x1 $y1 $x2 $y2 -fill brown]
}

for {set i 0} {$i < $natom} {incr i} {
    .c0 bind $r($i) <B1-Motion> {
    set ifmv 0
    set x1 [expr %x-$rad]
    set x2 [expr %x+$rad]
    set y1 [expr %y-$rad]
    set y2 [expr %y+$rad]
    .c0 coords current $x1 $y1 $x2 $y2
    # Now find which object the mouse is pointing..
    for {set j 1} {$j <= $natom} {incr j} {
        if {[.c0 gettags $j]!=""} then {
        set ii [expr $j-1]
        set x($ii) [expr (%x-$ox)/$scl]
        set y($ii) [expr ($ly-%y-$oy)/$scl]
        }
    }
    }
}

proc reset {} {
    global natom x y xini yini vxini vyini fx fy
    global scl rad ox oy ly xf yf vxf vyf r ang
    for {set i 0} {$i < $natom} {incr i} {
    set  x($i)  $xini($i)
    set  y($i)  $yini($i)
    set vx($i) $vxini($i)
    set vy($i) $vyini($i)
    set fx($i) 0
    set fy($i) 0

    set x1 [expr $x($i)*$scl-$rad+$ox]
    set x2 [expr $x($i)*$scl+$rad+$ox]
    set y1 [expr $ly-$y($i)*$scl-$rad-$oy]
    set y2 [expr $ly-$y($i)*$scl+$rad-$oy]
    set  xf($i) [format "%5.2f" [expr  $x($i)*$ang]]
    set  yf($i) [format "%5.2f" [expr  $y($i)*$ang]]
    set vxf($i) [format "%5.2f" [expr $vx($i)]]
    set vyf($i) [format "%5.2f" [expr $vy($i)]]
    .c0 coords $r($i) $x1 $y1 $x2 $y2
    }
    update
}

proc morsev {r} {
    set rang [expr $r*1.0e10]
    set ep 0.2703
    set al 1.1646
    set ro 3.2530
    set ev 1.6021892e-19
    set vval [expr $ep*(exp(-1.0*$al*($rang-ro))-2.0*exp(-$al*($rang-$ro)))]
    return $vval
}
proc morsevp {r} {
    set rang [expr $r*1.0e10]
    set ep 0.2703
    set al 1.1646
    set ro 3.2530
    set ev 1.6021892e-19
    set vpval [expr -2.0*$al*$ep*(exp(-2.0*$al*($rang-$ro))-exp(-$al*($rang-$ro))) *$ev*1.0e10]
    return $vpval
}

proc calcforce {} {
    global natom x y fx fy
    for {set i 0} {$i < $natom} {incr i} {
    set fx($i) 0
    set fy($i) 0
    for {set j 0} {$j < $natom} {incr j} {
        if {$i != $j} then {
        set xx [expr $x($i)-$x($j)]
        set yy [expr $y($i)-$y($j)]
        set rr [expr sqrt($xx*$xx+$yy*$yy)]
        set ff [morsevp $rr]
        set fx($i) [expr $fx($i)-$ff/$rr*$xx]
        set fy($i) [expr $fy($i)-$ff/$rr*$yy]
        }
    }
    }
}

while {1} {
    #global ifmv ifib
    if {$ifmv == 1} then {
    for {set i 0} {$i < $natom} {incr i} {
        set  x($i) [expr $x($i)+$dt*$vx($i)+($dt*$dt/2.0)*$fx($i)/$wm]
        set  y($i) [expr $y($i)+$dt*$vy($i)+($dt*$dt/2.0)*$fy($i)/$wm]
        set vx($i) [expr $vx($i)+$dt/2.0*$fx($i)/$wm]
        set vy($i) [expr $vy($i)+$dt/2.0*$fy($i)/$wm]
        #if {$ifib == 1} then {
        #}
    }
    calcforce
    for {set i 0} {$i < $natom} {incr i} {
        set vx($i) [expr $vx($i)+$dt/2.0*$fx($i)/$wm]
        set vy($i) [expr $vy($i)+$dt/2.0*$fy($i)/$wm]
        set x1 [expr $x($i)*$scl-$rad+$ox]
        set x2 [expr $x($i)*$scl+$rad+$ox]
        set y1 [expr $ly-$y($i)*$scl-$rad-$oy]
        set y2 [expr $ly-$y($i)*$scl+$rad-$oy]
        set  xf($i) [format "%5.2f" [expr  $x($i)*$ang]]
        set  yf($i) [format "%5.2f" [expr  $y($i)*$ang]]
        set vxf($i) [format "%5.2f" [expr $vx($i)]]
        set vyf($i) [format "%5.2f" [expr $vy($i)]]
        .c0 coords $r($i) $x1 $y1 $x2 $y2
    }
    update
    } else { update }
}