select.m

% Clear workspace and command window
close all;
clear all;
clc;

param = 5; % To select wi

w = [2*pi/5, 4*pi/5];
selection = 2; % Select w1
wc = w(selection);

wi = wc / param;

% Calculate Kp for |G(jw1)| = 0dB
% Solving the equations explained in the Method part:
Kp = wc^2/sqrt(wc^2 + wi^2);
Ki = wi * Kp;

% Create transfer functions
s = tf('s');
PI_controller = Kp * (1 + wi/s);
Plant = 1 / s;
G = PI_controller * Plant;

% Verify the magnitude at w1
[mag, ph] = bode(G, wc);
mag_db = 20*log10(mag)
ph

% Bode plot
omega = logspace(-2, 2, 1000);

bp = figure;
[mag, phase, wout] = bode(G, omega);

% Convert data
mag = squeeze(mag);
phase = squeeze(phase);
mag_db = 20*log10(mag);

% Generate Bode plots
subplot(2,1,1);
semilogx(wout, mag_db, 'LineWidth', 1.5); hold on;
plot([wout(1) wout(end)], [0 0], 'k');
grid on;
ylabel('Magnitude (dB)');
title(sprintf('Bode of G(jw) for w = %.2f rad/s', wc));

% Highlight frequencies...
[~, idx1] = min(abs(wout - wc));
plot([wc wc], ylim, 'r--', 'LineWidth', 1.2);
plot(wc, mag_db(idx1), 'ro', 'MarkerFaceColor', 'r', 'MarkerSize',8);
text(wc, mag_db(idx1), sprintf('%.2f dB', mag_db(idx1)), 'VerticalAlignment', 'bottom');

subplot(2,1,2);
semilogx(wout, phase, 'LineWidth', 1.5); hold on;
grid on;
ylabel('Phase (degrees)');
xlabel('Frequency (rad/s)');
phase_margin = 180 + phase(idx1)
title(sprintf('Phase margin = %.2f degrees', phase_margin));

plot([wc wc], ylim, 'r--', 'LineWidth', 1.2);
plot(wc, phase(idx1), 'ro', 'MarkerFaceColor', 'r', 'MarkerSize',8);
text(wc, phase(idx1), sprintf('%.1fdeg', phase(idx1)), 'VerticalAlignment', 'bottom');

plot_name = sprintf('plots/ode_plot_%d', selection);
saveas(bp, plot_name, 'png');

% Display transfer function, poles, and zeros
disp('Open loop:');
G
[zeros, poles, gain] = zpkdata(G, 'v');
disp('Zeros:');
disp(zeros);
disp('Poles:');
disp(poles);
disp('Gain:');
disp(gain);

wc
Kp
Ki