% This Matlab script file makes the two plots shown in class for
% the flexible satelite example.  Save this file as 'example.m'
% and type example at the matlab prompt.  Tested with matlab version
% 4.2 and 5.0.
% Date: November 1998
% Author: S A Bortoff

w=logspace(-5,1,200);                % range of frequencies

nump=0.1;                            % The plant numerator is defined here.
denp=[1 0];                          % But the plant denoninator is
denp=conv(denp,[1/0.01 1]);          % put together using the
denp=conv(denp,[1/10 1]);            % conv command.

numc=[100 1];                        % The controller numerator is 
numc=conv(numc,[1/0.02 1]);          % also put together using conv.
numc=10*numc;                        % The gain is 10.
denc=[1000 1];                       % The controller denominator.
denc=conv(denc,[1/0.2 1]);

num=conv(numc,nump);                 % The numerator and denominator of
den=conv(denc,denp);                 % C P (s).

num=[0 0 0 num];                     % This is so length(num)=length(den)

%
% Sensitivity function...
%

nums=den;
dens=num+den;                        % Works because length(num)=length(den).

% 
% Complementary sensitivity function...
% 

numt=num;
dent=num+den;

% 
% Uncertainty bound...
%

lm=(1+(w./0.03).*(w./0.03))/30;      % Note the .* which is the element - by -
lminv=1./lm;                         % element multiplication.
                                     % Thus, lm and lminv are vectors of the
                                     % same length as w.

% 
% Performance specification...
%

for j=1:length(w),
 if ( w(j) < 7.35*(10)^(-5) ),
  p(j)=12000;
 else
  p(j) = 10^(-8);                    % 10^(-8) is approximately 0, used because
 end;                                % if you use 0, loglog complains.
end;
pinv=1./p;                           % Again, note the ./ 

%
% At this point, lm, lminv, p and pinv are defined.  These are vectors.  
% of the same length as w.  Next, we are going to make the bode plots... 
%

[mags,phases]=bode(nums,dens,w);     % mag and phase of sensitivity
[magt,phaset]=bode(numt,dent,w);     % mag and phase of comp. sensitivity

%
% Make the plots.  Note use of hold command.  Change colors as desired.
% See help plot and help loglog.  First plot T and S...
%

figure(1);
hold off;
loglog(w,mags,'c','red');               % First plot the sensitivity function
text(10^0, 3, '|S(jw)|','c','red');     % along with a label.
hold on;
loglog(w,magt,'c','cyan');              % Plot the comp. sensitivity function
text(10^(-3), 3, '|T(jw)|','c','cyan'); % along with a label.

%
% Next we plot the weights lm(w) inverse and p(w) inverse
%

loglog(w,lminv,'--','c','red');
text(10^(-1),10^1,'lminv(w)','c','red');
loglog(w,pinv,'--','c','cyan');
text(1.2*10^(-5),4*10^(-4),'pinv(w)','c','cyan');

%
% Axes labels and title - important !
%

title('|S(jw)|, |T(jw)|, pinv(w), lminv(w) for Satelite Example')
xlabel('Frequency (log scale)');
ylabel('Magnitude')

%
% Finally, fix the axes...
%

axis([w(1) w(length(w)) 10^(-5) 10^2]); 
grid;

% If you want to print this out, save it as a postscript file
% using the print command (help print).

temp=mags.*p'+magt.*lm';           % Note the use of .* and ' (transpose)
figure(2);                         % Put plot on figure 2.
loglog(w,temp,'c','b','linewidth',2);
title('|S(jw)|*p(w) + |T(jw)| * lm(w) for Satelite Example')
xlabel('Frequency');
ylabel('Magnitude')
grid;                              % put grid lines

%
% To print, type help print :-)
%