Excel ASCII Art

ASCII_artAside from loading an image on a picture object, any image (supported formats “*.bmp;*.jpg;*.jpeg;*.gif”) can be represented in cells, like in ASCII art or colouring cells interior. For ASCII art, from 16M color (long) it should be reduced to 256 gray scale color, and finally converted to a 16 text scale (using in decreasing darkness the following characters: M # @ H X $ % + / ; : = – , .). It would be better if used a four reduced scale with Unicode symbols (9617 9618 9619 and 160 for whitest) Following code achieves both, colouring interior functionality is deactivated in code because it’s very sloooooow.
Option Explicit

Private Type tRGB
    R As Byte
    G As Byte
    B As Byte
End Type

Private Type OPENFILENAME
    lStructSize         As Long
    hwndOwner           As Long
    hInstance           As Long
    lpstrFilter         As Long
    lpstrCustomFilter   As Long
    nMaxCustFilter      As Long
    nFilterIndex        As Long
    lpstrFile           As Long
    nMaxFile            As Long
    lpstrFileTitle      As Long
    nMaxFileTitle       As Long
    lpstrInitialDir     As Long
    lpstrTitle          As Long
    flags               As Long
    nFileOffset         As Integer
    nFileExtension      As Integer
    lpstrDefExt         As Long
    lCustData           As Long
    lpfnHook            As Long
    lpTemplateName      As Long
End Type

Private Declare Function GetOpenFileName Lib "comdlg32.dll" Alias "GetOpenFileNameW" (pOpenfilename As OPENFILENAME) As Long
Private Declare Function GetDIBits Lib "gdi32" (ByVal aHDC As Long, _
                                                ByVal hBitmap As Long, _
                                                ByVal nStartScan As Long, _
                                                ByVal nNumScans As Long, _
                                                lpBits As Any, _
                                                lpBI As Any, _
                                                ByVal wUsage As Long) As Long
Private Declare Function GetDC Lib "user32" (ByVal hwnd As Long) As Long
Private Declare Function ReleaseDC Lib "user32" (ByVal hwnd As Long, _
                                                 ByVal hdc As Long) As Long

Public Sub RunAsciiArt2()

End Sub

Public Sub RunAsciiArt()
    Dim fileName As String
    Dim pic     As StdPicture
    Dim lgPos   As Long
    Dim bi()    As Long
    Dim pix()   As Long
    Dim hdc     As Long
    Dim ret()   As Byte
    Dim y       As Long
    Dim x       As Long
    Dim pal()   As Byte
    Dim Unipal()    As Integer
    Dim cRGB    As tRGB
    Dim bConversion As Long
    Dim lgColor     As Long
    Dim aOutput()   As Variant
    Dim rgOutput    As Excel.Range

    'On Error GoTo ERRLABEL

    fileName = GetFile()

    If VBA.Trim$(fileName)  vbNullString Then 'Better... if FileExists
        Set pic = LoadPicture(fileName)
        hdc = GetDC(Application.hwnd)

        ReDim bi(1063): bi(0) = 40
        GetDIBits hdc, pic.Handle, 0, 0, ByVal 0&, bi(0), 0

        ReDim pix(bi(1) - 1, Abs(bi(2)) - 1)

        If bi(2) > 0 Then bi(2) = -bi(2): bi(3) = &H200001: bi(4) = 0

        GetDIBits hdc, pic.Handle, 0, Abs(bi(2)), pix(0, 0), bi(0), 0
        ReleaseDC Application.hwnd, hdc

        ReDim ret((bi(1) + 2) * Abs(bi(2)) * 2)

        ReDim pal(15)
        pal(0) = &H4D:  pal(1) = &H23:  pal(2) = &H40:  pal(3) = &H48
        pal(4) = &H58:  pal(5) = &H24:  pal(6) = &H25:  pal(7) = &H2B
        pal(8) = &H2F:  pal(9) = &H3B:  pal(10) = &H3A: pal(11) = &H3D
        pal(12) = &H2D: pal(13) = &H2C: pal(14) = &H2E: pal(15) = &H20

        ReDim Unipal(3)
        Unipal(0) = 9617: Unipal(1) = 9618: Unipal(2) = 9619: Unipal(3) = 160

        With wsCanvas
            Set rgOutput = .Range(.Cells(1, 1), .Cells(UBound(pix, 2) + 1, UBound(pix, 1) + 1))
            ReDim aOutput(1 To rgOutput.Rows.Count, 1 To rgOutput.Columns.Count)
        End With

        With Application
            .ScreenUpdating = False
            .Calculation = xlCalculationManual
        End With

        For y = 0 To UBound(pix, 2)
            For x = 0 To UBound(pix, 1)
                'https://stackoverflow.com/questions/687261/converting-rgb-to-grayscale-intensity
                'https://en.wikipedia.org/wiki/Grayscale
                bConversion = (pix(x, y) And &HFF) * 0.2989 + _
                              (pix(x, y) \ &H100 And &HFF) * 0.587 + _
                              (pix(x, y) \ &H10000 And &HFF) * 0.114

                ret(lgPos) = pal(bConversion \ &H10)

                'ASCII art output:
                'aOutput(y + 1, x + 1) = VBA.Chr(ret(lgPos))            'For ASCII palette
                aOutput(y + 1, x + 1) = CharW(Unipal(bConversion \ 64)) 'For Block Unicode palette

                'Change from Pic BGR to RGB
                'lgColor = pix(x, y)
                'cRGB.R = lgColor And &HFF&
                'cRGB.G = (lgColor \ &H100&) And &HFF&
                'cRGB.B = lgColor \ &H10000
                'rgOutput(y + 1, x + 1).Interior.Color = VBA.RGB(cRGB.B, cRGB.G, cRGB.R)

                lgPos = lgPos + 2
            Next x

            ret(lgPos) = 13
            ret(lgPos + 2) = 10
            lgPos = lgPos + 4
        Next y

        'Open fileName & ".txt" For Binary As #1
        'Put #1, , ret
        'Close #1
        Erase ret()

        'ASCII art:
        rgOutput.Value2 = aOutput()

        With Application
            .ScreenUpdating = True
            .Calculation = xlCalculationAutomatic
        End With

        'Fit to selection
        With rgOutput
            .Columns.ColumnWidth = 1.43
            .Rows.RowHeight = 11.25
            .Select
            ActiveWindow.Zoom = True
        End With
    End If

    Exit Sub

ERRLABEL:
    MsgBox "Error", vbCritical
End Sub

Private Function GetFile() As String
    Dim ofn As OPENFILENAME
    Dim Out As String
    Dim i As Long

    ofn.nMaxFile = 260
    Out = String(260, vbNullChar)
    ofn.hwndOwner = Application.hwnd 'Application.hwnd
    ofn.lpstrTitle = StrPtr("Open image")
    ofn.lpstrFile = StrPtr(Out)
    ofn.lStructSize = Len(ofn)
    ofn.lpstrFilter = StrPtr("Supported image format" & vbNullChar & "*.bmp;*.jpg;*.jpeg;*.gif" & vbNullChar)
    If GetOpenFileName(ofn) Then
        i = InStr(1, Out, vbNullChar, vbBinaryCompare)
        If i Then GetFile = Left$(Out, i - 1)
    End If

End Function

Public Function CharW(ByVal dec As Long) As String
'https://en.wikipedia.org/wiki/Box-drawing_character
' Box drawings chars: https://www.unicode.org/charts/PDF/U2500.pdf --> from 9472 to 9599
' Blocks chars: https://www.unicode.org/charts/PDF/U2580.pdf       --> from 9600 to 9631
' https://www.vertex42.com/ExcelTips/unicode-symbols.html

    CharW = ChrW(dec)
End Function
 

3D imagery

From gis.stackexchange.com, there is this interesting thread about UAV and 3D mapping. Most of the information come from the answer of user ragnvald. This post is a combination of the other anwers with the useful information. What software to use for making digital elevation models from UAV aerial imagery? To solve the problem one needs to transform 2D images of 3D structures from different angles/perspectives into a solid model. This was earlier a manual job, but software allows for automated processes. Products from such processes can include Digital Elevation Models as well as Point Clouds, Digital Surface Models, Textured Digital Surface Models, Orthorectified Imagery and Classified Point Clouds. This is known as Structure from Motion analysis. Processing can be done locally, through online services or a combination of the two. The software is basically the same, although one will find that local processing allows for more specialized settings. Local processing nwiDU
  • pix4D is an overall good tool providing exceptional ease of use for inexperienced to more experienced users. In addition to provide a desktop tool it also integrates with a cloud based service. A youtube video introduces the software in a good way. It is a bit in the high end with re to price. you can get a trial version for 15 days. Here is the link for how to get a trial pix4d version. Moreover you would partially be able to control the processing steps using different parameters. Here is the link for the processing steps.
  • Agisoft PhotoScan is available for Windows, OSX and Linux. Agisoft is flexible on the platform side. Having tried it out it produces good terrain models and ortophotos. A feature which lacks on Pix4D is the ability to set up batch jobs. You can also script the processing using a python API.
  • Drone2Map from ESRI can create orthomosaics, 3D meshes and more. This is not an integrated ArcGIS desktop tool, but effectively a stand alone 64 bit application. On their webpage they state that: “Drone2Map for ArcGIS is powered by Pix4D.” Opposed to Pix4D the information provided to the user while calculating is poor. From what I understand this is Pix4D under the hood under a different licensing model.
  • OpenDroneMap runs on command line Ubuntu and is a project which will let you create Digital Elevation Models as well as the restof the products mentioned in the beginning of this answer. The latest addition to this project is a web front end. You can find the web front end, the main project and more code on github.
  • ERDAS Imagine have an UAV add-in module. Combined it’s more expensive than Agisoft Photoscan, but if you’re already an ERDAS user it integrates really well.
  • 3Dsurvey in addition to imagery analysis also offers what looks like good point cloud editing tools.
  • 3DFier. Delft University is doing a lot of research into CityGML, one of the upcoming standard formats for visualising 3D building information. They have a viewer and also this tool, to turn LIDAR info into CityGML. So, opensource, and supports a standard.
  • openMVG This package is used for inital structure from motion, camera calibration, feature detection and image orientation.
  • openMVS This package is used to generate a dense point cloud, dense textured mesh and a colorized mesh from initial OpenMVG results. Run OpenMVG output through OpenMVS to generate dense results and products. For DJI drones and sensors, they have a free imagery processing tool RAPID for DJI
Online services NBKCM Processing of aerial imagery is resource demanding and will require some serious hardware. Online services are therefore pay-as-you-go or tied to a license.
    • Dronemapper is an online service where you can upload imagery and have it processed.
    • Dronedeploy was initially made as a planning and processing framework for DJI drones. It now supports creating maps and 3D models for any drone imagery. One can also do analysis based on the imagery. Produces good maps. Their app is also a marketplace where you can install apps for free.
    • MicaSense MicaSense Atlas is a cloud-based data platform for processing, storage, management, presentation, and analytics of multispectral data captured with professional multispectral cameras like the MicaSense RedEdge and Parrot Sequoia.
    • Maps Made Easy is another provider of online processing and data management.
Combined services pix4D also provides an online processing service for users with a desktop license. Apps Apps seem to be more and more popular.
  • They usually serve several purposes.
  • Help the user to plan surveys
  • Take pictures according to the plan
  • Upload survey plans to the drone
  • Provide updated information from the drone during a survey
  • Help the user to adjust plans during a survey based on earlier flights
  • Facilitate for upload of data so that data can be processed using an online service.
Under this section we can recommend this two apps:
  • DroneDeploy has an app which has an extraordinary update rate. The app is very flexible and also has an option to use plugins to focus the survey effort. Recent (July 2017) updates fixed issues with image spacing, but at the same time put limitations on the flight speed.
  • Pix4D Capture provides a stable work tool for collecting aerial images. It is probably the best available. On the downside the app has some shortcomings (unclear exposure settings, accurate control of speed, manual and not flexible way of controlling flight direction, lacking upload of survey areas (only android has this) and more) which should be fixed. It is rarely updated.
Obsolete
  • VisualSFM is a software package which in combination with CSVS might be a way to go. A youtube movie on the www.flightriot.com webpage gives a practical example of potential endpoint products. The resulting products are not referenced and as such not useful for GIS work. Visual SFM seems to be a project which lacks momentum. As far as I can see there has not been any development on the software for a couple of years.
For any fast development, and preconsiderations, can recall to already computed DEM models. The big semi-global DEMs out there are: These are relatively raw data products, with some holes and artifacts. SRTM, at least, has been around quite a while now, and has been smoothed out and refined into national-scale production-ready data products by a number of organizations. The National Elevation Dataset is the official US DEM: http://ned.usgs.gov/

Neuronal networking in VBA+Excel

Not for a long time, I’m getting more and more engaged about neuronal networks. And as usual, I try to tame it with the help of Excel as the main IDE. Considering that they are neither the best application, nor the best suited language to do this kind of investigations, it seems there is a hard work to do here. I have come to some examples over the Net of people doing some incredible things, showing myself again how, when you dare to get something done, you get it done (… for sure). Lets start with some old articles and developments on the subject:
  • an spanish article (2002) about the creation of workbooks for ANN.
  • Phil Brierleys’s site, with TiberiusXL app (2003), about ANN and GAs.
  • An article from Pf. Cliff T. Ragsdale and Christopher W. Zobel (2010), tittled “A Simple Approach to Implementing and Training Neural Networks in Excel”, that I could not find open on the Net, but would love to read as I get my hands on it.
  • A not very descriptive article about NN, with some basic code.
  • Here is another article, it states Excel use, but can not see any clue of where it does (seems to be C++ related).
  • The site of Angshuman Saha, has some developments on data mining, with 3 files:
    • CTree – Builds Classification Tree (using a version of C4.5 algorithm).
    • NNClass – Builds Classification Models (using feedforward-backpropagation Neural Network).
    • NNPred – Builds Prediction Models (using feedforward-backpropagation Neural Network).
  • An StackOverflow thread about the topic.
In GitHub there are some projects too: If you need some explanation about Backpropagation, there is a profuse one by Matt Mazur here, althought his code has some errors feeding h1 and h2 values (net_h1= w1×i1 + w3×i3 + b1 (analog mistake for net_h2)). Also very illustrative on this topic the videos from Ryan Harris (I get there from Scott Turner‘s videos). And here come the two most incredible ones. As a VBA developer, having done some exceptional Works myself, I must surrender to this two, really impressive. They amuse me more as anything related to ANN seem to be magic, and that magic done in Excel looks astonishing:
  • Richard Maddison‘s creations of ANN. If you need further explanations, look in his website. Look at this one for example…
  • David Bots‘s creations of ANN. A sample for this autor:
One of the clearest explanations on how to implement ANN in Excel comes from the quatmacro blog, that has some links to Excellaneous‘s site where you can download some Excel related files… and some of them linked to FoxesTeam. Some more resources, via external add-ins:
  • XLSTAT-R engine, software for statistical analysis in Excel, for Excel
  • http://www.wardsystems.com/manuals/genehunter/index.html?examples_neural_net_for_bob_s_deli.htm
  • https://www.cheshireeng.com/Neuralyst/nrldemo.htm
  • http://www.deepexcel.net/
  • And here is a site that links to a plethora of add-ins for ANN in Excel.
  • And not only in Excel, there are script sub-sets to deal with ANN.
  • The Water Systems research group at the University of Adelaide School of Engineering has developed an add-in to implement ANNs for water resources modeling applications (such as flow forecasting, water quality forecasting, water treatment process modeling).
This posts is really informative also. Finally, I would recomend to take a look to this post analysis of top trending software for NN development, and to the master class on YouTube from Andrew Y. Ng (consider that the Youtube channel does not have all the Course videos that are show on the Andrew Ng course at Coursera). Other courses worth it are the CS231n lectures on Stanford and the MIT 6.034 course. You’ll not get as easy as a pick and go, but you will gain much confidence on how things operate.
Aleksey Teslenko on this Youtube video shows how to easily perform a NN on Excel, and recomended this book (only $4 on Amazon is a good choice). There I have linked two posts that show a little on the computations below (https://stevenmiller888.github.io/mind-how-to-build-a-neural-network/ and http://stevenmiller888.github.io/mind-how-to-build-a-neural-network-part-2/).

VBA bypass worksheet password protection

It’s quite common to run into worksheets that have been protected, by original author, so users can not modify formulas or select ranges. But it’s also more than often that oneself is in the bind to make modifications, so a bypass is needed. This is a popular shared code, that I’d modified a little, just to tidy and understand it a bit… and to extend functionality to the whole Workbook. It does not find original password, it just bypass the control protection using another one, 12 characters long, with the same hash as the unknown password. First 11 characters are either [A] or [B], and the last one is the real “free” one (can be any letter/number/special character -aside system reserved chars-). Not a lot of combinations to try: (127-32)·2^11 ~ 194.560
Public Sub sWorksheet_CrackUnprotect()
  Dim oWsh As Excel.Worksheet
  Dim aData(0 To 1) As String * 1
  Dim b1 As Byte, b2 As Byte, b3 As Byte
  Dim b4 As Byte, b5 As Byte, b6 As Byte
  Dim b7 As Byte, b8 As Byte, b9 As Byte
  Dim b10 As Byte, b11 As Byte, b12 As Long
  Dim strPassword As String

  On Error Resume Next
' Try with found Password (initial password = vbNullString)
  For Each oWsh In ThisWorkbook.Worksheets 
    With oWsh
      .Activate

      .Unprotect strPassword
      If .ProtectContents = False Then GoTo ByPass

      aData(0) = VBA.Chr(65)
      aData(1) = VBA.Chr(66)
      For b1 = 0 To 1: For b2 = 0 To 1: For b3 = 0 To 1
      For b4 = 0 To 1: For b5 = 0 To 1: For b6 = 0 To 1
      For b7 = 0 To 1: For b8 = 0 To 1: For b9 = 0 To 1
      For b10 = 0 To 1: For b11 = 0 To 1: For b12 = 32 To 126
        .Unprotect aData(b1) & aData(b2) & aData(b3) & _
                   aData(b4) & aData(b5) & aData(b6) & _
                   aData(b7) & aData(b8) & aData(b9) & _
                   aData(b10) & aData(b11) & VBA.Chr(b12)
        If .ProtectContents = False Then
          strPassword = aData(b1) & aData(b2) & aData(b3) & _
                        aData(b4) & aData(b5) & aData(b6) & _
                        aData(b7) & aData(b8) & aData(b9) & _
                        aData(b10) & aData(b11) & VBA.Chr(b12)
           GoTo Iterate
        End If
      Next: Next: Next
      Next: Next: Next
      Next: Next: Next
      Next: Next: Next
Bypass:
    End With
  Next oWsh

  On Error GoTo 0
End Sub
[/sourcecode]	

VBA Excel presentations

Ever wonder how can be performed an animation of a shape in Excel as it’s easily done in PowerPoint? Following code comes from BeyondExcel site, an excepcional place to learn new Excel tricks. Basically it performs the strech grow and spin effects for any shape, very cute. From here on, any other effects are easily achievable (show from direction, fade,…).
Public Sub Workbook_Open()
'   Description:Runs when workbook opens

    Dim n As Integer
    Dim oShp As Excel.Shape

    'Worksheets("Data").Activate
    Set oShp = ActiveSheet.Shapes(1) 'Selection.ShapeRange.Item(1)
    ActiveSheet.Range("A1").Select
    oShp.LockAspectRatio = False
    n = 5
    #If VBA7 Then
        n = 10 'is way faster
    #End If
    GrowShape oShp, n
    SpinShape oShp, n

End Sub

Public Function SpinShape(ByRef oShp As Excel.Shape, _
                          ByRef Step As Integer) As Boolean
'   Description:Expands a shape into view

'   Parameters: oShp       The shape to animate
'               Step        Larger #s animate faster
'                           Steps should divide 90 evenly

'   Example:    SpinShape ActiveSheet.Shapes("Logo"), 10

    Const PI As Double = 3.14159265358979

    Dim sng01 As Single: sng01 = PI / 180    '1 Degree in Radians

    Dim sgCenterX As Single     'Shape's center X coordinate
    Dim sgCenterY As Single     'Shape's center Y coordiante
    Dim sgWidth As Single       'Shape's width
    Dim sgHeight As Single      'Shape's height
    Dim lgRotate As Long        'Generic Counter for the loop

    With oShp
        .LockAspectRatio = False
       'Remember shape's original dimensions
        sgCenterX = .Width / 2 + .Left
        sgCenterY = .Height / 2 + .Top
        sgWidth = .Width
        sgHeight = .Height
        .Visible = True
       'Animation Loop
        For lgRotate = 0 To 360 Step Step
            .Width = sgWidth * Abs(Cos(lgRotate * sng01))
            .Left = sgCenterX - .Width / 2
            If lgRotate = 90 Or lgRotate = 270 Then .Flip msoFlipHorizontal
            DoEvents
        Next lgRotate
       'Restore shape's original dimensions
        .Width = sgWidth
        .Height = sgHeight
        .Left = sgCenterX - .Width / 2
        .Top = sgCenterY - .Height / 2
    End With

End Function

Public Function GrowShape(ByRef oShp As Excel.Shape, _
                          ByRef Step As Integer) As Boolean
'   Description:Expands a shape into view

'   Parameters: oShp       The shape to animate
'               Step        Larger #s animate faster

'   Example:    GrowShape ActiveSheet.Shapes("Logo"), 10

'   Note:       For best results, shape should be hidden before calling this routine

    Dim sgCenterX As Single    'Shape's center X coordinate
    Dim sgCenterY As Single    'Shape's center Y coordiante
    Dim sgWidth As Single      'Shape's width
    Dim sgHeight As Single     'Shape's height
    Dim lgAngle As Long        'Generic Counter for the loop

    With oShp
        ' Remember shape's original dimensions
        sgCenterX = .Width / 2 + .Left
        sgCenterY = .Height / 2 + .Top
        sgWidth = .Width
        sgHeight = .Height
        .Visible = True
        ' Animation Loop
        For lgAngle = 0 To VBA.CLng(sgWidth) Step Step
            .Width = lgAngle
            .Height = lgAngle * sgHeight / sgWidth
            .Left = sgCenterX - .Width / 2
            .Top = sgCenterY - .Height / 2
            DoEvents
        Next lgAngle
        ' Restore shape's original dimensions
        .Width = sgWidth
        .Height = sgHeight
        .Left = sgCenterX - .Width / 2
        .Top = sgCenterY - .Height / 2
    End With
End Function

[/sourcecode]	

VBA Linear regression vs Ramer-Douglas-Peucker_algorithm

It’s a bit annoying that wherever I look for it, linear regression is always beeing formulated as the minimization function of the vertical distances sum. This covers the casuistic where data is prominently in an horizontal distribution, so it has the best chance to assure that if range is below π/4, it will get the best line that fits the data. But, how about if data is spread in the half sector above π/4. Then the minimization function should look for the horizontal distance. Lets take a look for an ill conditioned dataset as the one pictured below: illConditionedDataOr figure a random set of points inside a vertical oriented ellipse (there is also some code to draw the ellipse perimeter). This kind of data set may never get a chance to be fitted if looking for vertical distance minimization. To solve this I’ve come to this code that solves the lack of criteria. It still can not do anything against three of the Anscombe’s quartet, but, you know, there should not be tried a linear regression against any ill conditioned data. And finally, but not the least, some code to compute Ramer-Douglas-Peucker (RDP) approximation to a dataset, very handy to reduce the number of points. I’d already developed an algorithm, not so CPU demanding, that also achieves similar reduction, but RPD is interesting as it’s not constrained by minimal length criteria… but it takes a life to compute some big sets of data.
Option Explicit

Public Const PI As Double = 3.14159265358979
Public Const PI_RAD As Double = PI / 180
Public Const EULER As Double = 2.71828182845905
Public Const g_BASE As Long = 0

Private Type tXYZ
    X As Double
    Y As Double
    Z As Double
End Type

Public Sub sLinearRegression()
    Dim oPoint() As tXYZ
    Dim lgItem As Long
    Dim aData As Variant
    Dim dbR²ToX As Double
    Dim dbR²ToY As Double
    Dim aRegressionV() As Double
    Dim aRegressionH() As Double

    aData = Selection.Value2
    ReDim oPoint(LBound(aData, 1) To UBound(aData, 1))
    For lgItem = LBound(aData, 1) To UBound(aData, 1)
        oPoint(lgItem).X = aData(lgItem, 1)
        oPoint(lgItem).Y = aData(lgItem, 2)
    Next lgItem

    aRegressionV() = fLinearRegression(oPoint())

    Erase oPoint()
    Erase aData
End Sub

Public Function fLinearRegressionRange(ByRef rgRangeData As Excel.Range) As Variant
    Dim oPoint() As tXYZ
    Dim lgItem As Long
    Dim aData As Variant

    aData = rgRangeData.Value2
    ReDim oPoint(LBound(aData, 1) To UBound(aData, 1))
    For lgItem = LBound(aData, 1) To UBound(aData, 1)
        oPoint(lgItem).X = aData(lgItem, 1)
        oPoint(lgItem).Y = aData(lgItem, 2)
    Next lgItem

    fLinearRegressionRange = fLinearRegression(oPoint())

    Erase oPoint()
    Erase aData
End Function

Public Function fLinearRegression(ByRef oPoint() As tXYZ) As Double()
' http://mathworld.wolfram.com/CorrelationCoefficient.html
' Check the Anscombe's quartet (ill conditioned data)
    'Y = ß0 + (ß1 * X)
    Dim lgItem As Long
    Dim lgPoints As Long
    Dim Xmean As Double
    Dim Ymean As Double
    Dim SSxx As Double
    Dim SSyy As Double
    Dim SSxy As Double
    'Dim SSR As Double
    'Dim SSE As Double
    Dim aReturnV(1 To 7) As Double
    Dim aReturnH(1 To 7) As Double

    For lgItem = LBound(oPoint) To UBound(oPoint)
        Xmean = Xmean + oPoint(lgItem).X
        Ymean = Ymean + oPoint(lgItem).Y
    Next lgItem
    lgPoints = (UBound(oPoint) - LBound(oPoint) + 1)
    Xmean = Xmean / lgPoints
    Ymean = Ymean / lgPoints

    SSxx = -(Xmean * Xmean) * lgPoints
    SSyy = -(Ymean * Ymean) * lgPoints
    SSxy = -(Xmean * Ymean) * lgPoints
    For lgItem = LBound(oPoint) To UBound(oPoint)
        ' For Correlation Coefficient
        SSxx = SSxx + (oPoint(lgItem).X * oPoint(lgItem).X)
        SSyy = SSyy + (oPoint(lgItem).Y * oPoint(lgItem).Y)
        SSxy = SSxy + (oPoint(lgItem).X * oPoint(lgItem).Y)
    Next lgItem

    ' Find SSR (sum of squared residuals, the quantity to minimize) and SSE (sum of squared errors)
    aReturnV(1) = SSxy / SSxx                    ' Slope --> ß1
    aReturnV(2) = Ymean - (aReturnV(1) * Xmean)  ' Origin --> ß0
    aReturnV(3) = SSxy ^ 2 / (SSxx * SSyy)       ' Correlation Coefficient --> R²
    aReturnV(4) = Xmean
    aReturnV(5) = Ymean
    aReturnV(6) = SSyy * (1 - aReturnV(3))   ' SSR = SSyy * (1 - R²)
    aReturnV(7) = SSyy - aReturnV(6)         ' SSE = ß1 * SSxy = ß1² * SSxx = SSyy - SSR

    ' Rotate the minimizer function 90º, to check ill conditioned data
    aReturnH(1) = SSxy / SSyy                    ' Slope --> ß1'
    aReturnH(2) = Xmean - (aReturnH(1) * Ymean)  ' Origin --> ß0'

    ' Get inverse correlation:
    aReturnH(1) = 1 / aReturnH(1)                ' ß1 = 1/ß1'
    aReturnH(2) = -aReturnH(2) * aReturnH(1)     ' ß0 = -ß0'* ß1

    aReturnH(3) = SSxy ^ 2 / (SSxx * SSyy)       ' Correlation Coefficient --> R²
    aReturnH(4) = Xmean
    aReturnH(5) = Ymean
    aReturnH(6) = SSxx * (1 - aReturnH(3))   ' SSR = SSxx * (1 - R²)
    aReturnH(7) = SSxx - aReturnH(6)         ' SSE = ß1 * SSxy = ß1² * SSyy = SSxx - SSR

    If aReturnH(6) > aReturnV(6) Then
        fLinearRegression = aReturnV()
    Else
        fLinearRegression = aReturnH()
    End If

    Erase aReturnV()
    Erase aReturnH()
End Function

Public Function fEllipseXL(ByRef aSemiAxis As Double, _
                           ByRef bSemiAxis As Double, _
                           ByRef Ø As Double, _
                           ByRef Xc As Double, _
                           ByRef Yc As Double, _
                           Optional ByRef ß As Double, _
                           Optional ByRef bFull As Boolean) As Double()
    Dim oCenter As tXYZ
    Dim oPoint() As tXYZ
    Dim aData() As Double
    Dim lgPoint As Long

    With oCenter
        .X = Xc
        .Y = Yc
    End With

    'oPoint() = fEllipse(aSemiAxis, bSemiAxis, Ø, oCenter, ß, bFull)
    oPoint() = fEllipsePolar(aSemiAxis, bSemiAxis, Ø, oCenter, ß, bFull)
    ReDim aData(LBound(oPoint) To UBound(oPoint), 1 To 2)
    For lgPoint = LBound(oPoint) To UBound(oPoint)
        aData(lgPoint, 1) = oPoint(lgPoint).X
        aData(lgPoint, 2) = oPoint(lgPoint).Y
    Next lgPoint

    fEllipseXL = aData()
End Function

Public Function fEllipsePolar(ByRef aSemiAxis As Double, _
                              ByRef bSemiAxis As Double, _
                              ByRef Ø As Double, _
                              ByRef oCenter As tXYZ, _
                              Optional ByRef ß As Double = 0, _
                              Optional ByRef bFull As Boolean = False) As tXYZ()
' For polar angles from center

    Dim lgAngle As Long
    Dim lgPoint As Long
    Dim dbRadius As Double
    Dim dbAngle As Double
    Dim oPoint() As tXYZ
    Dim oPointTmp As tXYZ
    Dim dbMajorAxis As Double
    Dim eccentricity As Double

    If bFull Then
        ' Fill initial values
        lgPoint = g_BASE - 1

        ' Select case major semiaxis
        If aSemiAxis >= bSemiAxis Then
            eccentricity = VBA.Sqr((aSemiAxis * aSemiAxis) - (bSemiAxis * bSemiAxis)) / aSemiAxis

            For lgAngle = 0 To 360 Step 10
                lgPoint = lgPoint + 1
                ReDim Preserve oPoint(0 To lgPoint)
                dbAngle = VBA.CDbl(lgAngle * PI_RAD)
                dbRadius = bSemiAxis / VBA.Sqr(1 - (eccentricity * VBA.Cos(dbAngle) ^ 2))
                oPoint(lgPoint).X = dbRadius * VBA.Cos(dbAngle)
                oPoint(lgPoint).Y = dbRadius * VBA.Sin(dbAngle)
            Next lgAngle

        ElseIf aSemiAxis < bSemiAxis Then
            eccentricity = VBA.Sqr((bSemiAxis * bSemiAxis) - (aSemiAxis * aSemiAxis)) / bSemiAxis

            For lgAngle = 0 To 360 Step 10
                lgPoint = lgPoint + 1
                ReDim Preserve oPoint(0 To lgPoint)
                dbAngle = VBA.CDbl(lgAngle * PI_RAD)
                dbRadius = aSemiAxis / VBA.Sqr(1 - (eccentricity * VBA.Cos(dbAngle) ^ 2))
                oPoint(lgPoint).X = dbRadius * VBA.Cos(dbAngle)
                oPoint(lgPoint).Y = dbRadius * VBA.Sin(dbAngle)
            Next lgAngle
        End If

        For lgPoint = LBound(oPoint) To UBound(oPoint)
            ' For a rotated ellipse:
            oPointTmp.X = (oPoint(lgPoint).X * Cos(Ø)) - (oPoint(lgPoint).Y * Sin(Ø)) + oCenter.X
            oPointTmp.Y = (oPoint(lgPoint).X * Sin(Ø)) + (oPoint(lgPoint).Y * Cos(Ø)) + oCenter.Y
            oPoint(lgPoint) = oPointTmp
        Next lgPoint
    Else
        ReDim oPoint(0)
    End If

    fEllipsePolar = oPoint()

End Function

Public Function fEllipse(ByRef aSemiAxis As Double, _
                         ByRef bSemiAxis As Double, _
                         ByRef Ø As Double, _
                         ByRef oCenter As tXYZ, _
                         Optional ByRef ß As Double = 0, _
                         Optional ByRef bFull As Boolean = False) As tXYZ()
' https://en.wikipedia.org/wiki/Ellipse
' Cannonical equation A·X² + B·X·Y + C·Y² + D·X + E·Y + F = 0
' canonical implicit equation: Xcan²/a² + Ycan²/b² = 1
' Xcan = (X - oCenter.X) * Cos(Ø) + (Y - oCenter.Y) * Sin(Ø)
' Ycan = -(X - oCenter.X) * Sin(Ø) + (Y - oCenter.Y) * Cos(Ø)
' Parametric form: (x,y)=(a*cos t,b*sin t) 0 <= t < 2*pi
' A = aSemiAxis² * (Sin(Ø))² + bSemiAxis² * (Cos(Ø))²
' B = 2 * (bSemiAxis² - aSemiAxis²) * Sin(Ø) * Cos(Ø)
' C = aSemiAxis² * (Cos(Ø))² + bSemiAxis² * (Sin(Ø))²
' D = -2·A·oCenter.X - B·oCenter.Y
' E = -B·oCenter.X - 2·C·oCenter.Y
' F = A·oCenter.X² + B·oCenter.X·oCenter.Y + C·oCenter.Y² - A²·B²
' ð = ((A·C) - (B²/4))·F + (B·E·D/4) - ((C·D²)/4) - ((A·E²)/4)

    Dim ð As Double 'discriminant for non degenerated ellipse
    Dim A As Double
    Dim B As Double
    Dim C As Double
    Dim D As Double
    Dim E As Double
    Dim F As Double
    Dim lgSegment As Long
    Dim lgPoint As Long
    Dim oPoint() As tXYZ
    Dim oPointTmp As tXYZ

    A = (aSemiAxis * Sin(Ø)) ^ 2 + (bSemiAxis * Cos(Ø)) ^ 2
    B = 2 * (bSemiAxis ^ 2 - aSemiAxis ^ 2) * Sin(Ø) * Cos(Ø)
    C = (aSemiAxis * Cos(Ø)) ^ 2 + (bSemiAxis * Sin(Ø)) ^ 2
    D = -2 * A * oCenter.X - B * oCenter.Y
    E = -B * oCenter.X - 2 * C * oCenter.Y
    F = (A * oCenter.X ^ 2) + (B * oCenter.X * oCenter.Y) + (C * oCenter.Y ^ 2) - (A * B) ^ 2

    ð = ((A * C) - (B ^ 2 / 4)) * F + (B * E * D / 4) - ((C * D ^ 2) / 4) - ((A * E ^ 2) / 4)

    If ð <= 0 Then ' not degenerated ellipse
        If bFull Then
            ReDim oPoint(0 To 40)
            ' Fill initial values
            lgPoint = -1
            For lgSegment = 0 To 19
                lgPoint = lgPoint + 1
                oPoint(lgPoint).X = -aSemiAxis + (aSemiAxis * lgSegment / 10)
                oPoint(lgPoint).Y = bSemiAxis * VBA.Sqr(1 - (oPoint(lgPoint).X / aSemiAxis) ^ 2)
            Next lgSegment

            lgPoint = lgPoint + 1
            oPoint(lgPoint).X = aSemiAxis
            oPoint(lgPoint).Y = 0

            For lgSegment = 19 To 1 Step -1
                lgPoint = lgPoint + 1
                oPoint(lgPoint).X = -aSemiAxis + (aSemiAxis * lgSegment / 10)
                oPoint(lgPoint).Y = -bSemiAxis * VBA.Sqr((1 - (oPoint(lgPoint).X / aSemiAxis) ^ 2))
            Next lgSegment

            lgPoint = lgPoint + 1
            oPoint(lgPoint) = oPoint(LBound(oPoint))

            For lgPoint = 0 To 40
                'For a rotated ellipse:
                oPointTmp.X = (oPoint(lgPoint).X * Cos(Ø)) - (oPoint(lgPoint).Y * Sin(Ø)) + oCenter.X
                oPointTmp.Y = (oPoint(lgPoint).X * Sin(Ø)) + (oPoint(lgPoint).Y * Cos(Ø)) + oCenter.Y
                oPoint(lgPoint) = oPointTmp
            Next lgPoint
        Else
            ReDim oPoint(0)
        End If

        fEllipse = oPoint()

    'ElseIf ð = 0 Then 'point ellipse
    '    ReDim oPoint(0)
    '    oPoint(0) = oCenter
    '    fEllipse = oPoint()

    End If

End Function

Private Function fDistanceToLine(ByVal px As Double, ByVal py As Double, _
                                 ByVal X1 As Double, ByVal Y1 As Double, _
                                 ByVal X2 As Double, ByVal Y2 As Double, _
                                 Optional ByRef t As Double) As Double
' Calculate the distance between the point and the segment.
' http://vb-helper.com/howto_distance_point_to_line.html
' https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line
' http://mathworld.wolfram.com/Point-LineDistance2-Dimensional.html
    Dim dX As Double
    Dim dY As Double

    dX = X2 - X1
    dY = Y2 - Y1
    If (dX = 0 And dY = 0) Then
    ' It's a point not a line segment.
        dX = px - X1
        dY = py - Y1
        near_x = X1
        near_y = Y1
        fDistanceToLine = VBA.Sqr(dX * dX + dY * dY)
        Exit Function
    End If

    ' Calculate the t that minimizes the distance.
    t = ((px - X1) * dX + (py - Y1) * dY) / (dX * dX + dY * dY)

    ' See if this represents one of the segment's end points or a point in the middle.
    If t  1 Then
        dX = px - X2
        dY = py - Y2
        near_x = X2
        near_y = Y2
    Else
        near_x = X1 + t * dX
        near_y = Y1 + t * dY
        dX = px - near_x
        dY = py - near_y
    End If

    fDistanceToLine = VBA.Sqr(dX * dX + dY * dY)
End Function

Public Function fDistance2DPoints(ByRef oPointA As tXYZ, ByRef oPointB As tXYZ) As Double
    Dim dX As Double
    Dim dY As Double

    dX = oPointA.X - oPointB.X
    dY = oPointA.Y - oPointB.Y

    fDistance2DPoints = VBA.Sqr((dX * dX) + (dY * dY))
End Function

Public Sub sDouglasPeucker()
    Dim lgItem As Long
    Dim lgR As Long
    Dim lgC As Long
    Dim aData As Variant
    Dim oPoint() As tXYZ
    Dim oPtFilter() As tXYZ
    Dim Threshold As Double

    aData = Selection.Value2
    ReDim oPoint(LBound(aData, 1) To UBound(aData, 1))
    For lgItem = LBound(aData, 1) To UBound(aData, 1)
        oPoint(lgItem).X = aData(lgItem, 1)
        oPoint(lgItem).Y = aData(lgItem, 2)
    Next lgItem

    Threshold = VBA.Val(VBA.InputBox("Threshold value", , 1))
    oPtFilter() = fDouglasPeucker(oPoint(), Threshold)

    'Show in Worksheet:
    With ActiveSheet
        lgR = Selection.Row - 1
        lgC = Selection.Column
        For lgItem = LBound(oPtFilter) To UBound(oPtFilter)
            lgR = lgR + 1
            .Cells(lgR, lgC + 2).Value2 = oPtFilter(lgItem).X
            .Cells(lgR, lgC + 3).Value2 = oPtFilter(lgItem).Y
        Next lgItem
    End With

    Erase oPoint()
    Erase aData
End Sub

Public Function fDouglasPeucker(ByRef oPoint() As tXYZ, _
                                Optional ByVal Threshold As Double = 0) As tXYZ()
' https://en.wikipedia.org/wiki/Ramer%E2%80%93Douglas%E2%80%93Peucker_algorithm
' Ramer–Douglas–Peucker algorithm (RDP)
' ToDo: Solve "crossing" segments (one/several vertex get trapped)...
' ToDo: if no more than one vertex get inside a segment, then if that vertex below threshold, can be deleted

    Dim oLine() As tXYZ
    Dim lgPoint As Long
    Dim lgSegment As Long
    Dim lgItem As Long
    Dim lgMove As Long
    Dim dbDistance As Double
    Dim dbDistanceMax As Double
    Dim dbCheck As Double
    Dim lgBreaker As Long
    Dim bClosest As Boolean
    Dim dbVectorMultiplier As Double
    Dim bBreaker As Boolean
    Dim bIterate As Boolean

    If Not (Not oPoint()) Then
        If Threshold = 0 Then
            fDouglasPeucker = oPoint()
            Exit Function
        End If

        ' Get K.P. for each point
        oPoint(LBound(oPoint)).Z = 0
        For lgPoint = (LBound(oPoint) + 1) To UBound(oPoint)
            oPoint(lgPoint).Z = oPoint(lgPoint - 1).Z + fDistance2DPoints(oPoint(lgPoint), oPoint(lgPoint - 1))
        Next lgPoint

        lgSegment = g_BASE + 1
        ReDim Preserve oLine(g_BASE + 0 To lgSegment)
        oLine(g_BASE + 0) = oPoint(LBound(oPoint))
        oLine(g_BASE + 1) = oPoint(UBound(oPoint))

        Do
            bIterate = False

            lgPoint = LBound(oPoint) ' initialze with first item
            For lgSegment = LBound(oLine) To UBound(oLine) - 1
                dbDistanceMax = Threshold
                bBreaker = False

                For lgPoint = (LBound(oPoint) + 1) To (UBound(oPoint) - 1) 'Avoid extremes, as they are already in the final set
                    dbDistance = fDistanceToLine(oPoint(lgPoint).X, oPoint(lgPoint).Y, _
                                                 oLine(lgSegment).X, oLine(lgSegment).Y, _
                                                 oLine(lgSegment + 1).X, oLine(lgSegment + 1).Y, _
                                                 dbVectorMultiplier)
                    ' First condition to apply only if is inside range: oLine(lgSegment).KP <= oNearestOnLine.KP = 0 And dbVectorMultiplier <= 1) Then
                        ' Second condition is to apply only if point is the closest to any other segment...
                        bClosest = True ' Initialize
                        For lgItem = LBound(oLine) To (UBound(oLine) - 1)
                            If lgItem  lgSegment Then ' avoid itself...
                            dbCheck = fDistanceToLine(oPoint(lgPoint).X, oPoint(lgPoint).Y, _
                                                      oLine(lgItem).X, oLine(lgItem).Y, _
                                                      oLine(lgItem + 1).X, oLine(lgItem + 1).Y, _
                                                      dbVectorMultiplier)
                                If (dbVectorMultiplier >= 0 And dbVectorMultiplier <= 1) Then ' Only if inside range...
                                    If dbCheck < dbDistance Then
                                        bClosest = False
                                        Exit For
                                    End If
                                End If
                            End If
                        Next lgItem

                        If bClosest Then
                            If dbDistanceMax  UBound(oPoint) Then Exit For ' only if points sorted...
            Next lgSegment
        Loop While bIterate

        fDouglasPeucker = oLine()
    End If
End Function

Public Function fDouglasPeuckerRange(ByRef rgRangeData As Excel.Range, _
                                     Optional ByVal Threshold As Double = 0) As Variant
    Dim lgItem As Long
    Dim aData As Variant
    Dim oPoint() As tXYZ
    Dim oPtFilter() As tXYZ

    aData = Selection.Value2
    ReDim Preserve oPoint(LBound(aData, 1) To UBound(aData, 1))
    For lgItem = LBound(aData, 1) To UBound(aData, 1)
        oPoint(lgItem).X = aData(lgItem, 1)
        oPoint(lgItem).Y = aData(lgItem, 2)
    Next lgItem

    Threshold = 1
    oPtFilter() = fDouglasPeucker(oPoint(), Threshold)
    ReDim aData(LBound(oPtFilter) To UBound(oPtFilter), 1 To 2)
    For lgItem = LBound(oPtFilter) To UBound(oPtFilter)
        aData(lgItem, 1) = oPtFilter(lgItem).X
        aData(lgItem, 2) = oPtFilter(lgItem).Y
    Next lgItem
    fDouglasPeuckerRange = aData

    Erase oPoint()
    Erase aData
End Function
[/sourcecode]	

VBA Procedure Builder with multiparameter

There are two special VBA events associated to macros or even shapes that could be build on automation. These are the Application.OnKey(Procedure), and the expression.OnAction (beeing expression an Excel.Shape, an Excel.FormsControl or even a CommandBar control). For such events any procedure in the Workbook can be called (has to be declared as Public). And they can even hold parameter(s), none/one/several. It’s very frustrating to deal with such a mess of [‘] and [“] characters, so it’s really come handy to have a call builder procedure. Just paste this code to any module:
Public Function fBuildCaller(ByVal bWorkbookName As Boolean, _
                             ByVal ProcName As String, _
                             ParamArray Args() As Variant) As Variant
' Function to build procedure with variable number of arguments
' Take care that if bWorkBookName = True, will be permanent linked to the Workbook
    Dim oItem As Variant
    Dim oSubItem As Variant
    Dim strDebug As String
    Dim lgRetVal As Long

    For Each oItem In Args
        If IsArray(oItem) Then
            For Each oSubItem In oItem
                strDebug = strDebug & " """ & oSubItem & ""","
            Next oSubItem
        Else
            strDebug = strDebug & " """ & oItem & ""","
        End If
    Next oItem
            
    If bWorkbookName Then
        lgRetVal = VBA.MsgBox("If bWorkbookName is set to True, will be permanently linked to Workbook, go with it?", _
                              vbYesNo + vbExclamation, "I N F O")
        If lgRetVal = vbNo Then bWorkbookName = False
    End If
    If strDebug = vbNullString Then
        strDebug = VBA.IIf(bWorkbookName, "'" & ThisWorkbook.Name & "'!", "") & _
                   "'" & ProcName & "'"
    Else
        strDebug = VBA.IIf(bWorkbookName, "'" & ThisWorkbook.Name & "'!", "") & _
                   "'" & ProcName & VBA.Mid$(strDebug, 1, Len(strDebug) - 1) & "'"
    End If
    fBuildCaller = strDebug
End Function

Public Function fBuildCaller2(ByVal ProcName As String, _
                              ParamArray Args() As Variant) As Variant
' Only working for PopUpMenus... better use fBuildCaller
' Has the advantage that "OnAction" is not linked to the WorkBook name
    Dim oItem As Variant
    Dim strDebug As String

    For Each oItem In Args
        strDebug = strDebug & Chr(34) & oItem + Chr(34) & ","
    Next

    If strDebug = vbNullString Then
        strDebug = ProcName
    Else
        strDebug = ProcName & "(" & VBA.Mid$(strDebug, 1, Len(strDebug) - 1) & ")"
    End If
    fBuildCaller2 = strDebug
End Function
And that’s it. Note1: as stated, parameter bWorkbookName is noticeable important. If the worksheets changes its name or it’s beeing copy/pasted, the item -with the event declaration-, the link to the original will be carried wherever it the item goes, and so Excel will complain on opening. Better of to set it to “False”. Note2: In fBuildCaller Arg variable can handle even an array, but it has to be exploded from 2D to 1D, as a linear array. Watch how the oSubItem works to deal with arrays. Note3: fBuildCaller2 is a similar procedure, but seems it’s only working with CommandBars, it has as advantange the better readibility of the procedure, but has its disadvantages too.

VBA Gauss-Jordan implementation

VBA has no implementation for array inversion, neither equations solver. So it comes very handy a Gauss-Jordan solver:
Public Function fGaussJordan(ByRef mArray() As Double) As Double()
    Dim lgR As Long
    Dim lgC As Long
    Dim lgPivot As Long
    Dim lgR_Homogenize As Long
    Dim dbTmp As Double
    Dim lgRetVal As Long
    Dim mArrayTmp() As Double
    Dim Nm As Integer
    
    On Error GoTo ErrControl
    
    Nm = UBound(mArray, 1) - LBound(mArray, 1) + 1
    ReDim mArrayTmp(LBound(mArray, 1) To UBound(mArray, 1), LBound(mArray, 2) To UBound(mArray, 2))
    
    ' Swap rows (if needed)
    If (mArray(0, 0) = 0) Then
        For lgR = LBound(mArray, 1) To UBound(mArray, 1)
            If (mArray(lgR, 0)  0) Then
                For lgC = LBound(mArray, 2) To UBound(mArray, 2)
                    mArrayTmp(0, lgC) = mArray(0, lgC)
                    mArray(0, lgC) = mArray(lgR, lgC)
                    mArray(lgR, lgC) = mArrayTmp(0, lgC)
                Next lgC
            End If
        Next lgR
    End If
    
    For lgPivot = LBound(mArray, 1) To UBound(mArray, 1)
        dbTmp = mArray(lgPivot, lgPivot)
        For lgC = LBound(mArray, 2) To UBound(mArray, 2)
            mArray(lgPivot, lgC) = mArray(lgPivot, lgC) / dbTmp
        Next lgC
        For lgR = LBound(mArray, 1) To UBound(mArray, 1)
            If (lgR = lgPivot) Then GoTo NextRow
            dbTmp = mArray(lgR, lgPivot)
            For lgR_Homogenize = LBound(mArray, 2) To UBound(mArray, 2)
                mArray(lgR, lgR_Homogenize) = mArray(lgR, lgR_Homogenize) - (dbTmp * mArray(lgPivot, lgR_Homogenize))
            Next lgR_Homogenize
NextRow:
        Next lgR
    Next lgPivot
    
    'Print solution
    ReDim mArrayTmp(LBound(mArray, 1) To UBound(mArray, 1))
    For lgR = LBound(mArray, 1) To UBound(mArray, 1)
        mArrayTmp(lgR) = mArray(lgR, Nm)
        'Debug.Print VBA.Format(mArray(lgR, Nm), "##,##0.00")
    Next lgR
    fGaussJordan = mArrayTmp()

ExitProc:
    Exit Function

ErrControl:
    lgRetVal = VBA.MsgBox("System has no solution", vbCritical)
End Function
[/sourcecode]	

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
 

VBA array functions (MatLabish/Pythonish implementation)

Nor VBA nor Visual Basic have natively implemented most of the array functions that Python or MatLab (https://www.mathworks.com/help/matlab/functionlist.html) have. But they can be coded to get similar functionality. Following are a bunch of functions (code not finished, or even not just started -for those that have a ‘!!!!! at the beginning of the description-), just to get fast creation and operation over matrices in VBA. Note: This is a work on progress, so it’ll grow in the future with new functions.
Option Explicit
'!!!!!!!!!!!!!!!!!!!
Public Const g_Base As Long = 0
'!!!!!!!!!!!!!!!!!!!
Public Function fNewArray(ByVal strText As String, _
                          Optional ByVal strColSeparator As String = " ", _
                          Optional ByVal strRowSeparator As String = ";", _
                          Optional ByVal strNewLine As String = "\") As Variant
'To reference all the elements in the mth row we type A(:m,).
'To reference all the elements in the nth column we type A(:,n).
'To reference all the elements in the mth to nth column we type A(:,m:n).
'To reference all the elements in the mth to nth row we type A(m:n,:).
'ToDo --> a([2,3,2,3],:)
     '--> Row Vector []
     '--> Column vector {}
    Dim mArray As Variant
     
    'Dim lgDim As Long
    Dim lgElement As Long
    Dim lgElements As Long
    'strEnclosing As String = "[]"
    Dim aVector() As String
    Dim aElement() As String
    Dim lgVector As Long
    Dim strVector As String
    Dim aCreator() As String
    
    strText = VBA.Trim$(strText)
    'Join lines...
    strText = VBA.Replace(strText, vbNewLine, "")
    strText = VBA.Replace(strText, strNewLine, "")
    
    aCreator() = VBA.Split(strText, ":")
    If LBound(aCreator)  UBound(aCreator) Then
    ' array = [first : second : ... : last]
        If (2 = (UBound(aCreator) - LBound(aCreator) + 1)) Then
            ReDim mArray(g_Base + 0 To g_Base + aCreator()(LBound(aCreator)), _
                         g_Base + 0 To g_Base + aCreator()(UBound(aCreator)))
        ElseIf (3 = (UBound(aCreator) - LBound(aCreator) + 1)) Then
            ReDim mArray(g_Base + 0 To g_Base + aCreator()(LBound(aCreator) + 0) - 1, _
                         g_Base + 0 To g_Base + aCreator()(LBound(aCreator) + 1) - 1, _
                         g_Base + 0 To g_Base + aCreator()(UBound(aCreator)) - 1)
        End If
    
    Else
        If strText Like "[[]*" Then
         '--> Row Vector []
            If strText Like "*]" Then
                strText = VBA.Mid$(strText, 2, VBA.Len(strText) - 2)
                strText = VBA.Trim$(strText)
                
                'Get vectors
                aVector = VBA.Split(strText, strRowSeparator)
                strVector = VBA.Trim$(aVector()(LBound(aVector)))
                
                'Avoid repeated separators
                Do While VBA.InStr(1, strVector, (strColSeparator & strColSeparator)) > 0
                    strVector = VBA.Replace(strVector, strColSeparator & strColSeparator, strColSeparator)
                Loop
                aElement = VBA.Split(strVector, strColSeparator)
                ReDim mArray(g_Base + LBound(aVector) To g_Base + UBound(aVector), _
                             g_Base + LBound(aElement) To g_Base + UBound(aElement))
                
                For lgVector = LBound(aVector) To UBound(aVector)
                    strVector = VBA.Trim$(aVector(lgVector))
                    
                    'Avoid repeated separators
                    Do While VBA.InStr(1, strVector, (strColSeparator & strColSeparator)) > 0
                        strVector = VBA.Replace(strVector, strColSeparator & strColSeparator, strColSeparator)
                    Loop
                    aElement = VBA.Split(strVector, strColSeparator)
                    For lgElement = LBound(aElement) To UBound(aElement)
                        mArray(g_Base + lgVector, g_Base + lgElement) = VBA.Val(aElement(lgElement))
                    Next lgElement
                Next lgVector
            End If
        
        ElseIf strText Like "{*" Then
         '--> Column vector {}
            If strText Like "*}" Then
                strText = VBA.Mid$(strText, 2, VBA.Len(strText) - 1)
            
                'Join lines...
                strText = VBA.Replace(strText, strNewLine, "")
                
                'Get Columns
                aVector = VBA.Split(strText, strColSeparator)
                strVector = VBA.Trim$(aVector()(LBound(aVector)))
                
                'Avoid repeated separators
                Do While VBA.InStr(1, strVector, (strRowSeparator & strRowSeparator)) > 0
                    strVector = VBA.Replace(strVector, strRowSeparator & strRowSeparator, strRowSeparator)
                Loop
                aElement = VBA.Split(strVector, strRowSeparator)
                ReDim mArray(g_Base + LBound(aElement) To g_Base + UBound(aElement), _
                             g_Base + LBound(aVector) To g_Base + UBound(aVector))
                
                For lgVector = LBound(aVector) To UBound(aVector)
                    strVector = VBA.Trim$(aVector(lgVector))
                    
                    'Avoid repeated separators
                    Do While VBA.InStr(1, strVector, (strRowSeparator & strRowSeparator)) > 0
                        strVector = VBA.Replace(strVector, strRowSeparator & strRowSeparator, strRowSeparator)
                    Loop
                    aElement = VBA.Split(strVector, strRowSeparator)
                    For lgElement = LBound(aElement) To UBound(aElement)
                        mArray(g_Base + lgVector, g_Base + lgElement) = VBA.Val(aElement(lgElement))
                    Next lgElement
                Next lgVector
            End If
        End If
    End If
    
    fNewArray = mArray
    Erase aCreator()
End Function

Public Function fNewArrayStr(ByVal strText As String, _
                             Optional ByVal strColSeparator As String = " ", _
                             Optional ByVal strRowSeparator As String = ";", _
                             Optional ByVal strNewLine As String = "\") As String()
    Dim vArray As Variant
    Dim aStr() As String
    Dim lgC As Long
    Dim lgR As Long
    
    vArray = fNewArray(strText, strColSeparator, strRowSeparator, strNewLine)
    If IsArray(vArray) Then
        ReDim aStr(LBound(vArray, 1) To UBound(vArray, 1), _
                   LBound(vArray, 2) To UBound(vArray, 2))
        For lgR = LBound(vArray, 1) To UBound(vArray, 1)
            For lgC = LBound(vArray, 2) To UBound(vArray, 2)
                aStr(lgR, lgC) = VBA.CStr(vArray(lgR, lgC))
            Next lgC
        Next lgR
        fNewArrayStr = aStr()
    
        Erase aStr()
    End If
    Erase vArray
End Function
Public Function fNewArrayDbl(ByVal strText As String, _
                             Optional ByVal strColSeparator As String = " ", _
                             Optional ByVal strRowSeparator As String = ";", _
                             Optional ByVal strNewLine As String = "\") As Double()
    Dim vArray As Variant
    Dim aDbl() As Double
    Dim lgC As Long
    Dim lgR As Long
    
    vArray = fNewArray(strText, strColSeparator, strRowSeparator, strNewLine)
    If IsArray(vArray) Then
        ReDim aDbl(LBound(vArray, 1) To UBound(vArray, 1), _
                   LBound(vArray, 2) To UBound(vArray, 2))
        For lgR = LBound(vArray, 1) To UBound(vArray, 1)
            For lgC = LBound(vArray, 2) To UBound(vArray, 2)
                aDbl(lgR, lgC) = VBA.Val(vArray(lgR, lgC))
            Next lgC
        Next lgR
        fNewArrayDbl = aDbl()
    
        Erase aDbl()
    End If
    Erase vArray
End Function
Public Function fNewArraySng(ByVal strText As String, _
                             Optional ByVal strColSeparator As String = " ", _
                             Optional ByVal strRowSeparator As String = ";", _
                             Optional ByVal strNewLine As String = "\") As Single()
    Dim vArray As Variant
    Dim aSng() As Single
    Dim lgC As Long
    Dim lgR As Long
    
    vArray = fNewArray(strText, strColSeparator, strRowSeparator, strNewLine)
    If IsArray(vArray) Then
        ReDim aSng(LBound(vArray, 1) To UBound(vArray, 1), _
                   LBound(vArray, 2) To UBound(vArray, 2))
        For lgR = LBound(vArray, 1) To UBound(vArray, 1)
            For lgC = LBound(vArray, 2) To UBound(vArray, 2)
                aSng(lgR, lgC) = VBA.CSng(VBA.Val(vArray(lgR, lgC)))
            Next lgC
        Next lgR
        fNewArraySng = aSng()
    
        Erase aSng()
    End If
    Erase vArray
End Function
Public Function fNewArrayLng(ByVal strText As String, _
                             Optional ByVal strColSeparator As String = " ", _
                             Optional ByVal strRowSeparator As String = ";", _
                             Optional ByVal strNewLine As String = "\") As Long()
    Dim vArray As Variant
    Dim aLng() As Long
    Dim lgC As Long
    Dim lgR As Long
    
    vArray = fNewArray(strText, strColSeparator, strRowSeparator, strNewLine)
    If IsArray(vArray) Then
        ReDim aLng(LBound(vArray, 1) To UBound(vArray, 1), _
                   LBound(vArray, 2) To UBound(vArray, 2))
        For lgR = LBound(vArray, 1) To UBound(vArray, 1)
            For lgC = LBound(vArray, 2) To UBound(vArray, 2)
                aLng(lgR, lgC) = VBA.CLng(VBA.Val(vArray(lgR, lgC)))
            Next lgC
        Next lgR
        fNewArrayLng = aLng()
    
        Erase aLng()
    End If
    Erase vArray
End Function
Public Function fNewArrayInt(ByVal strText As String, _
                             Optional ByVal strColSeparator As String = " ", _
                             Optional ByVal strRowSeparator As String = ";", _
                             Optional ByVal strNewLine As String = "\") As Integer()
    Dim vArray As Variant
    Dim aInt() As Integer
    Dim lgC As Long
    Dim lgR As Long
    
    vArray = fNewArray(strText, strColSeparator, strRowSeparator, strNewLine)
    If IsArray(vArray) Then
        ReDim aInt(LBound(vArray, 1) To UBound(vArray, 1), _
                   LBound(vArray, 2) To UBound(vArray, 2))
        For lgR = LBound(vArray, 1) To UBound(vArray, 1)
            For lgC = LBound(vArray, 2) To UBound(vArray, 2)
                aInt(lgR, lgC) = VBA.CInt(VBA.Val(vArray(lgR, lgC)))
            Next lgC
        Next lgR
        fNewArrayInt = aInt()
    
        Erase aInt()
    End If
    Erase vArray
End Function
Public Function fNewArrayBool(ByVal strText As String, _
                              Optional ByVal strColSeparator As String = " ", _
                              Optional ByVal strRowSeparator As String = ";", _
                              Optional ByVal strNewLine As String = "\") As Boolean()
    Dim vArray As Variant
    Dim aBool() As Boolean
    Dim lgC As Long
    Dim lgR As Long
    
    vArray = fNewArray(strText, strColSeparator, strRowSeparator, strNewLine)
    If IsArray(vArray) Then
        ReDim aBool(LBound(vArray, 1) To UBound(vArray, 1), _
                   LBound(vArray, 2) To UBound(vArray, 2))
        For lgR = LBound(vArray, 1) To UBound(vArray, 1)
            For lgC = LBound(vArray, 2) To UBound(vArray, 2)
                aBool(lgR, lgC) = VBA.CBool(VBA.Val(vArray(lgR, lgC)))
            Next lgC
        Next lgR
        fNewArrayBool = aBool()
    
        Erase aBool()
    End If
    Erase vArray
End Function
Public Function fNewArrayByte(ByVal strText As String, _
                              Optional ByVal strColSeparator As String = " ", _
                              Optional ByVal strRowSeparator As String = ";", _
                              Optional ByVal strNewLine As String = "\") As Byte()
    Dim vArray As Variant
    Dim aByte() As Byte
    Dim lgC As Long
    Dim lgR As Long
    
    vArray = fNewArray(strText, strColSeparator, strRowSeparator, strNewLine)
    If IsArray(vArray) Then
        ReDim aByte(LBound(vArray, 1) To UBound(vArray, 1), _
                   LBound(vArray, 2) To UBound(vArray, 2))
        For lgR = LBound(vArray, 1) To UBound(vArray, 1)
            For lgC = LBound(vArray, 2) To UBound(vArray, 2)
                aByte(lgR, lgC) = VBA.CByte(VBA.Val(vArray(lgR, lgC)))
            Next lgC
        Next lgR
        fNewArrayByte = aByte()
    
        Erase aByte()
    End If
    Erase vArray
End Function

Public Function fNewVector(ByVal strText As String, _
                           Optional ByVal strColSeparator As String = " ", _
                           Optional ByVal strRowSeparator As String = ";", _
                           Optional ByVal strNewLine As String = "\") As Variant
' vector = [first : step : last]
' To create a vector v with the first element f, last element l, and the difference between elements is any real number n
    Dim mVector As Variant
    
    Dim aCreator() As String
    Dim aElement() As String
    Dim dbFirst As Double
    Dim dbLast As Double
    Dim dbStep As Double
    Dim lgElement As Long
    Dim lgElements As Long
    Dim lgStep As Long
    Dim strVector As String
    
    strText = VBA.Trim$(strText)
    
    aCreator() = VBA.Split(strText, ":")
    If LBound(aCreator)  UBound(aCreator) Then
        dbFirst = VBA.Val(aCreator()(LBound(aCreator) + 0))
        dbLast = VBA.Val(aCreator()(LBound(aCreator) + 2))
        dbStep = VBA.Val(aCreator()(LBound(aCreator) + 1))
        If dbStep = 0 Then dbStep = 1
        lgElements = VBA.CLng((dbLast - dbFirst) / dbStep)
        ReDim mVector(g_Base + 0 To g_Base + lgElements - 1)
        lgStep = 0
        For lgElement = LBound(mVector) To UBound(mVector)
            mVector(lgElement) = dbFirst + (lgStep * dbStep)
            lgStep = lgStep + 1
        Next lgElement
        If mVector(UBound(mVector))  dbLast Then
            ReDim Preserve mVector(LBound(mVector) To UBound(mVector) + 1)
            mVector(UBound(mVector)) = dbLast
        End If
        fNewVector = mVector
    
    Else
        'Join lines...
        strText = VBA.Replace(strText, strNewLine, "")
        
        If VBA.InStr(1, strText, strRowSeparator) = 0 Then
        'row vector
            'Avoid repeated separators
            Do While VBA.InStr(1, strVector, (strColSeparator & strColSeparator)) > 0
                strVector = VBA.Replace(strVector, strColSeparator & strColSeparator, strColSeparator)
            Loop
            
            aElement() = VBA.Split(strText, strColSeparator)
            
            ReDim mVector(g_Base + LBound(aElement) To g_Base + UBound(aElement))
            For lgElement = LBound(aElement) To UBound(aElement)
                mVector(lgElement) = VBA.Val(aElement(lgElement))
            Next lgElement
            
        Else
        'column vector
            'Avoid repeated separators... not likelly on column vectors
            'Do While VBA.InStr(1, strVector, (strRowSeparator & strRowSeparator)) > 0
            '    strVector = VBA.Replace(strVector, strRowSeparator & strRowSeparator, strRowSeparator)
            'Loop
            
            aElement() = VBA.Split(strText, strRowSeparator)
            
            ReDim mVector(g_Base + LBound(aElement) To g_Base + UBound(aElement), g_Base)
            For lgElement = LBound(aElement) To UBound(aElement)
                mVector(lgElement, g_Base) = VBA.Val(aElement(lgElement))
            Next lgElement
        End If

        fNewVector = mVector
    End If
End Function
Public Function fNewVectorStr(ByVal strText As String, _
                              Optional ByVal strColSeparator As String = " ", _
                              Optional ByVal strRowSeparator As String = ";", _
                              Optional ByVal strNewLine As String = "\") As String()
    Dim vArray As Variant
    Dim aStr() As String
    Dim lgC As Long
    Dim lgR As Long
    
    vArray = fNewVector(strText, strColSeparator, strRowSeparator, strNewLine)
    If IsArray(vArray) Then
        If fNdims(vArray) = 1 Then
            ReDim aStr(LBound(vArray, 1) To UBound(vArray, 1))
            For lgR = LBound(vArray, 1) To UBound(vArray, 1)
                aStr(lgR) = VBA.CStr(VBA.Val(vArray(lgR)))
            Next lgR
        
        ElseIf fNdims(vArray) = 2 Then
            ReDim aStr(LBound(vArray, 1) To UBound(vArray, 1), _
                       LBound(vArray, 2) To UBound(vArray, 2))
            For lgR = LBound(vArray, 1) To UBound(vArray, 1)
                For lgC = LBound(vArray, 2) To UBound(vArray, 2)
                    aStr(lgR, lgC) = VBA.CStr(VBA.Val(vArray(lgR, lgC)))
                Next lgC
            Next lgR
        End If
        fNewVectorStr = aStr()
    
        Erase aStr()
    End If
    Erase vArray
End Function
Public Function fNewVectorDbl(ByVal strText As String, _
                              Optional ByVal strColSeparator As String = " ", _
                              Optional ByVal strRowSeparator As String = ";", _
                              Optional ByVal strNewLine As String = "\") As Double()
    Dim vArray As Variant
    Dim aDbl() As Double
    Dim lgC As Long
    Dim lgR As Long
    
    vArray = fNewVector(strText, strColSeparator, strRowSeparator, strNewLine)
    If IsArray(vArray) Then
        If fNdims(vArray) = 1 Then
            ReDim aDbl(LBound(vArray, 1) To UBound(vArray, 1))
            For lgR = LBound(vArray, 1) To UBound(vArray, 1)
                aDbl(lgR) = VBA.CDbl(VBA.Val(vArray(lgR)))
            Next lgR
        
        ElseIf fNdims(vArray) = 2 Then
            ReDim aDbl(LBound(vArray, 1) To UBound(vArray, 1), _
                       LBound(vArray, 2) To UBound(vArray, 2))
            For lgR = LBound(vArray, 1) To UBound(vArray, 1)
                For lgC = LBound(vArray, 2) To UBound(vArray, 2)
                    aDbl(lgR, lgC) = VBA.CDbl(VBA.Val(vArray(lgR, lgC)))
                Next lgC
            Next lgR
        End If
        fNewVectorDbl = aDbl()
    
        Erase aDbl()
    End If
    Erase vArray
End Function
Public Function fNewVectorSng(ByVal strText As String, _
                              Optional ByVal strColSeparator As String = " ", _
                              Optional ByVal strRowSeparator As String = ";", _
                              Optional ByVal strNewLine As String = "\") As Single()
    Dim vArray As Variant
    Dim aSng() As Single
    Dim lgC As Long
    Dim lgR As Long
    
    vArray = fNewVector(strText, strColSeparator, strRowSeparator, strNewLine)
    If IsArray(vArray) Then
        If fNdims(vArray) = 1 Then
            ReDim aSng(LBound(vArray, 1) To UBound(vArray, 1))
            For lgR = LBound(vArray, 1) To UBound(vArray, 1)
                aSng(lgR) = VBA.CSng(VBA.Val(vArray(lgR)))
            Next lgR
        
        ElseIf fNdims(vArray) = 2 Then
            ReDim aSng(LBound(vArray, 1) To UBound(vArray, 1), _
                       LBound(vArray, 2) To UBound(vArray, 2))
            For lgR = LBound(vArray, 1) To UBound(vArray, 1)
                For lgC = LBound(vArray, 2) To UBound(vArray, 2)
                    aSng(lgR, lgC) = VBA.CSng(VBA.Val(vArray(lgR, lgC)))
                Next lgC
            Next lgR
        End If
        fNewVectorSng = aSng()
    
        Erase aSng()
    End If
    Erase vArray
End Function
Public Function fNewVectorLng(ByVal strText As String, _
                              Optional ByVal strColSeparator As String = " ", _
                              Optional ByVal strRowSeparator As String = ";", _
                              Optional ByVal strNewLine As String = "\") As Long()
    Dim vArray As Variant
    Dim aLng() As Long
    Dim lgC As Long
    Dim lgR As Long
    
    vArray = fNewVector(strText, strColSeparator, strRowSeparator, strNewLine)
    If IsArray(vArray) Then
        If fNdims(vArray) = 1 Then
            ReDim aLng(LBound(vArray, 1) To UBound(vArray, 1))
            For lgR = LBound(vArray, 1) To UBound(vArray, 1)
                aLng(lgR) = VBA.CLng(VBA.Val(vArray(lgR)))
            Next lgR
        
        ElseIf fNdims(vArray) = 2 Then
            ReDim aLng(LBound(vArray, 1) To UBound(vArray, 1), _
                       LBound(vArray, 2) To UBound(vArray, 2))
            For lgR = LBound(vArray, 1) To UBound(vArray, 1)
                For lgC = LBound(vArray, 2) To UBound(vArray, 2)
                    aLng(lgR, lgC) = VBA.CLng(VBA.Val(vArray(lgR, lgC)))
                Next lgC
            Next lgR
        End If
        fNewVectorLng = aLng()
    
        Erase aLng()
    End If
    Erase vArray
End Function
Public Function fNewVectorInt(ByVal strText As String, _
                              Optional ByVal strColSeparator As String = " ", _
                              Optional ByVal strRowSeparator As String = ";", _
                              Optional ByVal strNewLine As String = "\") As Integer()
    Dim vArray As Variant
    Dim aInt() As Integer
    Dim lgC As Long
    Dim lgR As Long
    
    vArray = fNewVector(strText, strColSeparator, strRowSeparator, strNewLine)
    If IsArray(vArray) Then
        If fNdims(vArray) = 1 Then
            ReDim aInt(LBound(vArray, 1) To UBound(vArray, 1))
            For lgR = LBound(vArray, 1) To UBound(vArray, 1)
                aInt(lgR) = VBA.CInt(VBA.Val(vArray(lgR)))
            Next lgR
        
        ElseIf fNdims(vArray) = 2 Then
            ReDim aInt(LBound(vArray, 1) To UBound(vArray, 1), _
                       LBound(vArray, 2) To UBound(vArray, 2))
            For lgR = LBound(vArray, 1) To UBound(vArray, 1)
                For lgC = LBound(vArray, 2) To UBound(vArray, 2)
                    aInt(lgR, lgC) = VBA.CInt(VBA.Val(vArray(lgR, lgC)))
                Next lgC
            Next lgR
        End If
        fNewVectorInt = aInt()
    
        Erase aInt()
    End If
    Erase vArray
End Function
Public Function fNewVectorBool(ByVal strText As String, _
                               Optional ByVal strColSeparator As String = " ", _
                               Optional ByVal strRowSeparator As String = ";", _
                               Optional ByVal strNewLine As String = "\") As Boolean()
    Dim vArray As Variant
    Dim aBool() As Boolean
    Dim lgC As Long
    Dim lgR As Long
    
    vArray = fNewVector(strText, strColSeparator, strRowSeparator, strNewLine)
    If IsArray(vArray) Then
        If fNdims(vArray) = 1 Then
            ReDim aBool(LBound(vArray, 1) To UBound(vArray, 1))
            For lgR = LBound(vArray, 1) To UBound(vArray, 1)
                aBool(lgR) = VBA.CBool(VBA.Val(vArray(lgR)))
            Next lgR
        
        ElseIf fNdims(vArray) = 2 Then
            ReDim aBool(LBound(vArray, 1) To UBound(vArray, 1), _
                       LBound(vArray, 2) To UBound(vArray, 2))
            For lgR = LBound(vArray, 1) To UBound(vArray, 1)
                For lgC = LBound(vArray, 2) To UBound(vArray, 2)
                    aBool(lgR, lgC) = VBA.CBool(VBA.Val(vArray(lgR, lgC)))
                Next lgC
            Next lgR
        End If
        fNewVectorBool = aBool()
    
        Erase aBool()
    End If
    Erase vArray
End Function
Public Function fNewVectorByte(ByVal strText As String, _
                               Optional ByVal strColSeparator As String = " ", _
                               Optional ByVal strRowSeparator As String = ";", _
                               Optional ByVal strNewLine As String = "\") As Byte()
    Dim vArray As Variant
    Dim aByte() As Byte
    Dim lgC As Long
    Dim lgR As Long
    
    vArray = fNewVector(strText, strColSeparator, strRowSeparator, strNewLine)
    If IsArray(vArray) Then
        If fNdims(vArray) = 1 Then
            ReDim aByte(LBound(vArray, 1) To UBound(vArray, 1))
            For lgR = LBound(vArray, 1) To UBound(vArray, 1)
                aByte(lgR) = VBA.CByte(VBA.Val(vArray(lgR)))
            Next lgR
        
        ElseIf fNdims(vArray) = 2 Then
            ReDim aByte(LBound(vArray, 1) To UBound(vArray, 1), _
                       LBound(vArray, 2) To UBound(vArray, 2))
            For lgR = LBound(vArray, 1) To UBound(vArray, 1)
                For lgC = LBound(vArray, 2) To UBound(vArray, 2)
                    aByte(lgR, lgC) = VBA.CByte(VBA.Val(vArray(lgR, lgC)))
                Next lgC
            Next lgR
        End If
        fNewVectorByte = aByte()
    
        Erase aByte()
    End If
    Erase vArray
End Function

Public Function fLength(ByRef mArray As Variant) As Long
' length     Length of vector or largest array dimension
    Dim nDim As Long
    Dim lgDim As Long
    
    If IsArray(mArray) Then
        On Error GoTo ExitProc
        lgDim = 0
        Do
            lgDim = lgDim + 1
            If nDim  0)
        Loop
    End If

ExitProc:
    On Error GoTo 0
    fNdims = (lgDim - 1)
End Function

Public Function fNumEl(ByRef mArray As Variant) As Long
' numel      Number of array elements
    Dim lgDim As Long
    Dim lgElements As Long
    
    If IsArray(mArray) Then
        On Error GoTo ExitProc
        lgDim = 0
        Do
            lgDim = lgDim + 1
            lgElements = lgElements * (UBound(mArray, lgDim) - LBound(mArray, lgDim) + 1)
        Loop
    End If

ExitProc:
    On Error GoTo 0
    fNumEl = lgElements
End Function

Public Function IsColumn(ByRef mArray As Variant) As Boolean
' iscolumn   Determines whether input is column vector
    If IsArray(mArray) Then
        IsColumn = (UBound(mArray, 2) - LBound(mArray, 2) = 1)
    End If
End Function

'Public Function IsEmpty(ByRef mArray As Variant) As Boolean
'' isempty    Determines whether array is empty
'    If IsArray(mArray) Then
'        'IsEmpty = True
'    End If
'End Function
'Public Function IsMatrix(ByRef mArray As Variant) As Boolean
'' ismatrix   Determines whether input is matrix
'    If IsArray(mArray) Then
'        'IsMatrix = True
'    End If
'End Function

Public Function IsRow(ByRef mArray As Variant) As Boolean
' isrow      Determines whether input is row vector
    If IsArray(mArray) Then
        IsRow = (UBound(mArray, 1) - LBound(mArray, 1) = 1)
    End If
End Function

Public Function IsScalar(ByRef mArray As Variant) As Boolean
' isscalar   Determines whether input is scalar
    If IsArray(mArray) Then
        IsScalar = (Not IsArray(mArray))
    End If
End Function

Public Function IsVector(ByRef mArray As Variant) As Boolean
' isvector   Determines whether input is vector
    If IsArray(mArray) Then
        IsVector = (UBound(mArray, 1) = LBound(mArray, 1)) Or (UBound(mArray, 2) = LBound(mArray, 2))
    End If
End Function

Public Function fBlkDiag(ByVal mDiagonal As Variant, _
                         Optional ByVal lgDiagonal As Long = 0) As Variant
' blkdiag    Constructs block diagonal matrix from input arguments
'            placing the elements of vector mDiagonal on the lgDiagonal_th diagonal.
'     lgDiagonal=0 represents the main diagonal
'     lgDiagonal>0 is above the main diagonal
'     lgDiagonal<0 is below the main diagonal
    Dim mArray As Variant
    Dim lgElement As Long
    
    If IsArray(mDiagonal) Then
        ReDim mArray(LBound(mDiagonal) To UBound(mDiagonal), LBound(mDiagonal) To UBound(mDiagonal))
        If lgDiagonal = 0 Then
            For lgElement = LBound(mDiagonal) To UBound(mDiagonal)
                mArray(lgElement, lgElement) = mDiagonal(lgElement)
            Next lgElement
        ElseIf lgDiagonal  0 Then
            For lgElement = LBound(mDiagonal) To UBound(mDiagonal)
                mArray(lgElement, lgElement + lgDiagonal) = mDiagonal(lgElement)
            Next lgElement
        End If
        fBlkDiag = mArray
    End If
End Function

Public Function fCircShift(ByRef mArray As Variant, _
                           ByVal mShifter As Variant, _
                           Optional ByVal dimCirculate As Long = 0) As Boolean
' circshift  Shifts array circularly
' Y = circshift(A,K) circularly shifts the elements in array A by K positions.
' If K is an integer, then circshift shifts along the first dimension of A whose size does not equal 1.
' If K is a vector of integers, then each element of K indicates the shift amount in the corresponding dimension of A.
    Dim lgR As Long
    Dim lgC As Long
    Dim lgShift As Long
    Dim mArrayTmp As Variant
    
    If IsArray(mArray) Then
        If dimCirculate = 0 Then
            'Copy array
            mArrayTmp = mArray
            
            If IsArray(mShifter) Then
                For lgR = LBound(mArray, 1) To UBound(mArray, 1)
                    If lgR - mShifter(g_Base + 0)  dbThreshold, mArray(lgR, lgC), dbThreshold)
            Next lgC
        Next lgR
        
        fThreshold = mThreshold
    End If
End Function

Public Function fRound(ByVal mArray As Variant, _
                       Optional ByVal lgDigits As Long = 0) As Variant
    Dim lgR As Long
    Dim lgC As Long
    Dim mRound As Variant

    If IsArray(mArray) Then
        ReDim mThreshold(LBound(mArray, 1) To UBound(mArray, 1), LBound(mArray, 2) To UBound(mArray, 2))
        For lgR = LBound(mArray, 1) To UBound(mArray, 1)
            For lgC = LBound(mArray, 2) To UBound(mArray, 2)
                mRound(lgR, lgC) = VBA.Round(mArray(lgR, lgC), lgDigits)
            Next lgC
        Next lgR
        fRound = mRound
    End If
End Function

Public Function fMagnitude(ByVal mArray As Variant, _
                           Optional ByVal lgOrder As Long = 2) As Double
    Dim lgR As Long
    Dim lgC As Long
    Dim dbMagnitude As Long

    If IsArray(mArray) Then
        For lgR = LBound(mArray, 1) To UBound(mArray, 1)
            For lgC = LBound(mArray, 2) To UBound(mArray, 2)
                dbMagnitude = (mArray(lgR, lgC) * mArray(lgR, lgC))
            Next lgC
        Next lgR
        If lgOrder = 2 Then
            fMagnitude = VBA.Sqr(dbMagnitude)
        Else
            fMagnitude = dbMagnitude ^ (1 / lgOrder)
        End If
    End If
End Function
Also, the Gauss-Jordan reduction method is coded as:
Option Explicit

Dim mArray() As Double

Private Sub UserForm_Initialize()
    'call sSolve
End Sub

Private Sub cbSolve_Click()
    Call sSolve
End Sub

Private Sub sSolve()
    Call sArray_Load
    Call sGaussJordan(mArray())
End Sub

Private Sub sArray_Load()
    On Error GoTo ErrLec
    
    Dim lgR As Long
    Dim lgC As Long
    Dim lgRetVal As Long
    Dim Nm As Long
    
    If VarType(Selection) = vbObject Then
        mArray() = fNewArrayDbl(Me.txtSystem.Text, " ", ";", "\")
        Me.sbNum.Value = UBound(mArray, 1) - LBound(mArray, 1) + 1
        Nm = Me.sbNum.Value
    Else
        If Selection.Rows.Count > 1 And _
           Selection.Rows.Count > 1 Then
            Me.sbNum.Value = Selection.Rows.Count
            Nm = Me.sbNum.Value
            
            ReDim mArray(g_Base To (Nm - 1 + g_Base), g_Base To Nm + g_Base)
            For lgR = g_Base To Nm - 1 + g_Base
                For lgC = g_Base To Nm + g_Base
                    'If Cuadric.TextMatrix(lgR + 1, lgC) = "" Then Cuadric.TextMatrix(lgR + 1, lgC) = 0
                    'mArray(lgR, lgC) = Cuadric.TextMatrix(lgR + 1, lgC)
                    mArray(lgR, lgC) = Selection.Cells(lgR + 1, lgC + 1).Value2
                Next
            Next
        Else
            mArray() = fNewArrayDbl(Me.txtSystem.Text, " ", ";", "\")
            Me.sbNum.Value = UBound(mArray, 1) - LBound(mArray, 1) + 1
            Nm = Me.sbNum.Value
        End If
    End If
    
ExitProc:
    Exit Sub

ErrLec:
    lgRetVal = VBA.MsgBox("Error introducing data (" & Err.Description & ")", vbExclamation)
End Sub

Private Sub sbNum_Change()
    txtEquations.Value = VBA.Str(sbNum.Value)
End Sub

Public Sub sGaussJordan(ByRef mArray() As Double)
'Based on code found: http://mvb6.blogspot.com/2017/08/metodo-de-gauss-jordan-vb-60.html
    Dim lgR As Long
    Dim lgC As Long
    Dim lgPivot As Long
    Dim lgR_Homogenize As Long
    Dim dbTmp As Double
    Dim lgRetVal As Long
    Dim mArrayTmp() As Double
    Dim Nm As Integer
    
    On Error GoTo ErrControl
    
    Nm = UBound(mArray, 1) - LBound(mArray, 1) + 1
    ReDim mArrayTmp(LBound(mArray, 1) To UBound(mArray, 1), LBound(mArray, 2) To UBound(mArray, 2))
    
    ' Swap rows (if needed)
    If (mArray(0, 0) = 0) Then
        For lgR = LBound(mArray, 1) To UBound(mArray, 1)
            If (mArray(lgR, 0)  0) Then
                For lgC = LBound(mArray, 2) To UBound(mArray, 2)
                    mArrayTmp(0, lgC) = mArray(0, lgC)
                    mArray(0, lgC) = mArray(lgR, lgC)
                    mArray(lgR, lgC) = mArrayTmp(0, lgC)
                Next lgC
            End If
        Next lgR
    End If
    
    For lgPivot = LBound(mArray, 1) To UBound(mArray, 1)
        dbTmp = mArray(lgPivot, lgPivot)
        For lgC = LBound(mArray, 2) To UBound(mArray, 2)
            mArray(lgPivot, lgC) = mArray(lgPivot, lgC) / dbTmp
        Next lgC
        For lgR = LBound(mArray, 1) To UBound(mArray, 1)
            If (lgR = lgPivot) Then GoTo Es
            dbTmp = mArray(lgR, lgPivot)
            For lgR_Homogenize = LBound(mArray, 2) To UBound(mArray, 2)
                mArray(lgR, lgR_Homogenize) = mArray(lgR, lgR_Homogenize) - (dbTmp * mArray(lgPivot, lgR_Homogenize))
            Next lgR_Homogenize
Es:
        Next lgR
    Next lgPivot
    
    'Print solution
    For lgR = LBound(mArray, 1) To UBound(mArray, 1)
        Debug.Print mArray(lgR, Nm) 'vba.Format(mArray(lgR, Nm), "##,##0.00")
    Next lgR
    
ExitProc:
    Exit Sub

ErrControl:
    lgRetVal = VBA.MsgBox("System has no solution", vbCritical)
End Sub
Future implementations will be, for example, then capability to handle complex numbers, via the UDT tObject, with little to no change in the code (only replacing “) As Variant” with “) As tObject()” and ” As Variant” with “() As tObject”. Even NaN, NaT, Inf,… MatLab special reserved variables can be used along the code, setting TypeObj to 0 and giving .Text property the name of the reserved variable. This will be a possible group implementation to handle complex numbers inside VBA
Option Explicit

Public Enum eTypeObj
    eText = 0
    eNatural = 2
    eReal = 1
    eComplex = -1
End Enum
Public Type tObject
    Size As Long 'Total bytes
    
    TypeObj As Long
        'Text = 0
        '[R]Real = 1, [C]Complex = -1
        '[Z]Natural (Integers ±) = 2
    
    R As Double 'Real part
    I As Double 'Imaginary part
    
    Text As String
    'Name As String * 10
End Type

Public Function fNew(Optional ByVal TypeObj As Long = 0, _
                     Optional ByVal R As Double = 0, _
                     Optional ByVal c As Double = 0, _
                     Optional ByVal Text As String = "") As tObject
'Set new object
    With fNew
        .Size = 20 + Len(Text)
        .TypeObj = TypeObj '[Z]Natural = 2, [R]Real = 1, [C]Complex = -1, Text = 0
        
        .R = R
        .I = c
        
        .Text = Text
    End With
End Function
Public Function fComplex(ByRef dbReal As Double, _
                         ByRef dbImaginary As Double) As tObject
' complex   Create complex array
    fComplex = fNew(eComplex, dbReal, dbImaginary, "")
End Function
Public Function fAbs(ByRef oObject As tObject) As Double
' abs       Absolute value and complex magnitude
    With oObject
        If .TypeObj > 0 Then
            fAbs = VBA.Abs(.R)
        ElseIf .TypeObj = eComplex Then
            fAbs = VBA.Sqr(.R ^ 2 + .I ^ 2)
        End If
    End With
End Function
Public Function fAngle(ByRef oObject As tObject) As Double
' angle     Phase angle
    With oObject
        If VBA.Abs(.R)  0 Then
            fSignObj = fSign(oObject.R)
        ElseIf .TypeObj = eComplex Then
            fSignObj = fSign(oObject.R)
        End If
    End With
End Function
Public Function fUnwrap(ByRef oObject1 As tObject, _
                        ByRef oObject2 As tObject) As Double
' unwrap    Correct phase angles to produce smoother phase plots
'!!!!!
End Function
Public Function fReal(ByRef oObject As tObject) As Double
' real      Real part of complex number
    If oObject.TypeObj = eComplex Then fReal = oObject.R
End Function
Public Function fImag(ByRef oObject As tObject) As Double
' imag      Imaginary part of complex number
    If oObject.TypeObj = eComplex Then fImag = oObject.I
End Function

Public Function fComplexSum(ByRef oObject1 As tObject, _
                            ByRef oObject2 As tObject) As tObject
    If (VBA.Abs(oObject1.TypeObj) = eReal Or VBA.Abs(oObject2.TypeObj) = eReal) Then
        With fComplexSum
            .R = oObject1.R + oObject2.R
            If (oObject1.TypeObj = eComplex Or oObject2.TypeObj = eComplex) Then
                .TypeObj = eComplex
                .I = oObject1.I + oObject2.I
            Else
                .TypeObj = eReal
            End If
        End With
    End If
End Function

Public Function fComplexDiff(ByRef oObject1 As tObject, _
                             ByRef oObject2 As tObject) As tObject
    If (VBA.Abs(oObject1.TypeObj) = eReal Or VBA.Abs(oObject2.TypeObj) = eReal) Then
        With fComplexDiff
            .R = oObject1.R - oObject2.R
            If (oObject1.TypeObj = eComplex Or oObject2.TypeObj = eComplex) Then
                .TypeObj = eComplex
                .I = oObject1.I - oObject2.I
            Else
                .TypeObj = eReal
            End If
        End With
    End If
End Function

Public Function fComplexMult(ByRef oObject1 As tObject, _
                             ByRef oObject2 As tObject) As tObject
' z1·z2 = (a, b)·(c, d) = (a·c - b·d), (a·d - b·c)
    If (VBA.Abs(oObject1.TypeObj) = eReal Or VBA.Abs(oObject2.TypeObj) = eReal) Then
        With fComplexMult
            .R = oObject1.R * oObject2.R 'only the real part
            If (oObject1.TypeObj = eComplex Or oObject2.TypeObj = eComplex) Then
                .TypeObj = eComplex
                .I = oObject1.R * oObject2.I + oObject1.I * oObject2.R
                .R = .R - oObject1.I * oObject2.I
            Else
                .TypeObj = eReal
            End If
        End With
    End If
End Function

Public Function fComplexRec(ByRef oObject As tObject) As tObject
' 1/z = 1/(a, b) = (a, -b)/(a²+b²)
    Dim oReciproc As tObject
    Dim dbModule² As Double
    
    With oObject
        If .TypeObj = eComplex Then
            dbModule² = (.R ^ 2 + .I ^ 2)
            With fComplexRec
                .TypeObj = oObject.TypeObj
                .R = oObject.R / dbModule²
                .I = -oObject.I / dbModule²
            End With
        
        ElseIf .TypeObj = eReal Then
            With fComplexRec
                .TypeObj = oObject.TypeObj
                .R = oObject.R
            End With
        End If
    End With
End Function

Private Sub sComplexDiv()
    Dim oObject1 As tObject
    Dim oObject2 As tObject
    Dim oObject As tObject
    With oObject1
        .TypeObj = eComplex
        .R = 4
        .I = 3
    End With
    With oObject2
        .TypeObj = eComplex
        .R = 2
        .I = 1
    End With
    oObject = fComplexDiv(oObject1, oObject2)
Stop
End Sub

Public Function fComplexDiv(ByRef oObject1 As tObject, _
                            ByRef oObject2 As tObject) As tObject
' z1/z2 = (a, b)·[(c, d)/(c²+d²)] = (ac+bd , cb-da)/(c²+d²)
    'Dim oComplexRec As tObject
    Dim dbModule² As Double
    
    If (VBA.Abs(oObject1.TypeObj) = eReal Or VBA.Abs(oObject2.TypeObj) = eReal) Then
        With fComplexDiv
            If (oObject2.TypeObj = eComplex) Then
                dbModule² = (oObject2.R ^ 2 + oObject2.I ^ 2)
                .TypeObj = eComplex
                'oComplexRec = fComplexRec(oObject2)
                .R = (oObject1.R * oObject2.R + oObject1.I * oObject2.I) / dbModule²
                .I = (oObject1.I * oObject2.R - oObject1.R * oObject2.I) / dbModule²
            
            Else
                .TypeObj = eReal
                .R = oObject1.R / oObject2.R
                .I = oObject1.R / oObject2.R
            End If
        End With
    End If
End Function

Public Function fConj(ByRef oObject As tObject) As tObject
' conj      Complex conjugate
    With fConj
        If .TypeObj = eComplex Then
            .TypeObj = oObject.TypeObj
            .R = oObject.R
            .I = -oObject.I
        
        ElseIf .TypeObj = eReal Then
            .TypeObj = oObject.TypeObj
            .R = oObject.R
        End If
    End With
End Function