%  BODE PLOT OF GOL(s)

%  FILENAME: BODEGOL.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 of 
%       GOL(s), using the parameters specified by the process and tuning data.

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

  elseif m4 == 2
       junkwend = input('Enter value for highest frequency:  ') ;
       if isempty(junkwend)
           junkwend = wend ;
       end
       if junkwend <= wstart
          junkwend = wend;
          disp ('illegal entry, < lowest')
          disp ('press enter to continue')
          pause
       end
       wend = junkwend;


  elseif m4 == 3

     % Block to calculate GOL, 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 wstart and wend. Each element in the vector is then 
     % multiplied by jj to create the vector s. The vectors Gp, Gc and GOL 
     % are then calculated. Elementwise multiplication(.*), division(./) and 
     % powers(.^) are used in order that the expressions for Gp, Gc and GOL 
     % are carried out once for each element in the vectors making up
     % the expressions.
     
      jj = sqrt(-1) ;
      w = logspace(log10(wstart), log10(wend), 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)) ;
      
     % If block to prevent division by zero which would otherwise occur when
     % there is no integral mode.
      if ti == 0
         Gc = abs(kc)*(1 + td*s) ;
      else
         Gc = abs(kc)*(1 + (ti*s).^(-1) + td*s) ;
      end

      GOL = Gp.*Gc ;     % Calculate "open loop" frequency response
      
     % Calculate the vector of amplitude ratios. The function ABS calculates
     % the magnitude of each of the complex numbers in the vector GOL.
      ar = abs(GOL) ;

     % Calculate the vector of phase angles. The function ANGLE calculates
     % the phase angles of each of the complex numbers in the vector GOL.
      phase = angle(GOL) ;
      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('BODE PLOT OF GOL')
      xlabel('Frequency, w  (rad/time)')
      ylabel('Amplitude Ratio')

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

      figure (1)
      pause
      close all

      clear w s Gp Gc GOL ar phase
  
  else
      break
  end        % if m3
  
end    % while 1 (line 10)