%% Clear Workspace and Close figures
clear; close all; clc;

%% Intialize Laplace variable
s = zpk('s');

%% Path for functions, data and scripts
addpath('./mat/'); % Path for data

%% Colors for the figures
colors = colororder;

%% Frequency Vector [Hz]
freqs = logspace(0, 3, 1000);

%% Load the PSD of disturbances
load('uniaxial_disturbance_psd.mat', 'f', 'psd_ft', 'psd_xf');

%% Load Plants Dynamics
load('uniaxial_plants.mat', 'G_vc_light', 'G_md_light', 'G_pz_light', ...
                            'G_vc_mid',   'G_md_mid',   'G_pz_mid', ...
                            'G_vc_heavy', 'G_md_heavy', 'G_pz_heavy');

%% Load Damped Plants
load('uniaxial_damped_plants.mat', 'G_iff_vc_light', 'G_iff_md_light', 'G_iff_pz_light', ...
                                   'G_rdc_vc_light', 'G_rdc_md_light', 'G_rdc_pz_light', ...
                                   'G_dvf_vc_light', 'G_dvf_md_light', 'G_dvf_pz_light', ...
                                   'G_iff_vc_mid',   'G_iff_md_mid',   'G_iff_pz_mid', ...
                                   'G_rdc_vc_mid',   'G_rdc_md_mid',   'G_rdc_pz_mid', ...
                                   'G_dvf_vc_mid',   'G_dvf_md_mid',   'G_dvf_pz_mid', ...
                                   'G_iff_vc_heavy', 'G_iff_md_heavy', 'G_iff_pz_heavy', ...
                                   'G_rdc_vc_heavy', 'G_rdc_md_heavy', 'G_rdc_pz_heavy', ...
                                   'G_dvf_vc_heavy', 'G_dvf_md_heavy', 'G_dvf_pz_heavy');

% Damped Plant Dynamics
% <<ssec:uniaxial_position_control_damped_dynamics>>

% As was shown in Section ref:sec:uniaxial_active_damping, all three proposed active damping techniques yield similar damping plants.
% Therefore, /Integral Force Feedback/ will be used in this section to study the HAC-LAC performance.

% The obtained damped plants for the three nano-hexapod stiffnesses are shown in Figure ref:fig:uniaxial_hac_iff_damped_plants_masses.
% For $k_n = 0.01\,N/\mu m$ and $k_n = 1\,N/\mu m$, except from the nano-hexapod mode, the dynamics is quite simple and can be well approximated by a second order plant.
% However, this is not the case for the stiff nano-hexapod ($k_n = 100\,N/\mu m$) where two modes can be seen (Figure ref:fig:uniaxial_hac_iff_damped_plants_masses, right).

% This is due to the interaction between the micro-station (modelled modes at 70Hz, 140Hz and 320Hz) and the nano-hexapod.
% Such effect will be explained in Section ref:sec:uniaxial_support_compliance.


%% Damped plant - Robustness to change of sample's mass
figure;
tiledlayout(3, 3, 'TileSpacing', 'Compact', 'Padding', 'None');

ax1 = nexttile([2,1]);
hold on;
plot(freqs, abs(squeeze(freqresp(G_vc_light('d', 'f'), freqs, 'Hz'))), 'color', [colors(1,:), 0.5]);
plot(freqs, abs(squeeze(freqresp(G_iff_vc_light('d', 'f'), freqs, 'Hz'))), 'color', colors(1,:));
plot(freqs, abs(squeeze(freqresp(G_iff_vc_mid(  'd', 'f'), freqs, 'Hz'))), 'color', colors(2,:));
plot(freqs, abs(squeeze(freqresp(G_iff_vc_heavy('d', 'f'), freqs, 'Hz'))), 'color', colors(3,:));
loglog(10.^(0.4*cos([0:0.01:2*pi])+log10(100)), ...
       10.^(0.8*sin([0:0.01:2*pi]-pi/4)+log10(8e-8)), 'k--');
text(20, 4e-8, sprintf('Small\nInteraction'), 'horizontalalignment', 'center');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude [m/N]'); set(gca, 'XTickLabel',[]);
title('$k_n = 0.01\,N/\mu m$');
ylim([5e-10, 1e-3]);

