Excel - Cardinal Spline, Cubic Spline, and Monotone Cubic Spline array functions

MS Excel VBA code to be put in a module - array functions

In MS Excel, press Alt+F11 to open the Visual Basic Editor (Windows). Via top menu: Insert → Module. Copy the MS Excel VBA code of the functions into the module.

Generic helper function to read range data into single dimensional array.
Your input data can now be column organized (default table) or row organized (transposed default table).

Cardinal spline

Public Function CardinalSpline(DatapointsX, DatapointsY, ByVal numSegments As Integer, ByVal Tension As Single)

Dim xin As Variant, yin As Variant
    Dim ptsX(), ptsY()
    Dim iRow, i, op, r As Long
    Dim n As Integer
    Dim t, thi, tlo, x, xnow, y As Double
    Dim Outputs() As Variant
    
    'trick: transpose to get rid of second dimension of range
    xin = RangeToArray(DatapointsX)
    yin = RangeToArray(DatapointsY)
    n = UBound(xin)
	
    'sort in ascending X order
    Dim tmp
    If xin(1) > xin(n) Then
        'swap values
        For i = 1 To WorksheetFunction.RoundDown(n / 2, 0)
            tmp = xin(i)
            xin(i) = xin(n + 1 - i)
            xin(n + 1 - i) = tmp
            tmp = yin(i)
            yin(i) = yin(n + 1 - i)
            yin(n + 1 - i) = tmp
        Next
    End If
	
    ReDim ptsX(n + 2)
    ReDim ptsY(n + 2)
    For i = 1 To n
        ptsX(i) = xin(i)
        ptsY(i) = yin(i)
    Next
    ptsX(0) = ptsX(1)
    ptsX(n + 1) = ptsX(n)
    ptsX(n + 2) = ptsX(n)
    ptsY(0) = ptsY(1)
    ptsY(n + 1) = ptsY(n)
    ptsY(n + 2) = ptsY(n)

'spline calculation
    iRow = n + 1
    op = 0
    ReDim Outputs((iRow - 2) * numSegments, 1)
    For i = 1 To iRow - 2
        For n = 0 To numSegments - 1
            xnow = ptsX(i) + n / numSegments * (ptsX(i + 1) - ptsX(i))
            Outputs(op, 0) = xnow
            
            thi = 1
            tlo = 0
            r = 0
            Do
                r = r + 1
                t = (thi + tlo) / 2
                x = Tension * (-t ^ 3 + 2 * t ^ 2 - t) * ptsX(i - 1) + _
                    Tension * (-t ^ 3 + t ^ 2) * ptsX(i) + _
                              (2 * t ^ 3 - 3 * t ^ 2 + 1) * ptsX(i) + _
                    Tension * (t ^ 3 - 2 * t ^ 2 + t) * ptsX(i + 1) + _
                             (-2 * t ^ 3 + 3 * t ^ 2) * ptsX(i + 1) + _
                    Tension * (t ^ 3 - t ^ 2) * ptsX(i + 2)
                    
                If xnow > x Then tlo = t Else thi = t
            Loop Until (Abs(x - xnow) < 0.00001) Or (r > 10000)
			
            y = Tension * (-t ^ 3 + 2 * t ^ 2 - t) * ptsY(i - 1) + _
                Tension * (-t ^ 3 + t ^ 2) * ptsY(i) + _
                          (2 * t ^ 3 - 3 * t ^ 2 + 1) * ptsY(i) + _
                Tension * (t ^ 3 - 2 * t ^ 2 + t) * ptsY(i + 1) + _
                         (-2 * t ^ 3 + 3 * t ^ 2) * ptsY(i + 1) + _
                Tension * (t ^ 3 - t ^ 2) * ptsY(i + 2)
            
            If (n = 0) Then y = ptsY(i)
            
            Outputs(op, 1) = y
            op = op + 1
        Next n
    Next i
    Outputs(op, 0) = ptsX(iRow)
    Outputs(op, 1) = ptsY(iRow)
    CardinalSpline = Outputs
    
End Function

Cubic spline

