Das Programm stammt aus dem Quellcode des Freeware Programms XNUMBERS 5.6 - Multi Precision Floating Point Computing and Numerical Methods for EXCEL.
Ein Muss für jeden Excel-Anwender. Einfach suchen unter: Digilander
====================== ANFANG ======================
- 'Interpolates one point from sorted X,Y data pairs using Cubic Spline Interpolation
- Function Interpol_cspline(Xin As Range, Yin As Range, Xtarget As Double)
- Dim n%, ny% 'Xin and Yin counts
- Dim i% 'loop counter
- Dim ihi%, m%
- Dim a#, b#, p#, qn#, sig#, un#, h0#, h1#, h2#, Yout#
- ny = Yin.Cells.Count
- n = Xin.Cells.Count
- ' Next check to be sure that there are the came counts of Xin and Yin
- If n <> ny Then
- Interpol_cspline = CVErr(xlErrRef) 'Uneven counts of Xin and Yin"
- GoTo err_ret
- End If
- ReDim x(1 To n) As Double
- ReDim Y(1 To n) As Double
- 'ReDim u(1 To n - 1) As Single
- ReDim U(1 To n - 1) As Double 'mod. LV 31-8-02
- ReDim ypp(1 To n) As Double 'these are the 2nd derivative values
- 'populate the input arrays
- For i = 1 To n
- x(i) = Xin(i)
- If i > 1 Then
- If x(i) <= x(i - 1) Then 'Check if X values are not repeating and are sorted in ascending order
- Interpol_cspline = CVErr(xlErrValue) 'Error: Not sorted or repeating
- GoTo err_ret
- End If
- End If
- Y(i) = Yin(i)
- Next i
- ypp(1) = 0 'First knot boundary condition
- U(1) = 0 'First knot boundary condition
- For i = 2 To n - 1
- h0 = x(i + 0) - x(i - 1)
- h1 = x(i + 1) - x(i - 0)
- h2 = x(i + 1) - x(i - 1)
- sig = h0 / h2
- p = sig * ypp(i - 1) + 2#
- ypp(i) = (sig - 1) / p
- U(i) = (Y(i + 1) - Y(i)) / h1 - (Y(i) - Y(i - 1)) / h0
- U(i) = (6# * U(i) / h2 - sig * U(i - 1)) / p
- Next i
- qn = 0
- un = 0
- ypp(n) = (un - qn * U(n - 1)) / (qn * ypp(n - 1) + 1) 'Last knot boundary condition
- For i = n - 1 To 1 Step -1
- ypp(i) = ypp(i) * ypp(i + 1) + U(i) 'Backfill the 2nd derivatives
- Next i
- ''''''''''''''''''''''''''''''''''''''''
- 'Spline evaluation at target X point
- '''''''''''''''''''''''''''''''''''''''''
- 'Find correct interval using halving binary search
- i = 1
- ihi = n
- Do While (i <>
- m = (i + ihi) \ 2 'Calculate the midpoint
- If (Xtarget < ihi =" m" i =" m">
- Loop
- ' i = beginning of the correct interval
- h0 = x(i + 1) - x(i) 'Calc the width of the X interval
- a = (x(i + 1) - Xtarget) / h0
- b = (Xtarget - x(i)) / h0
- Yout = a * Y(i) + b * Y(i + 1) + ((a ^ 3 - a) * ypp(i) + (b ^ 3 - b) * ypp(i + 1)) * (h0 ^ 2) / 6#
- Interpol_cspline = Yout
- err_ret:
- End Function
Gruß Chiochip
Keine Kommentare:
Kommentar veröffentlichen