%% wartosci z projektu 1
z=tf('z')
s=tf('s')
kz=0.778
Tz=11.5
tauz=26.3
Ts=1.6
w0=0.4*kz
n1=7;
n2=9;
k=n1/n2
m1=n1
m2=4*n2
tt=2*n1+10
obiekt=tf([k],[1,m2,m1,1],'IODelay',tt)
dysk.Variable = 'z^-1';
dysk=c2d(obiekt,Ts)
C=tf([-0.25 1],[1],Ts)
%% Dane do równania diofantycznego
Num = dysk.num{1}
Den = dysk.den{1}
NumC=C.num{1}
dA=numel(Den)-1
dB=numel(Num)-2
dC=numel(NumC)-1
d=15
dF=15-1
dG=max(dA-1, dC-15)
lambda=0.5
b0=Num(4)
%% Elementy f i g
for i = 0:dF
 eval(['syms f', num2str(i)])
end
for i = 0:dG
 eval(['syms g', num2str(i)])
end
%% Równanie diofantyczne
a0=Den(1)
a1=Den(2)
a2=Den(3)
a3=Den(4)
c0=0
c1=0
c2=NumC(1)
c3=NumC(2)

eqns = [
 c2 == f1*a0 + a1
 c3 == f2*a0 + f1*a1 + a2
 0 == f3*a0 + f2*a1 + f1*a2 + a3
 0 == f4*a0 + f3*a1 + f2*a2 + f1*a3
 0 == f5*a0 + f4*a1 + f3*a2 + f2*a3
 0 == f6*a0 + f5*a1 + f4*a2 + f3*a3
 0 == f7*a0 + f6*a1 + f5*a2 + f4*a3
 0 == f8*a0 + f7*a1 + f6*a2 + f5*a3
 0 == f9*a0 + f8*a1 + f7*a2 + f6*a3
 0 == f10*a0 + f9*a1+ f8*a2 + f7*a3
 0 == f11*a0 + f10*a1 + f9*a2 + f8*a3
 0 == f12*a0 + f11*a1 + f10*a2 + f9*a3
 0 == f13*a0 + f12*a1 + f11*a2 + f10*a3
 0 == f14*a0 + f13*a1 + f12*a2 + f11*a3
 0 == f14*a1 + f13*a2 + f12*a3 + g0
 0 == f14*a2 + f13*a3 + g1
 0 == f14*a3 + g2];

sol=solve(eqns)
F(1) = 1;
for i = 1:dF
 F(i+1) = eval(['double(sol.f',num2str(i),')']);
end
for i = 0:dG
 G(i+1) = eval(['double(sol.g',num2str(i),')']);
end

%% Wielomiany
T=NumC
S=C
FF=fliplr(F)
GF=tf([FF],[1],Ts)
GB=tf(Num,1,Ts)
GR=(GF*GB+(lambda/b0)*C)
[NumR,DenR] = tfdata(GR)
celldisp(NumR)