function step2demo()
% STEP2DEMO -- 	Demonstrate relationships between pole locations and step response for a second
%				order system.
%	Click the mouse on the S plane to locate a complex pole.  Its conjugate appears automatically,
%	and the step response is drawn.  Click outside the graph to quit.
%
%	If the pole is close enough to the real axis, it is assumed real, and you are asked to
%	click a second pole.  Critical damping is supported, but not poles at the origin.
%
%	Pole locations are coarsely quantized, for simplicity.

%	E.W. Hansen, Engs 22, 97W
% 	Modified C.R. Sullivan 7/22/98 -- direct calculation of response.
%	Modified E.W. Hansen 1/30/99 -- 
%   Added performance measures, lines of constant wn and Mp in s plane
%   Modified 7/99 CR Sullivan Added manual/auto scale option, and intro dialog box.

button1=strcat('Fixed', '(for easy comparisons)')
button2=strcat('Variable', '(for best detailed view of a particular plot)')
ascale = menu('Choose a plot scaling option:',button1,button2);
ascale = ascale-1  % 0 for fixed, 1 for auto

menu(strvcat('An s-plane plot will be displayed.',...
'You can then click on the s-plane to place poles.',...
'To quit, click off the graph (but in the figure window).'),'OK');


% Initialize
MaxSigma = 3;
MaxOmega = 2;
PlotWidth = 2*MaxOmega;
corner = sqrt(MaxSigma^2 + MaxOmega^2);
Nt = 200;

th = linspace(pi/2, 3*pi/2, 51);	% Lines of constant wn
semicirc = exp(j*th);
Mp = [0.01, 0.1, 0.25];				% Lines of constant Mp
zeta = abs(log(Mp)) ./ sqrt(pi^2 + log(Mp).^2);
		
q = 0.05;							% Quantization on pole locations
OD = 1;								% Overdamped, critically damped, underdamped
CD = 2;
UD = 3;

h0 = [];							% Initialize object handles
h1 = [];						
h2 = [];

A = [-MaxSigma, -MaxSigma+PlotWidth, -MaxOmega, MaxOmega];

clf
subplot(1,2,1), 
	plot([-MaxSigma, MaxSigma/2], [0, 0], 'm', [0, 0], [-MaxOmega, MaxOmega], 'm'); grid
	hold on	
		
		for wn = [1, 2, 3],
			plot(wn*semicirc, 'k--');		% Constant wn lines
		end
		for z = zeta,
			th = pi-acos(z);
			plot([0, corner*cos(th)], [0, corner*sin(th)], 'k--')
		end		
	hold off
	axis(A), axis('square')
	title('CLICK ON A POLE LOCATION')

% Big event loop
while(1)
% GET POLE LOCATIONS OFF THE GRAPH
	subplot(1,2,1)
		[s, w] = ginput(1);			% Get a pole location off the graph
		s = q * round(s/q);			% Round the locations
		w = q * round(w/q);
		if (s<A(1) | s>A(2) | w<A(3) | w>A(4))	% Outside the frame?
			name = inputdlg('Enter text to add to the display:   ');
			if ~isempty(name),
 				text(-2.7,-1.8,name);
 			end
			return
		end
		
		if exist('h1'), set(h1, 'Visible', 'off'); end	% Clear previous marks
		if exist('h2'), set(h2, 'Visible', 'off'); end
		
		
	if w==0,					% Real pole?
		hold on
		h0 = plot(s, w, 'kx');
		hold off;
		title('SELECT A SECOND REAL POLE')
		[s2, w2] = ginput(1);
		s2 = q * round(s2/q);
		if s==s2,
			damping = CD;
		else
			damping = OD;
		end
		s = [s, s2];
		w = [0, 0];
		polelocation = [num2str(s(1)), ' , ', num2str(s(2))];	 
	else
		damping = UD;
		s = [s, s];				% Make a conjugate pair
		w = [w, -w];
		polelocation = [num2str(s(1)), ' ± j',num2str(abs(w(1)))];
	end

	p=s+j*w;					% Vector of pole locations
	a=-p(1);
	b=-p(2);
	wn=sqrt(a*b);				% Undamped natural frequency
	zeta=(a+b)/2/wn;			% Damping ratio
	if zeta<1,
		polelocation = [polelocation, ';   \zeta = ', num2str(zeta,3), ', \omega_n = ', num2str(wn,3)];
	end
	
	if exist('h0'), set(h0, 'Visible', 'off'); end
	hold on
	h1 = plot(s, w, 'kx'); title(['s = ', polelocation]);
	hold off
	
	if(damping==UD),
		hold on,
		h2 = plot([s(1), 0, s(2)], [w(1), 0, w(2)], 'r');
		hold off
	end
	xlabel('Real(s)'), ylabel('Imag(s)')
	
% COMPUTE AND PLOT THE RESPONSE	
	if(min(abs(s))==0 & min(abs(w))==0),	% System poles at the origin?
		title('SELECT SYSTEM POLES NOT AT THE ORIGIN')
	else
% 		p=s+j*w;
% 		a=-p(1);
% 		b=-p(2);
% 		wn=sqrt(a*b);
% 		zeta=(a+b)/2/wn;	
		
		if ascale		
		   if zeta~=0,								% Poles off imaginary axis
			tmax = 8/min(abs(s));				% 8 time constants for settling
		   else
			tmax = 8*2*pi/wn;					% or 8 full cycles
		   end
		else
		    tmax = 15;                          %fixed scale
		end
		
	
		t = linspace(0, round(tmax), Nt);		% Time vector for plotting		
		if zeta~=1,								% If we're not critically damped
% 			Response from the table #10.
			ystep =real((1 - b/(b-a)*exp(-a*t) + a /(b-a)*exp(-b*t))) ;

% 			Response from table #2
			yimp =real( a*b/(b-a) * (exp(-a*t) - exp(-b*t)) );
		else
			ystep = 1 - (1 + wn*t).*exp(-wn*t);	% Critically damped response
			yimp = wn^2*t.* exp(-wn*t);
		end
	
%		Compute performance measures, where appropriate
		if zeta<1
			[ypeak, ipeak] = max(ystep); 
			Mp = ypeak-1;						% Peak overshoot
			tp = t(ipeak);						% Peak time
			wd = abs(w(1));						% Ringing frequency
		else
			Mp = 0;	
		end
		n10 = max(find(ystep(1:ipeak)<=0.1));	% Risetime
		n90 = max(find(ystep(1:ipeak)<= 0.9));
		tr = t(n90) - t(n10);
		ns = max(find(abs(ystep-1) >= 0.01));	% Settling time
		ts = t(ns);
		
		
		subplot(1,2,2),
			plot(t, ystep, 'b', t, yimp, 'g'), grid
			axis('square')
			if not(ascale), axis([0,15, -0.2, 1.8]), end
			xlabel('time'), ylabel('responses')
			if(damping==OD),
				titlestring = ['OVERDAMPED,  t_r=', num2str(tr, 3)];
			else if (damping == CD),
				titlestring = ['CRITICALLY DAMPED,  t_r=', num2str(tr, 3)];
			else 
% 				titlestring = ['\zeta = ', num2str(zeta,3), ', \omega_n = ', num2str(wn,3)];
				if Mp>0,
					titlestring = ['PO=',num2str(Mp,2), ', '];
				else
					titlestring = [];
				end
				titlestring = [titlestring, 't_r=',num2str(tr,3), ', t_s=',num2str(ts,3),...
								', \omega_d=', num2str(wd,3)];
			end
			end
			title(titlestring)			
	end

end