% Define Loop Gain Function
s1 = 1+1j;
p1 = -2;
Tfunc = @(s) (s-s1).*(s-(conj(s1)))./(s-p1);

% Define Gamma Contour Function
r = 8;
center = [-1 0];
GammaFunc = @(theta) r*cos(-theta) +center(1) + 1j*( r*sin(-theta)+center(2));

% Get numerical Loop Gain
x = linspace(-10,10,101);
[sig, w] = meshgrid(x,x);
s = sig + 1j*w;
T = Tfunc(s);

% Plot |T(s)| and Phase(T(s)) in Complex Plane
figure(1)
subplot(2,1,1)
s = surf(sig,w,20*log10(abs(T)));
s.EdgeColor = 'none';
xlabel('Re\{s\}')
ylabel('Im\{s\}')
zlabel('||T(s)||_{dB}')

subplot(2,1,2)
s = surf(sig,w,angle(T)*180/pi);
s.EdgeColor = 'none';
xlabel('Re\{s\}')
ylabel('Im\{s\}')
zlabel('\angle T(s)')



% Evaluate Numerical T(Gamma)
theta = linspace(0,2*pi,100);
gamma = GammaFunc(theta);
gamma_on_surf = Tfunc(gamma);

% Plot results
figure(1)
subplot(2,1,1)
hold on;
plot3(real(gamma),imag(gamma),20*log10(abs(gamma_on_surf)),'-r','LineWidth',3)

subplot(2,1,2)
hold on;
plot3(real(gamma),imag(gamma),angle(gamma_on_surf)*180/pi,'-r','LineWidth',3)

figure(2)
subplot(2,2,1)
plot(theta,abs(gamma_on_surf),'LineWidth',3);
hold on;
xlabel('\theta')
ylabel('|T(\Gamma)|')

subplot(2,2,3)
plot(theta,angle(gamma_on_surf),'LineWidth',3);
hold on;
xlabel('\theta')
ylabel('\angle T(\Gamma)')

subplot(2,2,[2 4])
hold on;
plot(real(gamma_on_surf), imag(gamma_on_surf),'LineWidth',3);
plot(0, 0,'ok','LineWidth',3);
grid on
xlabel('Re\{T(\Gamma)\}')
ylabel('Im\{T(\Gamma)\}')
