unit Interpol;
{=================================================================}
{  Âåðñèÿ ôàéëà 3.0.000                                           }
{  îò 4 äåêàáðÿ 2002 ã.                                           }
{==========================} interface {==========================}

uses gltype,MathObj;

function Lagrange(var X1,Y1 :array of RealType;X:RealType;N,M:Integer):RealType;

PROCEDURE spline(var x,y,y2,u: array of RealType; n: integer; yp1,ypn: RealType);
function splint(var xa,ya,y2a: array of RealType; n: integer;x: RealType; VAR y: RealType):Integer;
{********************************************************************************************************}

                               IMPLEMENTATION

{********************************************************************************************************}


//**********************************************************************

//     ÈÍÒÅÐÏÎËßÖÈß ÏÎËÈÍÎÌÎÌ ËÀÃÐÀÍÆÀ

//**********************************************************************
function Lagrange(var X1,Y1 :array of RealType;X:RealType;N,M:Integer):RealType;
//
//    X - ÀÐÃÓÌÅÍÒ
//    X1 - ÌÀÑÑÈ ÇÍÀ×ÅÍÈÉ ÀÐÃÓÌÅÍÒÀ
//    Y1 - ÌÀÑÑÈ ÇÍÀ×ÅÍÈÉ ÔÓÍÊÖÈÈ
//    N - ÏÎÐßÄÎÊ ÏÎËÈÍÎÌÀ
//    M - HOMEP ÝËÅÌÅÍÒÀ, C KOTOPOÃO ÍÅÎÁÕÎÄÈÌÎ ÍÀ×ÀÒÜ ÈÍÒÅÐÏÎËßÖÈÞ
//
var i,j   : Integer;
    P1,P2 : RealType;
begin
      Result:=0.0;
      for j:=0 to N-1 do begin
       P1:=1.0;
       P2:=1.0;
       for i:=0 to N-1 do begin
        if i = j then continue;
        P1:=P1*(X-X1[I+M-1]);
        P2:=P2*(X1[J+M-1]-X1[I+M-1]);
       end;
       Result:=Result+P1/P2*Y1[J+M-1]
      end
end;

//**********************************************************************

//     ÈÍÒÅÐÏÎËßÖÈß ÑÏËÀÉÍÎÌ. ÏÎËÓ×ÅÍÈÅ ÊÓÁÈ×ÅÑÊÎÃÎ ÑÏËÀÉÍÀ

//**********************************************************************
PROCEDURE spline(var x,y,y2,u: array of RealType; n: integer; yp1,ypn: RealType);
VAR
   i,k: integer;
   p,qn,sig,un: RealType;
BEGIN
   IF (yp1 > 0.99e30) THEN BEGIN
      y2[0] := 0.0;
      u[0] := 0.0
   END ELSE BEGIN
      y2[0] := -0.5;
      u[0] := (3.0/(x[1]-x[0]))*((y[1]-y[0])/(x[1]-x[0])-yp1)
   END;

   FOR i := 1 to n-2 DO BEGIN
      sig := (x[i]-x[i-1])/(x[i+1]-x[i-1]);
      p := sig*y2[i-1]+2.0;
      y2[i] := (sig-1.0)/p;
      u[i] := (y[i+1]-y[i])/(x[i+1]-x[i])-(y[i]-y[i-1])/(x[i]-x[i-1]);
      u[i] := (6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p
   END;

   IF (ypn > 0.99e30) THEN BEGIN
      qn := 0.0;
      un := 0.0
   END ELSE BEGIN
      qn := 0.5;
      un := (3.0/(x[n-1]-x[n-2]))*(ypn-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]))
   END;

   y2[n-1] := (un-qn*u[n-2])/(qn*y2[n-2]+1.0);

   FOR k := n-2 DOWNTO 0 DO y2[k] := y2[k]*y2[k+1]+u[k]
END;
//**********************************************************************

//     ÈÍÒÅÐÏÎËßÖÈß ÑÏËÀÉÍÎÌ. ÂÛ×ÈÑËÅÍÈÅ ÔÓÍÊÖÈÈ ÊÓÁÈ×ÅÑÊÈÌ ÑÏËÀÉÍÎÌ

//**********************************************************************

function splint(var xa,ya,y2a: array of RealType; n: integer;x: RealType; VAR y: RealType):Integer;
VAR
   klo,khi,k: integer;
   h,b,a: RealType;
BEGIN
   klo := 0;
   khi := n-1;
   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.0) THEN Result:=1 else Result:=0;
   if Result <> 0 then exit;

   a := (xa[khi]-x)/h;
   b := (x-xa[klo])/h;
   y := a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0
END;

END.