'__Author: Dohyeon Bong @kaeru8714 2020'

Public m As Integer 'Number of Materials'
Public n As Integer 'Number of Points'
Public s As Integer 'Number of Resistance Force'
Public figOrigin(2) As Double
Function getLengthByCoord(x As Double, y As Double, A As Double, B As Double) As Double
    getLengthByCoord = Sqr((x - A) * (x - A) + (y - B) * (y - B)) 'calculate Euclidean distance'
End Function
Function getLengthByPoint(p As Long, q As Long) As Double
    Dim px As Double
    Dim qx As Double
    Dim py As Double
    Dim qy As Double
    px = Range("B10").Offset(p - 1, 0).Value
    qx = Range("B10").Offset(q - 1, 0).Value
    py = Range("C10").Offset(p - 1, 0).Value
    qy = Range("C10").Offset(q - 1, 0).Value
    getLengthByPoint = Sqr((px - qx) * (px - qx) + (py - qy) * (py - qy))
End Function
Function getOtherPointInMaterial(materialNo As Long, p As Long) As Long
    Dim matRange As Range
    Set matRange = Range("H10").Offset(materialNo - 1, 0)

    If matRange.Offset(0, 1) = p Then
        getOtherPointInMaterial = matRange.Offset(0, 2).Value
    End If
    If matRange.Offset(0, 2) = p Then
        getOtherPointInMaterial = matRange.Offset(0, 1).Value
    End If
End Function
Function getXdifference(p As Long, q As Long) As Double
    Dim px As Double
    Dim qx As Double

    px = Range("B10").Offset(p - 1, 0).Value
    qx = Range("B10").Offset(q - 1, 0).Value

    getXdifference = px - qx

End Function
Function getYdifference(p As Long, q As Long) As Double
    Dim py As Double
    Dim qy As Double

    py = Range("C10").Offset(p - 1, 0).Value
    qy = Range("C10").Offset(q - 1, 0).Value

    getYdifference = py - qy

End Function
Function getPointCoord(p As Integer)
    Dim Coord(2) As Double

    Coord(1) = Range("A" & (p + 9), "C" & (p + 9)).Cells(1, 2)
    Coord(2) = Range("A" & (p + 9), "C" & (p + 9)).Cells(1, 3)


    getPointCoord = Coord
End Function
Function getXEquilibrium(p As Long)

    Dim materials() As Long
    Dim materialRange As Range
    Dim matrixOrigin As Range
    Dim i As Long
    Dim j As Long
    Dim q As Long
    ReDim materials(m)

    Set materialRange = Range("I10", Range("J10").Offset(m - 1, 0))
    Set matrixOrigin = Range("O2")

    i = 0
    For Each k In materialRange
        If k.Value = p Then
            i = i + 1

            If k.Column = 9 Then 'Iｿｭ'
                materials(i) = k.Offset(0, -1).Value
            End If
            If k.Column = 10 Then 'Jｿｭ'
                materials(i) = k.Offset(0, -2).Value
            End If
        End If
    Next k

    For j = 1 To i
        q = getOtherPointInMaterial(materials(j), p)
        matrixOrigin.Offset(p - 1, materials(j) - 1).Value = getXdifference(p, q) / getLengthByPoint(p, q)
    Next

End Function
Function getYEquilibrium(p As Long)
    Dim materials() As Long
    Dim materialRange As Range
    Dim matrixOrigin As Range
    Dim i As Long
    Dim j As Long
    Dim q As Long
    ReDim materials(m)

    Set materialRange = Range("I10", Range("J10").Offset(m - 1, 0))
    Set matrixOrigin = Range("O2")

    i = 0
    For Each k In materialRange
        If k.Value = p Then
            i = i + 1

            If k.Column = 9 Then 'Iｿｭ'
                materials(i) = k.Offset(0, -1).Value

            End If
            If k.Column = 10 Then 'Jｿｭ'
                materials(i) = k.Offset(0, -2).Value
            End If
        End If
    Next k

    For j = 1 To i
        q = getOtherPointInMaterial(materials(j), p)
        matrixOrigin.Offset(p - 1 + n, materials(j) - 1).Value = getYdifference(p, q) / getLengthByPoint(p, q)
    Next
