%  Program "gpmargin.m" to demonstrate the basic Matlab functions 
%    that can be employed for determination of gain and phase margins
% Construct Loop Transfer Function

% Parameter values of model
%L = 0.01 ; %inertance in units of cm H2O s2/L
L = 0.01 ; %inertance in units of cm H2O s2/L
C = 0.1 ; %compliance in units of L/cm H2O
R = 1.0 ; %resistance in units of cm H2O s/L
lambda =  1

; %open-loop case

% set up transfer function
num = [1];
den = [L*C  (R*C)  lambda];
Hs = tf(num,den);


% Parameter values of second model
%num=[1];
%den=[1 3 2 0];
%Hs = tf(num, den);

f = 0.01:0.01:10;
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 (Bode diagrams)
input ('Display Gain & Phase Margins: Hit <return> to continue >>')
margin(mag,pha,w);
[GM,PM,wcM,wcP]=margin(mag,pha,w);