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
 

Leave a Reply

Your email address will not be published. Required fields are marked *