ax2 = nexttile([2,1]);
hold on;
plot(freqs, abs(squeeze(freqresp(G_md_light('d', 'f'), freqs, 'Hz'))), 'color', [colors(1,:), 0.5]);
plot(freqs, abs(squeeze(freqresp(G_iff_md_light('d', 'f'), freqs, 'Hz'))), 'color', colors(1,:));
plot(freqs, abs(squeeze(freqresp(G_iff_md_mid(  'd', 'f'), freqs, 'Hz'))), 'color', colors(2,:));
plot(freqs, abs(squeeze(freqresp(G_iff_md_heavy('d', 'f'), freqs, 'Hz'))), 'color', colors(3,:));
loglog(10.^(0.4*cos([0:0.01:2*pi])+log10(200)), ...
       10.^(0.8*sin([0:0.01:2*pi]-pi/4)+log10(2e-8)), 'k--');
text(40, 1e-8, sprintf('Small\nInteraction'), 'horizontalalignment', 'center');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
set(gca, 'XTickLabel',[]); set(gca, 'YTickLabel',[]);
title('$k_n = 1\,N/\mu m$');
ylim([5e-10, 1e-3]);

ax3 = nexttile([2,1]);
hold on;
plot(freqs, abs(squeeze(freqresp(G_pz_light('d', 'f'), freqs, 'Hz'))), 'color', [colors(1,:), 0.5], 'DisplayName', '$m_s = 1\,kg$, OL');
plot(freqs, abs(squeeze(freqresp(G_iff_pz_light('d', 'f'), freqs, 'Hz'))), 'color', colors(1,:), 'DisplayName', '$m_s = 1\,kg$, IFF');
plot(freqs, abs(squeeze(freqresp(G_iff_pz_mid(  'd', 'f'), freqs, 'Hz'))), 'color', colors(2,:), 'DisplayName', '$m_s = 25\,kg$, IFF');
plot(freqs, abs(squeeze(freqresp(G_iff_pz_heavy('d', 'f'), freqs, 'Hz'))), 'color', colors(3,:), 'DisplayName', '$m_s = 50\,kg$, IFF');
loglog(10.^(0.8*cos([0:0.01:2*pi])+log10(350)), ...
       10.^(1.2*sin([0:0.01:2*pi])+log10(8e-9)), 'k--', 'HandleVisibility', 'off');
text(200, 5e-7, sprintf('$\\mu$ Station\nCoupling'), 'horizontalalignment', 'center');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
set(gca, 'XTickLabel',[]); set(gca, 'YTickLabel',[]);
title('$k_n = 100\,N/\mu m$');
ylim([5e-10, 1e-3]);
ldg = legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 1);
ldg.ItemTokenSize = [20, 1];

ax1b = nexttile();
hold on;
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_vc_light('d', 'f'), freqs, 'Hz')))), 'color', [colors(1,:), 0.5]);
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_iff_vc_light('d', 'f'), freqs, 'Hz')))), 'color', colors(1,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_iff_vc_mid(  'd', 'f'), freqs, 'Hz')))), 'color', colors(2,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_iff_vc_heavy('d', 'f'), freqs, 'Hz')))), 'color', colors(3,:));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
xticks([1e0, 1e1, 1e2]);
yticks(-360:90:360);
ylim([-200, 20]);

ax2b = nexttile();
hold on;
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_md_light('d', 'f'), freqs, 'Hz')))), 'color', [colors(1,:), 0.5]);
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_iff_md_light('d', 'f'), freqs, 'Hz')))), 'color', colors(1,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_iff_md_mid(  'd', 'f'), freqs, 'Hz')))), 'color', colors(2,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_iff_md_heavy('d', 'f'), freqs, 'Hz')))), 'color', colors(3,:));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); set(gca, 'YTickLabel',[]);
xticks([1e0, 1e1, 1e2]);
yticks(-360:90:360);
ylim([-200, 20]);

ax3b = nexttile();
hold on;
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_pz_light('d', 'f'), freqs, 'Hz')))), 'color', [colors(1,:), 0.5]);
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_iff_pz_light('d', 'f'), freqs, 'Hz')))), 'color', colors(1,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_iff_pz_mid(  'd', 'f'), freqs, 'Hz')))), 'color', colors(2,:));
plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_iff_pz_heavy('d', 'f'), freqs, 'Hz')))), 'color', colors(3,:));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); set(gca, 'YTickLabel',[]);
xticks([1e0, 1e1, 1e2]);
yticks(-360:90:360);
ylim([-200, 20]);

