%Closed-Loop Responses from Open Loop Parameters
%Example of Open2cls with another pole (4 poles, 2 at 0 freq.)
%  frequency of the zero is implied by pole-1 frequency and ratio
 
r2deg = 180/pi; % for use in changing radians to degrees
% w -> Frequency points from 10^2 to 10^4, centered on 1000.
		w = logspace(2,4,201);
% t -> Time points from 0 to 10.  Modify to shift plots.
		t = linspace(0,.05,200); % Time points

r = 8; p1 = 2000; p2 = 4000; kdB = 193; typgr = 'b';
while (r >= 0)
% r is half the ratio of pole 1 freq to zero freq.
	r0 = input('Enter freq. ratio, r=(zero)/(2*pole1)), negative to quit: ');
	if isempty(r0),
		disp(r);
	else 
		r = r0;
	end % if
	if (r >= 0),
		p10 = input('Enter pole1 frequency in rad/sec: ');
		if isempty(p10),
			disp(p1);
		else
			p1 = p10;
		end % if
		z1 = p1/r;
		zero = z1
		p20 = input('Enter pole2 frequency in rad/sec: ');
		if isempty(p20),
			disp(p2);
		else
			p2 = p20;
		end % if
		z = [-z1]; % open loop zero
		p = [0;0;-p1;-p2]; % open loop poles
		kdB1 = input('Enter gain constant in dB: ');
		if isempty(kdB1),
			disp(kdB)
		else
			kdB = kdB1;
	 end;
		k = 10^(kdB/20);
		[nopen, dopen] = zp2tf(z,p,k); % open-loop num and den from poles, zeros, k
		den = nopen + dopen; % closed-loop denominator
		out = nopen; % numerator for output step response
		er = dopen; % numerator for error step response
		% integration of er by shifting right s^0 through s^4:
		rer = er*[0 1 0 0 0 ; 0 0 1 0 0 ; 0 0 0 1 0 ; 0 0 0 0 1 ; 0 0 0 0 0]; 

		disp('');
		disp('b=Bode plot (open loop)');
		disp('o=Step Output');
		disp('e=Step Error');
		disp('r=Ramp Error');
		disp('f=Frequency Response');
		typgr1 = input('Enter b, o, e, r, or f. ', 's');
		if ~isempty(typgr1), typgr = typgr1; end;
		clg;
		if strcmp(typgr,'b')
			bode(nopen,dopen,w); % open-loop transfer function
			title('Bode Plot (Open-Loop)')
 	elseif strcmp(typgr,'e')
			step(er,den,t); % error response to a step
			title('Error Response to Step')
			grid;
		elseif strcmp(typgr,'o')
			step(out,den,t); % output response to a step
			title('Output Response to Step')
			grid;
		elseif strcmp(typgr,'r')
			step(rer,den,t); % error response to a unit ramp
			title('Error Response to Ramp')
			grid;
		else
 		hout = freqs(out, den, w); % frequency response to output
 		magout =  20*log10(abs(hout)); % magnitude of hout
 		phout = r2deg*angle(hout); % phase of hout
 		her = freqs(er, den, w); % frequency response to error
 		mager = 20*log10(abs(her)); % magnitude of her
 		pher = r2deg*angle(her); % phase of her
			subplot(211), semilogx(w, magout, 'r', w, mager, '--r')
 		grid
 		title('Frequency Response Magnitude.  solid = output   dash = error')
 		xlabel('Omega / Omega-n')
 		ylabel('dB')
 		subplot(212), semilogx(w, phout, 'g', w, pher, '--g')
 		grid
 	 title('Frequency Response Phase.  solid = output   dash = error')
 		xlabel('Omega / Omega-n')
 		ylabel('degrees')
		end % if typgr
	end % if r
end % while
