1 文本格式
using System;
namespace Legalsoft.Truffer
{
/// <summary>
/// 三次样条插值
/// Cubic Spline Interpolation
/// Cubic spline interpolation object. Construct with x and y vectors, and
/// (optionally) values of the first derivative at the endpoints, then call
/// interp for interpolated values
/// </summary>
public class Spline_interp : Base_interp
{
private double[] y2 { get; set; }
public Spline_interp(double[] xv, double[] yv, double yp1 = 1.0e99, double ypn = 1.0e99) : base(xv, yv[0], 2)
{
this.y2 = new double[xv.Length];
sety2(xv, yv, yp1, ypn);
}
public Spline_interp(double[] xv, double yv, double yp1 = 1.0e99, double ypn = 1.0e99) : base(xv, yv, 2)
{
this.y2 = new double[xv.Length];
sety2(xv, y2, yp1, ypn);
}
/// <summary>
/// This routine stores an array y2[0..n - 1] with second derivatives of the
/// interpolating function at the tabulated points pointed to by xv, using
/// function values pointed to by yv.If yp1 and/or ypn are equal to 1.0E99 or
/// larger, the routine is signaled to set the corresponding boundary condition
/// for a natural spline, with zero second derivative on that boundary;
/// otherwise, they are the values of the first derivatives at the endpoints.
/// </summary>
/// <param name="xv"></param>
/// <param name="yv"></param>
/// <param name="yp1"></param>
/// <param name="ypn"></param>
public void sety2(double[] xv, double[] yv, double yp1, double ypn)
{
double[] u = new double[n - 1];
if (yp1 > 0.99e99)
{
y2[0] = u[0] = 0.0;
}
else
{
y2[0] = -0.5;
u[0] = (3.0 / (xv[1] - xv[0])) * ((yv[1] - yv[0]) / (xv[1] - xv[0]) - yp1);
}
for (int i = 1; i < n - 1; i++)
{
double sig = (xv[i] - xv[i - 1]) / (xv[i + 1] - xv[i - 1]);
double p = sig * y2[i - 1] + 2.0;
y2[i] = (sig - 1.0) / p;
u[i] = (yv[i + 1] - yv[i]) / (xv[i + 1] - xv[i]) - (yv[i] - yv[i - 1]) / (xv[i] - xv[i - 1]);
u[i] = (6.0 * u[i] / (xv[i + 1] - xv[i - 1]) - sig * u[i - 1]) / p;
}
double qn;
double un;
if (ypn > 0.99e99)
{
qn = un = 0.0;
}
else
{
qn = 0.5;
un = (3.0 / (xv[n - 1] - xv[n - 2])) * (ypn - (yv[n - 1] - yv[n - 2]) / (xv[n - 1] - xv[n - 2]));
}
y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.0);
for (int k = n - 2; k >= 0; k--)
{
y2[k] = y2[k] * y2[k + 1] + u[k];
}
}
/// <summary>
/// Given a value x, and using pointers to data xx and yy, this routine returns
/// an interpolated value y, and stores an error estimate dy. The returned
/// value is obtained by mm-point polynomial interpolation on the subrange
/// xx[jl..jl + mm - 1].
/// </summary>
/// <param name="jl"></param>
/// <param name="x"></param>
/// <returns></returns>
/// <exception cref="Exception"></exception>
public override double rawinterp(int jl, double x)
{
int klo = jl;
int khi = jl + 1;
double h = xx[khi] - xx[klo];
//if (h == 0.0)
if (Math.Abs(h) <= float.Epsilon)
{
throw new Exception("Bad input to routine splint");
}
double a = (xx[khi] - x) / h;
double b = (x - xx[klo]) / h;
double y = a * yy[klo] + b * yy[khi] + ((a * a * a - a) * y2[klo] + (b * b * b - b) * y2[khi]) * (h * h) / 6.0;
return y;
}
}
}
2 代码格式
using System;
namespace Legalsoft.Truffer
{
/// <summary>
/// 三次样条插值
/// Cubic Spline Interpolation
/// Cubic spline interpolation object. Construct with x and y vectors, and
/// (optionally) values of the first derivative at the endpoints, then call
/// interp for interpolated values
/// </summary>
public class Spline_interp : Base_interp
{
private double[] y2 { get; set; }
public Spline_interp(double[] xv, double[] yv, double yp1 = 1.0e99, double ypn = 1.0e99) : base(xv, yv[0], 2)
{
this.y2 = new double[xv.Length];
sety2(xv, yv, yp1, ypn);
}
public Spline_interp(double[] xv, double yv, double yp1 = 1.0e99, double ypn = 1.0e99) : base(xv, yv, 2)
{
this.y2 = new double[xv.Length];
sety2(xv, y2, yp1, ypn);
}
/// <summary>
/// This routine stores an array y2[0..n - 1] with second derivatives of the
/// interpolating function at the tabulated points pointed to by xv, using
/// function values pointed to by yv.If yp1 and/or ypn are equal to 1.0E99 or
/// larger, the routine is signaled to set the corresponding boundary condition
/// for a natural spline, with zero second derivative on that boundary;
/// otherwise, they are the values of the first derivatives at the endpoints.
/// </summary>
/// <param name="xv"></param>
/// <param name="yv"></param>
/// <param name="yp1"></param>
/// <param name="ypn"></param>
public void sety2(double[] xv, double[] yv, double yp1, double ypn)
{
double[] u = new double[n - 1];
if (yp1 > 0.99e99)
{
y2[0] = u[0] = 0.0;
}
else
{
y2[0] = -0.5;
u[0] = (3.0 / (xv[1] - xv[0])) * ((yv[1] - yv[0]) / (xv[1] - xv[0]) - yp1);
}
for (int i = 1; i < n - 1; i++)
{
double sig = (xv[i] - xv[i - 1]) / (xv[i + 1] - xv[i - 1]);
double p = sig * y2[i - 1] + 2.0;
y2[i] = (sig - 1.0) / p;
u[i] = (yv[i + 1] - yv[i]) / (xv[i + 1] - xv[i]) - (yv[i] - yv[i - 1]) / (xv[i] - xv[i - 1]);
u[i] = (6.0 * u[i] / (xv[i + 1] - xv[i - 1]) - sig * u[i - 1]) / p;
}
double qn;
double un;
if (ypn > 0.99e99)
{
qn = un = 0.0;
}
else
{
qn = 0.5;
un = (3.0 / (xv[n - 1] - xv[n - 2])) * (ypn - (yv[n - 1] - yv[n - 2]) / (xv[n - 1] - xv[n - 2]));
}
y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.0);
for (int k = n - 2; k >= 0; k--)
{
y2[k] = y2[k] * y2[k + 1] + u[k];
}
}
/// <summary>
/// Given a value x, and using pointers to data xx and yy, this routine returns
/// an interpolated value y, and stores an error estimate dy. The returned
/// value is obtained by mm-point polynomial interpolation on the subrange
/// xx[jl..jl + mm - 1].
/// </summary>
/// <param name="jl"></param>
/// <param name="x"></param>
/// <returns></returns>
/// <exception cref="Exception"></exception>
public override double rawinterp(int jl, double x)
{
int klo = jl;
int khi = jl + 1;
double h = xx[khi] - xx[klo];
//if (h == 0.0)
if (Math.Abs(h) <= float.Epsilon)
{
throw new Exception("Bad input to routine splint");
}
double a = (xx[khi] - x) / h;
double b = (x - xx[klo]) / h;
double y = a * yy[klo] + b * yy[khi] + ((a * a * a - a) * y2[klo] + (b * b * b - b) * y2[khi]) * (h * h) / 6.0;
return y;
}
}
}