Public Function CubicSpline(DatapointsX, DataPointsY, ByVal numSegments As Integer)
    'original code - see further below - shortened and optimized for Excel as array function
    
    Dim n As Long, i As Long, iRow, op, s As Long
    Dim xin As Variant, yin As Variant
    Dim p, sig, h, b, a As Double
    
    'trick: transpose to get rid of second dimension of range
    xin = RangeToArray(DatapointsX)
    yin = RangeToArray(DataPointsY)
    n = UBound(xin)
    
    ReDim u(n - 1) As Single
    ReDim yt(n) As Single
    yt(1) = 0
    u(1) = 0
    For i = 2 To n - 1
        sig = (xin(i) - xin(i - 1)) / (xin(i + 1) - xin(i - 1))
        p = sig * yt(i - 1) + 2
        yt(i) = (sig - 1) / p
        u(i) = (yin(i + 1) - yin(i)) / (xin(i + 1) - xin(i)) - (yin(i) - yin(i - 1)) / (xin(i) - xin(i - 1))
        u(i) = (6 * u(i) / (xin(i + 1) - xin(i - 1)) - sig * u(i - 1)) / p
    Next i
    yt(n) = 0
    For i = n - 1 To 1 Step -1
        yt(i) = yt(i) * yt(i + 1) + u(i)
    Next i
    
    ReDim Outputs((n - 1) * numSegments, 1) As Double
    op = 0
    For iRow = 1 To n - 1
        For s = 0 To numSegments - 1
            b = s / numSegments
            a = 1 - b
            h = xin(iRow + 1) - xin(iRow)
            Outputs(op, 0) = xin(iRow) + s / numSegments * h
            Outputs(op, 1) = a * yin(iRow) + b * yin(iRow + 1) + ((a ^ 3 - a) * yt(iRow) + (b ^ 3 - b) * yt(iRow + 1)) * (h ^ 2) / 6
            op = op + 1
        Next s
    Next
    Outputs(op, 0) = xin(n)
    Outputs(op, 1) = yin(n)
    CubicSpline = Outputs
End Function

Monotone cubic spline

Public Function MonotoneCubicSpline(DatapointsX, DataPointsY, ByVal numSegments As Integer)
    
    Dim n As Long, i As Long, iRow, op, s As Long
    Dim xin As Variant, yin As Variant
    Dim m, mNext, dx, dxNext, common, cl, m_, invDx, common_, h, c1, diff
    Dim cls()
    
    'trick: transpose to get rid of second dimension of range
    xin = RangeToArray(DatapointsX)
    yin = RangeToArray(DataPointsY)
    n = UBound(xin)
    
    'Get consecutive differences and slopes
    ReDim dxs(n)
    ReDim dys(n)
    ReDim ms(n)
    
    For i = 1 To n - 1
        dxs(i) = xin(i + 1) - xin(i)
        dys(i) = yin(i + 1) - yin(i)
        ms(i) = dys(i) / dxs(i)
    Next i
    
    'Get degree-1 coefficients
    ReDim c1s(1 To n + 1)
    c1s(1) = ms(1)
    For i = 1 To n - 1
        m = ms(i)
        mNext = ms(i + 1)
        If (m * mNext <= 0) Then
            c1s(i + 1) = 0
        Else
            dx = dxs(i)
            dxNext = dxs(i + 1)
            common = dx + dxNext
            c1s(i + 1) = 3 * common / ((common + dxNext) / m + (common + dx) / mNext)
        End If
    Next
    c1s(n) = ms(n - 1)
    
    'Get degree-2 and degree-3 coefficients
    ReDim c2s(n)
    ReDim c3s(n)
    For i = 1 To n - 1
        c1 = c1s(i)
        m_ = ms(i)
        invDx = 1 / dxs(i)
        common_ = c1 + c1s(i + 1) - m_ - m_
        c2s(i) = (m_ - c1 - common_) * invDx
        c3s(i) = common_ * invDx * invDx
    Next
    
    ReDim Outputs((n - 1) * numSegments, 1) As Double
    op = 0
    For iRow = 1 To n - 1
        For s = 0 To numSegments - 1
            h = xin(iRow + 1) - xin(iRow)
            Outputs(op, 0) = xin(iRow) + s / numSegments * h
            diff = s / numSegments * h
            Outputs(op, 1) = yin(iRow) + c1s(iRow) * diff + c2s(iRow) * diff * diff + c3s(iRow) * diff * diff * diff
            op = op + 1
        Next s
    Next
    Outputs(op, 0) = xin(n)
    Outputs(op, 1) = yin(n)
    MonotoneCubicSpline = Outputs