linkaxes([ax1,ax2,ax3,ax1b,ax2b,ax3b],'x');
xlim([1, 1e3]);

% Position Feedback Controller
% <<ssec:uniaxial_position_control_design>>

% The objective now is to design a position feedback controller for each of the three nano-hexapods that are robust to the change of sample's mass.

% The required feedback bandwidth was approximately determined in Section ref:sec:uniaxial_noise_budgeting:
% - $\approx 10\,\text{Hz}$ for the soft nano-hexapod ($k_n = 0.01\,N/\mu m$).
%   Near this frequency, the plants are equivalent to a mass line.
%   The gain of the mass line can vary up to a fact $\approx 5$ (suspended mass from $16\,kg$ up to $65\,kg$).
%   This means that the designed controller will need to have large gain margins to be robust to the change of sample's mass.
% - $\approx 50\,\text{Hz}$ for the relatively stiff nano-hexapod ($k_n = 1\,N/\mu m$).
%   Similarly to the soft nano-hexapod, the plants near the crossover frequency are equivalent to a mass line.
%   It will be probably easier to have a little bit more bandwidth in this configuration to be further away from the nano-hexapod suspension mode.
% - $\approx 100\,\text{Hz}$ for the stiff nano-hexapod ($k_n = 100\,N/\mu m$).
%   Contrary to the two first nano-hexapod stiffnesses, here the plants have more complex dynamics near the wanted crossover frequency.
%   The micro-station is not stiff enough to have a clear stiffness line at this frequency.
%   Therefore, there are both a change of phase and gain depending on the sample's mass.
%   This makes the robust design of the controller a little bit more complicated.


% Position feedback controllers are designed for each nano-hexapod such that it is stable for all considered sample masses with similar stability margins (see Nyquist plots in Figure ref:fig:uniaxial_nyquist_hac).
% These high authority controllers are generally composed of a lag at low frequency for disturbance rejection, a lead to increase the phase margin near the crossover frequency and a low pass filter to increase the robustness to high frequency dynamics.
% The controllers used for the three nano-hexapod are shown in Equation eqref:eq:uniaxial_hac_formulas, and the parameters are summarized in Table ref:tab:uniaxial_feedback_controller_parameters.

% \begin{subequations} \label{eq:uniaxial_hac_formulas}
% \begin{align}
% K_{\text{soft}}(s) &= g \cdot
%        \underbrace{\frac{s + \omega_0}{s + \omega_i}}_{\text{lag}} \cdot
%        \underbrace{\frac{1 + \frac{s}{\omega_c/\sqrt{a}}}{1 + \frac{s}{\omega_c \sqrt{a}}}}_{\text{lead}} \cdot
%        \underbrace{\frac{1}{1 + \frac{s}{\omega_l}}}_{\text{LPF}} \\
% K_{\text{mid}}(s) &= g \cdot
%        \underbrace{\left(\frac{s + \omega_0}{s + \omega_i}\right)^2}_{\text{2 lags}} \cdot
%        \underbrace{\frac{1 + \frac{s}{\omega_c/\sqrt{a}}}{1 + \frac{s}{\omega_c \sqrt{a}}}}_{\text{lead}} \cdot
%        \underbrace{\frac{1}{1 + \frac{s}{\omega_l}}}_{\text{LPF}} \\
% K_{\text{stiff}}(s) &= g \cdot
%        \underbrace{\left(\frac{1}{s + \omega_i}\right)^2}_{\text{2 lags}} \cdot
%        \underbrace{\left(\frac{1 + \frac{s}{\omega_c/\sqrt{a}}}{1 + \frac{s}{\omega_c \sqrt{a}}}\right)^2}_{\text{2 leads}} \cdot
%        \underbrace{\frac{1}{1 + \frac{s}{\omega_l}}}_{\text{LPF}}
% \end{align}
% \end{subequations}

