Wahrscheinlich ändert dein Makro irgendetwas ...
18.09.2017 16:30:56
Bernd
Option Explicit
Option Base 0
Public Function Cubspline(Method As Integer, xi As Double, _
xx As Object, yy As Object) As Double
' valid input for parameter "Method":
' 1 = cubic spline / 3 = own spline
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 Method = 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 Method = 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 Method = 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)
'Part one of Method 1
'Given arrays x(1:n) and y(1:n) containing a tabulated function,
'i.e., y i = f(xi), with x1
'and ypn for the first derivative of the interpolating 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)
'Part two of Method 1
'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
'| Method 3
'| 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.
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(Nmax) Then
'X outisde range. Linear interpolate
'Below min or max?
If X
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