% bodesketch.m - example of how to sketch a bode plot
% RMM, 29 Jul 06

aminit;					% initialize settings

% Define the transfer function we will use
a = 0.2;				% first pole
b = 1;					% zero
omega0 = 5; zeta = 0.05;		% second order poles
k = 10 * a / b * omega0^2;		% set gain to 10

% Define the actual transfer funtion
s = zpk('s');
H = k * (s + b) / ((s + a) * (s^2 + 2 * zeta * omega0 * s + omega0^2));

% Define the asymptotes, starting with magnitude (easy)
asymg = [  10 10  2 2 2/(20^2)]; 
asyfr = [0.01 0.2 1 5 100];

% For phase, we generate two segments, each fairly complex
asyfr1(1) = 0.01;			% initial phase
asyph1(1) = 0;

asyfr1(2) = 0.02;			% pole starts to have an effect
asyph1(2) = 0;

asyfr1(3) = 0.1;			% zero starts to cancel pole
asyph1(3) = asyph1(2) - 45 * (log10(asyfr1(3)) - log10(asyfr1(2)));

asyfr1(4) = 2;				% pole stops having an effect
asyph1(4) = asyph1(3);			% only zero remains

asyfr1(5) = 5;				% zero alone has effect
asyph1(5) = asyph1(4) + 45 * (log10(asyfr1(5)) - log10(asyfr1(4)));

% Now start the second segment
asyfr2(1) = 5;				% second order system kicks in
asyph2(1) = asyph1(5) - 180;

asyfr2(2) = 10;				% zero kicks out
asyph2(2) = -180;

asyfr2(3) = 100;			% tail end
asyph2(3) = -180;

% Extract out the Bode plot
[mag, phase, W] = bode(H);
mag = mag(:); phase = phase(:);

% Generate an initial bode plot that we can draw on top of
figure(1); clf; axes('position', [0.13 0.4096 0.7750 0.5]);
[magh, phaseh] = ambode(H);

% Replot with the asymptotes, starting with the magnitude
subplot(magh); hold on;

h = loglog(asyfr, asymg, 'r--');
set(h, 'Linewidth', AM_thick_linewidth, 'MarkerSize', AM_data_markersize);

h = loglog(W, mag, 'b-');
set(h, 'Linewidth', AM_bode_linewidth);

% Add on circles at the breakpoints
h = loglog(asyfr, asymg, 'ro', 'LineWidth', AM_ref_linewidth, ...
  'MarkerSize', AM_data_markersize);

amaxis([0.01 100 0.01 100]);		% set axis limits
set(gca, 'Ytick', [0.01 1 100]);	% set up tick marks
grid on;				% turn on grid marks
ylabel('mag');

% Now annotate the plot to show the features we want
amtext(a, 0.1, 's=a');
amtext(b, 0.1, 's=b');
amtext(omega0, 0.1, 's=omega0');

lgh = legend('Exact', 'Approx');
legend(lgh, 'boxoff');

% Now plot the phase
subplot(phaseh); hold on;

h = semilogx(asyfr1, asyph1, 'r--', asyfr2, asyph2, 'r--');
set(h, 'Linewidth', AM_thick_linewidth);

h = semilogx(W, phase, 'b-', 'Linewidth', AM_bode_linewidth);

% Add on circles at the breakpoints
semilogx(asyfr1, asyph1, 'ro', asyfr2, asyph2, 'ro', ...
  'LineWidth', AM_ref_linewidth, 'MarkerSize', AM_data_markersize);

set(gca, 'Ytick', [-180 -90 0]);	% set up tick marks
axis([0.01 100 -200 20]);		% set axis limits
grid on;				% turn on grid marks
ylabel('phase');
xlabel('freq');

% Now annotate the plot to show the features we want
text(a/10, -90, 's=a/10');
text(b/10, -90, 's=b/10');
text(10*a, -90, 's=10a');
% text(omega0, -90, 's=omega0');
text(10*b, -135, 's=10b');

amprint('bodesketch.eps');
