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 Expand viewCollapse view
`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 Expand viewCollapse view
`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