%                   
%   BODEDEMO1.M  -  This demo gives an example of some frequency response plots 
%                        for  H(s) = (s+0.1)/[s(s+1)] = b(s)/a(s)               
%
%   This file uses some Maltab built-in functions to plot the frequency response 
%   for the above system (and its individual components) - Bode, Nyquist, and 
%   Nichols plots.
%
%   A related file, BODEDEMO2.M, does the same thing, but we let Matlab produce its
%   own default plots (eg. no left hand arguments are used in the Bode commands).  
%   Once you understand the process, this approach is clearly the easier path
%   to follow!  The current file, although more work, gives a better idea of the
%   details.
%
%   File prepared by J. R. White, UMass-Lowell (last update: Feb. 2020)
%

      clear all,  close all,   nfig = 0;
%
%   set frequency range (logarithmic spacing)              
      w = logspace(-2,2,160);        
%           
%   data for the s+.1 factor                     
      b = [1 .1]; a = 1;  sys1 = tf(b,a); 
      [mag1,phase1] = bode(sys1,w);                        
      mag1 = 20*log10(mag1);     % needed to convert to db        
%
%   data for the 1/s factor                                
      b = 1; a = [1 0];   sys2 = tf(b,a);               
      [mag2,phase2] = bode(sys2,w);                   
      mag2 = 20*log10(mag2);     % needed to convert to db    
%
%   data for the 1/(s+1) factor            
      b = 1; a = [1 1];  sys3 = tf(b,a);          
      [mag3,phase3] = bode(sys3,w);                      
      mag3 = 20*log10(mag3);     % needed to convert to db         
%
%   composite magnitude and phase                           
      mag = mag1+mag2+mag3;  phase = phase1+phase2+phase3;        
%
%   Note: The magnitude and phase of the transfer function are stored in 3-D arrays
%            eg.  mag(#outputs, # inputs, # frequency values)
%         where this example has only a single input and single output
%                            
%   Bode magnitude plots                                                 
      tm1(:,1) = mag1(1,1,:);     tm2(:,1) = mag2(1,1,:);  
      tm3(:,1) = mag3(1,1,:);     tm(:,1) = mag(1,1,:);         
      nfig = nfig+1;   figure(nfig)
      subplot(2,2,1),semilogx(w,tm1,'r','LineWidth',2),grid     
      title('H(s) = s+0.1')
      xlabel('frequency (rad/sec)'),ylabel('|H(jw)| in db')  
      v = axis;  v(3)= -40;  v(4) = 40;   axis(v);    
      subplot(2,2,2),semilogx(w,tm2,'b','LineWidth',2),grid
      title('H(s) = 1/s')
      xlabel('frequency (rad/sec)'),ylabel('|H(jw)| in db'),axis(v);       
      subplot(2,2,3),semilogx(w,tm3,'g','LineWidth',2),grid
      title('H(s) = 1/(s+1)')
      xlabel('frequency (rad/sec)'),ylabel('|H(jw)| in db'),axis(v);             
      subplot(2,2,4),semilogx(w,tm,'k','LineWidth',2),grid
      title('H(s) = (s+0.1)/[s(s+1)]')
      xlabel('frequency (rad/sec)'),ylabel('|H(jw)| in db'),axis(v);                     
%                        
%   Bode phase plots                                 
      tp1(:,1) = phase1(1,1,:);     tp2(:,1) = phase2(1,1,:);  
      tp3(:,1) = phase3(1,1,:);     tp(:,1) = phase(1,1,:);          
      nfig = nfig+1;   figure(nfig)          
      subplot(2,2,1),semilogx(w,tp1,'r','LineWidth',2),grid    
      title('H(s) = s+0.1')
      xlabel('frequency (rad/sec)'),ylabel('angle')
      v = axis;  v(3)= -100;  v(4) = 100;   axis(v);   
      subplot(2,2,2),semilogx(w,tp2,'b','LineWidth',2),grid    
      title('H(s) = 1/s')
      xlabel('frequency (rad/sec)'),ylabel('angle'),axis(v);        
      subplot(2,2,3),semilogx(w,tp3,'g','LineWidth',2),grid 
      title('H(s) = 1/(s+1)')
      xlabel('frequency (rad/sec)'),ylabel('angle'),axis(v);         
      subplot(2,2,4),semilogx(w,tp,'k','LineWidth',2),grid
      title('H(s) = (s+0.1)/[s(s+1)]')
      xlabel('frequency (rad/sec)'),ylabel('angle'),axis(v);         
%                                       
%   Nichols plots (log-mag vs angle)                                
      nfig = nfig+1;   figure(nfig)                            
      subplot(1,2,1),plot(tp,tm,'r.','LineWidth',2),grid
      title('H(s) = (s+0.1)/[s(s+1)]')
      xlabel('angle'),ylabel('|H(jw)| in db')
      gtext('low freq'),gtext('high freq')                        
%                                           
%   Nyquist plots (Im vs Re)                        
      b =  [1 .1]; a = [1 1 0];                  
      h = freqs(b,a,w);   hr = real(h);   hi = imag(h);              
      subplot(1,2,2),plot(hr,hi,'b.','LineWidth',2),grid  
      title('H(s) = (s+0.1)/[s(s+1)]')
      xlabel('Real'),ylabel('Imag')                  
      gtext('low freq'),gtext('high freq')            
%
%   end of example