% Script gets z transform of closed loop response,
%   Type-2 second order loop, Eq. (7.22.3), and generates 
%   inverse transform as time series, one value per sample

clear
WnTs = .5; % Natural frequency times sample period, default
zeta = .5; % default
fprintf('Initial Values: WnTs = %5.3f -> Fs/Fn = %4.2f; zeta = %4.2f.\n',WnTs,2*pi/WnTs,zeta);
while (WnTs ~= 0)
	WnTsT = input('Input WnTs or negative of zeta (0 ends script): ');
	if ~isempty(WnTsT); % if an entry has been made
		if (WnTsT < 0), 
			zeta = -WnTsT; % - zeta is entered; negative value differentiates it from WnTs
		else
			WnTs = WnTsT; 
		end;
	end;
	if (WnTs == 0); break; end; % enter 0 to end
	Fs2Fn = 2*pi/WnTs;
	cols = ceil(8/WnTs); % approximately number of computed points; give duration of 8/Wn
	duration = zeros(1,cols); % array of zeros to control number of points
	num = [1 -1 duration]; % numerator plus zeros at end; division will stop when it becoms fractional
	den = [1 (WnTs*(2*zeta+WnTs)-2) (1-2*zeta*WnTs)]; % denominator
	TR = deconv(num,den); % divide den into num
	x= 0:cols-1; % xaxis range of plot.
	x = x*WnTs; % proper time scaling (change sample number to time)
	clf; % clear previous graphs; for v.3.5, clg #########
	plot(x,TR,'*')
	xlabel('OmegaN*t')
	ylabel('Unit Step Error Response')
	title(['Step Error with Zeta = ',num2str(zeta,'%4.2f'),', WnTs = ',num2str(WnTs,'%5.3f'),' -> Fs/Fn = ',num2str(Fs2Fn,'%4.2f')])
	grid
	fprintf('WnTs = %5.3f -> Fs/Fn = %4.2f, zeta = %4.2f, alpha = 1.\n',WnTs,Fs2Fn,zeta)
end

