%Z-plane plot for type-2 third-order loop with G given by Eq. (7.18) and (7.19)
%  No hold present.
% Locus is expected to cross unit circle at z = -1 so gain is given there.

clear
r=10; % r = rpz is the (non-zero) pole frequency divided by zero frequency
WpT=10;  % T is sampling period; Wp is pole radian frequency
ShowUC = 1; % 2 to show white unit circle, 1 for black, 0 for none
while(r~=0)
	fprintf('rpz = %d.  ',r);
	rr=input('input rpz (return to retain, 0 to end)  ');
	if ~isempty(rr), r=rr; end; % retain previous r if no input
	if r>1
		fprintf('WpT = %d.  ',WpT);
		WW=input('input OmegaP*Ts (return to retain)  ');
		if WW, WpT=WW; end;
		WzT=WpT/r; % zero radian frequency
		fprintf('\nrpz = Fpole/Fzero = %d,\nFpole/Fsample = %d,\nFzero/Fsample = %d.\n',r,WpT/2/pi,WzT/2/pi)
		a=exp(-WpT);
		p=[1;1;a]; % vector of pole frequencies in z plane
		z=[1-(1-a)*WpT/((1-a)*(r-1)+WpT);0]; % vector of zero frequencies in z plane
		k=(1-a)*(1-1/r)+WpT/r; % multiplies numerator
		[num,den]=zp2tf(z,p,k); % conversion from 3 items above to open-loop G numerator and denominator vectors
			% first vector element is coefficeint of highest power of z, etc.
		kzneg1=-den*[-1;1;-1;1]/(num*[-1;1;-1;1]); % gain multiplier of G to give closed loop pole
		fprintf('At z = -1, OmegaN^2*Ts/OmegaZ = %d\n and OmegaN*Ts = %d\n\n',kzneg1,sqrt(kzneg1*WpT/r));
		clf; % clear previous graphs; for v.3.5, clg #########
		rlocus(num,den) % do Evans plot
		title(['Z-plane Evans Plot for: rpz = ',num2str(r),', OmegaPTs = ',num2str(WpT),', OmegaZTs = ',num2str(WzT)])
		grid
		xy = linspace(0,2*pi);
		if ShowUC==1
			hold on; 
			plot(cos(xy),sin(xy),'k-');  % black unit circle.  Use if window has white background.
			hold off;
		end % if for black unit circle
		if ShowUC==2
			hold on; 
			plot(cos(xy),sin(xy),'w-'); % white unit circle.  Use if window has black background.
			hold off;
		end % if for white unit circle
	end %r
end %while

