% PlotFreqResponseMagnitude.m
%
% By: Neil E Cotter
% 11/09/2020
%
% Script file to plot magnitude frequency response curve with labels: |H| vs w (i.e., omega).
%
% User must set input variables (or defaults below are used):
%
% k = overall gain factor
% pole_array = array of pole locations as complex numbers
% pole_array = array of zero locations as complex numbers (set to [] if no zeros used)
% max_w = max angular frequency to use for plot

if ~exist('pole_array')
  pole_array = [-3+4j,-3-4j]*1e4;
end

if ~exist('zero_array')
  zero_array = [];
end

if ~exist('k')
  k = prod(pole_array);
end

if ~exist('max_w')
  max_w = 20e4;
end


num_poles = length(pole_array);
num_zeros = length(zero_array);

hold off  % Start a new plot.

% Close figure 2, if it exists.
if ishandle(2)
  close(2)
end

% Create a new figure 2.
figure(2)

% Reset all figure properties, except Position, Units, PaperPosition, and
%  PaperUnits.
clf('reset')


hold on

% Set defaults of fonts and line sizes for axes and labels.
set(0,'DefaultAxesFontSize',18)
set(0,'DefaultLineLineWidth',2)
set(0,'DefaultTextFontSize',18)

w_vec = 0:max_w/1e3:(max_w*0.95);
jw_vec = w_vec * j;

num_w = length(w_vec);

gain_vec = zeros(1,num_w);

for w_index = 1:length(jw_vec)
  jw = jw_vec(w_index);
  gain_vec(w_index) = k;
  for zero_index = 1:num_zeros
    gain_vec(w_index) = gain_vec(w_index)*abs(jw - zero_array(zero_index));
  end
  for pole_index = 1:num_poles
    gain_vec(w_index) = gain_vec(w_index)/abs(jw - pole_array(pole_index));
  end
end

gain_max = max(gain_vec);

xmin = -max_w/20;  % Freq axis starts just left of zero.
xmax = max_w;
ymin = -0.1*gain_max;  % Gain axis starts just below zero.
ymax = 1.2*gain_max;

axis_spec = [xmin, xmax, ymin, ymax];  % Put the values in an array.
axis(axis_spec)  % Set the axes limits.


% Draw freq and gain axes.
plot([xmin,xmax*0.95],[0,0],'k')
plot([0,0],[ymin,ymax*0.9],'k')
text(xmax*(1 - 0.75/20),0,'$w$','FontUnits','points','FontWeight','normal', ...
'interpreter','latex','FontSize',18,'FontName','Times')
text(-0.5e4,ymax*0.95,'$|H|$','FontUnits','points','FontWeight','normal', ...
'interpreter','latex','FontSize',18,'FontName','Times')

%plot(w_vec,gain_vec,'b-')
plot(w_vec,gain_vec,'b-')

hold off

