function [mag, phase] = bode_d(num, den, W, delay)
% [mag, phase] = bode_d(num, den, W, delay)
%
% A poor man's bode for use on computers without the Control Systems
% Toolbox.  This version of bode can also add a system delay, which the
% MATLAB version cannot.  The inputs are the same as for the MATLAB 
% version of bode except 'num' and 'den' must be coefficients of
% decreasing powers of 's' and the 'delay' value adds a system delay in 
% seconds.
%
% Anil Aswani: 12-4-04
% regards to Prof. Grizzle for the plotting commands 

num_power = fliplr(0:(length(num) - 1));
den_power = fliplr(0:(length(den) - 1));

num_W = repmat((j*W), length(num), 1).^repmat(num_power', 1, length(W));
den_W = repmat((j*W), length(den), 1).^repmat(den_power', 1, length(W));

G_num = num*num_W;
G_den = den*den_W;

G = G_num./G_den;

mag = abs(G);
phase = (unwrap(angle(G)) - delay*W)*180/pi;

subplot(2,1,1);
h=semilogx(W,20*log10(mag));
ylabel('Magnitude (dB)');
grid on;
set(h,'LineWidth',2);

subplot(2,1,2);
h=semilogx(W,phase);
xlabel('Frequency (rad/sec)');
ylabel('Phase (deg)');
grid on;
set(h,'LineWidth',2);