% Plot_Phase_vs_log_w.m
%
% By: Neil E Cotter
% 11/13/2020
%
% Script file to plot magnitude frequency response |H| in dB vs log omega.
%
% User must set input variables (or defaults below are used):
%
% k = overall gain factor
% pole_array = array of pole locations as complex numbers
% zero_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(abs(pole_array));
end

if ~exist('min_w')
  min_w = 1;
end

if ~exist('max_w')
  max_w = 1e8;
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)

log_w_vec = log10(min_w):log10(max_w)/1e3:log10(max_w)*0.89;
jw_vec = power(10,log_w_vec) * j;

num_w = length(log_w_vec);

phase_vec = zeros(1,num_w);

for w_index = 1:length(jw_vec)
  jw = jw_vec(w_index);
%  phase_vec(w_index) = 0;
  for zero_index = 1:num_zeros
    dist_vec = jw - zero_array(zero_index);
    phase_vec(w_index) = phase_vec(w_index) + 180/pi * atan2(imag(dist_vec),real(dist_vec));
  end
  for pole_index = 1:num_poles
    dist_vec = jw - pole_array(pole_index);
    phase_vec(w_index) = phase_vec(w_index) - 180/pi * atan2(imag(dist_vec),real(dist_vec));
  end
end

phase_min = min(phase_vec);
phase_max = max(phase_vec);

% Set plot box boundaries.
xmin = log10(min_w);  % Freq axis starts just left of zero.
xmax = log10(max_w);
ymin = floor(phase_min);
ymax = ceil(phase_max + 15);

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


% Draw freq and gain axes.
loglog([xmin,xmax*0.89],[0,0],'k')
loglog([xmin,xmin],[ymin,ymin+(ymax-ymin)*0.95],'k')
text(xmin+(xmax-xmin)*(1 - 0.098),0,'$log\; w$','FontUnits','points','FontWeight','normal', ...
'interpreter','latex','FontSize',18,'FontName','Times')  % Uses LateX string.
text(-0.17,ymin+(ymax-ymin)*0.975,'$\angle H\;$','FontUnits','points','FontWeight','normal', ...
'interpreter','latex','FontSize',18,'FontName','Times')  % Uses LateX string.

%plot(w_vec,gain_vec,'b-')
loglog(log_w_vec,phase_vec,'b-')

% Put tick marks at powers of 10 for freq.
xticks_vec = ceil(xmin):floor(xmax);
xticks(xticks_vec)

xticklabels_vec = [];
for decade_index = 1:length(xticks)
  log_freq = xticks_vec(decade_index);

  if log_freq == 0
    cell_tick_str = {'1'};
  elseif log_freq == 1
    cell_tick_str = {'10'};
  elseif log_freq == 2
    cell_tick_str = {'100'};
  elseif log_freq == 3
    cell_tick_str = {'1k'};
  elseif log_freq == 4
    cell_tick_str = {'10k'};
  elseif log_freq == 5
    cell_tick_str = {'100k'};
  elseif log_freq == 6
    cell_tick_str = {'1M'};
  elseif log_freq == 7
    cell_tick_str = {'10M'};
  elseif log_freq == 8
    cell_tick_str = {'100M'};
  elseif log_freq == 9
    cell_tick_str = {'1G'};
  elseif log_freq == 10
    cell_tick_str = {'10G'};
  elseif log_freq == 11
    cell_tick_str = {'100G'};
  end

  xticklabels_vec = [xticklabels_vec, cell_tick_str];
end

xticklabels(xticklabels_vec)

% Obsolete tick mark code.
%XTick = [1,10,100,1000,10000,100000,1000000];
%xticks([0,1,2,3,4,5,6]);
%xticklabels({'1','10','100','1k','10k','100k','1M'})
%XTickLabels = cellstr(num2str(round(log10(xticks(:))), '10^%d'));
%xticklabels(XTickLabels)

grid on

hold off

