# Excel physics

This post is a recopilation of functions that can be used to compute simple physics in Excel. Basically, most of the physics phenomenon deals with parabollic shoot. When considering this, one should recall on the Laws of Newton, and more specifically the 2nd one, Momentum conservation. Focussing only on collisions for a dynamic particles system, a good starting reference is this video, following this info, with code avaliable at GitHub: https://www.youtube.com/watch?v=irbshkdVFao Here the author points out the three basic problems that are not usually considered with collisions:
• excessive computations (dt too small, that’s when no collisions at fixed time increments)
• miss collisions (dt too large, that’s when space travelled at one fixed time increment is larger than the collision range for two particles)
• lack of position precission when calculating the collisions at fixed spaced moments.
More equations and explanation is given at this site. Further on, it could be brought some attention to drag resistance, and even Magnus effects, in order to refine the precission of the movements. Following is some of the Physic functions already implemented:
```Option Explicit

Private Const PI As Double = 3.14159265358979
Public Const EPSILON As Double = 0.0000001

Private Type tXYZ
X As Double
Y As Double
Z As Double
End Type
Private Type tCollision
Shp1 As Long
Shp2 As Long
Time As Double
End Type

Public Function fXYZ(Optional ByVal X As Double = 0, _
Optional ByVal Y As Double = 0, _
Optional ByVal Z As Double = 0) As tXYZ
With fXYZ
.X = X
.Y = Y
.Z = Z
End With
End Function

Public Function fVector(Optional ByVal dbModule As Double = 0, _
Optional ByVal ß As Double = 0, _
Optional ByVal Ø As Double = 0) As tXYZ
With fVector
.X = (dbModule * Cos(ß) * Cos(Ø))
.Y = (dbModule * Cos(ß) * Sin(Ø))
.Z = (dbModule * Sin(ß))
End With
End Function

Public Function fVectorModule(ByRef oVector As tXYZ) As Double
With oVector
fVectorModule = VBA.Sqr(.X ^ 2 + .Y ^ 2 + .Z ^ 2)
End With
End Function

Public Function fToDouble(ByRef vVariable As Variant) As Double()
Dim aDouble() As Double
Dim lgElement As Long

If IsArray(vVariable) Then
ReDim aDouble(LBound(vVariable) To UBound(vVariable))
For lgElement = LBound(vVariable) To UBound(vVariable)
aDouble(lgElement) = VBA.CDbl(vVariable(lgElement))
Next lgElement
fToDouble = aDouble()
End If
End Function

'-------------------------------

Public Sub sShoot()
Dim aTime() As Double
Dim oPoint() As tXYZ

oPoint() = fShoot(aTime:=fToDouble(fNewVector("0:1:10)")), _
dbStrength:=10, _
ß:=45, _
Ø:=0, _
oForce:=fXYZ(0, 0, 0), _
dbGravity:=9.81, _
dbMass:=10, _
dbMediaDensity:=0.1, _
dbAreaX:=10, _
dbAreaY:=10, _
dbAreaZ:=10, _
lgShapeX:=msoShapeRectangle, _
lgShapeY:=msoShapeRectangle, _
lgShapeZ:=msoShapeRectangle)
End Sub

Public Function fShoot(ByRef aTime() As Double, _
ByVal dbStrength As Double, _
ByVal ß As Double, _
ByVal Ø As Double, _
ByRef oForce As tXYZ, _
Optional ByVal dbGravity As Double = 9.81, _
Optional ByVal dbMass As Double = 0, _
Optional ByVal dbMediaDensity As Double = 0, _
Optional ByVal dbAreaX As Double = 0, _
Optional ByVal dbAreaY As Double = 0, _
Optional ByVal dbAreaZ As Double = 0, _
Optional ByVal lgShapeX As Long = msoShapeRectangle, _
Optional ByVal lgShapeY As Long = msoShapeRectangle, _
Optional ByVal lgShapeZ As Long = msoShapeRectangle) As tXYZ()
'oForce As tXYZ
' ••••••    vector
' ¤¤¤¤¤¤    projection on XY plane
' ......    arc
' ß  arc planeXY to vector
' Ø  arc planeXZ to vector XY projection
'
'       |   •
'       |  •
'       | •.
'       |•_.______
'      / ¤ . ß
'     /....¤
'    /  Ø    ¤
'
' Initial speed components:
' Vo,x = (Vo · Cos(ß)) · Cos(Ø)
' Vo,y = (Vo · Cos(ß)) · Sin(Ø)
' Vo,z = (Vo · Sin(ß))

' Dragg = (1 / 2) · MediaDensity · Cd · Area · V²
' F = m · a --> a = F / m  //  a = dV/dt  --> dV = a · dt = (F / m) · dt

' Speeds in any instant:
' Vx = Vox - (ResistanceX · t / Mass) + (t · ForceX / Mass)
' Vy = Voy - (ResistanceY · t / Mass) + (t · ForceY / Mass)
' Vz = Voz - (ResistanceZ · t / Mass) + (t · ForceZ / Mass) - (dbGravity · t)

' Position at any instant:
' X = Vox · t - (1/2 · ResistanceX · t² / Mass) + (1/2 · t² · ForceX / Mass)
' Y = Voy · t - (1/2 · ResistanceY · t² / Mass) + (1/2 · t² · ForceY / Mass)
' Z = Voz · t - (1/2 · ResistanceZ · t² / Mass) + (1/2 · t² · ForceZ / Mass) - (1/2 · dbGravity · t²)

Dim lgTime As Long
Dim t² As Double
Dim oPoint() As tXYZ
Dim oDrag As tXYZ
Dim oVel() As tXYZ
Dim oVo As tXYZ
Dim Vo As Double
Dim Speed As Double

'!!!!!!!!!!!
Vo = dbStrength
'Speed = fVectorModule(oVel(lgTime))
'!!!!!!!!!!!

'Initial speed components:
oVo = fVector(dbStrength, ß, Ø)

If dbMediaDensity  0 Then
'Dragg = (1 / 2) · CD · Area · V²
'F = m · a --> a = F / m  //  a = dV/dt  --> dV = a · dt = (F / m) · dt
If dbAreaX  0 Then
oDrag.X = fDrag(oVel(lgTime).X, dbMediaDensity, dbAreaX, lgShapeX)
End If
If dbAreaY  0 Then
oDrag.Y = fDrag(oVel(lgTime).Y, dbMediaDensity, dbAreaY, lgShapeY)
End If
If dbAreaZ  0 Then
oDrag.Z = fDrag(oVel(lgTime).Z, dbMediaDensity, dbAreaZ, lgShapeZ)
End If
End If

'Speeds at any instant:
ReDim oVel(LBound(aTime) To UBound(aTime))
For lgTime = LBound(aTime) To UBound(aTime)
With oVel(lgTime)
.X = oVo.X + ((oForce.X - oDrag.X) / dbMass) * aTime(lgTime)
.Y = oVo.Y + ((oForce.Y - oDrag.Y) / dbMass) * aTime(lgTime)
.Z = oVo.Z + (((oForce.Z - oDrag.Z) / dbMass) - dbGravity) * aTime(lgTime)
End With
Next lgTime

'Position at any instant:
ReDim oPoint(LBound(aTime) To UBound(aTime))
For lgTime = LBound(aTime) To UBound(aTime)
With oPoint(lgTime)
t² = aTime(lgTime) * aTime(lgTime)
.X = oVo.X * aTime(lgTime) _
+ (((oForce.X - oDrag.X) / dbMass)) * t²
.Y = oVo.Y * aTime(lgTime) _
+ (((oForce.Y - oDrag.Y) / dbMass)) * t²
.Z = oVo.Z * aTime(lgTime) _
+ (((oForce.Z - oDrag.Z) / dbMass) - dbGravity) * t²
End With
Next lgTime

fShoot = oPoint()
End Function

Public Function fDrag(Optional ByVal dbVelocity As Double = 0, _
Optional ByVal dbMediaDensity As Double = 0, _
Optional ByVal dbArea As Double = 0, _
Optional ByVal lgShape As Long = 0) As Double
' For a body following an unidirectional path, will compute the drag force opposed to movement
' Dragg = (1 / 2) · MediaDensity · Cd · Area · V²
Dim dbCd As Double

Select Case lgShape
'Case Is = msoShape...: dbCd = ...
End Select

fDrag = (1 / 2) * dbMediaDensity * dbCd * dbArea * (dbVelocity ^ 2)
End Function

Public Function fHooke() 'Optional ByVal dbStrength As Double = 0, _
Optional ByVal dbElasticity As Double = 0) As Boolean
'For any object collisioning with another one, Hooke law will have an effect on the shape of both objects
End Function

Public Sub Animate()
Dim oShpFrm As Excel.Shape
Dim lgShp As Long
Dim lgShpEval As Long
Dim oShp1 As Excel.Shape
Dim oShp2 As Excel.Shape

Dim Ovl1R As Single
Dim Ovl2R As Single
Dim CCDist As Single

Dim TopBox As Single
Dim BottBox As Single
Dim LeftBox As Single
Dim RightBox As Single

Dim CenterShp1 As tXYZ
Dim CenterShp2 As tXYZ
Dim DimShp1 As tXYZ
Dim DimShp2 As tXYZ
Dim vectorShp1 As tXYZ
Dim vectorShp2 As tXYZ
Dim Velocity As tXYZ
Dim CCAng As Single
Dim Shp2_Speed As Single
Dim Shp1_Speed As Single
Dim Angle_Shp2 As Single
Dim Angle_Shp1 As Single
Dim DX As Single
Dim DY As Single

Dim Start As Single
Dim TimeStep As Double
Dim TimeEval As Double
Dim TimeCollision As Double

With ActiveSheet
For Each oShpFrm In .Shapes
oShpFrm.Delete
Next oShpFrm

Velocity.X = 10 '.Range("Hspeed").Value
Velocity.Y = 10 '.Range("Vspeed").Value

TimeStep = 0.01

' Get frame limits
Set oShpFrm = .Shapes.AddShape(Type:=msoShapeRectangle, _
Left:=20, _
Top:=20, _
Width:=400, _
Height:=400)
'oShpFrm.Name = "Frame"
With oShpFrm
TopBox = .Top
BottBox = TopBox + .Height
LeftBox = .Left
RightBox = LeftBox + .Width
End With

'Random shape creation and speed vector assignment
DimShp1.X = (50 * Rnd())
DimShp1.Y = DimShp1.X '(50 * Rnd())
CenterShp1.X = LeftBox + ((RightBox - LeftBox - DimShp1.X) * Rnd())
CenterShp1.Y = TopBox + ((BottBox - TopBox - DimShp1.Y) * Rnd())
Set oShp1 = .Shapes.AddShape(Type:=msoShapeOval, _
Left:=CenterShp1.X, _
Top:=CenterShp1.Y, _
Width:=DimShp1.X, _
Height:=DimShp1.Y)
'oShp1.Name = "Oval1"
With vectorShp1
.X = Velocity.X * (((RightBox - LeftBox) / 1000) * Rnd())
.Y = Velocity.Y * (((BottBox - TopBox) / 1000) * Rnd())
End With
DimShp2.X = (50 * Rnd())
DimShp2.Y = DimShp2.X '(50 * Rnd())
CenterShp2.X = LeftBox + ((RightBox - LeftBox - DimShp2.X) * Rnd())
CenterShp2.Y = TopBox + ((BottBox - TopBox - DimShp2.Y) * Rnd())
Set oShp2 = .Shapes.AddShape(Type:=msoShapeOval, _
Left:=CenterShp2.X, _
Top:=CenterShp2.Y, _
Width:=DimShp2.X, _
Height:=DimShp2.Y)
'oShp1.Name = "Oval2"
With vectorShp2
.X = Velocity.X * (((RightBox - LeftBox) / 1000) * Rnd())
.Y = Velocity.Y * (((BottBox - TopBox) / 1000) * Rnd())
End With

Ovl1R = (DimShp1.X + DimShp1.Y) / 4
Ovl2R = (DimShp2.X + DimShp2.Y) / 4

' Random initial movements:
With vectorShp1
.X = Velocity.X * (((RightBox - LeftBox) / 1000) * Rnd())
.Y = Velocity.Y * (((BottBox - TopBox) / 1000) * Rnd())
End With
With vectorShp2
.X = Velocity.X * (((RightBox - LeftBox) / 1000) * Rnd())
.Y = Velocity.Y * (((BottBox - TopBox) / 1000) * Rnd())
End With

Do
With oShp1
.IncrementLeft vectorShp1.X
.IncrementTop vectorShp1.Y
CenterShp1.X = .Left + (DimShp1.X / 2)
CenterShp1.Y = .Top + (DimShp1.Y / 2)
End With
With vectorShp1
If (CenterShp1.X  RightBox - (DimShp1.X / 2)) Then .X = -.X
If (CenterShp1.Y  BottBox - (DimShp1.Y / 2)) Then .Y = -.Y
End With
With oShp2
.IncrementLeft vectorShp2.X
.IncrementTop vectorShp2.Y
CenterShp2.X = .Left + (DimShp2.X / 2)
CenterShp2.Y = .Top + (DimShp2.Y / 2)
End With
With vectorShp2
If (CenterShp2.X  RightBox - (DimShp2.X / 2)) Then .X = -.X
If (CenterShp2.Y  BottBox - (DimShp2.Y / 2)) Then .Y = -.Y
End With

'Distance between shapes
DX = (CenterShp1.X - CenterShp2.X)
DY = (CenterShp1.Y - CenterShp2.Y)
CCDist = Sqr(DX ^ 2 + DY ^ 2)

If CCDist  TimeEval Then
TimeCollision = TimeEval
Erase oCollision()
ReDim Preserve oCollision(g_Base)
oCollision(g_Base).Shp1 = lgShp
oCollision(g_Base).Shp2 = lgShpEval

Else 'If TimeCollision = TimeEval Then
'More than two objects colliding at the same moment:
TimeCollision = TimeEval
lgCollision = lgCollision + 1
ReDim Preserve oCollision(g_Base To lgCollision)
oCollision(lgCollision).Shp1 = lgShp
oCollision(lgCollision).Shp2 = lgShpEval
End If
End If
End If
Next lgShpEval

' Check collision against frame walls:
TimeEval = (CenterShp(lgShp).X - (LeftBox + DimShp(lgShp).X / 2)) _
/ vectorShp(lgShp).X
If TimeEval > 0 Then ' negative times implies that they are getting separated
If TimeCollision > TimeEval Then
TimeCollision = TimeEval
Erase oCollision()
lgCollision = g_Base
ReDim Preserve oCollision(g_Base To lgCollision)
oCollision(lgCollision).Shp1 = lgShp
oCollision(lgCollision).Shp2 = -xlEdgeLeft '7
ElseIf TimeCollision = TimeEval Then
TimeCollision = TimeEval
lgCollision = lgCollision + 1
ReDim Preserve oCollision(g_Base To lgCollision)
oCollision(g_Base).Shp1 = lgShp
oCollision(g_Base).Shp2 = -xlEdgeLeft '7
End If
End If
TimeEval = (RightBox - (DimShp(lgShp).X / 2) - CenterShp(lgShp).X) _
/ vectorShp(lgShp).X
If TimeEval > 0 Then ' negative times implies that they are getting separated
If TimeCollision > TimeEval Then
TimeCollision = TimeEval
Erase oCollision()
lgCollision = g_Base
ReDim Preserve oCollision(g_Base To lgCollision)
oCollision(lgCollision).Shp1 = lgShp
oCollision(lgCollision).Shp2 = -xlEdgeLeft '7
ElseIf TimeCollision = TimeEval Then
TimeCollision = TimeEval
lgCollision = lgCollision + 1
ReDim Preserve oCollision(g_Base To lgCollision)
oCollision(g_Base).Shp1 = lgShp
oCollision(g_Base).Shp2 = -xlEdgeRight '10
End If
End If
TimeEval = (CenterShp(lgShp).Y - (TopBox + DimShp(lgShp).Y / 2)) _
/ vectorShp(lgShp).Y
If TimeEval > 0 Then ' negative times implies that they are getting separated
If TimeCollision > TimeEval Then
TimeCollision = TimeEval
Erase oCollision()
lgCollision = g_Base
ReDim Preserve oCollision(g_Base To lgCollision)
oCollision(lgCollision).Shp1 = lgShp
oCollision(lgCollision).Shp2 = -xlEdgeLeft '7
ElseIf TimeCollision = TimeEval Then
TimeCollision = TimeEval
lgCollision = lgCollision + 1
ReDim Preserve oCollision(g_Base To lgCollision)
oCollision(g_Base).Shp1 = lgShp
oCollision(g_Base).Shp2 = -xlEdgeTop '8
End If
End If
TimeEval = (BottBox - (DimShp(lgShp).X / 2) - CenterShp(lgShp).Y) _
/ vectorShp(lgShp).Y
If TimeEval > 0 Then ' negative times implies that they are getting separated
If TimeCollision > TimeEval Then
TimeCollision = TimeEval
Erase oCollision()
lgCollision = g_Base
ReDim Preserve oCollision(g_Base To lgCollision)
oCollision(lgCollision).Shp1 = lgShp
oCollision(lgCollision).Shp2 = -xlEdgeLeft '7
ElseIf TimeCollision = TimeEval Then
TimeCollision = TimeEval
lgCollision = lgCollision + 1
ReDim Preserve oCollision(g_Base To lgCollision)
oCollision(g_Base).Shp1 = lgShp
oCollision(g_Base).Shp2 = -xlEdgeBottom '9
End If
End If
Next lgShp

' No object collide with anything until TimeCollision, so:
For lgShp = LBound(oShp) To UBound(oShp)
With oShp(lgShp)
.IncrementLeft (vectorShp(lgShp).X * TimeCollision)
.IncrementTop (vectorShp(lgShp).Y * TimeCollision)
CenterShp(lgShp).X = .Left + (DimShp(lgShp).X / 2)
CenterShp(lgShp).Y = .Top + (DimShp(lgShp).Y / 2)
End With
Next lgShp

Stop
' First check collisions against walls
lgCounter = LBound(oCollision)
For lgCollision = LBound(oCollision) To UBound(oCollision)
If oCollision(lgCollision).Shp2 < 0 Then 'Wall collision
If oCollision(lgCollision).Shp2 = xlEdgeLeft Then
With vectorShp(oCollision(lgCollision).Shp1)
.X = -.X
End With
End If
If oCollision(lgCollision).Shp2 = xlEdgeRight Then
With vectorShp(oCollision(lgCollision).Shp1)
.X = -.X
End With
End If
If oCollision(lgCollision).Shp2 = xlEdgeBottom Then
With vectorShp(oCollision(lgCollision).Shp1)
.Y = -.Y
End With
End If
If oCollision(lgCollision).Shp2 = xlEdgeTop Then
With vectorShp(oCollision(lgCollision).Shp1)
.Y = -.Y
End With
End If

Else
lgCounter = lgCounter + 1 'Counter with other particles
ReDim Preserve PtrCollision(g_Base To lgCounter)
PtrCollision(lgCounter) = lgCollision

'If they are not repeated...
'                    bStack = True
'                    For lgPtr = LBound(PtrCollision) To UBound(PtrCollision)
'                        If oCollision(PtrCollision(lgPtr)).Shp1 = ...lgShp1 Then
'                            bStack = False
'                            Exit For
'                        End If
'                        If oCollision(PtrCollision(lgPtr)).Shp2 = ...lgShp1 Then
'                            bStack = False
'                            Exit For
'                        End If
'                    Next lgPtr
'                    If bStack Then
'                        ReDim Preserve PtrObj(g_Base To lgCounter)
'                        PtrObj(lgCounter) = oCollision(lgCollision).Shp1
'                    End If
'
'                    bStack = True
'                    For lgPtr = LBound(PtrCollision) To UBound(PtrCollision)
'                        If oCollision(PtrCollision(lgPtr)).Shp1 = ...lgShp2 Then
'                            bStack = False
'                            Exit For
'                        End If
'                        If oCollision(PtrCollision(lgPtr)).Shp2 = ...lgShp2 Then
'                            bStack = False
'                            Exit For
'                        End If
'                    Next lgPtr
'                    If bStack Then
'                        ReDim Preserve PtrObj(g_Base To lgCounter)
'                        PtrObj(lgCounter) = oCollision(lgCollision).Shp2
'                    End If

End If
Next lgCollision

Stop
' Then process collisions against other particles
If Not (Not PtrCollision()) Then
'Create XYZ systems of equations for the momentum (call Gauss-Jordan solver):
ReDim mCollision(LBound(PtrCollision) To UBound(PtrCollision), _
LBound(PtrCollision) To UBound(PtrCollision) + 1)

ReDim PtrObj(LBound(PtrCollision) To UBound(PtrCollision))
For lgCollision = LBound(PtrCollision) To UBound(PtrCollision)
Next lgCollision
' Sort elements by Id
'Call fQuickSort_ArrayLng(PtrObj())

'For X direction
For lgCollision = LBound(PtrCollision) To UBound(PtrCollision)
'............
'                    mMomentum(lgCollision).X = mMomentum(lgCollision).X _
'                                             + vectorShp(lgShp1).X * (DimShp(lgShp1).X + DimShp(lgShp1).Y) / 2 _
'                                             + vectorShp(lgShp2).X * (DimShp(lgShp2).X + DimShp(lgShp2).Y) / 2
'                    mMomentum(lgCollision).Y = mMomentum(lgCollision).Y _
'                                             + vectorShp(lgShp1).Y * (DimShp(lgShp1).Y + DimShp(lgShp1).Y) / 2 _
'                                             + vectorShp(lgShp2).Y * (DimShp(lgShp2).Y + DimShp(lgShp2).Y) / 2

'                    'Distance between shapes (lgShp1, lgShp2)
'                     DX = (CenterShp(lgShp1).X - CenterShp(lgShp2).X)
'                     DY = (CenterShp(lgShp1).Y - CenterShp(lgShp2).Y)
'
'                     If DX  0 Then CCAng = Atn(DY / DX) Else CCAng = Pi / 2

'                     With vectorShp(oCollision(lgCollision).Shp1)
'                         Angle_Shp1 = Atn(.Y / .X)
'                         Shp1_Speed = Sqr(.X ^ 2 + .Y ^ 2)
'                     End With
'                     With vectorShp(oCollision(lgCollision).Shp2)
'                         Angle_Shp2 = Atn(.Y / .X)
'                         Shp2_Speed = Sqr(.X ^ 2 + .Y ^ 2)
'                     End With
'
'                     Angle_Shp1 = CCAng * 2 - Angle_Shp1
'                     Angle_Shp2 = CCAng * 2 - Angle_Shp2
'
'                     With vectorShp(oCollision(lgCollision).Shp1)
'                         .X = -Shp1_Speed * Cos(Angle_Shp1)
'                         .Y = Shp1_Speed * Sin(Angle_Shp1)
'                     End With
'                     With vectorShp(oCollision(lgCollision).Shp2)
'                         .X = Shp2_Speed * Cos(Angle_Shp2)
'                         .Y = -Shp2_Speed * Sin(Angle_Shp2)
'                     End With
'............
Next lgCollision
'Call fGaussJordan(mMomentum)

'For Y direction...
'...
'For Z direction...
'...
End If

Start = VBA.Timer()
Do While VBA.Timer() < (Start + TimeStep) 'TimeCollision
DoEvents
Loop
Loop
End With
End Sub

'Public Function fGaussJordan()
'End Function
'Public Function fQuickSort_ArrayLng()
'End Function

Public Sub NCradle()
Dim Plength As Single
Dim StartAng As Single
Dim NumBalls As Long
Dim StartA(1 To 5, 1 To 4) As Single
Dim StringA As Variant
Dim BallA As Variant
Dim TimeStep As Double
Dim i As Long
Dim Step As Single
Dim Level As Double
Dim StartLevel As Double
Dim SRotn As Single
Dim Start As Double
Dim Start2 As Double
Dim NextAng As Single
Dim AngA() As Double
Dim Taccn As Double
Dim V_1 As Double
Dim V_2 As Double
Dim Vav As Double
Dim Omav As Double
Dim Period As Double
Dim NumSteps As Long

Const BallR As Single = 25
Const StringTop As Single = 100
Const StringLength As Single = 200
Const String1X As Single = 250

Const g As Double = 9.8

Plength = Range("Plength").Value / 1000
StartAng = Range("StartAngle").Value * PI / 180
NumBalls = Range("NumNC").Value
StartLevel = StringLength / 10 * (1 - Cos(StartAng))

Period = 2 * PI * (Plength / g) ^ 0.5 * (1 + Sin(StartAng / 2) ^ 2 / 4 + Sin(StartAng / 2) ^ 4 * 9 / 64)
Range("period").Value = Period
NumSteps = Period / 0.05
TimeStep = Period / NumSteps
ReDim AngA(0 To NumSteps, 1 To 3)

AngA(0, 3) = -g * (Sin(StartAng))
AngA(0, 2) = TimeStep * AngA(0, 3) / 2
AngA(0, 1) = StartAng + AngA(0, 2) * TimeStep

For i = 1 To NumSteps
Taccn = -g * (Sin(AngA((i - 1), 1)))
AngA(i, 3) = Taccn
V_1 = AngA(i - 1, 2)
V_2 = V_1 + TimeStep * (Taccn * 1.5 - AngA((i - 1), 3) / 2)
AngA(i, 2) = V_2
Vav = (V_1 + V_2) / 2
Omav = Vav / Plength
AngA(i, 1) = AngA(i - 1, 1) + Omav * TimeStep
Next i

StringA = Array("NcLine1", "NcLine2", "NcLine3", "NcLine4", "NcLine5")
BallA = Array("Ncradle1", "Ncradle2", "Ncradle3", "Ncradle4", "Ncradle5")

Do
NextAng = StartAng
For Step = 1 To NumSteps
Start = Timer
For i = 1 To 5
StartA(i, 1) = String1X + (i - 1) * 2 * BallR
StartA(i, 2) = StringTop
If ((Step  NumSteps * 3 / 4) And i  NumSteps / 4 And Step  NumBalls) Then
SRotn = 0
StartA(i, 3) = StartA(i, 1)
StartA(i, 4) = StartA(i, 2) + StringLength
Else
SRotn = NextAng
StartA(i, 3) = StartA(i, 1) + StringLength * Sin(SRotn)
StartA(i, 4) = StartA(i, 2) + StringLength * Cos(SRotn)

End If

ActiveSheet.Shapes(StringA(i - 1)).Delete
With ActiveSheet.Shapes.AddLine(StartA(i, 1), StartA(i, 2), StartA(i, 3), StartA(i, 4))
.Name = StringA(i - 1)
.Line.Weight = 2
End With

With ActiveSheet.Shapes(BallA(i - 1))
.Left = StartA(i, 3) - BallR
.Top = StartA(i, 4) - BallR
.Width = BallR * 2
.Height = .Width
End With
Next i

If Step = Round(NumSteps / 4, 0) + 1 Or Step = Round(NumSteps * 3 / 4, 0) + 1 Then Beep

NextAng = AngA(Step, 1)
Level = StringLength / 10 * (1 - Cos(NextAng))

Do While Timer < Start + TimeStep
DoEvents
Loop
Next Step
Loop
End Sub
```