function [h1,h2,h3,h4,results]=gof(P,C,w,AM_data_linewidth,str,scale)
%Plotting GOF				
if nargin < 4
    AM_data_linewidth = 1.5
end
if nargin < 5
    str = 'b-';
end
if nargin < 6 
    scale = [min(w) max(w) 1e-3 1e2];
end
Pw = squeeze(freqresp(P,w));
Cw = squeeze(freqresp(C,w));

L=Pw.*Cw;
S=1./(1+L);
PS=Pw.*S;
CS=Cw.*S;
T=L.*S;

% Gain crossover frequency
[xxx,iii]=min(abs(abs(L)-1));
wgc=w(iii);
% Maximum sensititivites
[ms,iii]=max(abs(S));
ws=w(iii);
[mt,jjj]=max(abs(T));
wt=w(jjj);
[mcs,jjj]=max(abs(CS));
wcs=w(jjj);
[mps,jjj]=max(abs(PS));
wps=w(jjj);nnn=size(w);nnn=nnn(2);
nn1=floor(nnn/2);
[xx,jjj]=min(abs(abs(CS(nn1:nnn))-1));
if jjj+nn1 < length(w)
    wcs1=w(jjj+nn1);
else
    wcs1 = inf;
end
results = [wgc ms mt mps mcs ws wt wcs wps wcs1];

%Plot results
subplot(224);
ylabel('T');
loglog(w,abs(T),str,'Linewidth',AM_data_linewidth);
h1 = gcf;
%set(gca,'Xtick',[0.01 0.1 1 10 100]);
%set(gca,'Ytick',[0.01 0.1 1 10 100]);
axis(scale);grid on;hold on;
xlabel('w');ylabel('T');hold on
subplot(222);
h2 = loglog(w,abs(PS),str,'Linewidth',AM_data_linewidth);
%set(gca,'Xtick',[0.01 1 100]);
%set(gca,'Ytick',[0.01 0.1 1 10 100]);
axis(scale);grid on;hold on;
ylabel('PS');hold on
subplot(223);
h3 = loglog(w,abs(CS),str,'Linewidth',AM_data_linewidth);
%set(gca,'Xtick',[0.01 1 100]);
%set(gca,'Ytick',[0.01 0.1 1 10 100]);
axis(scale);grid on;hold on;
xlabel('w');ylabel('CS');
subplot(221);
h4 = loglog(w,abs(S),str,'Linewidth',AM_data_linewidth);
%set(gca,'Xtick',[0.01 1 100]);
%set(gca,'Ytick',[0.01 0.1 1 10 100]);
axis(scale);grid on;hold on;
ylabel('S');




