% desfcn.m
% Describing function and its negative inverse for relay with hysteresis
%   kja 070710

aminit;

figure(1); clf; xsubplot(331); hold on;
amaxis([-1 1 -1.2 1.2]);		% set axes so that units scale
daspect(daspect);

% Add axes (no labels)
arrow([-1 0], [1 0], AM_normal_arrowsize, 'LineWidth', AM_ref_linewidth);
arrow([0 -1.2], [0 1.3], AM_normal_arrowsize, 'LineWidth', AM_ref_linewidth);
axis off;

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

% Plot the curve
plot([-0.8 0.3 0.3 0.8], [-0.7 -0.7 0.7 0.7], 'b-', ...
  [-0.8 -0.3 -0.3 0.8], [-0.7 -0.7 0.7 0.7], 'b-', ...
  'LineWidth', AM_data_linewidth);
arrow([-0.3 0.11], [-0.3, 0.1], AM_normal_arrowsize, ...
  'FaceColor', color, 'EdgeColor', color);
arrow([0.3 -0.11], [0.3, -0.1], AM_normal_arrowsize, ...
  'FaceColor', color, 'EdgeColor', color);

% Add axis labels
amtext(-0.2, 1.3, '{\ity}');
amtext(0.9, -0.2, '{\itu}');

amtext(0.4, 0.2, '{\itc}');
amtext(0.1, 0.95, '{\itb}');

amprint('desfcn-hysteresis.eps');

%   Describing function
b=1;c=0.5;
%a=logspace(-2,2);
a=0.5:0.001:100;
N=4*b./(pi*a).*(sqrt(1-(c./a).^2)-i*c./a);
Ninv=1./N;

%   Linear system
w=logspace(-2,2,1000);s=i*w;
G=2./(s+1).^4;
[xxx,iii]=min(abs(imag(G(300:1000))+pi*c/4/b));

%   Plotting
figure(2); clf; xsubplot(331);
ampolezero(tf([1], [1]), [-1 3 -2 2]); daspect(daspect); hold on;

plot(real(G), imag(G), '--', -real(Ninv), -imag(Ninv), '-', ...
  'LineWidth', AM_data_linewidth);
plot(real(G(299+iii)), imag(G(299+iii)), 'ro', ...
  'MarkerSize', AM_bode_markersize);
set(gca,'Xtick',[-2 -1 0 1 2]);

amprint('desfcn-relay.eps');

wosc=w(299+iii);
[xxx,jjj]=min(abs(real(Ninv)+real(G(299+iii))))
aosc=a(jjj);

