{1D and 2D spline routines
author: Nikolai Shokhirev http://www.shokhirev.com/nikolai.html
http://www.chem.arizona.edu/~shokhirn/nikolai.html
created: 2002.02.02
last modified: 2004.12.27
ŠNikolai V. Shokhirev, 2002-2004 }
unit uSpline;
interface
uses
uMatTypes, uDynArrays;
{
Pascal version by Nikolai V. Shokhirev.
Note: Added the estimation of derivatives (NVS)
Ref: W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P.
Numerical Recipes in Fortran http://www.nr.com/
Input
x[Lo..Hi],
y[Lo..Hi] - y[i] = f(x[i])
yp1 - left boundary condition
yp1 > 0.99e30 - natural spline: f"(x[Lo]) = 0
yp1 < -0.99e30 - 3 or 4-point estimation of f'(x[Lo]) - default
otherwise f'[x[Lo]) = fp1
ypn - right boundary condition
yp1 > 0.99e30 - natural spline: f"(x[Hi]) = 0
yp1 < -0.99e30 - 3 or 4-point estimation of f'(x[Hi]) - default
otherwise f'[x[Hi]) = fpn
Output
y2[Lo..Hi] - y2[i] = f"(x[i]) }
procedure spline(const x, y, y2: IFArr1D; yp1: TFloat = -1.0e30;
ypn: TFloat = -1.0e30);
{
Pascal version by Nikolai V. Shokhirev.
Ref: W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P.
Numerical Recipes in Fortran http://www.nr.com/
Input
xa[Lo..Hi],
ya[Lo..Hi] - ya[i] = f(xa[i])
y2a[Lo..Hi] - y2a[i] = f"(xa[i])
x - 1D
Output
y = f(x) }
procedure splint(const xa, ya, y2a: IFArr1D; x: TFloat; var y: TFloat);
{
Pascal version by Nikolai V. Shokhirev.
Note: Added the estimation of derivatives (accuracy increase)
Ref: W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P.
Numerical Recipes in Fortran http://www.nr.com/
Purpose
Precomputing of the 2nd-derivative table (over x2).
Note
This routine is called once before calling splin2.
Input
x1a[Lo1..Hi1], - not used, included for consistency with splin2
x2a[Lo2..Hi2],
ya[Lo1..Hi1,Lo2..Hi2] - y[i, j] = f(x1a[i], x2a[j])
Output
y2a[Lo1..Hi1, Lo2..Hi2] - y2a[i,j] = f"(x1a[i], x2a[j]) over x2 }
procedure splie2(const x1a, x2a: IFArr1D; const ya, y2a: IFArr2D);
{
Pascal version by Nikolai V. Shokhirev.
Note: Added the estimation of derivatives (accuracy increase)
Ref: W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P.
Numerical Recipes in Fortran http://www.nr.com/
Purpose
Calculation of the bicubic spline interpolation at (x1, x2)
Input
x1a[Lo1..Hi1],
x2a[1..n],
ya[Lo1..Hi1,Lo2..Hi2] - y[i, j] = f(x1a[i], x2a[j])
y2a[Lo1..Hi1, Lo2..Hi2] - y2a[i,j] = f"(x1a[i], x2a[j]) over x2
x1, x2 - 2D coordinates
Output
y = f(x1, x2) }
procedure splin2(const x1a, x2a: IFArr1D; const ya, y2a: IFArr2D;
x1, x2: TFloat; var y: TFloat);
{ Xarr[Index2] < x <= Xarr[Index2+1] }
function Index2(x: TFloat; const Xa: IFArr1D): TInt;
{ Xarr[Index3] is nearest to x }
function Index0(x: TFloat; const Xa: IFArr1D): TInt;
{ Pascal version by Nikolai V. Shokhirev.
Ref: G.E.Forsythe, M.A.Malcolm, C.B.Moler.
Computer methods for mathematical computations. Prentice-Hall, 1977.
PURPOSE:
Compute the coefficients b,c,d for a cubic interpolating spline so that the
interpolated value is given by
s[x] = y[k] + b[k]*(x-x[k]) + c[k]*(x-x[k])**2 + d[k]*(x-x[k])**3
when x[k] <= x <= x[k+1]
The end conditions match the third derivatives of the interpolated curve to
the third derivatives of the unique polynomials thru the first four and
last four points.
Use Seval or Seval3 to evaluate the spline.
INPUT:
x - abscissas of knots
y - ordinates of knots
OUTPUT:
b - linear coeff
c - quadratic coeff.
d - cubic coeff.
All arrays have the same limits }
procedure FMMSpline(const x, y, b, c, d: IFArr1D);
{ Pascal version by Nikolai V. Shokhirev.
Ref: G.E.Forsythe, M.A.Malcolm, C.B.Moler.
Computer methods for mathematical computations. Prentice-Hall, 1977.
PURPOSE:
Construct the natural spline thru a set of points
NOTES:
A natural spline has zero second derivative at both endpoints.
INPUT:
x, y - coordinates of knots
OUTPUT:
b, c, d - linear, quadratic and cubic coefficients }
procedure NaturalSpline(const x, y, b, c, d: IFArr1D);
{ Pascal version by Nikolai V. Shokhirev.
Ref: G.E.Forsythe, M.A.Malcolm, C.B.Moler.
Computer methods for mathematical computations. Prentice-Hall, 1977.
PURPOSE:
Evaluate the cubic spline function
Seval=y[i]+b[i]!(u-x[i])+c[i]*(u-x[i])**2+d[i]*(u-x[i])**3
where x[i] <= u < x[i+1]
NOTES:
if u<x[Lo], i = Lo is used;if u>x[Hi], i = Hi is used
INPUT:
u - abscissa at which the spline is to be evaluated
x - abscissas of knots
y - ordinates of knots
b, c, d - linear, quadratic, cubic coeff
IN-OUT
i - last used index; should be set to x.Lo1 before the 1st call }
function Seval(const u: TFloat; var i: TInt; const x, y, b, c, d: IFArr1D): TFloat;
{ Pascal version by Nikolai V. Shokhirev.
Ref: G.E.Forsythe, M.A.Malcolm, C.B.Moler.
Computer methods for mathematical computations. Prentice-Hall, 1977.
PURPOSE:
Evaluate the cubic spline function
Seval=y[i]+b[i]!(u-x[i])+c[i]*(u-x[i])**2+d[i]*(u-x[i])**3
where x[i] <= u < x[i+1]
NOTES:
if u<x[Lo], i=1 is used;if u>x[Hi], i=Dim is used
u - abscissa at which the spline is to be evaluated
x - abscissas of knots
y - ordinates of knots
b, c, d - linear, quadratic, cubic coeff
f, fp, fpp, fppp - function, 1st, 2nd, 3rd derivatives
IN-OUT
i - last used index; should be set to x.Lo1 before the 1st call }
procedure Seval3(const u: TFloat; var i: TInt; const x, y, b, c, d: IFArr1D;
var f, fp, fpp, fppp: TFloat);
{
PURPOSE:
Linear interpolation of 1D data
NOTES:
Use LinSpline1D for fast lookup of dense data
Use FMMSpline / Seval or spline / splint for higher accuracy
INPUT:
x - 1D point ( Dat1D.x1 )
RESULT:
f(x) - interpolated value ( Dat1D.x2 ); }
function LinSpline1D(x: TFloat; const Dat1D: IFData1D): TFloat;
{
PURPOSE:
Linear interpolation of 2D data
NOTES:
Use LinSpline2D for fast lookup of dense data
Use splie2 / spline for higher accuracy
INPUT:
(t1, t2) - 2D point ( Dat2D.x1, Dat2D.x2 )
RESULT:
f(t1, t2) - interpolated value ( Dat2D.z ); }
function LinSpline2D(t1, t2: TFloat; const Dat2D: IFData2D): TFloat;
implementation
uses
Math;
end.
Back
Generated by Lore's Source to HTML Converter(http://www.newty.de/lsc/index.html)