clear
clc
close all

% Define some somewhat arbitrary values for pole and zero locations:
mp = 0.95;
pp = pi/8;
mz = 1.1;
pz = 1.234*pp;

% Make 2 complex conjugate poles and 2 complex conjugate zeros:
p = [mp*(cos(pp)+1i*sin(pp)), mp*(cos(pp)-1i*sin(pp))];  % Pole locations
z = [mz*(cos(pz)+1i*sin(pz)), mz*(cos(pz)-1i*sin(pz))];  % Zero locations
den = poly(p);         % Form denominator polynomial
num = poly(z);         % Form numerator polynomial
dt = 1/100;            % Time step (not used)
sys = tf(num,den,dt);  % Form system (not used)

% Define the Z-Transform numerically for the 2 pole, 2 zero system:
dz = -1.2:0.01:1.2;
x = dz;
y = dz;
lx = length(x);
ly = length(y);
Z  = zeros(lx,ly);
for kx = 1:lx
    for ky = 1:ly
        z_step = (x(kx)+1i*y(ky));
        Z(kx,ky) = (num(1) + num(2)*z_step^-1 + num(3)*z_step^-2)/...
            (den(1) + den(2)*z_step^-1 + den(3)*z_step^-2);
    end
end

% Define z-Transform just around the unit circle:
xw = zeros(1001,1);
yw = zeros(1001,1);
zw = zeros(1001,1);
for kw = 1:1001
    w = 2*pi/1000*(kw-1);
    xw(kw) = cos(w);
    yw(kw) = sin(w);
    z_step = (xw(kw)+1i*yw(kw));
    zw(kw) = (num(1) + num(2)*z_step^-1 + num(3)*z_step^-2)/...
        (den(1) + den(2)*z_step^-1 + den(3)*z_step^-2);
end

%======================
% PLOT RESULTS BELOW
%======================

% Plot zgrid by itself to see frequency and damping contours:
figure;
zgrid;
title('Z-Plane: Contours of Constant Frequency and Damping')
ylabel('Imaginary Part, z');
xlabel('Real Part, z');
axis square

% Plot pole-zero plot for our example system:
figure;
plot(real(p),imag(p),'b x',real(z),imag(z),'b o',...
    'LineWidth',2,'MarkerSize',10);
zgrid;
title('Z-Plane: Example Pole-Zero Plot')
ylabel('Imaginary Part, z');
xlabel('Real Part, z');
axis square

% Plot z-surface, poles, and unit circle for our example system:
figure;
hz = 30;
[X,Y] = meshgrid(dz);
plot3(0,0,0,'k-');
hold on
plot3([real(p(1)) real(p(1))], [imag(p(1)) imag(p(1))], [0 hz],'k-x',...
    'LineWidth',2,'MarkerSize',10);
plot3([real(p(2)) real(p(2))], [imag(p(2)) imag(p(2))], [0 hz],'k-x',...
    'LineWidth',2,'MarkerSize',10);
mesh(Y,X,abs(Z));
plot3(xw, yw,abs(zw),'g-','LineWidth',2,'MarkerSize',10);
hold off
zlim([0 hz]);
xlim([-1.2 1.2]);
ylim([-1.2 1.2]);
ylabel('Imaginary Part, z');
xlabel('Real Part, z');
zgrid;