End Function

Sub Clear()
    Range("N1", Range("N1").End(xlToRight).End(xlToRight).End(xlDown)).ClearContents
    Range("J2", "L8").ClearContents
    Range("K10", Range("K10").End(xlDown).End(xlDown)).ClearContents
    Range("L10", Range("K10").End(xlDown).End(xlDown)).ClearContents

End Sub
Function init()
    Call Clear
    Dim matrixRange As Range
    Dim i As Integer

    Set matrixRange = Range(Range("O2"), Range("O2").Offset(2 * n - 1, 2 * n))


    For i = 1 To n
        Cells(2, 14).Offset(i - 1, 0).Value = "X" & i
        Cells(2, 14).Offset(i - 1 + n, 0).Value = "Y" & i
    Next

    For i = 1 To m
        Cells(1, 14).Offset(0, i).Value = "F" & i
    Next


    For Each k In matrixRange
        k.Value = 0
    Next



End Function
Function Inverse(m As Range) As Variant
        Dim minv As Variant
        minv = Application.MInverse(m)
        Inverse = minv
    End Function
Function Multify(m As Range, n As Range) As Variant
    Dim mulv As Variant
    mulv = Application.MMult(m, n)
    Multify = mulv
End Function

Function calculateInverseMatrix()
    Dim matrixRange As Range
    Dim inverseMatrixRange As Range
    Dim i As Integer

    Set matrixRange = Range(Range("O2"), Range("O2").Offset(2 * n - 1, 2 * n - 1))
    Set inverseMatrixRange = Range(Range("O" & (2 * n + 4)), Range("O" & (2 * n + 4)).Offset(2 * n - 1, 2 * n - 1))


    inverseMatrixRange.Value = Inverse(matrixRange)

End Function
Function calculateMaterialForce()
    On Error GoTo ErrorHandler
    Dim inverseMatrixRange As Range
    Dim constMatrixRange As Range
    Dim answerMatrixRange As Range
    Dim k As Long

    Set inverseMatrixRange = Range(Range("O" & (2 * n + 4)), Range("O" & (2 * n + 4)).Offset(2 * n - 1, 2 * n - 1))
    Set constMatrixRange = Range(Range("O2").Offset(0, 2 * n), Range("O2").Offset(2 * n - 1, 2 * n))
    Set answerMatrixRange = Range(Range("L10"), Range("L10").Offset(2 * n - 1, 0))
    Set resistanceForceRange = Range(Range("K2"), Range("K2").Offset(s - 1, 0))

    If IsError(inverseMatrixRange.Range("A1").Value) Then
        MsgBox "Determinant of matrix is zero.", , "Error!"
    End If

    answerMatrixRange.Value = Multify(inverseMatrixRange, constMatrixRange)

    For Each c In Range(Range("L10"), Range("L10").Offset(2 * n - 1, 0))
        c.Value = Round(c.Value, 5)
    Next

    Range(Range("K2"), Range("K2").Offset(s - 1, 0)).Value = answerMatrixRange.Range("A" & m + 1, "A" & m + s).Value
    answerMatrixRange.Range("A" & m + 1, "A" & m + s).ClearContents
ErrorHandler:



End Function
Function calculateMaterialLength()
    Dim i As Long
    Dim x As Double
    Dim y As Double
    Dim A As Double
    Dim B As Double
    For i = 1 To m
        x = getPointCoord(Range("H" & (i + 9), "J" & (i + 9)).Cells(1, 2).Value)(1)
        y = getPointCoord(Range("H" & (i + 9), "J" & (i + 9)).Cells(1, 2).Value)(2)
        A = getPointCoord(Range("H" & (i + 9), "J" & (i + 9)).Cells(1, 3).Value)(1)
        B = getPointCoord(Range("H" & (i + 9), "J" & (i + 9)).Cells(1, 3).Value)(2)
        Range("K" & (i + 9)).Value = getLengthByCoord(x, y, A, B)
    Next i
