%% 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
addpath('./src/'); % Path for functions

%% 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');

% Sensitivity to disturbances
% From the Uni-axial model, the transfer function from the disturbances ($f_s$, $x_f$ and $f_t$) to the displacement $d$ are computed.

% This is done for *two extreme sample masses* $m_s = 1\,\text{kg}$ and $m_s = 50\,\text{kg}$ and *three nano-hexapod stiffnesses*:
% - $k_n = 0.01\,N/\mu m$ that could represent a voice coil actuator with soft flexible guiding
% - $k_n = 1\,N/\mu m$ that could represent a voice coil actuator with a stiff flexible guiding or a mechanically amplified piezoelectric actuator
% - $k_n = 100\,N/\mu m$ that could represent a stiff piezoelectric stack actuator

% The obtained sensitivity to disturbances for the three nano-hexapod stiffnesses are shown in Figure ref:fig:uniaxial_sensitivity_disturbances_nano_hexapod_stiffnesses for the light sample (same conclusions can be drawn with the heavy one).

% #+begin_important
% From Figure ref:fig:uniaxial_sensitivity_disturbances_nano_hexapod_stiffnesses, following can be concluded for the *soft nano-hexapod*:
% - It is more sensitive to forces applied on the sample (cable forces for instance), which is expected due to the lower stiffness
% - Between the suspension mode of the nano-hexapod (here at 5Hz) and the first mode of the micro-station (here at 70Hz), the disturbances induced by the stage vibrations are filtered out.
% - Above the suspension mode of the nano-hexapod, the sample's motion is unaffected by the floor motion, and therefore the sensitivity to floor motion is almost $1$.
% #+end_important


%% Sensitivity to disturbances for three different nano-hexpod stiffnesses
figure;
tiledlayout(1, 3, 'TileSpacing', 'Compact', 'Padding', 'None');

ax1 = nexttile();
hold on;
plot(freqs, abs(squeeze(freqresp(G_vc_light('d', 'fs'), freqs, 'Hz'))));
plot(freqs, abs(squeeze(freqresp(G_md_light('d', 'fs'), freqs, 'Hz'))));
plot(freqs, abs(squeeze(freqresp(G_pz_light('d', 'fs'), freqs, 'Hz'))));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude $d/f_{s}$ [m/N]'); xlabel('Frequency [Hz]');
xticks([1e0, 1e1, 1e2]);

ax2 = nexttile();
hold on;
plot(freqs, abs(squeeze(freqresp(G_vc_light('d', 'ft'), freqs, 'Hz'))));
plot(freqs, abs(squeeze(freqresp(G_md_light('d', 'ft'), freqs, 'Hz'))));
plot(freqs, abs(squeeze(freqresp(G_pz_light('d', 'ft'), freqs, 'Hz'))));
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude $d/f_{t}$ [m/N]'); xlabel('Frequency [Hz]');
xticks([1e0, 1e1, 1e2]);

ax3 = nexttile();
hold on;
plot(freqs, abs(squeeze(freqresp(G_vc_light('d', 'xf'), freqs, 'Hz'))), 'DisplayName', '$k_n = 0.01\,N/\mu m$');
plot(freqs, abs(squeeze(freqresp(G_md_light('d', 'xf'), freqs, 'Hz'))), 'DisplayName', '$k_n = 1   \,N/\mu m$');
plot(freqs, abs(squeeze(freqresp(G_pz_light('d', 'xf'), freqs, 'Hz'))), 'DisplayName', '$k_n = 100 \,N/\mu m$');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Amplitude $d/x_{f}$ [m/m]'); xlabel('Frequency [Hz]');
xticks([1e0, 1e1, 1e2]);
legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 1);

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

% Open-Loop Dynamic Noise Budgeting
% Now, the power spectral density of the disturbances is taken into account to estimate the residual motion $d$ in each case.

