[Back to GRAPHICS SWAG index] [Back to Main SWAG index] [Original]

     Program Interpolation;
    {
                 ÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛ
                 ÛÛÛÝÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÜÞÛÛÛ±±
                 ÛÛÛÝÛÛ                                      ÛÛÞÛÛÛ±±
                 ÛÛÛÝÛÛ  Least squares method interpolation  ÛÛÞÛÛÛ±±
                 ÛÛÛÝÛÛ                                      ÛÛÞÛÛÛ±±
                 ÛÛÛÝÛÛ           Aleksandar Dlabac          ÛÛÞÛÛÛ±±
                 ÛÛÛÝÛÛ    (C) 1995. Dlabac Bros. Company    ÛÛÞÛÛÛ±±
                 ÛÛÛÝÛÛ    ------------------------------    ÛÛÞÛÛÛ±±
                 ÛÛÛÝÛÛ      adlabac@urcpg.urc.cg.ac.yu      ÛÛÞÛÛÛ±±
                 ÛÛÛÝÛÛ      adlabac@urcpg.pmf.cg.ac.yu      ÛÛÞÛÛÛ±±
                 ÛÛÛÝÛÛ                                      ÛÛÞÛÛÛ±±
                 ÛÛÛÝßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßßÞÛÛÛ±±
                 ÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛÛ±±
                   ±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±±
    }

    {
        Interpolation is a technique used to calculate estimated function of
    which we have some experimental points. Firstly, we have to choose order of
    interpolation. Order of interpolation represents the order of estimated
    function. For example, if order is 5 function will be in form of:

        F(x) = a5*x^5 + a4*x^4 + a3*x^3 + a2*x^2 + a1*x + a0

        This program (procedure Interpolate) returns coefficients a0..an, where
    n is order. We have to determine the order considering positions of points.
    }
       Uses Crt, Graph;

       Type Row    = array [0..10] of real;     { Maximal order               }
            Matrix = array [0..10] of Row;      { is 10.                      }
            Data   = array [1..200] of integer; { Data for interpolation      }
            Coeff  = array [0..10] of real;     { Coefficients (also max. 10) }

       Const Order = 5;  { Order of interpolated curve. Should be smaller than  }
                         { number of points. For orders greater than 6-7 system }
                         { (of equations) becomes unstable, so curve could be   }
                         { wrong. Also, for higher orders real type should be   }
                         { changed to extended, and $N switch should be on.     }

       Var Gd, Gm, I, J : integer;
           R            : real;
           InversionOK  : Boolean;
           Mat          : Matrix;
           X, Y         : Data;
           C, YM        : Coeff;

       Function Power (Base,Pow:integer) : real;
         Var Temp : real;
           Begin
             If (Base<0) and (Pow mod 2=1) then
               Temp:=-Exp (Pow*Ln (-Base))
                         else
               Temp:=Exp (Pow*Ln (Abs (Base)));
             Power:=Temp
           End;

       Procedure MatrixInversion (Var A:Matrix; N:integer);
         Var I, J, K : integer;
             Factor  : real;
             Temp    : Row;
             B       : Matrix;
           Begin
             InversionOK:=False;
             For I:=0 to N do
               For J:=0 to N do
                 If I=J then
                   B [I,J]:=1
                        else
                   B [I,J]:=0;
             For I:=0 to N do
               Begin
                 For J:=I+1 to N do
                   If Abs (A [I,I])<Abs (A [J,I]) then
                     Begin
                       Temp:=A [I];
                       A [I]:=A [J];
                       A [J]:=Temp;
                       Temp:=B [I];
                       B [I]:=B [J];
                       B [J]:=Temp
                     End;
                 If A [I,I]=0 then Exit;
                 Factor:=A [I,I];
                 For J:=N downto 0 do
                   Begin
                     B [I,J]:=B [I,J]/Factor;
                     A [I,J]:=A [I,J]/Factor
                   End;
                 For J:=I+1 to N do
                   Begin
                     Factor:=-A [J,I];
                     For K:=0 to N do
                       Begin
                         A [J,K]:=A [J,K]+A [I,K]*Factor;
                         B [J,K]:=B [J,K]+B [I,K]*Factor
                       End
                   End
               End;
             For I:=N downto 1 do
               Begin
                 For J:=I-1 downto 0 do
                   Begin
                     Factor:=-A [J,I];
                     For K:=0 to N do
                       Begin
                         A [J,K]:=A [J,K]+A [I,K]*Factor;
                         B [J,K]:=B [J,K]+B [I,K]*Factor
                       End
                   End
               End;
             A:=B;
             InversionOK:=True
           End;

       Procedure Interpolate (X,Y:Data; var A:Coeff);
         Var I, J, K : integer;
           Begin
             For I:=0 to Order do
               For J:=0 to Order do
                 Mat [I,J]:=0;
             For I:=0 to Order do
               For J:=0 to Order do
                 For K:=1 to 10 do
                   Mat [I,J]:=Mat [I,J]+Power (X [K],I+J);
             MatrixInversion (Mat,Order);
             For I:=0 to Order do
               Begin
                 YM [I]:=0;
                 A [I]:=0
               End;
             For I:=0 to Order do
               For J:=1 to 10 do
                 YM [I]:=YM [I]+Y [J]*Power (X [J],I);
             For I:=0 to Order do
               For J:=0 to Order do
                 A [I]:=A [I]+Mat [I,J]*YM [J];
           End;

         Begin
           DetectGraph (Gd,Gm);
           InitGraph (Gd,Gm,'');
           X [1]:=100;  Y [1]:=100;
           X [2]:=150;  Y [2]:=180;
           X [3]:=200;  Y [3]:=250;
           X [4]:=250;  Y [4]:=270;
           X [5]:=300;  Y [5]:=300;
           X [6]:=350;  Y [6]:=320;
           X [7]:=400;  Y [7]:=290;
           X [8]:=450;  Y [8]:=230;
           X [9]:=500;  Y [9]:=130;
           X [10]:=550; Y [10]:=90;
           Interpolate (X,Y,C);
           SetFillStyle (SolidFill,Red);
           For I:=1 to 10 do                        {< Draws red points (input }
             Bar (X [I]-1,Y [I]-1,X [I]+1,Y [I]+1); {< values - X & Y).        }
           SetColor (Green);
           For I:=10 to 55 do                       {< Draws calculated curve. }
             Begin
               R:=0;
               For J:=0 to Order do
                 R:=R+C [J]*Power (I*10,J);
               If I=10 then
                 MoveTo (I*10,Round (R))
                       else
                 LineTo (I*10,Round (R))
             End;
           Repeat Until ReadKey<>#0;
           CloseGraph
         End.

[Back to GRAPHICS SWAG index] [Back to Main SWAG index] [Original]