clear
clc
close all

% Define some somewhat arbitrary values for pole and zero locations:
mp = [0.95, .98];
pp = [pi/8, pi/5];
mz = 1.1;
pz = 1.234*pp(1);

% Make 4 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)
n = length(p);         % n = number of poles
m = length(z);         % m = number of zeros

% Define the Z-Transform numerically for the 4 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_step_num = 0; z_step_den = 0;
        for km = 1:m+1
            z_step_num = z_step_num + num(km)*z_step^-(km-1);
        end
        for kn = 1:n+1
            z_step_den = z_step_den + den(kn)*z_step^-(kn-1);
        end
        Z(kx,ky) = z_step_num/z_step_den;
    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));
    z_step_num = 0; z_step_den = 0;
    for km = 1:m+1
        z_step_num = z_step_num + num(km)*z_step^-(km-1);
    end
    for kn = 1:n+1
        z_step_den = z_step_den + den(kn)*z_step^-(kn-1);
    end
    zw(kw) = z_step_num/z_step_den;
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 = 50;
[X,Y] = meshgrid(dz);
plot3(0,0,0,'k-');
hold on
for k1 = 1:n
    plot3([real(p(k1)) real(p(k1))], [imag(p(k1)) imag(p(k1))], [0 hz],...
        'k-x','LineWidth',2,'MarkerSize',10);
end
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;

