% function [eigv] = OL2Ellipse (num,den,r1,r2,f,x0,y0);
%
% Function OL2Ellipse calculates whether an open-loop transfer function
%
% GOL(s) = num/den
%
% intersects ellipse with the following parameters:
% r1 - major axis, r2 - minor axis, f - angle of major axis, 
% x0 - displacement of center of ellipse on real axis, 
% y0 - displacement of center of ellipse on imaginary axis.
%
% The result are eigenvalues (eigv)
%
% The calculation is based on Shorten-Curran-Wulff theorem

function [eigv] = OL2Ellipse (num,den,r1,r2,f,x0,y0);

lennum = length(num);
lenden = length(den);

num_tmp = [zeros(1,lenden-lennum) num];

if (roots(num_tmp+den) > 0)						% Testing stability of the system
   disp('Closed-loop system is unstable')
else
   disp('Closed-loop system is stable')
end;

numv = num(lennum:-1:1);					% Flip both vectors (from the lowest to the highest power)
denv = den(lenden:-1:1);
indxn = 0:lennum-1;
indxd = 0:lenden-1;

i = sqrt(-1);									% Transformation from s (Laplace) to frequency domain (s = i*w)
multn = i.^indxn;
multd = i.^indxd;
numv = numv.*multn;
denv = denv.*multd;

numv = numv(lennum:-1:1);					% Flip both vectors back (from the highest to the lowest power)
denv = denv(lenden:-1:1);

numr = real(numv);							% Find real and imaginary components
numi = imag(numv);
denr = real(denv);
deni = imag(denv);

numReG1 = conv(numr,denr);
numReG2 = conv(numi,deni);
numReG = numReG1 + numReG2;

denReG1 = conv(denr,denr);
denReG2 = conv(deni,deni);
denReG = denReG1 + denReG2;
denImG = denReG;

numImG1 = conv(numr,deni);
numImG2 = conv(numi,denr);
numImG = -numImG1 + numImG2;

numReG =  [zeros(1,length(denReG)-length(numReG)) numReG] - x0*denReG;		% Shift real and imaginary components
numImG =  [zeros(1,length(denImG)-length(numImG)) numImG] - y0*denImG;		% by displacements x0 and y0

K1 = (r1-r2)/r1;
K2 = tan(f);
tmp = 1/(1+K2^2);
Atr = [1-K1*tmp -K1*K2*tmp; -K1*K2*tmp 1-K1*K2^2*tmp];			% Transformation matrix (shrinking in direction of
																					% ellipse major axis)
numRe1 = Atr(1,1)*numReG + Atr(1,2)*numImG;							% Transformed real and imaginary components
numIm1 = Atr(2,1)*numReG + Atr(2,2)*numImG;

tmp1 = numRe1(length(numRe1):-1:1);					% Flip vectors
tmp2 = numIm1(length(numIm1):-1:1);
tmp3 = denReG(length(denReG):-1:1);

while (tmp1(1) == 0 & tmp2(1) == 0 & tmp3(1) == 0)		% Delete redundant components
   tmp1 = tmp1(2:length(tmp1));
   tmp2 = tmp2(2:length(tmp2));
   tmp3 = tmp3(2:length(tmp3));   
end;

numRe1 = tmp1(length(tmp1):-1:1);					% Flip vectors (back)
numIm1 = tmp2(length(tmp2):-1:1);
denReG = tmp3(length(tmp3):-1:1);

num2 = conv(numRe1,numRe1) + conv(numIm1,numIm1);		% Numerator squared 
den2 = conv(denReG,denReG);									% Denominator squared

lennum = length(num2);
lenden = length(den2);

											
tmp1 = num2(lennum:-1:1);			% Flip vectors	
tmp2 = den2(lenden:-1:1);
indxn = 0:2*lennum-1;				
indxd = 0:2*lenden-1;

i = sqrt(-1);
multn = (-i).^indxn;					
multd = (-i).^indxd;
tmp3 = [];
tmp4 = [];

for k = 1:lennum,
   tmp3 = [tmp3 tmp1(k) 0];		% Make transformation w1 = w^2
end
for k = 1:lenden,
   tmp4 = [tmp4 tmp2(k) 0];
end

tmp3 = tmp3.*multn;					% Transformation from frequency to s-domain
tmp4 = tmp4.*multd;

num2 = tmp3(2*lennum-1:-1:1);		% Final transfer function in s-domain
den2 = tmp4(2*lenden-1:-1:1);

[A,b,c,d] = tf2ss (num2,den2);										% Transformation from transfer function to state-space
eigv = eig(A*(A-b*c/(d-r2^2)));										% Calculation of eigenvalues
[res,resi] = find(abs(imag(eigv))<1e-10 & real(eigv) < 0);		% Find negative real roots (if exist)

if ~isempty(res)
   disp('Ellipse has been touched or intersected')				% There are negative real roots
else
   disp('Ellipse has not been touched')
end;