End Function
Function calculateEquilibriumCondition()
    Dim d As Long

    d = m + s - 2 * n


    If Not d = 0 Then
        MsgBox "d = " & Str(m) & " +" & Str(s) & " - 2*" & Str(n) & " = " & Str(d), , "ERROR!"
        End
    End If

    MsgBox "d = " & Str(m) & " +" & Str(s) & " - 2*" & Str(n) & " = " & Str(d), , "Equilibrium Condition OK"
    MsgBox "2020 Dohyeon Bong @kaeru8714", , "Author"
End Function
Function calculateForceEquilibrium()
    Dim p As Long
    For p = 1 To n
        getXEquilibrium (p)
        getYEquilibrium (p)
    Next p
End Function


Sub main()
    m = Cells(2, 2).Value
    n = Cells(1, 2).Value
    s = WorksheetFunction.CountIf(Range("D10", "E" & (n + 9)), True)

    Call init
    Call calculateEquilibriumCondition
    Call calculateMaterialLength

    Dim i As Long
    Dim j As Long
    Dim k As Long
    Dim p As Long

    Call calculateForceEquilibrium

    Range("N1").Value = "Matrix"
    Range("N" & 2 * n + 3).Value = "InvMatrix"

    Dim xPoints() As Long
    Dim yPoints() As Long
    ReDim xPoints(n)
    ReDim yPoints(n)
    Dim fixRange As Range
    Set fixRange = Range("D10", Range("D10").Offset(n, 1))

    i = 0
    j = 0
    For Each c In fixRange
        If c.Value = True Then
            If c.Column = 4 Then
                i = i + 1
                xPoints(i) = c.Offset(0, -3).Value
            End If
            If c.Column = 5 Then
                j = j + 1
                yPoints(j) = c.Offset(0, -4).Value
            End If
        End If
    Next c


    For k = 1 To i
        Range("N1").Offset(0, m + k).Value = "R" & xPoints(k) & "x"
        Range("J2").Offset(k - 1, 0).Value = "R" & xPoints(k) & "x"
        Range("N1").Offset(xPoints(k), m + k).Value = -1
    Next k

    For k = 1 To j
        Range("N1").Offset(0, m + k + i).Value = "R" & yPoints(k) & "y"
        Range("J2").Offset(k + i - 1, 0).Value = "R" & yPoints(k) & "y"
        Range("N" & n + 1).Offset(yPoints(k), m + k + i).Value = -1
    Next k


    Dim loadRange As Range
    Set loadRange = Range("F10", Range("F10").Offset(n, 1))
    Dim lastOrder As Long
    lastOrder = m + k + i
    i = 0
    j = 0
    Dim xLoads() As Long
    Dim yLoads() As Long
    ReDim xLoads(n)
    ReDim yLoads(n)
    For Each c In loadRange
        If Not c.Value = 0 Then
            If c.Column = 6 Then
                i = i + 1
                xLoads(i) = c.Offset(0, -5).Value
            End If
            If c.Column = 7 Then
                j = j + 1
                yLoads(j) = c.Offset(0, -6).Value
            End If
        End If
    Next c

    Range("N1").Offset(0, lastOrder).Value = "const"
    For k = 1 To i
        Range("N1").Offset(xLoads(k), lastOrder).Value = Range("F10").Offset(xLoads(k) - 1, 0).Value
    Next k

    lastOrder = lastOrder + i
    For k = 1 To j
        Range("N" & n + 1).Offset(yLoads(k), lastOrder).Value = Range("F10").Offset(yLoads(k) - 1, 1).Value
    Next k

    Call calculateInverseMatrix
    Call calculateMaterialForce


End Sub

Sub getMN()
    Sheet1.Range("B1").Value = Range("A10").End(xlDown).Value
    Sheet1.Range("B2").Value = Range("H10").End(xlDown).Value
End Sub