function [magh, phaseh] = ambode(sys, varargin)
global AM_data_linewidth AM_ref_linewidth AM_bode_linewidth
global AM_bode_magloc AM_bode_phaseloc;

% Generate the Bode plot data using the MATLAB bode command
[mag, phase, omega] = bode(sys, varargin{:});
phase = unwrap(pi*phase/180) * 180/pi;

% Get the dimensions of the current axis, which will divide up
curpos = get(gca, 'Position');
height = curpos(4) * 0.45;
magbot = curpos(2) + curpos(4) * 0.55;
phabot = curpos(2);

% Turn off the background axes
axis off;

% Magnitude plot
magh = axes('position', [curpos(1), magbot, curpos(3), height]);
h = loglog(omega, squeeze(mag));
set(h, 'LineWidth', AM_bode_linewidth);

% Reset the axes to something reasonable
omega_min = min(omega);
omega_max = max(omega);
mag_min = 10^floor(log10(min(mag)));
mag_max = 10^ceil(log10(max(mag)));

set(gca, 'XTickLabel', {});
axis([omega_min, omega_max, mag_min, mag_max]);
ylabel('Gain');

% Set the tick labels at reasonable powers of 10
mag_incr = ceil((log10(mag_max) - log10(mag_min))/3);
set(gca, 'YTick', 10.^(log10(mag_min):mag_incr:log10(mag_max)));

% Set the frequency labels at reasonable powers of 10
omega_incr = ceil((log10(omega_max) - log10(omega_min))/4);
set(gca, 'XTick', 10.^(log10(omega_min):omega_incr:log10(omega_max)));

% Phase plot
phaseh = axes('position', [curpos(1), phabot, curpos(3), height]);
h = semilogx(omega, squeeze(phase));
set(h, 'LineWidth', AM_bode_linewidth);
ylabel('Phase [deg]'); xlabel('Frequency [rad/s]');

% Set the axis limits to something sensible
phase_min = floor(min(phase)/90) * 90; phase_max = ceil(max(phase)/90) * 90;
if (phase_max < 0) phase_max = 0; end;
if (phase_min > 0) phase_min = 0; end;
if (phase_min == phase_max) 
  phase_max = phase_max + 45;
  phase_min = phase_min - 45;
end
axis([min(omega), max(omega), phase_min, phase_max]);

% Set the tick labels to be at 90 degree increments
set(gca, 'YTick', [phase_min:90:phase_max]);
set(gca, 'XTick', 10.^(log10(omega_min):omega_incr:log10(omega_max)));