% #+name: tab:uniaxial_feedback_controller_parameters
% #+caption: Parameters used for the position feedback controllers
% #+attr_latex: :environment tabularx :width \linewidth :align lXXX
% #+attr_latex: :center t :booktabs t
% |        | *Soft*                                    | *Moderately stiff*                         | *Stiff*                                  |
% |--------+-------------------------------------------+--------------------------------------------+------------------------------------------|
% | *Gain* | $g = 4 \cdot 10^5$                        | $g = 3 \cdot 10^6$                         | $g = 6 \cdot 10^12$                      |
% | *Lead* | $a = 5$, $\omega_c = 20\,Hz$              | $a = 4$, $\omega_c = 70\,Hz$               | $a = 5$, $\omega_c = 100\,Hz$            |
% | *Lag*  | $\omega_0 = 5\,Hz$, $\omega_i = 0.01\,Hz$ | $\omega_0 = 20\,Hz$, $\omega_i = 0.01\,Hz$ | $\omega_i = 0.01\,Hz$   |
% | *LPF*  | $\omega_l = 200\,Hz$                      | $\omega_l = 300\,Hz$                       | $\omega_l = 500\,Hz$                     |

% The loop gains for the three nano-hexapod are shown in Figure ref:fig:uniaxial_loop_gain_hac.
% We can see that:
% - for the soft and moderately stiff nano-hexapod, the crossover frequency varies a lot with the sample mass.
%   This is due to the fact that the crossover frequency corresponds to the mass line of the plant.
% - for the stiff nano-hexapod, the obtained crossover frequency is not at high as what was estimated necessary.
%   The crossover frequency in that case is close to the stiffness line of the plant, which makes the robust design of the controller easier.

% Note that these controllers were quickly tuned by hand and not designed using any optimization methods.
% The goal is just to have a first estimation of the attainable performance.


%% High Authority Controller - Soft Nano-Hexapod
% Lead to increase phase margin
a  = 5;  % Amount of phase lead / width of the phase lead / high frequency gain
wc = 2*pi*20; % Frequency with the maximum phase lead [rad/s]
H_lead = (1 + s/(wc/sqrt(a)))/(1 + s/(wc*sqrt(a)));

% Lag at low frequency
H_lag = (s + 2*pi*5)/(s + 2*pi*0.01);

% Low Pass filter to increase robustness
H_lpf = 1/(1 + s/2/pi/200);

% High Authority Controller
K_hac_vc = 4e5    * ... % Gain
           H_lead * ... % Lead
           H_lag  * ... % Lag
           H_lpf;       % LPF
K_hac_vc.InputName = {'d'};
K_hac_vc.OutputName = {'f'};

%% High Authority Controller - Mid Stiffness Nano-Hexapod
% Lead to increase phase margin
a  = 4;  % Amount of phase lead / width of the phase lead / high frequency gain
wc = 2*pi*70; % Frequency with the maximum phase lead [rad/s]

H_lead = (1 + s/(wc/sqrt(a)))/(1 + s/(wc*sqrt(a)));

% Lag at low frequency
H_lag = ((s + 2*pi*15)/(s + 2*pi*0.01))^2;

% Low Pass filter to increase robustness
H_lpf = 1/(1 + s/2/pi/300);

% High Authority Controller
K_hac_md = 3e6    * ... % Gain
           H_lead * ... % Lead
           H_lag  * ... % Lag
           H_lpf;       % LPF
K_hac_md.InputName = {'d'};
K_hac_md.OutputName = {'f'};

%% High Authority Controller - Stiff Nano-Hexapod
% Lead to increase phase margin
a  = 5;  % Amount of phase lead / width of the phase lead / high frequency gain
wc = 2*pi*100; % Frequency with the maximum phase lead [rad/s]
H_lead = ((1 + s/(wc/sqrt(a)))/(1 + s/(wc*sqrt(a))))^2;

% Integrator
H_int = 1/(s + 2*pi*0.01)^2;

% Low Pass filter to increase robustness
H_lpf = 1/(1 + s/2/pi/500);

% High Authority Controller
K_hac_pz = 6e12   * ... % Gain
           H_lead * ... % Lead
           H_int  * ... % Lag
           H_lpf;       % LPF
K_hac_pz.InputName = {'d'};
K_hac_pz.OutputName = {'f'};

%% Save High Authority Controllers
save('./mat/uniaxial_high_authority_controllers.mat', ...
        'K_hac_vc', 'K_hac_md', 'K_hac_pz');

