Public Class Form1 Private Const Xmax = 10 Private Const Ymax = 10 Private Const DispR = 0.5 Private Const EPSI = 1.0# Private Const SIGMA = 1 Private Structure Point_Data Public X As Double : Public Y As Double End Structure Private N As Integer : Private R(,) As Point_Data : Private F() As Point_Data Private V(,) As Point_Data : Private E1 As Double : Private E2 As Double Private MG As Double : Private DT As Double : Private ID1 As Integer Private ID2 As Integer Dim img As New Bitmap(300, 300) '300x300サイズのImage Private Sub 初期設定() Dim i As Integer, k As Integer, NN As Integer Dim DX As Double, DY As Double N = Val(TextBox1.Text) : DT = Val(TextBox2.Text) ReDim R(1, N) : ReDim F(N) : ReDim V(1, N) E1 = -48 * EPSI * SIGMA ^ 12 : E2 = 24 * EPSI * SIGMA ^ 12 MG = 1 : ID1 = 0 : ID2 = 1 For i = 1 To N F(i).X = 0 : F(i).Y = 0 : V(0, i).X = 0 : V(0, i).Y = 0 Next Randomize() : i = 0 '分子の生成 Do While i < N R(0, i + 1).X = Rnd() * (Xmax - 2) + 1 : R(0, i + 1).Y = Rnd() * (Ymax - 2) + 1 k = 0 : NN = 0 Do '便宜的に100回繰り返しても近すぎる座標のとき諦める k = k + 1 : NN = NN + 1 DX = R(0, k).X - R(0, i + 1).X : DY = R(0, k).Y - R(0, i + 1).Y Loop Until DX * DX + DY * DY < 1 Or NN > 100 If k > i Or NN > 100 Then i = i + 1 Loop 表示() End Sub Private Sub Circle(g As Graphics, X As Double, Y As Double, R As Double) Dim IX As Integer, IY As Integer, IR As Integer IR = R * 60 + 0.5 IX = (X - R) * 30 + 0.5 : IY = (Y - R) * 30 + 0.5 g.FillEllipse(Brushes.Red, IX, IY, IR, IR) End Sub Private Sub 表示() Dim g As Graphics = Graphics.FromImage(img) g.Clear(Color.White) For i = 1 To N Circle(g, R(ID1, i).X, R(ID1, i).Y, DispR) Next g.Dispose() 'リソースを解放する PictureBox1.Image = img End Sub Private Sub 力計算() Dim DX As Double, DY As Double, DR As Double, FF As Double For i = 1 To N F(i).X = 0 : F(i).Y = 0 For j = 1 To N If i <> j Then DX = R(ID1, j).X - R(ID1, i).X : DY = R(ID1, j).Y - R(ID1, i).Y DR = Math.Sqrt(DX * DX + DY * DY) FF = 0 If DR < 1 Then DR = 1 FF = E1 / (DR ^ 13) + E2 / (DR ^ 7) F(i).X = F(i).X + FF * DX / DR F(i).Y = F(i).Y + FF * DY / DR End If Next Next End Sub Private Sub 速度計算() Dim T As Double, VX As Double, VY As Double For i = 1 To N V(ID2, i).X = V(ID1, i).X + DT * F(i).X / MG V(ID2, i).Y = V(ID1, i).Y + DT * F(i).Y / MG Next T = 0 ' 温度境界条件 For i = 1 To N VX = V(ID2, i).X : VY = V(ID2, i).Y T = T + VX * VX + VY * VY Next T = Math.Sqrt(T) * 0.1 If T > 0.00001 Then For i = 1 To N V(ID2, i).X = V(ID2, i).X / T V(ID2, i).Y = V(ID2, i).Y / T Next End If End Sub Private Sub 位置計算() For i = 1 To N R(ID2, i).X = R(ID1, i).X + DT * V(ID2, i).X R(ID2, i).Y = R(ID1, i).Y + DT * V(ID2, i).Y Next End Sub Private Sub 壁衝突() For i = 1 To N If R(ID2, i).X < (0.5 * SIGMA) Then R(ID2, i).X = (0.5 * SIGMA) V(ID2, i).X = Math.Abs(V(ID2, i).X) End If If R(ID2, i).X > Xmax - (0.5 * SIGMA) Then R(ID2, i).X = Xmax - (0.5 * SIGMA) V(ID2, i).X = -Math.Abs(V(ID2, i).X) End If If R(ID2, i).Y < (0.5 * SIGMA) Then R(ID2, i).Y = (0.5 * SIGMA) V(ID2, i).Y = Math.Abs(V(ID2, i).Y) End If If R(ID2, i).Y > Ymax - (0.5 * SIGMA) Then R(ID2, i).Y = Ymax - (0.5 * SIGMA) V(ID2, i).Y = -Math.Abs(V(ID2, i).Y) End If Next End Sub Private Sub Form1_Load(sender As System.Object, e As System.EventArgs) Handles MyBase.Load With PictureBox1 .Top = 10 : .Left = 10 : .Width = img.Width : .Height = img.Height End With End Sub Private Sub Button1_Click(sender As System.Object, e As System.EventArgs) Handles Button1.Click 初期設定() Timer1.Interval = Val(TextBox3.Text) Timer1.Enabled = True End Sub Private Sub Timer1_Tick(sender As System.Object, e As System.EventArgs) Handles Timer1.Tick Dim IDD As Integer 力計算() 速度計算() 位置計算() 壁衝突() 表示() IDD = ID1 : ID2 = ID1 : ID1 = IDD End Sub Private Sub Button2_Click(sender As System.Object, e As System.EventArgs) Handles Button2.Click Timer1.Enabled = False End Sub End Class