End Function

Regular function - Single Y value output only, given X value

In MS Excel, press Alt+F11 to open the Visual Basic Editor (Windows). Via top menu: Insert → Module. Copy the MS Excel VBA code of the functions into the module.

Cardinal spline

Public Function CardinalSpline_SP(DatapointsX, DatapointsY, ByVal tension As Single, atX)

Dim xin As Variant, yin As Variant
    Dim ptsX()
    Dim ptsY()
    Dim iRow, i, op, r As Long
    Dim n As Integer
    Dim t, thi, tlo, x, xnow, y As Double
    Dim klo, khi As Long
    Dim X1, X2
    
    'trick: transpose to get rid of second dimension of range
    xin = RangeToArray(DatapointsX)
    yin = RangeToArray(DatapointsY)
    n = UBound(xin)
	
    'sort in ascending X order
    Dim tmp
    If xin(1) > xin(n) Then
        'swap values
        For i = 1 To WorksheetFunction.RoundDown(n / 2, 0)
            tmp = xin(i)
            xin(i) = xin(n + 1 - i)
            xin(n + 1 - i) = tmp
            tmp = yin(i)
            yin(i) = yin(n + 1 - i)
            yin(n + 1 - i) = tmp
        Next
    End If

ReDim ptsX(n + 2)
    ReDim ptsY(n + 2)
    For i = 1 To n
        ptsX(i) = xin(i)
        ptsY(i) = yin(i)
    Next
    ptsX(0) = ptsX(1)
    ptsX(n + 1) = ptsX(n)
    ptsX(n + 2) = ptsX(n)
    ptsY(0) = ptsY(1)
    ptsY(n + 1) = ptsY(n)
    ptsY(n + 2) = ptsY(n)
    
    'spline calculation
    For i = 0 To n + 2
        If ptsX(i) > atX Then Exit For
    Next
    i = i - 1
    
    If atX = ptsX(i) Then
        CardinalSpline_SP = ptsY(i)
    Else
        thi = 1
        tlo = 0
        r = 0
        Do
            r = r + 1
            t = (thi + tlo) / 2
            x = tension * (-t ^ 3 + 2 * t ^ 2 - t) * ptsX(i - 1) + _
                tension * (-t ^ 3 + t ^ 2) * ptsX(i) + _
                          (2 * t ^ 3 - 3 * t ^ 2 + 1) * ptsX(i) + _
                tension * (t ^ 3 - 2 * t ^ 2 + t) * ptsX(i + 1) + _
                         (-2 * t ^ 3 + 3 * t ^ 2) * ptsX(i + 1) + _
                tension * (t ^ 3 - t ^ 2) * ptsX(i + 2)
            If atX > x Then tlo = t Else thi = t
        Loop Until (Abs(x - atX) < 0.00001) Or (r > 10000)
        y = tension * (-t ^ 3 + 2 * t ^ 2 - t) * ptsY(i - 1) + _
            tension * (-t ^ 3 + t ^ 2) * ptsY(i) + _
                      (2 * t ^ 3 - 3 * t ^ 2 + 1) * ptsY(i) + _
            tension * (t ^ 3 - 2 * t ^ 2 + t) * ptsY(i + 1) + _
                     (-2 * t ^ 3 + 3 * t ^ 2) * ptsY(i + 1) + _
            tension * (t ^ 3 - t ^ 2) * ptsY(i + 2)
        CardinalSpline_SP = y
    End If
End Function

Cubic spline

