diff --git a/Analysis/effect_ground_motion.m b/Analysis/effect_ground_motion.m
new file mode 100644
index 0000000..1754061
--- /dev/null
+++ b/Analysis/effect_ground_motion.m
@@ -0,0 +1,48 @@
+%%
+clear; close all; clc;
+
+%%
+load('./mat/Gd_ol_cl.mat', 'Gd_ol_20', 'Gd_cl_20');
+
+%%
+load('./mat/weight_Wxg.mat', 'Wxg')
+
+%%
+load('./mat/G_gm_to_dh.mat', 'G_gm_to_dh')
+load('./mat/psd_ground_motion.mat', 'psd_f', 'pxx')
+
+%%
+bodeFig({Gd_ol_20(1, 1), G_gm_to_dh})
+
+%% PSD
+freqs = logspace(-2, 2, 1000);
+
+gm_ol = abs(squeeze(freqresp(Wxg*Gd_ol_20(1, 1), freqs, 'Hz')));
+gm_cl = abs(squeeze(freqresp(Wxg*Gd_cl_20(1, 1), freqs, 'Hz')));
+
+figure;
+hold on;
+plot(freqs, gm_ol)
+plot(freqs, gm_cl)
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('Frequency [$Hz$]'); ylabel('PSD [$m/\sqrt{Hz}$]');
+hold off;
+
+%% CAS
+freqs = logspace(-1, 2, 1000);
+
+gm_ol = abs(squeeze(freqresp(Wxg*Gd_ol_20(1, 1), freqs, 'Hz')));
+gm_cl = abs(squeeze(freqresp(Wxg*Gd_cl_20(1, 1), freqs, 'Hz')));
+
+dw = freqs - [0, freqs(1:end-1)];
+
+figure;
+hold on;
+plot(freqs, cumsum(gm_ol'.*dw))
+plot(freqs, cumsum(gm_cl'.*dw))
+set(gca, 'XScale', 'log');
+% set(gca, 'YScale', 'log');
+xlabel('Frequency [Hz]'); ylabel('CAS [m]');
+hold off;
+
+
diff --git a/Analysis/simulation_analysis.m b/Analysis/simulation_analysis.m
new file mode 100644
index 0000000..5553436
--- /dev/null
+++ b/Analysis/simulation_analysis.m
@@ -0,0 +1,109 @@
+%%
+clear; close all; clc;
+
+%%
+exp_ol   = load('./data/exp_open_loop.mat', 'Dmeas');
+exp_cl   = load('./data/exp_close_loop_xyz.mat', 'Dmeas');
+exp_without_nass = load('./data/exp_whitout_nass.mat', 'Dmeas');
+load('./mat/sim_conf.mat', 'sim_conf');
+
+exp_without_nass.Dmeas.Data(:, 3) = exp_without_nass.Dmeas.Data(:, 3) - exp_without_nass.Dmeas.Data(end, 3);
+
+%% With and without NASS
+steady_i = ceil(length(exp_ol.Dmeas.Time)/2);
+
+figure;
+hold on;
+plot(exp_without_nass.Dmeas.Data(steady_i:end, 1),exp_without_nass.Dmeas.Data(steady_i:end, 2))
+plot(exp_cl.Dmeas.Data(steady_i:end, 1),exp_cl.Dmeas.Data(steady_i:end, 2))
+xlim([-1e-6, 1e-6]);
+ylim([-1e-6, 1e-6]);
+hold off;
+xlabel('Displacement - $x$ [m]'); ylabel('Displacement - $y$ [m]');
+
+exportFig('exp_w_wo_nass_xy', 'half-short')
+
+figure;
+hold on;
+plot(exp_without_nass.Dmeas.Time(steady_i:end), exp_without_nass.Dmeas.Data(steady_i:end, 2));
+plot(exp_without_nass.Dmeas.Time(steady_i:end), exp_cl.Dmeas.Data(steady_i:end, 2));
+legend({'$y$ - without NASS', '$y$ - with NASS'})
+hold off;
+xlabel('Time [s]'); ylabel('Displacement [m]');
+
+exportFig('exp_w_wo_nass_y', 'half-small')
+
+
+
+%% Compare OL and CL - Time
+figure;
+hold on;
+plot(exp_ol.Dmeas.Time, exp_ol.Dmeas.Data(:, 1));
+plot(exp_cl.Dmeas.Time, exp_cl.Dmeas.Data(:, 1));
+legend({'x - OL', 'x - CL'})
+hold off;
+xlabel('Time [s]'); ylabel('Displacement [m]');
+
+exportFig('exp_control_time_x', 'normal-normal')
+
+
+figure;
+hold on;
+plot(exp_ol.Dmeas.Time, exp_ol.Dmeas.Data(:, 2));
+plot(exp_cl.Dmeas.Time, exp_cl.Dmeas.Data(:, 2));
+legend({'y - OL', 'y - CL'})
+hold off;
+xlabel('Time [s]'); ylabel('Displacement [m]');
+
+exportFig('exp_control_time_y', 'normal-normal')
+
+
+figure;
+hold on;
+plot(exp_ol.Dmeas.Time, exp_ol.Dmeas.Data(:, 3));
+plot(exp_cl.Dmeas.Time, exp_cl.Dmeas.Data(:, 3));
+legend({'z - OL', 'z - CL'})
+hold off;
+xlabel('Time [s]'); ylabel('Displacement [m]');
+
+exportFig('exp_control_time_z', 'normal-normal')
+
+%%
+steady_i = ceil(length(exp_ol.Dmeas.Time)/2);
+% steady_i = 1;
+
+figure;
+hold on;
+plot(exp_ol.Dmeas.Data(steady_i:end, 1),exp_ol.Dmeas.Data(steady_i:end, 2))
+plot(exp_cl.Dmeas.Data(steady_i:end, 1),exp_cl.Dmeas.Data(steady_i:end, 2))
+legend({'OL', 'CL'})
+hold off;
+xlabel('Displacement - $x$ [s]'); ylabel('Displacement - $y$ [m]');
+
+
+figure;
+hold on;
+plot3(exp_ol.Dmeas.Data(steady_i:end, 1),exp_ol.Dmeas.Data(steady_i:end, 2),exp_ol.Dmeas.Data(steady_i:end, 3))
+plot3(exp_cl.Dmeas.Data(steady_i:end, 1),exp_cl.Dmeas.Data(steady_i:end, 2),exp_cl.Dmeas.Data(steady_i:end, 3))
+legend({'OL', 'CL'})
+hold off;
+xlabel('Displacement - $x$ [s]'); ylabel('Displacement - $y$ [m]'); zlabel('Displacement - $z$ [m]');
+
+
+%% Compare OL and CL - PSD
+han_windows_ol = hanning(ceil(length(exp_ol.Dmeas.Time(steady_i:end))/2));
+[psd_x_ol, freqs_x_ol] = pwelch(exp_ol.Dmeas.Data(steady_i:end, 1), han_windows_ol, 0, [], 1/sim_conf.Ts);
+
+han_windows_cl = hanning(ceil(length(exp_cl.Dmeas.Time(steady_i:end))/2));
+[psd_x, freqs_x] = pwelch(exp_cl.Dmeas.Data(steady_i:end, 1), han_windows_cl, 0, [], 1/sim_conf.Ts);
+
+figure;
+hold on;
+plot(freqs_x_ol, sqrt(psd_x_ol));
+plot(freqs_x, sqrt(psd_x));
+set(gca,'xscale','log'); set(gca,'yscale','log');
+xlabel('Frequency [Hz]'); ylabel('PSD [$m/\sqrt{Hz}$]');
+legend({'x - OL', 'x - CL'})
+hold off;
+
+exportFig('exp_psd_x', 'normal-normal')
diff --git a/Assemblage.slx b/Assemblage.slx
index b886624..d3701b7 100644
Binary files a/Assemblage.slx and b/Assemblage.slx differ
diff --git a/run_simulations.m b/run_simulations.m
new file mode 100644
index 0000000..a09705b
--- /dev/null
+++ b/run_simulations.m
@@ -0,0 +1,41 @@
+%% Open Loop simulation and save the final state
+steady_time = 10;
+
+initializeSimConf(struct('Tsim', steady_time, 'cl_time', steady_time));
+
+set_param('Assemblage',...
+          'SaveFinalState','on',...
+          'FinalStateName','myOperPoint',...
+          'SaveCompleteFinalSimState','on'...
+);
+
+sim('Assemblage')
+
+save('./data/myOperPoint.mat', 'myOperPoint');
+
+set_param('Assemblage',...
+          'SaveFinalState','off',...
+          'SaveCompleteFinalSimState','off'...
+);
+
+save('./data/exp_open_loop.mat', 'Dmeas');
+
+%% Close the Loop and start from steady state
+sim_time = 10;
+
+initializeSimConf(struct('Tsim', steady_time+sim_time, 'cl_time', steady_time));
+
+load('./data/myOperPoint.mat', 'myOperPoint');
+
+set_param('Assemblage',...
+          'LoadInitialState','on',...
+          'InitialState','myOperPoint'...
+);
+
+sim('Assemblage')
+
+set_param('Assemblage',...
+          'LoadInitialState','off' ...
+);
+
+save('./data/exp_close_loop_xyz.mat', 'Dmeas');