%FREQUENCY RESPONSE
%  Choose loop and parameters dynamically.
%    Compare 1st-order to Fig. 2.6,
%            2nd-order to Fig. 7.6-7.15
%            3rd-order to Fig. 7.20
 
r2deg = 180/pi; % for use in changing radians to degrees
w = logspace(-1,1,51); % Frequency range
Wn = 1 % Natural frequency normalized

go = 2;
while (go ~= 0)
 order = input('Enter order (1-3), 0 to quit: ','s');
	if strcmp(order,'0'),
  go = 0;
 elseif strcmp(order,'1'),
  disp('First Order');
  K = input('First Order.  Input K: ');
  out = [0 K]; % from Eq. (2.25b)
  er = [1 0]; % from Eq. (2.26b)
  den = [1 K];
 elseif strcmp(order,'3'),
  rn = input('Third Order, Type 2, zeta = 0.5.  Enter Wn/W+: ');
  out = [0 0 1 1]; % as in SS Added Pole , paragraph 6.M.5
  er = [rn 1 0 0]; % from Eq. (6.M.4)
  den = [rn 1 1 1];
 elseif (go~=0),
  alpha = input('Second Order. Input alpha: ');
  zeta = input('Input zeta: ');
  out = [2*alpha*zeta 1]; % from Eq. (6.4a)
  er = [1 2*zeta*(1-alpha) 0]; % from Eq. (6.4b)
  den = [1 2*zeta 1];
 end

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

 clf 
 subplot(211), semilogx(w, magout, 'r', w, mager, ':r')
 grid
 title('Frequency Response Magnitude.  solid = output   dotted = error')
 xlabel('Omega / Omega-n')
 ylabel('dB')
 subplot(212), semilogx(w, phout, 'g', w, pher, ':g')
 grid
  title('Frequency Response Phase.  solid = output   dotted = error')
 xlabel('Omega / Omega-n')
 ylabel('degrees')
end % while go
