%theta = [ 0.4 0.5 2.5 0.1];
theta = [ 1.4 0.8 0.9 0.1 ];

B = [ theta(3) theta(4) ]; %% B i A procesa
A = [ 1 theta(1) theta(2)];

% minimalni stupnjevi za slucaj s integralnim djelovanjem
% n = 2
% deg R, S, T = deg A = n
% deg Am = n
% deg Ao = n

s = tf('s');
z=tf('z',1);
%Am = s^2 + 2*s + 1;
Am = s^2 / 1.8^2 + 2*0.6901/1.8*s+1;
Am = c2d(Am,1,'matched')

Ao = z^2;

Az = Am*Ao

syms z
Az_coefs = Az.num{1}
syms r1 r2 s0 s1 s2
A_pol = poly2sym(A,z)
B_pol = poly2sym(B,z)

% u R ubacujemo integralno djelovanje
R_pol = (r2 + r1*z)*(z-1);
S_pol = s2 + s1*z + s0*z^2;

lijeva_strana = flip(coeffs(A_pol*R_pol + B_pol*S_pol,z))

% rijesi diofantsku jednadzbu
rjesenje = solve(lijeva_strana-Az_coefs)

r1 = rjesenje.r1
r2 = rjesenje.r2

s0 = rjesenje.s0
s1 = rjesenje.s1
s2 = rjesenje.s2

syms z
Ao = z^2;
T_pol = polyval(Am.num{1},1) / polyval(B, 1) * Ao

disp 'polinom R'
R = vpa(expand(eval(R_pol)),4)

disp 'polinom S'
S = vpa(expand(eval(S_pol)),4)

disp 'polinom T'

T = vpa(T_pol,4)