% Program "gpmargin.m" to demonstrate the basic Matlab functions 
% that can be employed for determination of gain and phase margins
% Construct Loop Transfer Function of the Pupil model
num=[0.16];					% DC Gain 
den=[0.001 0.03 0.3 1]; 	% Coefficients of the expanded denominator
Hs = tf(num, den);

f = 0.01:0.01:10;			% Compute Bode from 0.01 to 10 Hz
w  = 2*pi*f;
% Compute and plot Bode diagrams
input ('Display Bode Plots: Hit <return> to continue >>')
figure
bode(Hs,w)
[mag,pha]=bode(Hs,w); %save mag & phase of Bode plot to workspace

% Compute and plot Nichols chart
input ('Display Nichols Chart: Hit <return> to continue >>')
figure
nichols(Hs,w)
[mag,pha]=nichols(Hs,w); %save mag & phase of Nichols plot to workspace

% Compute and plot Nyquist diagram
input ('Display Nyquist Diagram: Hit <return> to continue >>')
figure
nyquist(Hs,w)
[nyr,nyi]=nyquist(Hs,w); %save real & imag parts of Nyquist plot to workspace

% Compute and plot gain and phase margins
input ('Display Gain & Phase Margins: Hit <return> to continue >>')
margin(mag,pha,w);
[GM,PM,wcM,wcP]=margin(mag,pha,w);