%% Compute Loop gain for Nyquist Plot
L_vc_light = squeeze(freqresp(K_hac_vc*G_iff_vc_light('d', 'f'), freqs, 'Hz'));
L_vc_mid   = squeeze(freqresp(K_hac_vc*G_iff_vc_mid(  'd', 'f'), freqs, 'Hz'));
L_vc_heavy = squeeze(freqresp(K_hac_vc*G_iff_vc_heavy('d', 'f'), freqs, 'Hz'));

L_md_light = squeeze(freqresp(K_hac_md*G_iff_md_light('d', 'f'), freqs, 'Hz'));
L_md_mid   = squeeze(freqresp(K_hac_md*G_iff_md_mid(  'd', 'f'), freqs, 'Hz'));
L_md_heavy = squeeze(freqresp(K_hac_md*G_iff_md_heavy('d', 'f'), freqs, 'Hz'));

L_pz_light = squeeze(freqresp(K_hac_pz*G_iff_pz_light('d', 'f'), freqs, 'Hz'));
L_pz_mid   = squeeze(freqresp(K_hac_pz*G_iff_pz_mid(  'd', 'f'), freqs, 'Hz'));
L_pz_heavy = squeeze(freqresp(K_hac_pz*G_iff_pz_heavy('d', 'f'), freqs, 'Hz'));

%% Nyquist Plot - Hight Authority Controller for all three nano-hexapod stiffnesses and all sample masses
figure;
hold on;
plot(real(L_vc_light), +imag(L_vc_light), '-', 'color', [colors(1,:), 0.5], 'DisplayName', '$k_n = 0.01\,N/\mu m$')
plot(real(L_vc_light), -imag(L_vc_light), '-', 'color', [colors(1,:), 0.5], 'HandleVisibility', 'off')
plot(real(L_vc_mid  ), +imag(L_vc_mid  ), '-', 'color', [colors(1,:), 0.5], 'HandleVisibility', 'off')
plot(real(L_vc_mid  ), -imag(L_vc_mid  ), '-', 'color', [colors(1,:), 0.5], 'HandleVisibility', 'off')
plot(real(L_vc_heavy), +imag(L_vc_heavy), '-', 'color', [colors(1,:), 0.5], 'HandleVisibility', 'off')
plot(real(L_vc_heavy), -imag(L_vc_heavy), '-', 'color', [colors(1,:), 0.5], 'HandleVisibility', 'off')

plot(real(L_md_light), +imag(L_md_light), '-', 'color', [colors(2,:), 0.5], 'DisplayName', '$k_n = 1\,N/\mu m$')
plot(real(L_md_light), -imag(L_md_light), '-', 'color', [colors(2,:), 0.5], 'HandleVisibility', 'off')
plot(real(L_md_mid  ), +imag(L_md_mid  ), '-', 'color', [colors(2,:), 0.5], 'HandleVisibility', 'off')
plot(real(L_md_mid  ), -imag(L_md_mid  ), '-', 'color', [colors(2,:), 0.5], 'HandleVisibility', 'off')
plot(real(L_md_heavy), +imag(L_md_heavy), '-', 'color', [colors(2,:), 0.5], 'HandleVisibility', 'off')
plot(real(L_md_heavy), -imag(L_md_heavy), '-', 'color', [colors(2,:), 0.5], 'HandleVisibility', 'off')

plot(real(L_pz_light), +imag(L_pz_light), '-', 'color', [colors(3,:), 0.5], 'DisplayName', '$k_n = 100\,N/\mu m$')
plot(real(L_pz_light), -imag(L_pz_light), '-', 'color', [colors(3,:), 0.5], 'HandleVisibility', 'off')
plot(real(L_pz_mid  ), +imag(L_pz_mid  ), '-', 'color', [colors(3,:), 0.5], 'HandleVisibility', 'off')
plot(real(L_pz_mid  ), -imag(L_pz_mid  ), '-', 'color', [colors(3,:), 0.5], 'HandleVisibility', 'off')
plot(real(L_pz_heavy), +imag(L_pz_heavy), '-', 'color', [colors(3,:), 0.5], 'HandleVisibility', 'off')
plot(real(L_pz_heavy), -imag(L_pz_heavy), '-', 'color', [colors(3,:), 0.5], 'HandleVisibility', 'off')
plot(-1, 0, 'kx', 'HandleVisibility', 'off');
hold off;
set(gca, 'XScale', 'lin'); set(gca, 'YScale', 'lin');
xlabel('Real'); ylabel('Imag');
legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 1);
xlim([-3.8, 0.2]); ylim([-2, 2]);
axis square;



