function secondrand()
% SECONDRAND -- Demonstrate relationships between pole locations and response of second-order
%	systems.
%
%	Click the mouse on the S plane to locate a complex pole.  Its conjugate appears automatically,
%	and two system behaviors are 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.
%
%   The two responses shown each time are randomly selected from all the possibilities for:
%   -Different 2nd order LTI system configurations.
%   -Different initial conditions.
%   -Differnt step, constant, or zero inputs.
%   However, the values range over a finite range to keep the outputs within a moderate range of 
%   amplitudes near one.
%
%	Pole locations are coarsely quantized, for simplicity and repeatablility.


%   C. Sullivan, 7/2004
%   Based on the program step2demo.m that is similar but gives response of one second-order system
%   to two particular inputs:  History of that:
%	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')
    time=clock;
     
    ht=text(-2.7,-3.8,num2str(time(3:6)+3),'fontsize',6);

% 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.
            C = 2*(rand-0.5);
            B = 2*(rand-0.5);
            AA = 2*(rand-0.5);
			y1 =real((C - AA*exp(-a*t) + B*exp(-b*t))) ;

% 			Response from table #2
            C = 2*(rand-0.5);
            B = 2*(rand-0.5);
            AA = 2*(rand-0.5);
			y2 =real((C - AA*exp(-a*t) + B*exp(-b*t))) ;
		else
			y1 = 1 - (1 + wn*t).*exp(-wn*t);	% Critically damped response
			y2 = wn^2*t.* exp(-wn*t);
		end
	
%		Compute performance measures, where appropriate
		if zeta<1
		%	[ypeak, ipeak] = max(y1); 
		%%	Mp = ypeak-1;						% Peak overshoot
		%	tp = t(ipeak);						% Peak time
			wd = abs(w(1));						% Ringing frequency
        end
        %else
		%	Mp = 0;	
        %end
		%n10 = max(find(y1(1:ipeak)<=0.1));	% Risetime
		%n90 = max(find(y1(1:ipeak)<= 0.9));
		%tr = t(n90) - t(n10);
		%ns = max(find(abs(y1-1) >= 0.01));	% Settling time
		%ts = t(ns);
		
		
		subplot(1,2,2),
			plot(t, y1, 'b', t, y2, 'g'), grid
			axis('square')
			%if not(ascale), axis([0,15, -0.2, 1.8]), end
			xlabel('time'), ylabel('responses')
				if zeta<1
				    titlestring = ['\omega_d=', num2str(wd,3)];
                else
                    titlestring = ' ';
                end
	
         title(titlestring)
       subplot(1,2,1)
         time=clock;
         delete(ht);
         ht=text(-2.7,-3.8,num2str(time(3:6)+3),'fontsize',6);
	end

end