Public Function CubicSpline_SP(DatapointsX, DatapointsY, ByVal atX As Single) As Double
    Dim n, i As Long
    Dim xin, yin As Variant
    Dim p, sig, h, b, a, k As Double
    Dim klo, khi As Long
    Dim X1, X2
    
    'trick: transpose to get rid of second dimension of range - with input data sitting in a regular table with records in rows
    xin = RangeToArray(DatapointsX)
    yin = RangeToArray(DatapointsY)
    
    n = UBound(xin)
    
    ReDim u(n - 1) As Single
    ReDim yt(n) As Single
    yt(1) = 0
    u(1) = 0
    For i = 2 To n - 1
        sig = (xin(i) - xin(i - 1)) / (xin(i + 1) - xin(i - 1))
        p = sig * yt(i - 1) + 2
        yt(i) = (sig - 1) / p
        u(i) = (yin(i + 1) - yin(i)) / (xin(i + 1) - xin(i)) - (yin(i) - yin(i - 1)) / (xin(i) - xin(i - 1))
        u(i) = (6 * u(i) / (xin(i + 1) - xin(i - 1)) - sig * u(i - 1)) / p
    Next i
    yt(n) = 0
    For i = n - 1 To 1 Step -1
        yt(i) = yt(i) * yt(i + 1) + u(i)
    Next i
    
    X1 = WorksheetFunction.XLookup(atX, xin, xin, , -1)
    X2 = WorksheetFunction.XLookup(atX, xin, xin, , 1)
    klo = WorksheetFunction.Match(X1, xin, 0)
    khi = WorksheetFunction.Match(X2, xin, 0)
    
    If (klo = khi) Then
        CubicSpline_SP = yin(klo)
    Else
        h = xin(khi) - xin(klo)
        a = (xin(khi) - atX) / h
        b = (atX - xin(klo)) / h
        CubicSpline_SP = a * yin(klo) + b * yin(khi) + ((a ^ 3 - a) * yt(klo) + (b ^ 3 - b) * yt(khi)) * (h ^ 2) / 6
    End If
End Function

Monotone cubic spline

Public Function MonotoneCubicSpline_SP(DatapointsX, DatapointsY, ByVal atX)
    
    Dim n As Long, i As Long, iRow, op, s As Long
    Dim xin As Variant, yin As Variant
    Dim m, mNext, dx, dxNext, common, c1, m_, invDx, common_, h, diff
    Dim c1s()
    Dim klo, khi As Long
    Dim X1, X2
    
    'trick: transpose to get rid of second dimension of range
    xin = RangeToArray(DatapointsX)
    yin = RangeToArray(DatapointsY)
    n = UBound(xin)
    
    'Get consecutive differences and slopes
    ReDim dxs(n)
    ReDim dys(n)
    ReDim ms(n)
    
    For i = 1 To n - 1
        dxs(i) = xin(i + 1) - xin(i)
        dys(i) = yin(i + 1) - yin(i)
        ms(i) = dys(i) / dxs(i)
    Next i
    
    'Get degree-1 coefficients
    ReDim c1s(1 To n + 1)
    c1s(1) = ms(1)
    For i = 1 To n - 1
        m = ms(i)
        mNext = ms(i + 1)
        If (m * mNext <= 0) Then
            c1s(i + 1) = 0
        Else
            dx = dxs(i)
            dxNext = dxs(i + 1)
            common = dx + dxNext
            c1s(i + 1) = 3 * common / ((common + dxNext) / m + (common + dx) / mNext)
        End If
    Next
    c1s(n) = ms(n - 1)
    
    'Get degree-2 and degree-3 coefficients
    ReDim c2s(n)
    ReDim c3s(n)
    For i = 1 To n - 1
        c1 = c1s(i)
        m_ = ms(i)
        invDx = 1 / dxs(i)
        common_ = c1 + c1s(i + 1) - m_ - m_
        c2s(i) = (m_ - c1 - common_) * invDx
        c3s(i) = common_ * invDx * invDx
    Next
    
    X1 = WorksheetFunction.XLookup(atX, xin, xin, , -1)
    X2 = WorksheetFunction.XLookup(atX, xin, xin, , 1)
    klo = WorksheetFunction.Match(X1, xin, 0)
    khi = WorksheetFunction.Match(X2, xin, 0)
    
    If (klo = khi) Then
        MonotoneCubicSpline_SP = yin(klo)
    Else
        If klo > khi Then 'X-axis values were in descending order
            Dim tmp
            tmp = klo
            klo = khi
            khi = tmp
        End If
        diff = (atX - xin(klo))
        MonotoneCubicSpline_SP = yin(klo) + c1s(klo) * diff + c2s(klo) * diff * diff + c3s(klo) * diff * diff * diff
    End If
End Function

Excel - Splines UDF

MS Excel VBA code to be put in a module - array functions

Regular function - Single Y value output only, given X value