% #+name: fig:uniaxial_nyquist_hac
% #+caption: Nyquist Plot - Hight Authority Controller for all three nano-hexapod stiffnesses (soft one in blue, moderately stiff in red and very stiff in yellow) and all sample masses (corresponding to the three curves of each color)
% #+RESULTS:
% [[file:figs/uniaxial_nyquist_hac.png]]


%% Loop Gain - High Authority Controller for all three nano-hexapod stiffnesses and all sample masses
figure;
tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None');

ax1 = nexttile([2,1]);
hold on;
plot(freqs, abs(L_vc_light), 'color', [colors(1,:), 0.5], 'DisplayName', '$k_n = 0.01\,N/\mu m$');
plot(freqs, abs(L_vc_mid),   'color', [colors(1,:), 0.5], 'HandleVisibility', 'off');
plot(freqs, abs(L_vc_heavy), 'color', [colors(1,:), 0.5], 'HandleVisibility', 'off');
plot(freqs, abs(L_md_light), 'color', [colors(2,:), 0.5], 'DisplayName', '$k_n = 1\,N/\mu m$');
plot(freqs, abs(L_md_mid),   'color', [colors(2,:), 0.5], 'HandleVisibility', 'off');
plot(freqs, abs(L_md_heavy), 'color', [colors(2,:), 0.5], 'HandleVisibility', 'off');
plot(freqs, abs(L_pz_light), 'color', [colors(3,:), 0.5], 'DisplayName', '$k_n = 100\,N/\mu m$');
plot(freqs, abs(L_pz_mid),   'color', [colors(3,:), 0.5], 'HandleVisibility', 'off');
plot(freqs, abs(L_pz_heavy), 'color', [colors(3,:), 0.5], 'HandleVisibility', 'off');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Loop Gain'); set(gca, 'XTickLabel',[]);
ylim([1e-3, 1e3]);
legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 1);

ax2 = nexttile;
hold on;
plot(freqs, 180/pi*unwrap(angle(L_vc_light)), 'color', [colors(1,:), 0.5]);
plot(freqs, 180/pi*unwrap(angle(L_vc_mid  )), 'color', [colors(1,:), 0.5]);
plot(freqs, 180/pi*unwrap(angle(L_vc_heavy)), 'color', [colors(1,:), 0.5]);
plot(freqs, 180/pi*unwrap(angle(L_md_light)), 'color', [colors(2,:), 0.5]);
plot(freqs, 180/pi*unwrap(angle(L_md_mid  )), 'color', [colors(2,:), 0.5]);
plot(freqs, 180/pi*unwrap(angle(L_md_heavy)), 'color', [colors(2,:), 0.5]);
plot(freqs, 180/pi*unwrap(angle(L_pz_light)), 'color', [colors(3,:), 0.5]);
plot(freqs, 180/pi*unwrap(angle(L_pz_mid  )), 'color', [colors(3,:), 0.5]);
plot(freqs, 180/pi*unwrap(angle(L_pz_heavy)), 'color', [colors(3,:), 0.5]);
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
hold off;
yticks(-360:90:360);
ylim([-270, 0]);

linkaxes([ax1,ax2],'x');
xlim([1, 500]);

% Closed-Loop Noise Budgeting
% <<ssec:uniaxial_position_control_cl_noise_budget>>

