%% wartosci z projektu 1
Ts=1.6
z=tf('z',Ts)
z.Variable = 'z^-1'
s=tf('s')
kz=0.778
Tz=11.5
tauz=26.3
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=c2d(obiekt,Ts)
dysk.Variable='z^-1'
Gdysk=dysk*z
C=tf([1 -0.25],1,Ts);
GC=C*z^-1
%% Dane do równania diofantycznego
Num = Gdysk.num{1}
Den = Gdysk.den{1}
NumC=C.num{1}
dA=numel(Den)-1
dB=numel(Num)-2
dC=numel(NumC)-1
d=16
dF=d-1
dG=max(dA-1, dC-15)
lambda=0.75
b0=Num(1)
%% 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=NumC(1)
c1=NumC(2)

eqns = [
 c1 == f1*a0 + a1
 0 == 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 == f15*a0 + f14*a1 + f13*a2 + f12*a3
 0 == f15*a1 + f14*a2 + f13*a3 + g0
 0 == f15*a2 + f14*a3 + g1
 0 == f15*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
GB=0.02417+0.02331*z^-1+1.248e-05*z^-2
GFf=tf(F,1,Ts)
GF=GFf*z^-16
GGg=tf(G,1,Ts)
GG=GGg*z^-2
GA=1-1.6712*z^-1+0.7323*z^-2-1.469e-18*z^-3
S=GG
T=GC
R=GB*GF+(0.5/b0)*GC
%% Liczniki
NumA=GA.num{1}
NumB=GB.num{1}
NumC=GC.num{1}
NumS=GG.num{1}
NumT=NumC
NumR=R.num{1}