%FREQUENCY RESPONSE
% derivative of output and integral of error
% reference Fig. 7.16 and 7.17
%  Choose loop and parameters dynamically.
 
r2deg = 180/pi; % for use in changing radians to degrees
w = logspace(-1,1,51); % Frequency range
Wn = 1 % Natural frequency normalized

alpha = 0;
while (alpha >= 0)
	alpha = input('Second Order. Input alpha (negative to quit): ');
	if (alpha >= 0)
		zeta = input('Input zeta: ');
 	outderiv = [2*alpha*zeta 1 0]; % from Eq. (6.4a)
 	erint = [0 1 2*zeta*(1-alpha)]; % from Eq. (6.4b)
 	den = [1 2*zeta 1];

 	hout = freqs(outderiv, den, w);
 	magout =  20*log10(abs(hout));
 	phout = r2deg*angle(hout);
 	her = freqs(erint, den, w);
 	mager = 20*log10(abs(her));
	 pher = r2deg*angle(her);

	 clf % For v. 4, clf #########
 	subplot(211), semilogx(w, magout, 'r', w, mager, '--r')
 	grid
 	title('Magnitude: solid = derivative of output   dashed = integral of error')
 	xlabel('Omega / Omega-n')
 	ylabel('dB')
 	subplot(212), semilogx(w, phout, 'b', w, pher, '--b')
 	grid
 	title('Phase:  solid = derivative of output   dashed = integral of error')
 	xlabel('Omega / Omega-n')
	 ylabel('degrees')
	end % if
end % while go