% The high authority position feedback controllers are then implemented and the closed-loop sensitivity to disturbances are computed.
% These are compared with the open-loop and damped plants cases in Figure ref:fig:uniaxial_sensitivity_dist_hac_lac for just one configuration (moderately stiff nano-hexapod with 25kg sample's mass).
% As expected, the sensitivity to disturbances is decreased in the controller bandwidth and slightly increase outside this bandwidth.


%% Compute Closed Loop Systems
G_hac_iff_vc_light = feedback(G_iff_vc_light, K_hac_vc, 'name', -1);
G_hac_iff_vc_mid   = feedback(G_iff_vc_mid  , K_hac_vc, 'name', -1);
G_hac_iff_vc_heavy = feedback(G_iff_vc_heavy, K_hac_vc, 'name', -1);

G_hac_iff_md_light = feedback(G_iff_md_light, K_hac_md, 'name', -1);
G_hac_iff_md_mid   = feedback(G_iff_md_mid  , K_hac_md, 'name', -1);
G_hac_iff_md_heavy = feedback(G_iff_md_heavy, K_hac_md, 'name', -1);

G_hac_iff_pz_light = feedback(G_iff_pz_light, K_hac_pz, 'name', -1);
G_hac_iff_pz_mid   = feedback(G_iff_pz_mid  , K_hac_pz, 'name', -1);
G_hac_iff_pz_heavy = feedback(G_iff_pz_heavy, K_hac_pz, 'name', -1);

%% Verify Stability
isstable(G_hac_iff_vc_light) && isstable(G_hac_iff_vc_mid) && isstable(G_hac_iff_vc_heavy)

isstable(G_hac_iff_md_light) && isstable(G_hac_iff_md_mid) && isstable(G_hac_iff_md_heavy)

isstable(G_hac_iff_pz_light) && isstable(G_hac_iff_pz_mid) && isstable(G_hac_iff_pz_heavy)

%% Change of sensitivity to disturbances with LAC and with HAC-LAC
figure;
tiledlayout(1, 3, 'TileSpacing', 'Compact', 'Padding', 'None');

ax1 = nexttile();
hold on;
plot(freqs, abs(squeeze(freqresp(G_md_mid(    'd', 'fs'), freqs, 'Hz'))), 'k-');
plot(freqs, abs(squeeze(freqresp(G_iff_md_mid('d', 'fs'), freqs, 'Hz'))), 'color', colors(1,:));
plot(freqs, abs(squeeze(freqresp(G_hac_iff_md_mid('d', 'fs'), freqs, 'Hz'))), 'color', colors(2,:));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xticks([1e0, 1e1, 1e2]);
ylabel('Amplitude $d/f_{s}$ [m/N]'); xlabel('Frequency [Hz]');

ax2 = nexttile();
hold on;
plot(freqs, abs(squeeze(freqresp(G_md_mid(    'd', 'ft'), freqs, 'Hz'))), 'k-');
plot(freqs, abs(squeeze(freqresp(G_iff_md_mid('d', 'ft'), freqs, 'Hz'))), 'color', colors(1,:));
plot(freqs, abs(squeeze(freqresp(G_hac_iff_md_mid('d', 'ft'), freqs, 'Hz'))), 'color', colors(2,:));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xticks([1e0, 1e1, 1e2]);
ylabel('Amplitude $d/f_{t}$ [m/N]'); xlabel('Frequency [Hz]');

ax3 = nexttile();
hold on;
plot(freqs, abs(squeeze(freqresp(G_md_mid(    'd', 'xf'), freqs, 'Hz'))), 'k-', 'DisplayName', 'OL');
plot(freqs, abs(squeeze(freqresp(G_iff_md_mid('d', 'xf'), freqs, 'Hz'))), 'color', colors(1,:), 'DisplayName', 'IFF');
plot(freqs, abs(squeeze(freqresp(G_hac_iff_md_mid('d', 'xf'), freqs, 'Hz'))), 'color', colors(2,:), 'DisplayName', 'HAC-IFF');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xticks([1e0, 1e1, 1e2]);
ylabel('Amplitude $d/x_{f}$ [m/m]'); xlabel('Frequency [Hz]');
legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 1);

linkaxes([ax1,ax2,ax3],'x');
xlim([1, 500]);



% #+name: fig:uniaxial_sensitivity_dist_hac_lac
% #+caption: Change of sensitivity to disturbances with LAC and with HAC-LAC. Nano-Hexapod with $k_n = 1\,N/\mu m$ and sample mass of $25\,kg$ are used.
% #+RESULTS:
% [[file:figs/uniaxial_sensitivity_dist_hac_lac.png]]

