[Solved]Cubic Spline interpolation: Porting VBA to OOo Basic

[johnroberts]

Hi

I am trying to migrate to OOo after years of Excel torture. I have a problem with a small VBA program that I've used many times successfully in Excel '97. It is a custom function that does cubic spline interpolation of data. For engineering applications this is a life saver. It is called as:
cubspline(<i>;<X_position_to interpolate_Y>;<LIST_of_X_input_values>;<LIST_of_corresponding_Y_input_values>)
where i=1 or 3 (usually 1) to choose between the typical and a non-overshoot constrained algorithm for the spline
The VBA code is contained in an XLS file courtesy of http://www.korf.co.uk/spline.xls. I am not very familiar with either VBA or OOo Basic. I tried using this on an OOo installation in SuSE 10.3 that is supposed to support VBA macros. When I try to use it, I get an OOo Basic error message 380 (don't know what that is ...). The original VBA code is:

Code: Select all

Attribute VB_Name = "Module1"
Option Explicit
Option Base 0


Function cubspline(Metode As Integer, xi As Double, xx As Object, yy As Object) As Double

Dim i As Integer
Dim yi As Double
Dim x() As Double
Dim y() As Double
Dim y2() As Double
Dim j As Integer

If Metode = 1 Then
  'Numerical Recipes are 1 based
  j = 0
Else
  'Others are 0 based
  j = -1
End If

For i = 1 To UBound(xx())
  If yy(i) <> "" Then
    j = j + 1
    ReDim Preserve x(j)
    ReDim Preserve y(j)
    x(j) = CDbl(xx(i))
    y(j) = CDbl(yy(i))
  End If
Next i

If Metode = 1 Then
  'NR cubic spline
  'Get y2
  ReDim y2(1 To UBound(x()))
  Call spline(x(), y(), UBound(x()), 10 ^ 30, 10 ^ 30, y2())
  'Get y
  Call splint(x(), y(), y2(), UBound(x()), xi, yi)
ElseIf Metode = 3 Then
  'Own cubic spline
  yi = SplineX3(xi, x(), y())
End If

'Return
cubspline = yi


End Function

Sub spline(x() As Double, y() As Double, N As Integer, yp1 As Double, ypn As Double, y2() As Double)

