%  BODE PLOT OF Gp(s)

%  FILENAME: BODEGP.M
%  DATE: October 1994
%  VERSION: 02
%  COPYRIGHT 1993, 1994 M.Gretzinger & T.Marlin
%            McMaster University   

% This m-file is called from STABLE.M. It sets up a menu to enable the user
%       to specify the start and end frequency and create a bode plot for
%       the process only, Gp(s), using the parameters specified by the 
%       process data.

while 1
  clc
  disp('*************************************************************')
  disp('*          BODE PLOT OF PROCESS, Gp(s)                      *')
  disp('*                                                           *')
  disp('*************************************************************')
  disp('')
  disp('SELECT THE APPROPRIATE MENU ITEM')
  disp('')
  disp('MODIFY...')
  disp('                            PRESENT VALUES')
  fprintf('1) Lowest Frequency         %4.5g\n',w1)
  fprintf('2) Highest Frequency        %4.2f\n',w2)
  disp('3) Create Bode plot')
  disp('4) Return to main menu')
  disp('')
        
  m5 = input('Enter the desired selection:   ');
  
  if isempty(m5), break, end
      
  m5 = round(m5);
  if (m5<1) | (m5>4)
      disp('')
      disp('Not a valid selection')
      disp('')
      disp('Press ENTER to continue')
      pause;
  end    
  
  if m5 == 1
       junkw1 = input('Enter value for lowest frequency:  ') ;
       if isempty(junkw1)
           junkw1 = w1 ;
       end
       if junkw1 >= w2
          disp ('entry illegal, > highest')
          disp ('press enter to continue')
          junkw1 = w1;
          pause
       end
       w1 = junkw1;

  elseif m5 == 2
       junkw2 = input('Enter value for highest frequency:  ') ;
       if isempty(junkw2)
           junkw2 = w2  ;
       end
       if junkw2 <= w1
          disp ('entry illegal, < lowest')
          disp ('press enter to continue')
          junkw2 = w2;
          pause
       end
       w2 = junkw2;

  elseif m5 == 3
     
     % Block to calculate Gp, phase angle and amplitude ratio over the 
     % chosen frequency range. A vector of frequencies, called w, is defined
     % using logspace. The vector of 500 frequencies are evenly spaced in the 
     % log domain between w1 and w2. Each element in the vector is then 
     % multiplied by jj to create the vector s. The vector Gp 
     % is then calculated. Elementwise multiplication(.*), division(./) and 
     % powers(.^) are used in order that the expression for Gp
     % is carried out once for each element of the vectors making up
     % the expression.

      jj = sqrt(-1) ;
      w = logspace(log10(w1), log10(w2), 500) ;
      s = jj*w ;

      Gp = abs(kp).*exp(-theta*s).*(taulead*s + 1)...
           ./((tau1*s + 1).*(tau2*s + 1).*(tau3^2*s.^2 + 2*xi*tau3*s + 1)) ;
      
     % Calculate the vector of amplitude ratios. The function ABS calculates
     % the magnitude of each of the complex numbers in the vector Gp.
      ar = abs(Gp) ;           
      
     % Calculate the vector of phase angles. The function ANGLE calculates
     % the phase angles of each of the complex numbers in the vector Gp.
      phase = angle(Gp) ;     
      phase = unwrap(phase) ;      % Unwraps the phase so it is continuous
                                   % above or below +/- 180 degree boundaries
     
     % Create Bode plot
      clf; axis ('auto')
      subplot(2,1,1), loglog(w, ar)
      title('PROCESS (GP) BODE PLOT')
      xlabel('Frequency, w  (rad/time)')
      ylabel('Amplitude Ratio')

      subplot(212), semilogx(w, phase*180/pi)
      xlabel('Frequency, w  (rad/time)')
      ylabel('Phase Angle (degrees)')

      figure(1)
      pause
      close all
      clear w s Gp ar phase
  
  else
      break
  end        % if m5
  
end    % while 1 (line 12)
