Macro function PINTERP interpolates X-Y data tables arranged in either column or row format. PINTERP is especially suitable for very large spreadsheet data tables, because it uses an efficient bi-section search to bracket the value to be interpolated, and uses the average of the two possible quadratic interpolation polynominals over the nearest 4 data points in the table instead of a cubic spline. This approach appears as accurate as cubic interpolation for equally spaced points, and better than quadratic interpolation for non equally spaced points, but does not require solution of the cubic spline matrix for every function call, which would quickly become a problem if pasted to a huge table sending different data sets on every call.
OpenOffice Calc Version (Not Excel Compatible - See below for Excel VBA Version):
Code: Select all
Function PINTERP(xtab, ytab, x, Optional Linear)
'
' PINTERP Version 1.0.0
' PINTERP Version 1.0.2 Minor comment changes
'
' OpenOffice.org Calc BASIC Spreadsheet Interpolation Macro
'
' Copyright (C) 2010, 2018 E Geyer
'
' This PINTERP function is free software: you can redistribute it
' and/or modify it under the terms of the GNU Lesser General Public
' License as published by the Free Software Foundation, either
' version 3 of the License, or (at your option) any later version.
'
' This program is distributed in the hope that it will be useful,
' but WITHOUT ANY WARRANTY; without even the implied warranty of
' MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
' GNU Lesser General Public License for more details.
'
' You should have received a copy of the GNU Lesser General Public License
' along with this program. If not, see <http://www.gnu.org/licenses/>.
'
' E Geyer
' geyerej@gmail.com
'
' Input:
' xtab: array of n x values
' ytab: array of n y values
' x: interpolant
' Linear: optional boolean; forces linear interpolation if >0
'
' Output:
' y: interpolated y at x
'
' Interpolation method:
' n=2: linear
' n=3: quadratic
' n>3: average of two possible quadratics through nearest 4 local points
' Linear: linear if optional Linear argument is true
'
' Cubic spline not used due to spreadsheet macro run time issues
' with large data sets.
'
Dim coltables, descending As Boolean
Dim nx, lo, hi, mid, nquad As Long
Dim quad, b1, b2 As Double
Dim delhilo, dellolo1, delhilo1, delhi1hi,delhi1lo As Double
'
' Default is more accurate quadratic interpolation
' Linear interpolation is for exceptional cases only
'
If IsMissing(Linear) Then Linear=0
'
' Determine if input tables are in columns or rows
'
coltables=(UBound(xtab,1)> UBound(xtab,2))
'
If (coltables) Then
'
' Data tables are arranged in columns
'
' Check Input Data
' Avoid loop counting data as it is inefficient for large tables
'
nx=UBound(xtab,1)
If (nx<2) Then
PINTERP="Error: nX<2"
Exit Function
End If
'
If (nx<>UBound(ytab,1)) Then
PINTERP="Error: nX<>nY"
Exit Function
End If
'
' Allow limited extrapolation
'
If (((xtab(1,1)-x(1,1))/(xtab(2,1)-xtab(1,1))) >0.20 or _
((x(1,1)-xtab(nx,1))/(xtab(nx,1)-xtab(nx-1,1)))>0.20) Then
PINTERP="Error: >20% extrapolation"
Exit Function
End If
'
' Locate interval containing x with bisection
' Run time is very important for huge spreadsheet tables:
' a) Do not test for table order(via descending Xor...) inside Loop
' b) Do not check equality inside Loop
' Run time is not important at all for small tables:
' a) Accept 1 extra iteration to avoid equality check in loop
'
descending=(xtab(1,1)>xtab(nx,1))
lo=1
hi=nX
If (descending) Then
Do While (lo<hi)
'do not use average to avoid addition overflow; note backslash integer divide
mid=lo+(hi-lo)\2
If (x(1,1)<xtab(mid,1)) Then
lo=mid+1
Else
hi=mid
End If
Loop
Else
Do While (lo<hi)
'note no averaging to prevent addition overflow; and backslash integer divide
mid=lo+(hi-lo)\2
If (x(1,1)>xtab(mid,1)) Then
lo=mid+1
Else
hi=mid
End If
Loop
End If
If (descending Xor (x(1,1)>xtab(1,1))) Then
lo=hi-1
Else
hi=lo+1
End If
'
' Check for divide by zero and save result for re-use.
'
delhilo=(xtab(hi,1)-xtab(lo,1))
If (delhilo=0) Then
PINTERP="Error: Equal X values"
Exit Function
End If
'
' Linear for 2 points
'
If nx<3 or Linear Then
PINTERP=ytab(lo,1)+(x(1,1)-xtab(lo,1))/delhilo*(ytab(hi,1)-ytab(lo,1))
Exit Function
End IF
'
' Compute possible Newton quadratic polynomials
'
nquad=0
quad=0
If lo>1 Then
'
dellolo1=(xtab(lo,1)-xtab(lo-1,1))
If (dellolo1=0) Then
PINTERP="Error: Equal X values"
Exit Function
End If
'
delhilo1=(xtab(hi,1)-xtab(lo-1,1))
If (delhilo1=0) Then
PINTERP="Error: Equal X values"
Exit Function
End If
'
b1=(ytab(lo,1)-ytab(lo-1,1))/dellolo1
b2=((ytab(hi,1)-ytab(lo,1))/delhilo-b1)/delhilo1
quad=ytab(lo-1,1)+b1*(x(1,1)-xtab(lo-1,1))+b2*(x(1,1)-xtab(lo-1,1))*(x(1,1)-xtab(lo,1))
nquad=nquad+1
End If
If hi<nx Then
'
delhi1hi=(xtab(hi+1,1)-xtab(hi,1))
If (delhi1hi=0) Then
PINTERP="Error: Equal X values"
Exit Function
End If
'
delhi1lo=(xtab(hi+1,1)-xtab(lo,1))
If (delhi1lo=0) Then
PINTERP="Error: Equal X values"
Exit Function
End If
'
b1=(ytab(hi,1)-ytab(lo,1))/delhilo
b2=((ytab(hi+1,1)-ytab(hi,1))/delhi1hi-b1)/delhi1lo
quad=quad+ytab(lo,1)+b1*(x(1,1)-xtab(lo,1))+b2*(x(1,1)-xtab(lo,1))*(x(1,1)-xtab(hi,1))
nquad=nquad+1
End If
'
PINTERP=quad/nquad
Else
'
' Data tables are arranged in rows
'
' Check Input Data
' Avoid loop counting data as it is inefficient for large tables
'
nx=UBound(xtab,2)
If (nx<2) Then
PINTERP="Error: nX<2"
Exit Function
End If
'
If (nx<>UBound(ytab,2)) Then
PINTERP="Error: nX<>nY"
Exit Function
End If
'
' Allow limited extrapolation
'
If (((xtab(1,1)-x(1,1))/(xtab(1,2)-xtab(1,1))) >0.20 or _
((x(1,1)-xtab(1,nx))/(xtab(1,nx)-xtab(1,nx-1)))>0.20) Then
PINTERP="Error: >20% extrapolation"
Exit Function
End If
'
' Locate interval containing x with bisection
' Run time is very important for huge spreadsheet tables:
' a) Do not test for table order(via descending Xor...) inside Loop
' b) Do not check equality inside Loop
' Run time is not important at all for small tables:
' a) Accept 1 extra iteration to avoid equality check in loop
'
descending=(xtab(1,1)>xtab(1,nx))
lo=1
hi=nX
If (descending) Then
Do While (lo<hi)
'do not use average to avoid addition overflow; note backslash integer divide
mid=lo+(hi-lo)\2
If (x(1,1)<xtab(1,mid)) Then
lo=mid+1
Else
hi=mid
End If
Loop
Else
Do While (lo<hi)
'note no averaging to prevent addition overflow; and backslash integer divide
mid=lo+(hi-lo)\2
If (x(1,1)>xtab(1,mid)) Then
lo=mid+1
Else
hi=mid
End If
Loop
End If
If (descending Xor (x(1,1)>xtab(1,1))) Then
lo=hi-1
Else
hi=lo+1
End If
'
' Check for divide by zero and save result for re-use.
'
delhilo=(xtab(1,hi)-xtab(1,lo))
If (delhilo=0) Then
PINTERP="Error: Equal X values"
Exit Function
End If
'
' Linear for 2 points
'
If nx<3 or Linear Then
PINTERP=ytab(1,lo)+(x(1,1)-xtab(1,lo))/delhilo*(ytab(1,hi)-ytab(1,lo))
Exit Function
End IF
'
' Compute possible Newton quadratic polynomials
'
nquad=0
quad=0
If lo>1 Then
'
dellolo1=(xtab(1,lo)-xtab(1,lo-1))
If (dellolo1=0) Then
PINTERP="Error: Equal X values"
Exit Function
End If
'
delhilo1=(xtab(1,hi)-xtab(1,lo-1))
If (delhilo1=0) Then
PINTERP="Error: Equal X values"
Exit Function
End If
'
b1=(ytab(1,lo)-ytab(1,lo-1))/dellolo1
b2=((ytab(1,hi)-ytab(1,lo))/delhilo-b1)/delhilo1
quad=ytab(1,lo-1)+b1*(x(1,1)-xtab(1,lo-1))+b2*(x(1,1)-xtab(1,lo-1))*(x(1,1)-xtab(1,lo))
nquad=nquad+1
End If
If hi<nx Then
'
delhi1hi=(xtab(1,hi+1)-xtab(1,hi))
If (delhi1hi=0) Then
PINTERP="Error: Equal X values"
Exit Function
End If
'
delhi1lo=(xtab(1,hi+1)-xtab(1,lo))
If (delhi1lo=0) Then
PINTERP="Error: Equal X values"
Exit Function
End If
'
b1=(ytab(1,hi)-ytab(1,lo))/delhilo
b2=((ytab(1,hi+1)-ytab(1,hi))/delhi1hi-b1)/delhi1lo
quad=quad+ytab(1,lo)+b1*(x(1,1)-xtab(1,lo))+b2*(x(1,1)-xtab(1,lo))*(x(1,1)-xtab(1,hi))
nquad=nquad+1
End If
'
PINTERP=quad/nquad
End If
End Function
Excel VBA Version:
Code: Select all
Function PINTERP(xtab As Range, ytab As Range, x, Optional Linear)
' PINTERP Version 1.0.0: Original OpenOffice Macro
' PINTERP Version 1.0.1: Excel VBA conversion
' PINTERP Version 1.0.2: Minor comment changes
' Excel VBA Spreadsheet Interpolation Macro
'
' Copyright (C) 2010, 2012, 2018 E Geyer
'
' This PINTERP function is free software: you can redistribute it
' and/or modify it under the terms of the GNU Lesser General Public
' License as published by the Free Software Foundation, either
' version 3 of the License, or (at your option) any later version.
'
' This program is distributed in the hope that it will be useful,
' but WITHOUT ANY WARRANTY; without even the implied warranty of
' MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
' GNU Lesser General Public License for more details.
'
' You should have received a copy of the GNU Lesser General Public License
' along with this program. If not, see <http://www.gnu.org/licenses/>.
'
' E Geyer
' geyerej@gmail.com
'
' Input:
' xtab: array of n x values
' ytab: array of n y values
' x: interpolant
' Linear: optional boolean; forces linear interpolation if >0
'
' Output:
' y: interpolated y at x
'
' Interpolation method:
' n=2: linear
' n=3: quadratic
' n>3: average of two possible quadratics through nearest 4 local points
' Linear: linear if optional Linear argument is true
'
' Cubic spline not used due to spreadsheet macro run time issues
' with large data sets.
'
Dim descending As Boolean
Dim nx, lo, hi, mid, nquad As Long
Dim quad, b1, b2 As Double
Dim delhilo, dellolo1, delhilo1, delhi1hi, delhi1lo As Double
'
' Default is more accurate quadratic interpolation
' Linear interpolation is for exceptional cases only
'
If IsMissing(Linear) Then Linear = 0
'
' Check Input Data
' Avoid loop counting data as it is inefficient for large tables
'
nx = xtab.Count
If (nx < 2) Then
PINTERP = "Error: nX<2"
Exit Function
End If
'
If (nx <> ytab.Count) Then
PINTERP = "Error: nX<>nY"
Exit Function
End If
'
' Allow limited extrapolation
'
If (((xtab(1) - x) / (xtab(2) - xtab(1))) > 0.2 Or _
((x - xtab(nx)) / (xtab(nx) - xtab(nx - 1))) > 0.2) Then
PINTERP = "Error: >20% extrapolation"
Exit Function
End If
'
' Locate interval containing x with bisection
' Run time is very important for huge spreadsheet tables:
' a) Do not test for table order(via descending Xor...) inside Loop
' b) Do not check equality inside Loop
' Run time is not important at all for small tables:
' a) Accept 1 extra iteration to avoid equality check in loop
'
descending = (xtab(1) > xtab(nx))
lo = 1
hi = nx
If (descending) Then
Do While (lo < hi)
'do not use average to avoid addition overflow; note backslash integer divide
mid = lo + (hi - lo) \ 2
If (x < xtab(mid)) Then
lo = mid + 1
Else
hi = mid
End If
Loop
Else
Do While (lo < hi)
'note no averaging to prevent addition overflow; and backslash integer divide
mid = lo + (hi - lo) \ 2
If (x > xtab(mid)) Then
lo = mid + 1
Else
hi = mid
End If
Loop
End If
If (descending Xor (x > xtab(1))) Then
lo = hi - 1
Else
hi = lo + 1
End If
'
' Check for divide by zero and save result for re-use.
'
delhilo = (xtab(hi) - xtab(lo))
If (delhilo = 0) Then
PINTERP = "Error: Equal X values"
Exit Function
End If
'
' Linear for 2 points
'
If nx < 3 Or Linear Then
PINTERP = ytab(lo) + (x - xtab(lo)) / delhilo * (ytab(hi) - ytab(lo))
Exit Function
End If
'
' Compute possible Newton quadratic polynomials
'
nquad = 0
quad = 0
If lo > 1 Then
'
dellolo1 = (xtab(lo) - xtab(lo - 1))
If (dellolo1 = 0) Then
PINTERP = "Error: Equal X values"
Exit Function
End If
'
delhilo1 = (xtab(hi) - xtab(lo - 1))
If (delhilo1 = 0) Then
PINTERP = "Error: Equal X values"
Exit Function
End If
'
b1 = (ytab(lo) - ytab(lo - 1)) / dellolo1
b2 = ((ytab(hi) - ytab(lo)) / delhilo - b1) / delhilo1
quad = ytab(lo - 1) + b1 * (x - xtab(lo - 1)) + b2 * (x - xtab(lo - 1)) * (x - xtab(lo))
nquad = nquad + 1
End If
If hi < nx Then
'
delhi1hi = (xtab(hi + 1) - xtab(hi))
If (delhi1hi = 0) Then
PINTERP = "Error: Equal X values"
Exit Function
End If
'
delhi1lo = (xtab(hi + 1) - xtab(lo))
If (delhi1lo = 0) Then
PINTERP = "Error: Equal X values"
Exit Function
End If
'
b1 = (ytab(hi) - ytab(lo)) / delhilo
b2 = ((ytab(hi + 1) - ytab(hi)) / delhi1hi - b1) / delhi1lo
quad = quad + ytab(lo) + b1 * (x - xtab(lo)) + b2 * (x - xtab(lo)) * (x - xtab(hi))
nquad = nquad + 1
End If
'
PINTERP = quad / nquad
End Function
