unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs;
procedure SPLINE(X:array of real;var Y:array of real;
                      N:integer; YP1, YPN:real; VAR Y2:array of real);
procedure SPLINT(XA, YA, Y2A:array of real; N:integer; X:real;var Y:real);

implementation
procedure SPLINE(X:array of real;var Y:array of real;
                     N:integer; YP1, YPN:real; var Y2:array of real);
var
    U:array[0..100] of real;
    AAA,SIG,BBB,CCC,P,QN,UN:real;
    I,K:integer;
begin
    If YP1 > 9.9E+29 Then
    begin
        Y2[1]:=0;
        U[1]:=0;
    end
    Else
    begin
        Y2[1]:=-0.5;
        AAA:=(Y[2] - Y[1]) / (X[2] - X[1]);
        U[1]:=(3 / (X[2] - X[1])) * (AAA - YP1);
    end;
    For I:=2 To N - 1 do
    begin
        SIG:=(X[I] - X[I - 1]) / (X[I + 1] - X[I - 1]);
        P:=SIG * Y2[I - 1] + 2;
        Y2[I]:=(SIG - 1) / P;
        AAA:=(Y[I + 1] - Y[I]) / (X[I + 1] - X[I]);
        BBB:=(Y[I] - Y[I - 1]) / (X[I] - X[I - 1]);
        CCC:=X[I + 1] - X[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 - (Y[N] - Y[N - 1]) / (X[N] - X[N - 1]);
        UN:=(3 / (X[N] - X[N - 1])) * AAA;
    end;
    Y2[N]:=(UN - QN * U[N - 1]) / (QN * Y2[N - 1] + 1);
    For K:=N - 1 DownTo 1 do
        Y2[K]:=Y2[K] * Y2[K + 1] + U[K];
end;

procedure SPLINT(XA, YA, Y2A:array of real; N:integer; X:real;var Y:real);
label 1;
var
    K,KLO,KHI:integer;
    H,A,B,AAA,BBB,Q,QQ:real;
begin
    KLO:=1;
    KHI:=N;
1:  If KHI - KLO > 1 Then
    begin
        K:=(KHI + KLO) Div 2;
        If XA[K] > X Then
            KHI:=K
        Else
            KLO:=K;
        GoTo 1;
    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;
end.