% descfcn.m - describing function plots
% RMM, 3 Jan 08

aminit;

% Linear system
[num, den] = pade(2, 2);
L = tf([4], [1 2 1]) * tf(num, den);
w = logspace(-3, 4, 1000);

% Nonlinearity
a = 0:0.01:1;
Na = atan(4*a);

% Generate the plot axes
figure(1); clf; subplot(221); 
ampolezero(tf([1], [1]), [-5 5 -5 5]); daspect(daspect); hold on;

% Plot the describing function curve
plot(-1.25-3.75*a, 0.9*Na - 1.5, 'Linewidth', AM_data_linewidth, ...
  'Color', [0 0.5 0]);
amtext(-5.5, -1, '-1/N(a)');

% Plot the frequency response 
H = freqresp(L, w); H = H(:);
plot(real(H), imag(H), 'b');
amtext(3.5, -2.5, 'G(iw)');

% Add a circle at the intersection (manually determined)
H = freqresp(L, .775);
plot(real(H), imag(H), 'ro', 'markersize', AM_bode_markersize);

% set up the color
color = 'b';
if (AMPRINT_FLAG==2) color = 'k'; end;

% Put some arrows on the diagrams
H1 = freqresp(L, .6); H2 = freqresp(L, .601);
arrow([real(H1) imag(H1)], [real(H2) imag(H2)], ... 
  AM_normal_arrowsize, 'FaceColor', color, 'EdgeColor', color);

% set up the color
color = [0 0.5 0];
if (AMPRINT_FLAG==2) color = 'k'; end;

% Put some arrows on the diagrams
mid = 18;
H1 = -1.25-3.75*a(mid) + i*(0.9*Na(mid) - 1.5);
H2 = -1.25-3.75*a(mid+1) + i*(0.9*Na(mid+1) - 1.5);
arrow([real(H1) imag(H1)], [real(H2) imag(H2)], ... 
  AM_small_arrowsize, 'FaceColor', color, 'EdgeColor', color);

amprint('desfcn-nyquist.eps');