% The Cumulative Amplitude Spectrum of the relative motion $d$ due to both the floor motion $x_f$ and the stage vibrations $f_t$ are shown in Figure ref:fig:uniaxial_cas_d_disturbances_stiffnesses for the three nano-hexapod stiffnesses.

% It is shown that the effect of the floor motion is much less than the stage vibrations, except for the soft nano-hexapod below 5Hz.


%% Cumulative Amplitude Spectrum of the relative motion d, due to both the floor motion and the stage vibrations
figure;
tiledlayout(1, 1, 'TileSpacing', 'Compact', 'Padding', 'None');

ax1 = nexttile();
hold on;
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_vc_light('d', 'ft'), f, 'Hz'))).^2)))), '-', 'color', colors(1,:), 'DisplayName', '$f_t$');
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_md_light('d', 'ft'), f, 'Hz'))).^2)))), '-', 'color', colors(2,:), 'DisplayName', '$f_t$');
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_ft.*abs(squeeze(freqresp(G_pz_light('d', 'ft'), f, 'Hz'))).^2)))), '-', 'color', colors(3,:), 'DisplayName', '$f_t$');
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_xf.*abs(squeeze(freqresp(G_vc_light('d', 'xf'), f, 'Hz'))).^2)))), '--', 'color', colors(1,:), 'DisplayName', '$x_f$ ($k_n = 0.01\,N/\mu m$)');
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_xf.*abs(squeeze(freqresp(G_md_light('d', 'xf'), f, 'Hz'))).^2)))), '--', 'color', colors(2,:), 'DisplayName', '$x_f$ ($k_n = 1   \,N/\mu m$)');
plot(f, sqrt(flip(-cumtrapz(flip(f), flip(psd_xf.*abs(squeeze(freqresp(G_pz_light('d', 'xf'), f, 'Hz'))).^2)))), '--', 'color', colors(3,:), 'DisplayName', '$x_f$ ($k_n = 100 \,N/\mu m$)');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Cumulative Ampl. Spectrum [m]'); xlabel('Frequency [Hz]');
legend('location', 'southwest', 'FontSize', 8, 'NumColumns', 2);

xlim([1, 500]);
ylim([1e-12, 1e-6])



% #+name: fig:uniaxial_cas_d_disturbances_stiffnesses
% #+caption: Cumulative Amplitude Spectrum of the relative motion d, due to both the floor motion and the stage vibrations (light sample)
% #+RESULTS:
% [[file:figs/uniaxial_cas_d_disturbances_stiffnesses.png]]

% The total cumulative amplitude spectrum for the three nano-hexapod stiffnesses and for the two sample's masses are shown in Figure ref:fig:uniaxial_cas_d_disturbances_payload_masses.
% The conclusion is that the sample's mass has little effect on the cumulative amplitude spectrum of the relative motion $d$.


%% Cumulative Amplitude Spectrum of the relative motion d due to all disturbances, for two sample masses
figure;
tiledlayout(1, 1, 'TileSpacing', 'Compact', 'Padding', 'None');

ax1 = nexttile();
hold on;
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', colors(1,:), 'DisplayName', '$m_s = 1\,kg$');
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', colors(2,:), 'DisplayName', '$m_s = 1\,kg$');
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', colors(3,:), 'DisplayName', '$m_s = 1\,kg$');
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', colors(1,:), 'DisplayName', '$m_s = 50\,kg$ ($k_n = 0.01\,N/\mu m$)');
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', colors(2,:), 'DisplayName', '$m_s = 50\,kg$ ($k_n = 1\,N/\mu m$)');
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', colors(3,:), 'DisplayName', '$m_s = 50\,kg$ ($k_n = 100\,N/\mu m$)');
hold off;
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
ylabel('Cumulative Ampl. Spectrum [m]'); xlabel('Frequency [Hz]');
legend('location', 'southwest', 'FontSize', 8, 'NumColumns', 2);

xlim([1, 500]);
ylim([1e-11, 3e-6])