PROGRAM DemCheb;
{Predvadi Chebyshevovu interpolaci funkce od Xmin do Xmax
funkce f(x) = 1 pro x < delt, f(x) = 0 pro x >= delt
pro ctyri zadane rady Cheb. polynomu
zapisuje vysledku do chebd1.dat
}
const
  np = 60;
TYPE
   Double = Real;
   RealArrayNP = ARRAY [1..np] OF real;
var
   C,cint,cder : RealArrayNP;
   xmax,xmin,delt,eps,sumeps,dx1,yf,ych,y1f,y1ch,yif,yich,xh : Real;
   nrad,mrad,ndel,i : integer;
   f1 : text;
FUNCTION func(y: real): real;
BEGIN
        if (y < delt) then func :=1
        else func :=  0;
        END;
PROCEDURE chebft(a,b: real;
               VAR c: RealArrayNP;
                   n: integer);
{Chebyshevuv fit funkce FUNC v intervalu [A,B] stupne maximalne N,
 procedura pocita koeficieny C[1..N] takove, ze
  FUNC ~= (suma od 1 do N z C[K]*T[K-1](y)) - C[1]/2,
      kde y = x-0.5*[b+a]/0.5[b-a]
 Tato procedura ma byt  uzivana s neprilis velkymi N (30 - 50),
 pole C by melo byt zkraceno na mensi hodnote M tak ze C[M+1]
 a dalsi jsou zanedbatelne.}
CONST
   pi = 3.141592653589793;
VAR
   k,j: integer;
   y,fac,bpa,bma: real;
   sum: double;
   f: ^RealArrayNP;
LABEL 1;
BEGIN
   new(f);
   bma := 0.5*(b-a);
   bpa := 0.5*(b+a);
   FOR k := 1 TO n DO BEGIN
      y := cos(pi*(k-0.5)/n);
      f^[k] := func(y*bma+bpa)
   END;
   fac := 2.0/n;
   FOR j := 1 TO n DO BEGIN
      sum := 0.0;
      FOR k := 1 TO n DO
         sum := sum+f^[k]*cos(pi*(j-1)*(k-0.5)/n);
      c[j] := fac*sum
   END;
   dispose(f)
END;
FUNCTION chebev(a,b: real;
              VAR c: RealArrayNP;
                  m: integer;
                  x: real): real;
{Pocita hodnotu funkce v bodu x pomoci Chebyshevova polynomu
 v bode y urceneho z bodu x v intevalu A, B.
 Koeficienty Chebyshevova polynomu byly spocitany pomoci CHEBFT.}
VAR
   d,dd,sv,y,y2: real;
   j: integer;
BEGIN
   IF (x-a)*(x-b) > 0.0 THEN BEGIN
      writeln('pause in CHEBEV - x not in range.');
      readln
   END;
   d := 0.0;
   dd := 0.0;
   y := (2.0*x-a-b)/(b-a);
   y2 := 2.0*y;
   FOR j := m DOWNTO 2 DO BEGIN
      sv := d;
      d := y2*d-dd+c[j];
      dd := sv
   END;
   chebev := y*d-dd+0.5*c[1]
END;

VAR
   j:integer;
   NR: array[1..4] of integer;
   ypdat: array[1..4] of real;
LABEL 1;
BEGIN
   ASSIGN(f1,'chebd1.dat');
   REWRITE(f1);
   Writeln('ZADEJ xmin,xmax,nr[1-4],delt');
   Readln(xmin,xmax,nr[1],nr[2],nr[3],nr[4],delt);
 1:   Writeln(' zadej ndel');
      Readln(ndel);
      dx1 := (xmax-xmin)/(ndel - 1);
      writeln('    x         yt       ych     y1t   y1ch   yich');
   FOR i := 1 to ndel DO BEGIN
       xh := xmin+(i-1)*dx1;
       yf := func(xh);
       for j:=1 to 4 do begin
            nrad := nr[j];
            Chebft(xmin,xmax,C,nrad);
            ypdat[j]:= chebev(xmin,xmax,c,nrad,xh);
            end;
       writeln(xh,yf,ypdat[1],ypdat[2],ypdat[3],ypdat[4]);
       writeln(f1,xh:10:6,' ',yf:13:-4,' ',ypdat[1]:13:-4,' ',
                  ypdat[2]:13:-4,' ',ypdat[3]:13:-4,' ',ypdat[4]:13:-4)
       END;
   CLOSE(f1)
END.