'Given arrays x(1:n) and y(1:n) containing a tabulated function, i.e., y i = f(xi), with
'x1<x2< :::<xN , and given values yp1 and ypn for the first derivative of the inter-
'polating function at points 1 and n, respectively, this routine returns an array y2(1:n) of
'length n which contains the second derivatives of the interpolating function at the tabulated
'points xi. If yp1 and/or ypn are equal to 1 * 10^30 or larger, the routine is signaled to set
'the corresponding boundary condition for a natural spline, with zero second derivative on
'that boundary.
'Parameter: NMAX is the largest anticipated value of n.

      Dim Nmax As Integer
      Nmax = 500
      
      Dim i As Integer
      Dim k As Integer
      
      Dim p As Double
      Dim qn As Double
      Dim sig As Double
      Dim un As Double
      ReDim u(Nmax) As Double
      
      'The lower boundary condition is set either to be natural
      If (yp1 > 9.9E+29) Then
        y2(1) = 0#
        u(1) = 0#
      Else
        'or else to have a specicied first derivative.
        y2(1) = -0.5
        u(1) = (3# / (x(2) - x(1))) * ((y(2) - y(1)) / (x(2) - x(1)) - yp1)
      End If
      
      'This is the decomposition loop of the tridiagonal
      'algorithm. y2 and u are used for temporary
      'storage of the decomposed factors.
      
      For i = 2 To N - 1
        sig = (x(i) - x(i - 1)) / (x(i + 1) - x(i - 1))
        p = sig * y2(i - 1) + 2#
        y2(i) = (sig - 1#) / p
        u(i) = (6# * ((y(i + 1) - y(i)) / (x(i + 1) - x(i)) - (y(i) - y(i - 1)) / (x(i) - x(i - 1))) / (x(i + 1) - x(i - 1)) - sig * u(i - 1)) / p
      Next i
      
      'The upper boundary condition is set either to be natural
      If (ypn > 9.9E+29) Then
        qn = 0#
        un = 0#
      Else
        'or else to have a specified first derivative.
        qn = 0.5
        un = (3# / (x(N) - x(N - 1))) * (ypn - (y(N) - y(N - 1)) / (x(N) - x(N - 1)))
      End If
      y2(N) = (un - qn * u(N - 1)) / (qn * y2(N - 1) + 1#)
      
      'This is the backsubstitution loop of the tridiagonal algorithm.
      For k = N - 1 To 1 Step -1
        y2(k) = y2(k) * y2(k + 1) + u(k)
      Next k
      
End Sub
      
Sub splint(xa() As Double, ya() As Double, y2a() As Double, N As Integer, x As Double, y As Double)
      
'Given the arrays xa(1:n) and ya(1:n) of length n, which tabulate a function (with the
'xai 's in order), and given the array y2a(1:n), which is the output from spline above,
'and given a value of x, this routine returns a cubic-spline interpolated value y.
      
      Dim k As Integer
      Dim khi As Integer
      Dim klo As Integer
      
      Dim A As Double
      Dim B As Double
      Dim h As Double
      
      'We will the right place in the table by means of bisection.
      klo = 1
      khi = N
      
      While (khi - klo > 1)
        k = (khi + klo) / 2
        If (xa(k) > x) Then
          khi = k
        Else
          klo = k
        End If
      Wend
      
      'klo and khi now bracket the input value of x.
      h = xa(khi) - xa(klo)
      If (h = 0) Then MsgBox ("bad xa input in splint")
      
      'Cubic spline polynomial is now evaluated.
      A = (xa(khi) - x) / h
      B = (x - xa(klo)) / h
      y = A * ya(klo) + B * ya(khi) + ((A ^ 3 - A) * y2a(klo) + (B ^ 3 - B) * y2a(khi)) * (h ^ 2) / 6#
      
End Sub

Public Function SplineX3(x As Double, xx() As Double, yy() As Double) As Double
'|-------------------------------------------------------------------------------
'| Function returns y value for a corresponding x value, based on cubic spline.
'| Will never oscillates or overshoot. No need to solve matrix.
'| Also calculate constants for cubic in case needed (for integration).
'|
'| xx(0 to No_of_lines) is x values
'|    * Must be unique (no two consequetive ones the same)
'|    * Must be in ascending order
'|    * No of lines = Number of points - 1
'| yy(0 to No_of_lines) is y values
'|
'| Uses function dxx to prevent div by zero.
'|
'| Developer: C Kruger, Guildford, UK
'| Date: December 2001
'|-------------------------------------------------------------------------------

Dim i As Integer
Dim j As Integer
Dim Nmax As Integer
Dim Num As Integer

'1st and 2nd derivative for left and right ends of line
Dim gxx(0 To 1) As Double
Dim ggxx(0 To 1) As Double

'Constants for cubic equations
Dim A As Double     'Also for linear extrapolation
Dim B As Double     'Also for linear extrapolation
Dim C As Double
Dim D As Double

'Number of lines = points - 1
Nmax = UBound(xx())

'(1a) Find LineNumber or segment. Linear extrapolate if outside range.
Num = 0
If x < xx(0) Or x > xx(Nmax) Then
  'X outisde range. Linear interpolate
  'Below min or max?
  If x < xx(0) Then Num = 1 Else Num = Nmax
  B = (yy(Num) - yy(Num - 1)) / dxx(xx(Num), xx(Num - 1))
  A = yy(Num) - B * xx(Num)
  SplineX3 = A + B * x
  Exit Function
  
'(1b) Find LineNumber or segment. Linear extrapolate if outside range.
Else
  'X in range. Get line.
  For i = 1 To Nmax
    If x <= xx(i) Then
      Num = i
      Exit For
    End If
  Next i
End If

'(2) Calc first derivative (slope) for intermediate points
For j = 0 To 1          'Two points around line
  i = Num - 1 + j
  If i = 0 Or i = Nmax Then
    'Set very large slope at ends
    gxx(j) = 10 ^ 30
  ElseIf (yy(i + 1) - yy(i) = 0) Or (yy(i) - yy(i - 1) = 0) Then
    'Only check for 0 dy. dx assumed NEVER equals 0 !
    gxx(j) = 0
  ElseIf ((xx(i + 1) - xx(i)) / (yy(i + 1) - yy(i)) + (xx(i) - xx(i - 1)) / (yy(i) - yy(i - 1))) = 0 Then
    'Pos PLUS neg slope is 0. Prevent div by zero.
    gxx(j) = 0
  ElseIf (yy(i + 1) - yy(i)) * (yy(i) - yy(i - 1)) < 0 Then
    'Pos AND neg slope, assume slope = 0 to prevent overshoot
    gxx(j) = 0
  Else
    'Calculate an average slope for point based on connecting lines
    gxx(j) = 2 / (dxx(xx(i + 1), xx(i)) / (yy(i + 1) - yy(i)) + dxx(xx(i), xx(i - 1)) / (yy(i) - yy(i - 1)))
  End If
Next j

'(3) Reset first derivative (slope) at first and last point
If Num = 1 Then
  'First point has 0 2nd derivative
  gxx(0) = 3 / 2 * (yy(Num) - yy(Num - 1)) / dxx(xx(Num), xx(Num - 1)) - gxx(1) / 2
End If
If Num = Nmax Then
  'Last point has 0 2nd derivative
  gxx(1) = 3 / 2 * (yy(Num) - yy(Num - 1)) / dxx(xx(Num), xx(Num - 1)) - gxx(0) / 2
End If

'(4) Calc second derivative at points
ggxx(0) = -2 * (gxx(1) + 2 * gxx(0)) / dxx(xx(Num), xx(Num - 1)) + 6 * (yy(Num) - yy(Num - 1)) / dxx(xx(Num), xx(Num - 1)) ^ 2
ggxx(1) = 2 * (2 * gxx(1) + gxx(0)) / dxx(xx(Num), xx(Num - 1)) - 6 * (yy(Num) - yy(Num - 1)) / dxx(xx(Num), xx(Num - 1)) ^ 2

'(5) Calc constants for cubic
D = 1 / 6 * (ggxx(1) - ggxx(0)) / dxx(xx(Num), xx(Num - 1))
C = 1 / 2 * (xx(Num) * ggxx(0) - xx(Num - 1) * ggxx(1)) / dxx(xx(Num), xx(Num - 1))
B = (yy(Num) - yy(Num - 1) - C * (xx(Num) ^ 2 - xx(Num - 1) ^ 2) - D * (xx(Num) ^ 3 - xx(Num - 1) ^ 3)) / dxx(xx(Num), xx(Num - 1))
A = yy(Num - 1) - B * xx(Num - 1) - C * xx(Num - 1) ^ 2 - D * xx(Num - 1) ^ 3

'Return function
SplineX3 = A + B * x + C * x ^ 2 + D * x ^ 3

''Alternative method based on Numerical Recipes.
''Shorter but does not calc cubic constants A, B, C, D
'i = Num
'A = (xx(i) - x) / (xx(i) - xx(i - 1))
'B = 1 - A
'Cy = 1 / 6 * (A ^ 3 - A) * (6 * (yy(i) - yy(i - 1)) - 2 * (gxx(i) + 2 * gxx(i - 1)) * (xx(i) - xx(i - 1)))
'Dy = 1 / 6 * (B ^ 3 - B) * (2 * (2 * gxx(i) + gxx(i - 1)) * (xx(i) - xx(i - 1)) - 6 * (yy(i) - yy(i - 1)))
''Return function
'SplineX3 = A * yy(i - 1) + B * yy(i) + Cy + Dy

End Function

Public Function dxx(x1 As Double, x0 As Double) As Double
'Calc Xi - Xi-1 to prevent div by zero

dxx = x1 - x0
If dxx = 0 Then dxx = 10 ^ 30

End Function

and after OOo imports the code it is interpreted as (as seen from the OOo Basic Editor):

Code: Select all

Attribute VB_Name = "Module1"
Option Explicit
Option Base 0


Function cubspline(Metode As Integer, xi As Double, xx As Object, yy As Object) As Double

Dim i As Integer
Dim yi As Double
Dim x() As Double
Dim y() As Double
Dim y2() As Double
Dim j As Integer

If Metode = 1 Then
  'Numerical Recipes are 1 based
  j = 0
Else
  'Others are 0 based
  j = -1
End If

For i = 1 To UBound(xx())
  If yy(i) <> "" Then
    j = j + 1
    ReDim Preserve x(j)
    ReDim Preserve y(j)
    x(j) = CDbl(xx(i))
    y(j) = CDbl(yy(i))
  End If
Next i

If Metode = 1 Then
  'NR cubic spline
  'Get y2
  ReDim y2(1 To UBound(x()))
  Call spline(x(), y(), UBound(x()), 10 ^ 30, 10 ^ 30, y2())
  'Get y
  Call splint(x(), y(), y2(), UBound(x()), xi, yi)
ElseIf Metode = 3 Then
  'Own cubic spline
  yi = SplineX3(xi, x(), y())
End If

'Return
cubspline = yi


End Function

Sub spline(x() As Double, y() As Double, N As Integer, yp1 As Double, ypn As Double, y2() As Double)

'Given arrays x(1:n) and y(1:n) containing a tabulated function, i.e., y i = f(xi), with
'x1<x2< :::<xN , and given values yp1 and ypn for the first derivative of the inter-
'polating function at points 1 and n, respectively, this routine returns an array y2(1:n) of
'length n which contains the second derivatives of the interpolating function at the tabulated
'points xi. If yp1 and/or ypn are equal to 1 * 10^30 or larger, the routine is signaled to set
'the corresponding boundary condition for a natural spline, with zero second derivative on
'that boundary.
'Parameter: NMAX is the largest anticipated value of n.

      Dim Nmax As Integer
      Nmax = 500
      
      Dim i As Integer
      Dim k As Integer
      
      Dim p As Double
      Dim qn As Double
      Dim sig As Double
      Dim un As Double
      ReDim u(Nmax) As Double
      
      'The lower boundary condition is set either to be natural
      If (yp1 > 9.9E+29) Then
        y2(1) = 0#
        u(1) = 0#
      Else
        'or else to have a specicied first derivative.
        y2(1) = -0.5
        u(1) = (3# / (x(2) - x(1))) * ((y(2) - y(1)) / (x(2) - x(1)) - yp1)
      End If
      
      'This is the decomposition loop of the tridiagonal
      'algorithm. y2 and u are used for temporary
      'storage of the decomposed factors.
      
      For i = 2 To N - 1
        sig = (x(i) - x(i - 1)) / (x(i + 1) - x(i - 1))
        p = sig * y2(i - 1) + 2#
        y2(i) = (sig - 1#) / p
        u(i) = (6# * ((y(i + 1) - y(i)) / (x(i + 1) - x(i)) - (y(i) - y(i - 1)) / (x(i) - x(i - 1))) / (x(i + 1) - x(i - 1)) - sig * u(i - 1)) / p
      Next i
      
      'The upper boundary condition is set either to be natural
      If (ypn > 9.9E+29) Then
        qn = 0#
        un = 0#
      Else
        'or else to have a specified first derivative.
        qn = 0.5
        un = (3# / (x(N) - x(N - 1))) * (ypn - (y(N) - y(N - 1)) / (x(N) - x(N - 1)))
      End If
      y2(N) = (un - qn * u(N - 1)) / (qn * y2(N - 1) + 1#)
      
      'This is the backsubstitution loop of the tridiagonal algorithm.
      For k = N - 1 To 1 Step -1
        y2(k) = y2(k) * y2(k + 1) + u(k)
      Next k
      
End Sub
      
Sub splint(xa() As Double, ya() As Double, y2a() As Double, N As Integer, x As Double, y As Double)
      
'Given the arrays xa(1:n) and ya(1:n) of length n, which tabulate a function (with the
'xai 's in order), and given the array y2a(1:n), which is the output from spline above,
'and given a value of x, this routine returns a cubic-spline interpolated value y.
      
      Dim k As Integer
      Dim khi As Integer
      Dim klo As Integer
      
      Dim A As Double
      Dim B As Double
      Dim h As Double
      
      'We will the right place in the table by means of bisection.
      klo = 1
      khi = N
      
      While (khi - klo > 1)
        k = (khi + klo) / 2
        If (xa(k) > x) Then
          khi = k
        Else
          klo = k
        End If
      Wend
      
      'klo and khi now bracket the input value of x.
      h = xa(khi) - xa(klo)
      If (h = 0) Then MsgBox ("bad xa input in splint")
      
      'Cubic spline polynomial is now evaluated.
      A = (xa(khi) - x) / h
      B = (x - xa(klo)) / h
      y = A * ya(klo) + B * ya(khi) + ((A ^ 3 - A) * y2a(klo) + (B ^ 3 - B) * y2a(khi)) * (h ^ 2) / 6#
      
End Sub

Public Function SplineX3(x As Double, xx() As Double, yy() As Double) As Double
'|-------------------------------------------------------------------------------
'| Function returns y value for a corresponding x value, based on cubic spline.
'| Will never oscillates or overshoot. No need to solve matrix.
'| Also calculate constants for cubic in case needed (for integration).
'|
'| xx(0 to No_of_lines) is x values
'|    * Must be unique (no two consequetive ones the same)
'|    * Must be in ascending order
'|    * No of lines = Number of points - 1
'| yy(0 to No_of_lines) is y values
'|
'| Uses function dxx to prevent div by zero.
'|
'| Developer: C Kruger, Guildford, UK
'| Date: December 2001
'|-------------------------------------------------------------------------------

Dim i As Integer
Dim j As Integer
Dim Nmax As Integer
Dim Num As Integer

'1st and 2nd derivative for left and right ends of line
Dim gxx(0 To 1) As Double
Dim ggxx(0 To 1) As Double

'Constants for cubic equations
Dim A As Double     'Also for linear extrapolation
Dim B As Double     'Also for linear extrapolation
Dim C As Double
Dim D As Double

'Number of lines = points - 1
Nmax = UBound(xx())

'(1a) Find LineNumber or segment. Linear extrapolate if outside Dim oSheet as Object
oSheet = ThisComponent.CurrentController.ActiveSheet
oSheet.getCellRangeByName($1).
Num = 0
If x < xx(0) Or x > xx(Nmax) Then
  'X outisde Dim oSheet as Object
oSheet = ThisComponent.CurrentController.ActiveSheet
oSheet.getCellRangeByName($1). Linear interpolate
  'Below min or max?
  If x < xx(0) Then Num = 1 Else Num = Nmax
  B = (yy(Num) - yy(Num - 1)) / dxx(xx(Num), xx(Num - 1))
  A = yy(Num) - B * xx(Num)
  SplineX3 = A + B * x
  Exit Function
  
'(1b) Find LineNumber or segment. Linear extrapolate if outside Dim oSheet as Object
oSheet = ThisComponent.CurrentController.ActiveSheet
oSheet.getCellRangeByName($1).
Else
  'X in Dim oSheet as Object
oSheet = ThisComponent.CurrentController.ActiveSheet
oSheet.getCellRangeByName($1). Get line.
  For i = 1 To Nmax
    If x <= xx(i) Then
      Num = i
      Exit For
    End If
  Next i
End If

'(2) Calc first derivative (slope) for intermediate points
For j = 0 To 1          'Two points around line
  i = Num - 1 + j
  If i = 0 Or i = Nmax Then
    'Set very large slope at ends
    gxx(j) = 10 ^ 30
  ElseIf (yy(i + 1) - yy(i) = 0) Or (yy(i) - yy(i - 1) = 0) Then
    'Only check for 0 dy. dx assumed NEVER equals 0 !
    gxx(j) = 0
  ElseIf ((xx(i + 1) - xx(i)) / (yy(i + 1) - yy(i)) + (xx(i) - xx(i - 1)) / (yy(i) - yy(i - 1))) = 0 Then
    'Pos PLUS neg slope is 0. Prevent div by zero.
    gxx(j) = 0
  ElseIf (yy(i + 1) - yy(i)) * (yy(i) - yy(i - 1)) < 0 Then
    'Pos AND neg slope, assume slope = 0 to prevent overshoot
    gxx(j) = 0
  Else
    'Calculate an average slope for point based on connecting lines
    gxx(j) = 2 / (dxx(xx(i + 1), xx(i)) / (yy(i + 1) - yy(i)) + dxx(xx(i), xx(i - 1)) / (yy(i) - yy(i - 1)))
  End If
Next j

'(3) Reset first derivative (slope) at first and last point
If Num = 1 Then
  'First point has 0 2nd derivative
  gxx(0) = 3 / 2 * (yy(Num) - yy(Num - 1)) / dxx(xx(Num), xx(Num - 1)) - gxx(1) / 2
End If
If Num = Nmax Then
  'Last point has 0 2nd derivative
  gxx(1) = 3 / 2 * (yy(Num) - yy(Num - 1)) / dxx(xx(Num), xx(Num - 1)) - gxx(0) / 2
End If

'(4) Calc second derivative at points
ggxx(0) = -2 * (gxx(1) + 2 * gxx(0)) / dxx(xx(Num), xx(Num - 1)) + 6 * (yy(Num) - yy(Num - 1)) / dxx(xx(Num), xx(Num - 1)) ^ 2
ggxx(1) = 2 * (2 * gxx(1) + gxx(0)) / dxx(xx(Num), xx(Num - 1)) - 6 * (yy(Num) - yy(Num - 1)) / dxx(xx(Num), xx(Num - 1)) ^ 2

'(5) Calc constants for cubic
D = 1 / 6 * (ggxx(1) - ggxx(0)) / dxx(xx(Num), xx(Num - 1))
C = 1 / 2 * (xx(Num) * ggxx(0) - xx(Num - 1) * ggxx(1)) / dxx(xx(Num), xx(Num - 1))
B = (yy(Num) - yy(Num - 1) - C * (xx(Num) ^ 2 - xx(Num - 1) ^ 2) - D * (xx(Num) ^ 3 - xx(Num - 1) ^ 3)) / dxx(xx(Num), xx(Num - 1))
A = yy(Num - 1) - B * xx(Num - 1) - C * xx(Num - 1) ^ 2 - D * xx(Num - 1) ^ 3

'Return function
SplineX3 = A + B * x + C * x ^ 2 + D * x ^ 3

''Alternative method based on Numerical Recipes.
''Shorter but does not calc cubic constants A, B, C, D
'i = Num
'A = (xx(i) - x) / (xx(i) - xx(i - 1))
'B = 1 - A
'Cy = 1 / 6 * (A ^ 3 - A) * (6 * (yy(i) - yy(i - 1)) - 2 * (gxx(i) + 2 * gxx(i - 1)) * (xx(i) - xx(i - 1)))
'Dy = 1 / 6 * (B ^ 3 - B) * (2 * (2 * gxx(i) + gxx(i - 1)) * (xx(i) - xx(i - 1)) - 6 * (yy(i) - yy(i - 1)))
''Return function
'SplineX3 = A * yy(i - 1) + B * yy(i) + Cy + Dy

End Function

Public Function dxx(x1 As Double, x0 As Double) As Double
'Calc Xi - Xi-1 to prevent div by zero

dxx = x1 - x0
If dxx = 0 Then dxx = 10 ^ 30

End Function

I am trying to figure out what goes wrong with this. If anyone can help, pointing out what is wrong - or alternatively - pointing to a similar program in OOo Basic that does work.


Many Thanks in advance
Best Regards
JR

[Hagar Delest]

Hi,

I move your thread in the Macro forum.

Note that VBA support is limited so there might be functions that are not supported. But wait for some macro gurus.

[Villeroy]

Hi

I am trying to migrate to OOo after years of Excel torture. I have a problem with a small VBA program that I've used many times successfully in Excel '97.

Welcome to the other spreadsheet hell!
What makes you think that Calc is "better" than Excel '97? I don't think so. If you definitively want to do this type of analysis in a spreadsheet, I would suggest Gnumeric. It is small and fast, which implies that it has not that many features but it is well known for it's wide range of analysis tools and Excel compatibility (regarding most aspects of the formula language and naming converntions).
Since version 1.7.11 it has cubic splines built-in.

[TerryE]

JR, You've got me a little confused. The programme is essentially just Basic, and GWBasic at that, so it should work unchanged in OOo. The only issue is to do with the calling parameters and how it is called. Am I right in assuming that you want to call this as a spreadsheet function? You shouldn't need to use UNO API calls such as getCellRangeByName at all, and in fact if you are using this as a spreadsheet function you shouldn't.

Onr of the main differences between VBA and OOoBasic is that if you refer to, say, B3:B20 as a parameter in such a function, in VBA this passed by reference as a range object. In OOo this is passed as

• A variant scalar if this refers to a single cell
• An array (1 to 1, 1 to N) if this refers to a row
• An array (1 to M, 1 to 1) if this refers to a column
• An array (1 to M, 1 to N) if this refers to a table

So you need to make the initialising code which sets x and y handle these three cases.

[johnroberts]

Hi guys

Am I right in assuming that you want to call this as a spreadsheet function?

Assumption correct... The idea being of the interpolation scheme function to be available just like any other implicit function. I am not that familiar with Basic (any variant of Basic...); some C/C++ and FORTRAN only

I would suggest Gnumeric...Since version 1.7.11 it has cubic splines built-in

. I have Gnumeric running on Puppy LiveCDs... What is the function's name?

- Do you see something that is obviously wrong with the "equivalent" OOo Basic code?
- What is this error 380?
- Is there a way to run it in OOo in debug mode to see what happens?

[Villeroy]

I have Gnumeric running on Puppy LiveCDs... What is the function's name?

I have a pretty old version of Gnumeric. This is what I found via Google: http://lists.altlinux.org/pipermail/sis ... 38336.html

...
libgnomeoffice - Library for writing gnome office programs
* Tue Jul 31 2007 Alexey Morsov <swi на altlinux> 0.4.2-alt1
- Avoid crash if libxml2 returns ERROR or NONE when guessing encoding.
- do not leak the mime type
- fixed signedness issue spotted with help of sparse
- needed for gnumeric 1.7.11
* Sat Jun 02 2007 Alexey Morsov <swi на altlinux> 0.4.0-alt1
- new version for Gnumeric 1.7.10
- added new cubic spline support
- new convenience function
- fixes
...

Have a look in menu "Tools". Gnumeric implements much of its analysis via wizards.

[TerryE]

Do you see something that is obviously wrong with the "equivalent" OOo Basic code?

Yes, all of the changes are irrelevant. As I said the original VBA is a closer starting point and you need to rewrite the section that I discussed.

[Villeroy]

This module is derived from the original VBA and returns nothing useful, but it runs.

Code: Select all

Option Explicit

Function cubspline(Metode, xi, xx, yy)

Dim i As Integer
Dim yi As Double
Dim x() As Double
Dim y() As Double
Dim y2() As Double
Dim j As Integer

If Metode = 1 Then
  'Numerical Recipes are 1 based
  j = 0
Else
  'Others are 0 based
  j = -1
End If

For i = 1 To UBound(xx())
  If yy(i,1) <> "" Then
    j = j + 1
    ReDim Preserve x(j)
    ReDim Preserve y(j)
    x(j) = CDbl(xx(i,1))
    y(j) = CDbl(yy(i,1))
  End If
Next i

If Metode = 1 Then
  'NR cubic spline
  'Get y2
  ReDim y2(1 To UBound(x()))
  Call spline(x(), y(), UBound(x()), 10 ^ 30, 10 ^ 30, y2())
  'Get y
  Call splint(x(), y(), y2(), UBound(x()), xi, yi)
ElseIf Metode = 3 Then
  'Own cubic spline
  yi = SplineX3(xi, x(), y())
End If

'Return
cubspline = yi


End Function

Sub spline(x() As Double, y() As Double, N As Integer, yp1 As Double, ypn As Double, y2() As Double)

'Given arrays x(1:n) and y(1:n) containing a tabulated function, i.e., y i = f(xi), with
'x1<x2< :::<xN , and given values yp1 and ypn for the first derivative of the inter-
'polating function at points 1 and n, respectively, this routine returns an array y2(1:n) of
'length n which contains the second derivatives of the interpolating function at the tabulated
'points xi. If yp1 and/or ypn are equal to 1 * 10^30 or larger, the routine is signaled to set
'the corresponding boundary condition for a natural spline, with zero second derivative on
'that boundary.
'Parameter: NMAX is the largest anticipated value of n.

      Dim Nmax As Integer
      Nmax = 500
     
      Dim i As Integer
      Dim k As Integer
     
      Dim p As Double
      Dim qn As Double
      Dim sig As Double
      Dim un As Double
      Dim u(Nmax) As Double
     
      'The lower boundary condition is set either to be natural
      If (yp1 > 9.9E+29) Then
        y2(1) = 0#
        u(1) = 0#
      Else
        'or else to have a specicied first derivative.
        y2(1) = -0.5
        u(1) = (3# / (x(2) - x(1))) * ((y(2) - y(1)) / (x(2) - x(1)) - yp1)
      End If
     
      'This is the decomposition loop of the tridiagonal
      'algorithm. y2 and u are used for temporary
      'storage of the decomposed factors.
     
      For i = 2 To N - 1
        sig = (x(i) - x(i - 1)) / (x(i + 1) - x(i - 1))
        p = sig * y2(i - 1) + 2#
        y2(i) = (sig - 1#) / p
        u(i) = (6# * ((y(i + 1) - y(i)) / (x(i + 1) - x(i)) - (y(i) - y(i - 1)) / (x(i) - x(i - 1))) / (x(i + 1) - x(i - 1)) - sig * u(i - 1)) / p
      Next i
     
      'The upper boundary condition is set either to be natural
      If (ypn > 9.9E+29) Then
        qn = 0#
        un = 0#
      Else
        'or else to have a specified first derivative.
        qn = 0.5
        un = (3# / (x(N) - x(N - 1))) * (ypn - (y(N) - y(N - 1)) / (x(N) - x(N - 1)))
      End If
      y2(N) = (un - qn * u(N - 1)) / (qn * y2(N - 1) + 1#)
     
      'This is the backsubstitution loop of the tridiagonal algorithm.
      For k = N - 1 To 1 Step -1
        y2(k) = y2(k) * y2(k + 1) + u(k)
      Next k
     
End Sub
     
Sub splint(xa() As Double, ya() As Double, y2a() As Double, N As Integer, x As Double, y As Double)
     
'Given the arrays xa(1:n) and ya(1:n) of length n, which tabulate a function (with the
'xai 's in order), and given the array y2a(1:n), which is the output from spline above,
'and given a value of x, this routine returns a cubic-spline interpolated value y.
     
      Dim k As Integer
      Dim khi As Integer
      Dim klo As Integer
     
      Dim A As Double
      Dim B As Double
      Dim h As Double
     
      'We will the right place in the table by means of bisection.
      klo = 1
      khi = N
     
      While (khi - klo > 1)
        k = (khi + klo) / 2
        If (xa(k) > x) Then
          khi = k
        Else
          klo = k
        End If
      Wend
     
      'klo and khi now bracket the input value of x.
      h = xa(khi) - xa(klo)
      If (h = 0) Then MsgBox ("bad xa input in splint")
     
      'Cubic spline polynomial is now evaluated.
      A = (xa(khi) - x) / h
      B = (x - xa(klo)) / h
      y = A * ya(klo) + B * ya(khi) + ((A ^ 3 - A) * y2a(klo) + (B ^ 3 - B) * y2a(khi)) * (h ^ 2) / 6#
     
End Sub

Public Function SplineX3(x As Double, xx() As Double, yy() As Double) As Double
'|-------------------------------------------------------------------------------
'| Function returns y value for a corresponding x value, based on cubic spline.
'| Will never oscillates or overshoot. No need to solve matrix.
'| Also calculate constants for cubic in case needed (for integration).
'|
'| xx(0 to No_of_lines) is x values
'|    * Must be unique (no two consequetive ones the same)
'|    * Must be in ascending order
'|    * No of lines = Number of points - 1
'| yy(0 to No_of_lines) is y values
'|
'| Uses function dxx to prevent div by zero.
'|
'| Developer: C Kruger, Guildford, UK
'| Date: December 2001
'|-------------------------------------------------------------------------------

Dim i As Integer
Dim j As Integer
Dim Nmax As Integer
Dim Num As Integer

'1st and 2nd derivative for left and right ends of line
Dim gxx(0 To 1) As Double
Dim ggxx(0 To 1) As Double

'Constants for cubic equations
Dim A As Double     'Also for linear extrapolation
Dim B As Double     'Also for linear extrapolation
Dim C As Double
Dim D As Double

'Number of lines = points - 1
Nmax = UBound(xx())

'(1a) Find LineNumber or segment. Linear extrapolate if outside range.
Num = 0
If x < xx(0) Or x > xx(Nmax) Then
  'X outisde range. Linear interpolate
  'Below min or max?
  If x < xx(0) Then Num = 1 Else Num = Nmax
  B = (yy(Num) - yy(Num - 1)) / dxx(xx(Num), xx(Num - 1))
  A = yy(Num) - B * xx(Num)
  SplineX3 = A + B * x
  Exit Function
 
'(1b) Find LineNumber or segment. Linear extrapolate if outside range.
Else
  'X in range. Get line.
  For i = 1 To Nmax
    If x <= xx(i) Then
      Num = i
      Exit For
    End If
  Next i
End If

'(2) Calc first derivative (slope) for intermediate points
For j = 0 To 1          'Two points around line
  i = Num - 1 + j
  If i = 0 Or i = Nmax Then
    'Set very large slope at ends
    gxx(j) = 10 ^ 30
  ElseIf (yy(i + 1) - yy(i) = 0) Or (yy(i) - yy(i - 1) = 0) Then
    'Only check for 0 dy. dx assumed NEVER equals 0 !
    gxx(j) = 0
  ElseIf ((xx(i + 1) - xx(i)) / (yy(i + 1) - yy(i)) + (xx(i) - xx(i - 1)) / (yy(i) - yy(i - 1))) = 0 Then
    'Pos PLUS neg slope is 0. Prevent div by zero.
    gxx(j) = 0
  ElseIf (yy(i + 1) - yy(i)) * (yy(i) - yy(i - 1)) < 0 Then
    'Pos AND neg slope, assume slope = 0 to prevent overshoot
    gxx(j) = 0
  Else
    'Calculate an average slope for point based on connecting lines
    gxx(j) = 2 / (dxx(xx(i + 1), xx(i)) / (yy(i + 1) - yy(i)) + dxx(xx(i), xx(i - 1)) / (yy(i) - yy(i - 1)))
  End If
Next j

'(3) Reset first derivative (slope) at first and last point
If Num = 1 Then
  'First point has 0 2nd derivative
  gxx(0) = 3 / 2 * (yy(Num) - yy(Num - 1)) / dxx(xx(Num), xx(Num - 1)) - gxx(1) / 2
End If
If Num = Nmax Then
  'Last point has 0 2nd derivative
  gxx(1) = 3 / 2 * (yy(Num) - yy(Num - 1)) / dxx(xx(Num), xx(Num - 1)) - gxx(0) / 2
End If

'(4) Calc second derivative at points
ggxx(0) = -2 * (gxx(1) + 2 * gxx(0)) / dxx(xx(Num), xx(Num - 1)) + 6 * (yy(Num) - yy(Num - 1)) / dxx(xx(Num), xx(Num - 1)) ^ 2
ggxx(1) = 2 * (2 * gxx(1) + gxx(0)) / dxx(xx(Num), xx(Num - 1)) - 6 * (yy(Num) - yy(Num - 1)) / dxx(xx(Num), xx(Num - 1)) ^ 2

'(5) Calc constants for cubic
D = 1 / 6 * (ggxx(1) - ggxx(0)) / dxx(xx(Num), xx(Num - 1))
C = 1 / 2 * (xx(Num) * ggxx(0) - xx(Num - 1) * ggxx(1)) / dxx(xx(Num), xx(Num - 1))
B = (yy(Num) - yy(Num - 1) - C * (xx(Num) ^ 2 - xx(Num - 1) ^ 2) - D * (xx(Num) ^ 3 - xx(Num - 1) ^ 3)) / dxx(xx(Num), xx(Num - 1))
A = yy(Num - 1) - B * xx(Num - 1) - C * xx(Num - 1) ^ 2 - D * xx(Num - 1) ^ 3

'Return function
SplineX3 = A + B * x + C * x ^ 2 + D * x ^ 3

''Alternative method based on Numerical Recipes.
''Shorter but does not calc cubic constants A, B, C, D
'i = Num
'A = (xx(i) - x) / (xx(i) - xx(i - 1))
'B = 1 - A
'Cy = 1 / 6 * (A ^ 3 - A) * (6 * (yy(i) - yy(i - 1)) - 2 * (gxx(i) + 2 * gxx(i - 1)) * (xx(i) - xx(i - 1)))
'Dy = 1 / 6 * (B ^ 3 - B) * (2 * (2 * gxx(i) + gxx(i - 1)) * (xx(i) - xx(i - 1)) - 6 * (yy(i) - yy(i - 1)))
''Return function
'SplineX3 = A * yy(i - 1) + B * yy(i) + Cy + Dy

End Function

Public Function dxx(x1 As Double, x0 As Double) As Double
'Calc Xi - Xi-1 to prevent div by zero

dxx = x1 - x0
If dxx = 0 Then dxx = 10 ^ 30

End Function

[johnroberts]

The version of Gnumeric in Puppy is also very old...1.6.3. "Cubic Spline support" is not very enlightening in the sense that it may actually mean having cubic spline smoothing as an option when creating X-Y plot graphs (this is the most widely used case...). I googled persistently on "Gnumeric" and "cubic spline interpolation" and found a couple of references on "Time Series Analysis Functions plugin" for Gnumeric. If this is a standard item of the current version or not, or an additional extra library, or something else that the user must install manually, eludes me. Looking for it led to some package with "extras" in Sourceforge but I don't know if this is it. It gave me the impression that it's comprised of Python scripts... Not having a recent Gnumeric, i cannot test it...

Yes, all of the changes are irrelevant. As I said the original VBA is a closer starting point and you need to rewrite the section that I discussed.

This module is derived from the original VBA and returns nothing useful, but it runs.

That's what I was afraid of. That it actually turned into something completely bogus... . (...hmm...and they claim the version of OOo in SuSE is compiled with VBA support...go figure...)

Welcome to the other spreadsheet hell!

A better one, instead of the "Error 380" MessageBox I've been getting, would be "Error 666: Successfully misguided user"

The programme is essentially just Basic, and GWBasic at that

I think I have an old GWBasic book stashed somewhere, together with the old DOS books. Will have a look...

One of the main differences between VBA and OOoBasic is that if you refer to, say, B3:B20 as a parameter in such a function, in VBA this passed by reference as a range object. In OOo this is passed as

* A variant scalar if this refers to a single cell
* An array (1 to 1, 1 to N) if this refers to a row
* An array (1 to M, 1 to 1) if this refers to a column
* An array (1 to M, 1 to N) if this refers to a table

The input <X_position_to interpolate_Y>, would (should) always be a single cell. And the <LIST_of_X_input_values>, <LIST_of_corresponding_Y_input_values> are always 1-dimensional arrays, usually columns (wouldn't hurt if they could be columns or rows - regardless, to leave open the option of formatting to the user...). From a math POV, only requirements would be that they should have the same dimension and that values in the list of X ordinates are monotonic...

Where could I find a good book on programming OOo Basic? I had a look in Amazon and my local technical bookstore, but I found only the "usual" "The complete guide to xyz..." type books. I have a rather bitter experience from those...They are typically 1000 pages, covering all and nothing with only 3-4 pages to programming script languages, like

"...You can customize your wonderful XYZ software with macros and scripts..."

and maybe a "Hello world" sample and that is all...

What would you recommend?

[Villeroy]

Since the functions do not require API-calls, you have all the required documentation in the online help and here http://wiki.services.openoffice.org/wik ... nd_Dialogs
The Basic-IDE is basic indeed. Learn about stop-marks and variable watching: http://wiki.services.openoffice.org/wik ... d_Debugger

Debug one instance of the formula, so the code won't be called a thousand times.
You need to know that cell functions pass over arrays (xx and yy in your code) as 2-D structures. This is why I changed "x(j) = CDbl(xx(i))" to "x(j) = CDbl(xx(i,1))" since the former assumed a flat array.
Like almost all VBA code I ever read, this code is of very poor quality. In particular it is extremely inefficient. It should be rewritten as array function.

[johnroberts]

Thanks for the input! Will look into it... What is not clear to me yet (regarding OOo Basic), is how a function, such as the one under discussion, receives data from the spreadsheet. The calculation part should be quite straightforward...

Like almost all VBA code I ever read, this code is of very poor quality.

I could say that this is not my code, but that would be mean and ungrateful to the people that wrote it. Even if rough in appearance, (in Excel at least...) it does what it was supposed to do, well. This is really an essential tool when using spreadsheets in Engineering applications. Will revert for advice after doing some reading...

[Villeroy]

It is your code as soon as you are running it against your data. If your code was constructed by someone else, it's all matter of trust. I could make it run without errors but returning wrong results (zeros).
Create a test spreadsheet with some data.
Copy my version of the code.
Tools>Macros>Organize>Basic...
Pick your document's macro container, library "Standard" and hit [New...] in order to add a new module in library "Standard" of your test doc.
Select all default content of the new module (Ctrl+A)
Paste the code from the clipboard (Ctrl+V)
Save the document
Test the function as you did in Excel.

This macro language ("StarBasic") is much more primitive than VBA, in particular it is not that optimized. When you run the function against some 100 x/y-pairs, you can watch the sheet recalculate in seconds, due to the fact that the very same array of numbers is processed more than once for each of the 100 formulas. I did not start debugging that shit.

Look what I've found:
http://quantlib.org/
http://quantlib.org/docs.shtml
http://quantlibaddin.org/

[johnroberts]

It is your code as soon as you are running it against your data. If your code was constructed by someone else, it's all matter of trust. I could make it run without errors but returning wrong results (zeros).

...Cannot claim ownership over the code, even though data is mine...By using it I am accepting fully both it's strong points and shortcomings.

This macro language ("StarBasic") is much more primitive than VBA, in particular it is not that optimized. When you run the function against some 100 x/y-pairs, you can watch the sheet recalculate in seconds, due to the fact that the very same array of numbers is processed more than once for each of the 100 formulas.

100 pairs is quite a number. If n=100 is the number of points, theoretical size of the problem, in terms of the linear system involved (if you're solving by brute force) is 4x(n-1), so every time you use the formula, you tackle a 396x396 linear system (albeit indirectly...)

These links you provided, look interesting. A tool like that, would be a very useful framework in porting large projects to Calc.

...Copy my version of the code
...I did not start debugging that shit

I lost you here... Which code do you mean? You meant to post some code?

In the meantime, trying to see if things can be simplified somewhat, I found another snippet of code doing the same thing, from a different source. I tested it in Excel and numerically it behaves almost the same (same results on same data, numerical differences were insignificant). This one, does not have the extra, overshoot avoidance algorithm (I can live without it, straight vanilla is good enough for now...), so it is much shorter (only 100 lines). Here it is:

Code: Select all

Attribute VB_Name = "modSpline"
Option Base 1

'********************  Cubic_Spline by SRS1 Software ****************
'
'
' Version 1.01


Function cubic_spline(input_column As Range, _
                      output_column As Range, _
                      x As Range)

'Purpose:   Given a data set consisting of a list of x values
'           and y values, this function will smoothly interpolate
'           a resulting output (y) value from a given input (x) value


' This counts how many points are in "input" and "output" set of data
Dim input_count As Integer
Dim output_count As Integer

input_count = input_column.Rows.Count
output_count = output_column.Rows.Count

' Next check to be sure that "input" # points = "output" # points
If input_count <> output_count Then
    cubic_spline = "Something's messed up!  The number of indeces number of output_columnues don't match!"
    GoTo out
End If
 
ReDim xin(input_count) As Single
ReDim yin(input_count) As Single

Dim c As Integer

For c = 1 To input_count
xin(c) = input_column(c)
yin(c) = output_column(c)
Next c

'''''''''''''''''''''''''''''''''''''''
' values are populated
'''''''''''''''''''''''''''''''''''''''
Dim n As Integer 'n=input_count
Dim i, k As Integer 'these are loop counting integers
Dim p, qn, sig, un As Single
ReDim u(input_count - 1) As Single
ReDim yt(input_count) As Single 'these are the 2nd deriv values

n = input_count
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
    
qn = 0
un = 0

yt(n) = (un - qn * u(n - 1)) / (qn * yt(n - 1) + 1)

For k = n - 1 To 1 Step -1
    yt(k) = yt(k) * yt(k + 1) + u(k)
Next k


''''''''''''''''''''
'now eval spline at one point
'''''''''''''''''''''
Dim klo, khi As Integer
Dim h, b, a As Single

' first find correct interval
klo = 1
khi = n
Do
k = khi - klo
If xin(k) > x Then
khi = k
Else
klo = k
End If
k = khi - klo
Loop While k > 1
h = xin(khi) - xin(klo)
a = (xin(khi) - x) / h
b = (x - xin(klo)) / h
y = a * yin(klo) + b * yin(khi) + ((a ^ 3 - a) * yt(klo) + (b ^ 3 - b) * yt(khi)) * (h ^ 2) / 6


cubic_spline = y

out:
End Function

On this code, locating the essential stuff is easier. All the calculation tidbits should be Ok irrespective of VBA or other type of Basic, I guess... To put my finger on it, I would ask:

- this piece of code that evaluates the other dimension of the 1-dimensional arrays holding the input, how should it be given in OOo Basic?

Code: Select all

input_count = input_column.Rows.Count
output_count = output_column.Rows.Count

- if above line is OK, then all it takes is the correct function definition, because this

Code: Select all

Function cubic_spline(input_column As Range, _
                      output_column As Range, _
                      x As Range)

does not seem to work.

I tried an on line VBA to OpenOffice Basic converter found here and with input as above, the result is

Code: Select all

Function cubic_spline(input_column As Dim oSheet as Object
oSheet = ThisComponent.CurrentController.ActiveSheet
oSheet.getCellRangeByName($1), output_column As Dim oSheet as Object
oSheet = ThisComponent.CurrentController.ActiveSheet
oSheet.getCellRangeByName($1), x As Dim oSheet as Object
oSheet = ThisComponent.CurrentController.ActiveSheet
oSheet.getCellRangeByName($1))

which also does not seem to work...(the As Dim expression, is it valid?)

What should the correct definition be like?

And how could arguments be passed if getCellRangeByName is not used? In Excel the usual way to call the function - say for data in rows 1-10, A column for X, B column for Y and and target X value in E1, would be:

Code: Select all

cubic_spline(A1:A10;B1:B10;E1)

-

[Villeroy]

See attachment. It is the same code I posted Friday. Now I used it in a proper way. I do not claim any merits nor ownership nor do I even understand the implementation details. I just solved syntactical problems, adding one restriction: All vectors need to be vertical.

[johnroberts]

Code: Select all

 7	0.12183678	0,12183678
12	0.20792038	0,20792038
17	0.29236933	0,29236933
22	0.37460716	0,37460716
27	0.45399027	0,45399027
32	0.52991923	0,52991923
37	0.60181493	0,60181493
42	0.66913051	0,66913051
47	0.73135363	0,73135363
52	0.78801053	0,78801053
57	0.83867086	0,83867086
62	0.88294595	0,88294595
67	0.92051040	0,92051040
72	0.95103521	0,95103521
77	0.97444895	0,97444895
82	0.98997300	0,98997300

First column is angle in degrees and using as a test function sin(), namely sin(angle*pi()/180), the second column contains the interpolated values in Excel and the third column the interpolated values in OOo Calc (input data in both cases were 5deg.-85deg. in 5 deg. intervals).

I'd say we call it a success!!! Both thumbs up!!

... adding one restriction: All vectors need to be vertical

...a minor one; it is only logical for them to be vertical; w/o this restriction the user may have a bit more latitude in formating his worksheet, but this is trivial; important is that it should work.

I just solved syntactical problems

It would be beneficial for anybody reading this to pay attention to these syntactical differences, between the original code and the working code inside cubic_spline.ods

Many, many thanks

[johnroberts]

The differences between the two codes are in the following section only:

========================================================================
Attribute VB_Name = "Module1"
Option Explicit
Option Base 0


Function cubspline(Metode As Integer, xi As Double, xx As Object, yy As Object) As Double

Dim i As Integer
Dim yi As Double
Dim x() As Double
Dim y() As Double
Dim y2() As Double
Dim j As Integer

If Metode = 1 Then
'Numerical Recipes are 1 based
j = 0
Else
'Others are 0 based
j = -1
End If

For i = 1 To UBound(xx())
If yy(i,1) <> "" Then
j = j + 1
ReDim Preserve x(j)
ReDim Preserve y(j)
x(j) = CDbl(xx(i,1))
y(j) = CDbl(yy(i,1))
End If
Next i
===================================================================================

and using these cues, the modified, shorter version of the code should be:

Code: Select all

Function cubic_spline(input_column , output_column , x) As Double

'Purpose:   Given a data set consisting of a list of x values
'           and y values, this function will smoothly interpolate
'           a resulting output (y) value from a given input (x) value


' This counts how many points are in "input" and "output" set of data
Dim input_count As Integer
Dim output_count As Integer

input_count = Ubound(input_column())
output_count = Ubound(output_column())

' Next check to be sure that "input" # points = "output" # points
If input_count <> output_count Then
    cubic_spline = "Something's messed up!  The number of indeces number of output_columnues don't match!"
    GoTo out
End If
 
ReDim xin(input_count) As Double
ReDim yin(input_count) As Double

Dim c As Integer

For c = 1 To input_count
xin(c) = input_column(c,1)
yin(c) = output_column(c,1)
Next c

'''''''''''''''''''''''''''''''''''''''
' values are populated
'''''''''''''''''''''''''''''''''''''''
Dim n As Integer 'n=input_count
Dim i, k As Integer 'these are loop counting integers
Dim p, qn, sig, un As Double
ReDim u(input_count - 1) As Double
ReDim yt(input_count) As Double 'these are the 2nd deriv values

n = input_count
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
    
qn = 0
un = 0

yt(n) = (un - qn * u(n - 1)) / (qn * yt(n - 1) + 1)

For k = n - 1 To 1 Step -1
    yt(k) = yt(k) * yt(k + 1) + u(k)
Next k


''''''''''''''''''''
'now eval spline at one point
'''''''''''''''''''''
Dim klo, khi As Integer
Dim h, b, a As Double

' first find correct interval
klo = 1
khi = n
Do
k = khi - klo
If xin(k) > x Then
khi = k
Else
klo = k
End If
k = khi - klo
Loop While k > 1
h = xin(khi) - xin(klo)
a = (xin(khi) - x) / h
b = (x - xin(klo)) / h
y = a * yin(klo) + b * yin(khi) + ((a ^ 3 - a) * yt(klo) + (b ^ 3 - b) * yt(khi)) * (h ^ 2) / 6


cubic_spline = y

out:
End Function