Броуновское движение (отталкивание рандомных частиц друг от друга и заданных границ) - VB

Option Explicit
 'X направлен вниз (совпадает с Top, Height)
 'Y -------- влево (----------- Left, Width)
 
 Const kv As Double = 0.999 'потеря энергии при ударе
 Const dV As Double = 0.001 'потеря скорости в очередном Timer1_Timer()
 
 Const N As Long = 15
 Const D As Double = 40
 Const TI As Long = 10
 Const Z As Double = 1.000001
 
 Dim meWidth As Double, meHeight As Double
 Dim meWidthOld As Double, meHeightOld As Double
 Dim a_VX(N) As Double
 Dim a_VY(N) As Double
 Dim a_V(N) As Double
 Dim a_X(N) As Double, a_Y(N) As Double
 Dim TN As Double, TN0 As Double, tik As Long
 Dim II As Integer, JJ As Integer
 Dim vXij As Double, vYij As Double
 Dim xI As Double, yI As Double, xRel As Double, yRel As Double, xJ As Double, yJ As Double
 
 Private Sub Form_Load()
 Dim x1 As Double, y1 As Double, x2 As Double, y2 As Double, X As Double, Y As Double, xn As Double, yn As Double
 Dim I As Integer, J As Integer
 Timer1.Enabled = False
 Me.BorderStyle = 1
 Me.Left = 0
 Me.Top = 0
 Me.Height = 9480
 Me.Width = 18090
 Me.ScaleMode = 3
 Me.BackColor = vbGreen / 1.5
 Shape1(0).Height = D
 Shape1(0).Width = D
 Shape1(0).Shape = 3 '2
 Shape1(0).BackStyle = 1
 Shape1(0).BackColor = vbYellow
 Shape1(0).BorderColor = vbYellow
 a_X(0) = Me.ScaleHeight / 2
 Shape1(0).Top = a_X(0)
 a_Y(0) = 0
 Shape1(0).Left = a_Y(0)
 a_VX(0) = 0
 a_VY(0) = 10
 a_V(0) = Sqr(a_VX(0) * a_VX(0) + a_VY(0) * a_VY(0))
 Y = Me.ScaleWidth / 2 - 3 * D
 Do
   X = Me.ScaleHeight / 2 - J * D * Z / 2
   For J = 0 To J
     I = I + 1
     Load Shape1(I)
     With Shape1(I)
      .Visible = True
      .Shape = 3
      .Width = D
      .Height = D
      a_X(I) = X
      .Top = a_X(I)
      a_Y(I) = Y
      .Left = a_Y(I)
     End With
     
     If I > N Then Exit For
     X = X + D * Z
   Next J
   Y = Y + D * Z * Sqr(3) / 2
 Loop While I < N
 
 For I = N + 1 To N + 6
   Load Shape1(I)
   With Shape1(I)
     .Visible = True
     .Shape = 4
     .BackColor = vbGreen / 2
     .BorderColor = vbBlack
     If I <= N + 2 Then
       .Width = D * 4
       .Height = Me.ScaleHeight - 2 * D
       .Top = D
       .Left = -D * 3.5 - (I = N + 2) * (Me.ScaleWidth + D * 3)
     Else
       .Width = Me.ScaleWidth / 2 - D * 1.5
       .Height = D * 4
       If I <= N + 4 Then
         .Top = -3.5 * D
         .Left = D - (I = N + 4) * (Me.ScaleWidth / 2 - D / 2)
       Else
         .Top = -D * 3.5 + (Me.ScaleHeight + D * 3)
         .Left = D - (I = N + 6) * (Me.ScaleWidth / 2 - D / 2)
       End If
     End If
   End With
 
 Next I
 
 near_Bac
 Timer1_Timer
 Timer1.Enabled = True
 End Sub
 
 Sub near_Bac()
 Dim I As Integer, J As Integer, T As Double
 Dim V As Double, xn As Double, yn As Double
 Dim dX As Double, dY As Double
 
 TN = 1.79769313486231E+308
 For I = 0 To N - 1
   For J = I + 1 To N
     vXij = a_VX(J) - a_VX(I)
     vYij = a_VY(J) - a_VY(I)
     If vXij <> 0 Or vYij <> 0 Then
       T = 1.79769313486231E+308
       dX = a_X(I) - a_X(J)
       dY = a_Y(I) - a_Y(J)
       V = Sqr(vXij * vXij + vYij * vYij)
       yn = (dX * vYij - vXij * dY) / V
       If Abs(yn) <= D Then
         xn = (dX * vXij + vYij * dY) / V
         If xn >= 0 Then T = (xn - Sqr(D * D - yn * yn)) / V
       End If
       If T < TN Then
         TN = T: II = I: JJ = J
       End If
     End If
   Next J
 Next I
 For I = 0 To N
   If a_VX(I) > 0 Then
     T = (Me.ScaleHeight - a_X(I) - D * 1.5) / a_VX(I): If T < TN Then II = I: JJ = -1: TN = T
   ElseIf a_VX(I) < 0 Then
     T = (D / 2 - a_X(I)) / a_VX(I): If T < TN Then II = I: JJ = -1: TN = T
   End If
   If a_VY(I) > 0 Then
     T = (Me.ScaleWidth - a_Y(I) - D * 1.5) / a_VY(I): If T < TN Then II = I: JJ = -2: TN = T
   ElseIf a_VY(I) < 0 Then
     T = (D / 2 - a_Y(I)) / a_VY(I): If T < TN Then II = I: JJ = -2: TN = T
   End If
 Next I
 If TN < 0 Then TN = 0
 TN0 = TN
 tik = 0
 End Sub
 
 Private Sub Timer1_Timer()
 Dim I As Long
 Dim K As Double, dX As Double, dY As Double
 Dim VIR As Double, VJR As Double
 Do
   If TN < 1 Then
       For I = 0 To N
         a_Y(I) = a_Y(I) + a_VY(I) * TN: Shape1(I).Left = a_Y(I)
         a_X(I) = a_X(I) + a_VX(I) * TN: Shape1(I).Top = a_X(I)
         If a_V(I) > 0 Then
           K = (a_V(I) - TN0 * dV) / a_V(I): If K < 0 Then K = 0
           a_V(I) = a_V(I) * K
           a_VX(I) = a_VX(I) * K
           a_VY(I) = a_VY(I) * K
         End If
       Next I
       Select Case JJ
         Case -1: a_VX(II) = -a_VX(II) * kv: a_VY(II) = a_VY(II) * kv
         Case -2: a_VY(II) = -a_VY(II) * kv: a_VX(II) = a_VX(II) * kv
         Case Else
           dX = (a_X(II) - a_X(JJ)) / D
           dY = (a_Y(II) - a_Y(JJ)) / D
           vXij = a_VX(JJ) - a_VX(II)
           vYij = a_VY(JJ) - a_VY(II)
           VIR = (dY * vYij + dX * vXij)
           VJR = (dX * vYij - vXij * dY)
           a_VX(JJ) = (a_VX(II) - VJR * dY) * kv
           a_VY(JJ) = (a_VY(II) + VJR * dX) * kv
           a_VX(II) = (a_VX(II) + VIR * dX) * kv
           a_VY(II) = (a_VY(II) + VIR * dY) * kv
           a_V(II) = Sqr(a_VX(II) * a_VX(II) + a_VY(II) * a_VY(II))
           a_V(JJ) = Sqr(a_VX(JJ) * a_VX(JJ) + a_VY(JJ) * a_VY(JJ))
       End Select
       near_Bac
   Else
     For I = 0 To N
       a_Y(I) = a_Y(I) + a_VY(I): Shape1(I).Left = a_Y(I)
       a_X(I) = a_X(I) + a_VX(I): Shape1(I).Top = a_X(I)
     Next I
     TN = TN - 1
   End If
   If TN < 1 Then
     Timer1.Interval = TN * TI
     Me.Caption = Timer1.Interval
   Else
     Timer1.Interval = TI
   End If
 Loop While Timer1.Interval = 0
 End Sub