% The cumulative amplitude spectrum of the motion $d$ is computed for all nano-hexapod configurations, all sample masses and in the open-loop (OL), damped (IFF) and position controlled (HAC-IFF) cases.
% The results are shown in Figure ref:fig:uniaxial_cas_hac_lac.
% Obtained root mean square values of the distance $d$ are better for the soft nano-hexapod ($\approx 25\,nm$ to $\approx 35\,nm$ depending on the sample's mass) than for the stiffer nano-hexapod (from $\approx 30\,nm$ to $\approx 70\,nm$).


%% Cumulative Amplitude Spectrum for all three nano-hexapod stiffnesses - Comparison of OL, IFF and HAC-LAC cases
figure;
tiledlayout(1, 3, 'TileSpacing', 'Compact', 'Padding', 'None');

ax1 = nexttile();
hold on;
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_vc_mid('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_vc_mid('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [0,0,0,0.5], 'DisplayName', 'OL');
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_vc_light('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_vc_light('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [0,0,0,0.5], 'HandleVisibility', 'off');
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_vc_heavy('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_vc_heavy('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [0,0,0,0.5], 'HandleVisibility', 'off');
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_iff_vc_mid('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_iff_vc_mid('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [colors(1,:), 0.5], 'DisplayName', 'IFF');
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_iff_vc_light('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_iff_vc_light('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [colors(1,:), 0.5], 'HandleVisibility', 'off');
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_iff_vc_heavy('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_iff_vc_heavy('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [colors(1,:), 0.5], 'HandleVisibility', 'off');
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_hac_iff_vc_mid('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_hac_iff_vc_mid('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [colors(2,:), 0.5], 'DisplayName', 'HAC-IFF');
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_hac_iff_vc_light('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_hac_iff_vc_light('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [colors(2,:), 0.5], 'HandleVisibility', 'off');
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_hac_iff_vc_heavy('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_hac_iff_vc_heavy('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [colors(2,:), 0.5], 'HandleVisibility', 'off');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xticks([1e0, 1e1, 1e2]);
ylabel('CAS of $d$ [m]'); xlabel('Frequency [Hz]');
legend('location', 'southwest', 'FontSize', 8, 'NumColumns', 1);
title('$k_n = 0.01\,N/\mu m$');

ax2 = nexttile();
hold on;
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_md_mid('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_md_mid('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [0,0,0,0.5]);
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_md_light('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_md_light('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [0,0,0,0.5]);
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_md_heavy('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_md_heavy('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [0,0,0,0.5]);
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_iff_md_mid('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_iff_md_mid('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [colors(1,:), 0.5]);
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_iff_md_light('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_iff_md_light('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [colors(1,:), 0.5]);
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_iff_md_heavy('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_iff_md_heavy('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [colors(1,:), 0.5]);
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_hac_iff_md_mid('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_hac_iff_md_mid('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [colors(2,:), 0.5]);
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_hac_iff_md_light('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_hac_iff_md_light('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [colors(2,:), 0.5]);
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_hac_iff_md_heavy('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_hac_iff_md_heavy('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [colors(2,:), 0.5]);
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xticks([1e0, 1e1, 1e2]);
xlabel('Frequency [Hz]'); set(gca, 'YTickLabel',[]);
title('$k_n = 1\,N/\mu m$');

ax3 = nexttile();
hold on;
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_pz_mid('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_pz_mid('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [0,0,0,0.5]);
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_pz_light('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_pz_light('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [0,0,0,0.5]);
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_pz_heavy('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_pz_heavy('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [0,0,0,0.5]);
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_iff_pz_mid('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_iff_pz_mid('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [colors(1,:), 0.5]);
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_iff_pz_light('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_iff_pz_light('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [colors(1,:), 0.5]);
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_iff_pz_heavy('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_iff_pz_heavy('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [colors(1,:), 0.5]);
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_hac_iff_pz_mid('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_hac_iff_pz_mid('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [colors(2,:), 0.5]);
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_hac_iff_pz_light('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_hac_iff_pz_light('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [colors(2,:), 0.5]);
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_hac_iff_pz_heavy('d', 'ft'), f, 'Hz'))).^2 + ...
                                     psd_xf.*abs(squeeze(freqresp(G_hac_iff_pz_heavy('d', 'xf'), f, 'Hz'))).^2)))), '-', ...
     'color', [colors(2,:), 0.5]);
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
xticks([1e0, 1e1, 1e2]);
xlabel('Frequency [Hz]'); set(gca, 'YTickLabel',[]);
title('$k_n = 100\,N/\mu m$');

linkaxes([ax1,ax2,ax3],'xy');
xlim([1, 500]);
ylim([2e-10, 3e-6])
