fndwmal.m

pk=2;                           % plant Bode gain
pnum=[0 0 0 1];                 % coeffs of the plant numerator
pden1=[1 0];                    % plant denominator factor
pden2=[1/3 1];                  % another plant denominator factor
pden3=[1/4 1];                  % last plant denominator factor

pden=conv( pden1, conv(pden2, pden3));

ck=2;

num=ck*pk*pnum;             % numerator with new Bode gain
den=pden;                   % denominator
% --- compute phase margin

[mag, pha, w]=bode(num, den);
[Gm, Pm, Wcg,Wcp]=margin(mag,pha,w);

% --- display results and prompt user for lead angle phimax

str1   =' Phase margin = ';
str    = [str1 num2str(Pm)] % display the string

phimax = input(' Please enter phimax  --- > ');
% --- compute al from phimax

leadrad = pi/180*phimax; % phimax in radians
sn  = sin(leadrad);      % its sinus
al  = sqrt ((1+sn)/(1-sn));
% --- wm is the frequency where the maximum lead angle is produced
%     and is certainly larger than Wcp.

% We are going to plot the phase margin versus wm
nwm = 50;                % number of wm
awm = linspace(Wcp, 10*Wcp, 50); % array of wm-s.

apm = zeros(awm);        % make space for array of phase margin

% --- next loop prepares the data for the graph of phase
%     margin as a function of wm

for l=1:nwm,
    wm = awm(l);
    cnum = [1/(wm/al) 1];% numerator 1 + s/Z (Z=wm/al)
    cden = [1/(wm*al) 1];% denominator + s/P (P=wm*al)
    num  = ck*pk*conv(cnum, pnum);
    den  =    conv(cden, pden);
    % --- Compute and save phase margin for this design

    [mag, pha, w]=bode(num, den); [Gm, Pm, Wcg,Wcp]=margin(mag,pha,w);
    apm(l)= Pm;          % array of phase margins, relative to wm
end
plot(awm, apm);
xlabel('wn (rad/sec)');
ylabel('Phase margin (degrees)');
grid

str = ['phi-max = ' num2str(phimax)];
title(str);
pause;
% --- inform designer of maximum phase margin and prompt him for wm

[maxpm, ipm]=max(apm);
str1 =['Phase margin = ' num2str(maxpm)];
str2 =['  at wm = ' num2str(awm(ipm))];
str = [str1 str2]

wm = input('Please enter wm --- > ');

% --- test result


cnum = [1/(wm/al) 1];  % numerator 1 + s/Z (Z=wm/al)
cden = [1/(wm*al) 1];  % denominator + s/P (P=wm*al)
num  = ck*conv(cnum, pnum);
den  =    conv(cden, pden);

[mag,pha,w] = bode(num,den); [Gm, Pm, Wcg, Wcp]=margin(mag,pha,w);
% --- Display result and bode diagram
str  =['Phase margin = ' num2str(maxpm)]
pause;
bode(num,den);
title(' Bode diagram - open loop ');