%  **** EXAMPLE 4.16 SERIES REACTORS FREQ. RESPONSE  ***
%  ****    FROM MARLIN 1994
%  COPYRIGHT 1994 McMASTER UNIVERSITY
%  WRITTEN BY T. MARLIN
%  FILENAME: EXAM4_16.M
%  LAST UPDATED: NOVEMBER 1994
%  VERSION:  1.0
%  DESCRIPTION OF FILE:
%  this MATLAB M-file calculates and plots the
%  frequency response for two series chemical reactors
%  in example 4.16 
%  ******************************************************

clear;                  % clear workspace
clf  ;                  % clear graph
axis('auto');           % auto scale
clc ;                   %  clear the screen
%close all

%  *******************************************
%   PROCESS DATA
%  *******************************************

v = 1.05    ;
f = 0.085   ;
k = 0.040   ;
cas= 0.414  ;   % initial steady-state outlet
A = 0.925   ;   % amplitude of the sine

%  ********************************************
%  parameters in the linear model
%  ********************************************

kp = (f/(f+v*k))^2  ;
tau= (v/(f+v*k))    ;

%  *********************************************
%    SIMULATION DATA
%    note that frequency range is equally spaced
%    in log scale (not linear)
%  *********************************************

wstart = .001  ;   % the smallest frequency
wend   = 1000   ;   %  the highest frequency
wtimes = 1200  ;   %  number of points in frequency range

omega = logspace ( log10(wstart), log10(wend), wtimes);
jj = sqrt(-1) ;   %  define the complex variable

%  *********************************************
%  put calculations here
%  *********************************************

for kk = 1:wtimes

     % **********************
     % this is calculation using
     %  equations from the text
     % ***********************
  
     ar(kk) = kp/(1 + omega(kk)^2 * tau^2 );
     phase(kk) = atan ( -2*omega(kk)*tau/ ...
                (1-tau^2 * omega(kk)^2))   ;

     % **************************
     %  this is calculation using
     %  complex algebra , s = jj*omega
     % **************************
     
     s = jj*omega(kk)          ;
     G(kk) = kp /( ( tau*s + 1)^2) ; 
     AR(kk) = abs (G(kk));         ;
     PHASE(kk) = angle (G(kk))     ;

end  % for cnt

%  convert from radians to degrees for Bode plot

     phase = phase *360/(2*pi);% function unwrap
                               % corrects for "jumps"
                               % at +- n*90 degrees
     PHASE = unwrap(PHASE) *360/(2*pi);
% ***************************************************
%  plot the results in Bode plot
% ***************************************************

subplot (2,1,1), loglog(omega,ar, omega, AR, '--')
xlabel ('frequency, rad/time ')
ylabel ('amplitude ratio')
title('Example 4.16, solid line is eqn(4-107), dashed from complex algebra')

subplot (2,1,2), semilogx(omega,phase, omega, PHASE,'--')
xlabel ('frequency, rad/time ')
ylabel ('phase angle,  degrees')

figure (1) ; 
pause ; 
close all;