spline.pas

一个详细功能齐全的档案管理信息系统 跟大家分享一下。这是老师给的光盘里带的

unit Spline;

interface

uses
  Windows, Messages, SysUtils, Classes;

type
  TSpline = class(TComponent)
  private
    { Private declarations }
    XA, YA : array of Double;   //输入点坐标
    Y2A    : array of Double;   //二阶导数
    num    : integer;           //输入点数

    procedure Switch(i,j:Integer);
  protected
    { Protected declarations }
  public
    { Public declarations }
    procedure Spline(YP1,YPN : Double);     //
    procedure Splint(X:Double;var Y:Double);//
    procedure Sort(dir:Integer); //1:从小道大  -1：从大道小
    procedure Clear;
    procedure AddXY(x,y:Double);
    function  AssignParam(x,y: array of Double; N: Integer): Integer;
    procedure SplineCurve(x : array of Double;var Y : array of Double);
    procedure SplinePoint(x : Double;var y : Double);
  published
    { Published declarations }
  end;

procedure Register;

implementation

procedure Register;
begin
  RegisterComponents('System', [TSpline]);
end;

{ TSpline }

procedure TSpline.AddXY(x, y: Double);   //增加点
begin
  Inc(num);
  SetLength(XA,num);
  SetLength(YA,num);
  SetLength(Y2A,num);
  XA[High(XA)] := x;
  YA[High(YA)] := y;
end;
//设点
function TSpline.AssignParam(x,y: array of double; N: Integer): Integer;
var
  i : Integer;
  yp1,ypn : Double;
begin
  num := N;
  SetLength(XA,num);
  SetLength(YA,num);
  SetLength(Y2A,num);

  for i:=Low(XA) to High(XA) do
  begin
    XA[i] := x[i];
    YA[i] := y[i];
    Y2A[i] := 0;
  end;

  Sort(1);  //排序（从小到大）

  yp1 := 0;
  ypn := 0;

  Spline(yp1,ypn);
  result := 1;
end;
//清除点
procedure TSpline.Clear;
begin
  num := 0;
  SetLength(XA,num);
  SetLength(YA,num);
  SetLength(Y2A,num);
end;
//排序
procedure TSpline.Sort(dir : Integer);
var
  i,j:Integer;
begin
  for i:=Low(XA) to High(XA)-1 do
  begin
    for j:=i+1 to High(XA) do
    begin
      if dir = 1 then     //从小到大
      begin
        if XA[j] < XA[i] then
          Switch(i,j);
      end
      else if dir = -1 then  //从大到小
      begin
        if XA[j] > XA[i] then
          Switch(i,j);
      end;
    end;
  end;
end;

procedure TSpline.Spline(YP1,YPN : Double);
var
    U:array of Double;
    AAA,SIG,BBB,CCC,P,QN,UN:Double;
    I,K:integer;
begin
  try
    SetLength(U,num);
    If YP1 > 9.9E+29 Then
    begin
        Y2A[0]:=0;
        U[0]:=0;
    end
    Else
    begin
        Y2A[0]:=-0.5;
        AAA:=(YA[1] - YA[0]) / (XA[1] - XA[0]);
        U[0]:=(3 / (XA[1] - XA[0])) * (AAA - YP1);
    end;
    For I:=1 To num - 2 do
    begin
        SIG:=(XA[I] - XA[I - 1]) / (XA[I + 1] - XA[I - 1]);
        P:=SIG * Y2A[I - 1] + 2;
        Y2A[I]:=(SIG - 1) / P;
        AAA:=(YA[I + 1] - YA[I]) / (XA[I + 1] - XA[I]);
        BBB:=(YA[I] - YA[I - 1]) / (XA[I] - XA[I - 1]);
        CCC:=XA[I + 1] - XA[I - 1];
        U[I]:=(6 * (AAA - BBB) / CCC - SIG * U[I - 1]) / P;
    end;
    If YPN > 9.9E+29 Then
    begin
        QN:=0;
        UN:=0;
    end
    Else
    begin
        QN:=0.5;
        AAA:=YPN - (YA[num-1] - YA[num - 2]) / (XA[num-1] - XA[num - 2]);
        UN:=(3 / (XA[num-1] - XA[num - 2])) * AAA;
    end;
    Y2A[num-1]:=(UN - QN * U[num - 3]) / (QN * Y2A[num - 2] + 1);
    For K:=num - 2 DownTo 0 do
        Y2A[K]:=Y2A[K] * Y2A[K + 1] + U[K];
  except
    for I := Low(Y2A) to High(Y2A) do
      Y2A[I] := 0;
  end;
end;

procedure TSpline.SplineCurve(x: array of Double;var Y: array of Double);
var
  i : Integer;
  yt : Double;
begin
  if High(Y) <> High(X) then
  begin
    for i:=Low(Y) to High(Y) do
      Y[i] := 0;
    exit;
  end;
  for i:=Low(X) to High(X) do
  begin
    Splint(X[i],yt);
    Y[i] := yt;
  end;
end;

procedure TSpline.SplinePoint(x: Double; var y: Double);
begin
  Splint(x,y);
end;

procedure TSpline.Splint(X: Double; var Y: Double);
var
  K,KLO,KHI:integer;
  H,A,B,AAA,BBB,Q,QQ:real;
begin
  KLO:=1;
  KHI:=num;
  while (KHI - KLO > 1) do
  begin
    K:=(KHI + KLO) div 2;
    if XA[K] > X then
        KHI:=K
    else
        KLO:=K;
  end;
  H:=XA[KHI] - XA[KLO];
  If H = 0 Then
  begin
//      ShowMessage('  PAUSE  "BAD  XA  INPUT"');
      Exit;
  end;
  A:=(XA[KHI] - X) / H;
  B:=(X - XA[KLO]) / H;
  AAA:=A * YA[KLO] + B * YA[KHI];
  if A = 0 then
    Q:= 0
  else
    if A > 0 then
      Q:= exp(3*ln(A))
    else
      Q:= -exp(3*ln(-A));
  if B = 0 then
    QQ:= 0
  else
    if B > 0 then
      QQ:= exp(3*ln(B))
    else
      QQ:= -EXP(3*ln(-B));
  BBB:=(Q - A) * Y2A[KLO] + (QQ - B) * Y2A[KHI];
  Y:=AAA + BBB * (H * H) / 6;
end;

procedure TSpline.Switch(i, j: Integer);
var
  temp : Double;
begin
  if (i>=Low(XA)) and (i<=High(XA)) and
     (j>=Low(XA)) and (j<=High(XA)) then
  begin
    temp := XA[i];
    XA[i] := XA[j];
    XA[j] := temp;
    temp := YA[i];
    YA[i] := YA[j];
    YA[j] := temp;
  end;
end;

end.