diff --git a/Analysis/analysis_perturbations.m b/Analysis/analysis_perturbations.m
index 2712117..5659f3d 100644
--- a/Analysis/analysis_perturbations.m
+++ b/Analysis/analysis_perturbations.m
@@ -5,7 +5,7 @@ clear; close all; clc;
 load('./mat/G_xg_to_d.mat', 'G_xg_to_d');
 
 %% Load shape of the perturbation
-load('./mat/weight_Wxg.mat', 'Wxg');
+load('./mat/perturbations.mat', 'Wxg');
 
 %% Effect of the perturbation on the output
 freqs = logspace(-1, 3, 1000);
diff --git a/Analysis/effect_ground_motion.m b/Analysis/effect_ground_motion.m
index 1754061..120341b 100644
--- a/Analysis/effect_ground_motion.m
+++ b/Analysis/effect_ground_motion.m
@@ -5,7 +5,7 @@ clear; close all; clc;
 load('./mat/Gd_ol_cl.mat', 'Gd_ol_20', 'Gd_cl_20');
 
 %%
-load('./mat/weight_Wxg.mat', 'Wxg')
+load('./mat/perturbations.mat', 'Wxg')
 
 %%
 load('./mat/G_gm_to_dh.mat', 'G_gm_to_dh')
@@ -44,5 +44,3 @@ set(gca, 'XScale', 'log');
 % set(gca, 'YScale', 'log');
 xlabel('Frequency [Hz]'); ylabel('CAS [m]');
 hold off;
-
-
diff --git a/Analysis/ground_motion_ol_cl.m b/Analysis/ground_motion_ol_cl.m
index e620911..19ad8f0 100644
--- a/Analysis/ground_motion_ol_cl.m
+++ b/Analysis/ground_motion_ol_cl.m
@@ -123,7 +123,7 @@ exportFig('gm_control_psd_y', 'normal-normal')
 
 %% Compare OL and CL - PSD
 load('./mat/G_xg_to_d.mat', 'G_xg_to_d');
-load('./mat/weight_Wxg.mat', 'Wxg');
+load('./mat/perturbations.mat', 'Wxg');
 load('./mat/T_S.mat', 'S', 'T');
 
 freqs = logspace(-1, 3, 1000);
diff --git a/Assemblage.slx b/Assemblage.slx
index 3163c73..d771139 100644
Binary files a/Assemblage.slx and b/Assemblage.slx differ
diff --git a/active_damping/iff_control.m b/active_damping/iff_control.m
index a16c2ff..898ab60 100644
--- a/active_damping/iff_control.m
+++ b/active_damping/iff_control.m
@@ -19,4 +19,4 @@ K_iff_heavy_vc = 20/s*tf(eye(6));
 K_iff_heavy_pz = 535/s*tf(eye(6));
 
 %%
-save('./mat/K_iff_crit.mat', 'K_iff_light_vc', 'K_iff_light_pz', 'K_iff_heavy_vc', 'K_iff_heavy_pz');
\ No newline at end of file
+save('./mat/K_iff_crit.mat', 'K_iff_light_vc', 'K_iff_light_pz', 'K_iff_heavy_vc', 'K_iff_heavy_pz');
diff --git a/active_damping/iff_identification.m b/active_damping/iff_identification.m
index c2db8db..4a01dac 100644
--- a/active_damping/iff_identification.m
+++ b/active_damping/iff_identification.m
@@ -9,12 +9,12 @@ initializeSample(struct('mass', 1));
 
 initializeNanoHexapod(struct('actuator', 'lorentz'));
 K_iff = K_iff_light_vc; %#ok
-save('./mat/K_iff.mat', 'K_iff');
+save('./mat/controllers.mat', 'K_iff');
 G_iff_light_vc = identifyPlant();
 
 initializeNanoHexapod(struct('actuator', 'piezo'));
 K_iff = K_iff_light_pz; %#ok
-save('./mat/K_iff.mat', 'K_iff');
+save('./mat/controllers.mat', 'K_iff');
 G_iff_light_pz = identifyPlant();
 
 %% Heavy Sample
@@ -22,12 +22,12 @@ initializeSample(struct('mass', 50));
 
 initializeNanoHexapod(struct('actuator', 'lorentz'));
 K_iff = K_iff_heavy_vc; %#ok
-save('./mat/K_iff.mat', 'K_iff');
+save('./mat/controllers.mat', 'K_iff');
 G_iff_heavy_vc = identifyPlant();
 
 initializeNanoHexapod(struct('actuator', 'piezo'));
 K_iff = K_iff_heavy_pz;
-save('./mat/K_iff.mat', 'K_iff');
+save('./mat/controllers.mat', 'K_iff', '-append');
 G_iff_heavy_pz = identifyPlant();
 
 %% Save the obtained transfer functions
diff --git a/demonstration/Micro_Station_Displacement.slx b/demonstration/Micro_Station_Displacement.slx
index 44d4e0c..f585019 100644
Binary files a/demonstration/Micro_Station_Displacement.slx and b/demonstration/Micro_Station_Displacement.slx differ
diff --git a/demonstration/error_NASS.m b/demonstration/error_NASS.m
new file mode 100644
index 0000000..f8cec8d
--- /dev/null
+++ b/demonstration/error_NASS.m
@@ -0,0 +1,50 @@
+%% Plot all 6 errors expressed in the NASS base
+figure;
+
+%% Tx
+subaxis(2, 3, 1);
+hold on;
+plot(error_nass.Time, error_nass.Data(:, 1), 'k-', 'DisplayName', '$\epsilon_x$');
+legend();
+hold off;
+xlabel('Time (s)'); ylabel('Position (m)');
+
+%% Ty
+subaxis(2, 3, 2);
+hold on;
+plot(error_nass.Time, error_nass.Data(:, 2), 'k-', 'DisplayName', '$\epsilon_y$');
+legend();
+hold off;
+xlabel('Time (s)'); ylabel('Position (m)');
+
+%% Tz
+subaxis(2, 3, 3);
+hold on;
+plot(error_nass.Time, error_nass.Data(:, 3), 'k-', 'DisplayName', '$\epsilon_z$');
+legend();
+hold off;
+xlabel('Time (s)'); ylabel('Position (m)');
+
+%% Rx
+subaxis(2, 3, 4);
+hold on;
+plot(error_nass.Time, error_nass.Data(:, 4), 'k-', 'DisplayName', '$\epsilon_{\theta_x}$');
+legend();
+hold off;
+xlabel('Time (s)'); ylabel('Rotation (rad)');
+
+%% Ry
+subaxis(2, 3, 5);
+hold on;
+plot(error_nass.Time, error_nass.Data(:, 5), 'k-', 'DisplayName', '$\epsilon_{\theta_y}$');
+legend();
+hold off;
+xlabel('Time (s)'); ylabel('Rotation (rad)');
+
+%% Rz
+subaxis(2, 3, 6);
+hold on;
+plot(error_nass.Time, error_nass.Data(:, 6), 'k-', 'DisplayName', '$\epsilon_{\theta_z}$');
+legend();
+hold off;
+xlabel('Time (s)'); ylabel('Rotation (rad)');
diff --git a/demonstration/sample_pos_init.m b/demonstration/sample_pos_init.m
index 636d42b..e958e61 100644
--- a/demonstration/sample_pos_init.m
+++ b/demonstration/sample_pos_init.m
@@ -3,7 +3,7 @@ clear; close all; clc;
 
 %% Initialize simulation configuration
 opts_sim = struct(...
-    'Tsim', 2 ...
+    'Tsim', 1 ...
 );
 
 initializeSimConf(opts_sim);
@@ -14,19 +14,22 @@ load('./mat/sim_conf.mat', 'sim_conf')
 time_vector = 0:sim_conf.Ts:sim_conf.Tsim;
 
 % Translation Stage
-ty = 0*ones(length(time_vector), 1);
+ty = 0.05*ones(length(time_vector), 1);
 
 % Tilt Stage
 ry = 2*pi*(3/360)*ones(length(time_vector), 1);
+% ry = 2*pi*(3/360)*sin(2*pi*time_vector);
 
 % Spindle
 rz = 2*pi*1*(time_vector);
+% rz = 2*pi*(190/360)*ones(length(time_vector), 1);
 
 % Micro Hexapod
 u_hexa = zeros(length(time_vector), 6);
 
 % Gravity Compensator system
 mass = zeros(length(time_vector), 2);
+mass(:, 2) = pi;
 
 opts_inputs = struct(...
     'ty', ty, ...
diff --git a/demonstration/setpoint_vs_position.m b/demonstration/setpoint_vs_position.m
new file mode 100644
index 0000000..917fe85
--- /dev/null
+++ b/demonstration/setpoint_vs_position.m
@@ -0,0 +1,56 @@
+%%
+figure;
+
+%% Tx
+subaxis(2, 3, 1);
+hold on;
+plot(pos.Time, pos.Data(:, 1), 'k-');
+plot(setpoint.Time, setpoint.Data(:, 1), 'k--');
+legend({'x - pos', 'x - setpoint'});
+hold off;
+xlabel('Time (s)'); ylabel('Position (m)');
+
+%% Ty
+subaxis(2, 3, 2);
+hold on;
+plot(pos.Time, pos.Data(:, 2), 'k-');
+plot(setpoint.Time, setpoint.Data(:, 2), 'k--');
+legend({'y - pos', 'y - setpoint'});
+hold off;
+xlabel('Time (s)'); ylabel('Position (m)');
+
+%% Tz
+subaxis(2, 3, 3);
+hold on;
+plot(pos.Time, pos.Data(:, 3), 'k-');
+plot(setpoint.Time, setpoint.Data(:, 3), 'k--');
+legend({'z - pos', 'z - setpoint'});
+hold off;
+xlabel('Time (s)'); ylabel('Position (m)');
+
+%% Rx
+subaxis(2, 3, 4);
+hold on;
+plot(pos.Time, pos.Data(:, 4), 'k-');
+plot(setpoint.Time, setpoint.Data(:, 4), 'k--');
+legend({'$\theta_x$ - pos', '$\theta_x$ - setpoint'});
+hold off;
+xlabel('Time (s)'); ylabel('Rotation (rad)');
+
+%% Ry
+subaxis(2, 3, 5);
+hold on;
+plot(pos.Time, pos.Data(:, 5), 'k-');
+plot(setpoint.Time, setpoint.Data(:, 5), 'k--');
+legend({'$\theta_y$ - pos', '$\theta_y$ - setpoint'});
+hold off;
+xlabel('Time (s)'); ylabel('Rotation (rad)');
+
+%% Rz
+subaxis(2, 3, 6);
+hold on;
+plot(pos.Time, pos.Data(:, 6), 'k-');
+plot(setpoint.Time, setpoint.Data(:, 6), 'k--');
+legend({'$\theta_z$ - pos', '$\theta_z$ - setpoint'});
+hold off;
+xlabel('Time (s)'); ylabel('Rotation (rad)');
diff --git a/identification/id_flexible_rigid.m b/identification/id_flexible_rigid.m
new file mode 100644
index 0000000..249bfed
--- /dev/null
+++ b/identification/id_flexible_rigid.m
@@ -0,0 +1,55 @@
+%% Script Description
+% Determine if we take into account the flexibilities,
+% does that changes a lot
+
+%%
+clear; close all; clc;
+
+%% Initialize all the stage by default
+run init_data.m
+
+%% Options for Linearized
+options = linearizeOptions;
+options.SampleTime = 0;
+
+%% Name of the Simulink File
+mdl = 'simscape_id_micro_station';
+
+%% Micro-Hexapod
+% Input/Output definition
+io(1) = linio([mdl, '/Micro-Station/Fm_ext'],1,'openinput');
+io(2) = linio([mdl, '/Micro-Station/Fg_ext'],1,'openinput');
+io(3) = linio([mdl, '/Micro-Station/Dm_inertial'],1,'output');
+io(4) = linio([mdl, '/Micro-Station/Ty_inertial'],1,'output');
+io(5) = linio([mdl, '/Micro-Station/Ry_inertial'],1,'output');
+io(6) = linio([mdl, '/Micro-Station/Dg_inertial'],1,'output');
+
+
+%% Run the linearization
+initializeTy();
+
+G_ms_flexible = linearize(mdl, io, 0);
+
+% Input/Output names
+G_ms_flexible.InputName  = {'Fmx', 'Fmy', 'Fmz',...
+                   'Fgx', 'Fgy', 'Fgz'};
+G_ms_flexible.OutputName = {'Dmx', 'Dmy', 'Dmz', ...
+                   'Tyx', 'Tyy', 'Tyz', ...
+                   'Ryx', 'Ryy', 'Ryz', ...
+                   'Dgx', 'Dgy', 'Dgz'};
+
+%% Run the linearization
+initializeTy(struct('rigid', true));
+
+G_ms_ty_rigid = linearize(mdl, io, 0);
+
+% Input/Output names
+G_ms_ty_rigid.InputName  = {'Fmx', 'Fmy', 'Fmz',...
+                            'Fgx', 'Fgy', 'Fgz'};
+G_ms_ty_rigid.OutputName = {'Dmx', 'Dmy', 'Dmz', ...
+                            'Tyx', 'Tyy', 'Tyz', ...
+                            'Ryx', 'Ryy', 'Ryz', ...
+                            'Dgx', 'Dgy', 'Dgz'};
+
+%% Save the obtained transfer functions
+save('./mat/id_micro_station_flexibility.mat', 'G_ms_flexible', 'G_ms_ty_rigid');
diff --git a/identification/id_flexible_rigid_plots.m b/identification/id_flexible_rigid_plots.m
new file mode 100644
index 0000000..6e8307e
--- /dev/null
+++ b/identification/id_flexible_rigid_plots.m
@@ -0,0 +1,99 @@
+%% Script Description
+% Determine if we take into account the flexibilities,
+% does that changes a lot
+
+%%
+clear; close all; clc;
+
+%% Load Configuration file
+load('./mat/config.mat', 'save_fig', 'freqs');
+
+%% Load the obtained transfer functions
+load('./mat/id_micro_station_flexibility.mat', 'G_ms_flexible', 'G_ms_ty_rigid');
+
+%% Get Measurement Object
+load('2018_01_12.mat', 'm_object');
+
+% Get Measurements Data
+opts = struct('freq_min', 10, 'est_backend', 'idfrd');
+meas_sys = getDynamicTFs(m_object, 'marble', 'hexa', {{'tx', 'tx'},{'ty', 'ty'},{'tz', 'tz'}}, opts);
+
+%%
+dir = 'y';
+
+figure;
+% Amplitude
+ax1 = subaxis(2,1,1);
+hold on;
+plot(freqs, abs(squeeze(freqresp(G_ms_flexible(['Dg' dir], ['Fg' dir]), freqs, 'Hz'))));
+plot(freqs, abs(squeeze(freqresp(G_ms_ty_rigid(['Dg' dir], ['Fg' dir]), freqs, 'Hz'))), '--');
+set(gca,'xscale','log'); set(gca,'yscale','log');
+ylabel('Amplitude [m/N]');
+set(gca, 'XTickLabel',[]);
+legend({'Flexible', 'Ty - Rigid'});
+hold off;
+% Phase
+ax2 = subaxis(2,1,2);
+hold on;
+plot(freqs, 180/pi*angle(squeeze(freqresp(G_ms_flexible(['Dg' dir], ['Fg' dir]), freqs, 'Hz'))));
+plot(freqs, 180/pi*angle(squeeze(freqresp(G_ms_ty_rigid(['Dg' dir], ['Fg' dir]), freqs, 'Hz'))), '--');
+set(gca,'xscale','log');
+ylim([-180, 180]);
+yticks([-180, -90, 0, 90, 180]);
+xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+hold off;
+linkaxes([ax1,ax2],'x');
+
+%%
+dir = 'y';
+
+figure;
+% Amplitude
+ax1 = subaxis(2,1,1);
+hold on;
+plot(freqs, abs(squeeze(freqresp(G_ms_flexible(['Dm' dir], ['Fm' dir]), freqs, 'Hz'))));
+plot(freqs, abs(squeeze(freqresp(G_ms_ty_rigid(['Dm' dir], ['Fm' dir]), freqs, 'Hz'))), '--');
+set(gca,'xscale','log'); set(gca,'yscale','log');
+ylabel('Amplitude [m/N]');
+set(gca, 'XTickLabel',[]);
+legend({'Flexible', 'Ty - Rigid'});
+hold off;
+% Phase
+ax2 = subaxis(2,1,2);
+hold on;
+plot(freqs, 180/pi*angle(squeeze(freqresp(G_ms_flexible(['Dm' dir], ['Fm' dir]), freqs, 'Hz'))));
+plot(freqs, 180/pi*angle(squeeze(freqresp(G_ms_ty_rigid(['Dm' dir], ['Fm' dir]), freqs, 'Hz'))), '--');
+set(gca,'xscale','log');
+ylim([-180, 180]);
+yticks([-180, -90, 0, 90, 180]);
+xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+hold off;
+linkaxes([ax1,ax2],'x');
+
+%%
+dir = 'z';
+
+figure;
+% Amplitude
+ax1 = subaxis(2,1,1);
+hold on;
+plot(freqs, abs(squeeze(freqresp(G_ms_flexible(['Dg' dir], ['Fm' dir]), freqs, 'Hz'))));
+plot(freqs, abs(squeeze(freqresp(G_ms_ty_rigid(['Dg' dir], ['Fm' dir]), freqs, 'Hz'))), '--');
+plot(freqs, abs(squeeze(freqresp(meas_sys(['Dm' dir], ['Fh' dir]), freqs, 'Hz'))), '.');
+set(gca,'xscale','log'); set(gca,'yscale','log');
+ylabel('Amplitude [m/N]');
+set(gca, 'XTickLabel',[]);
+legend({'Flexible', 'Ty - Rigid'});
+hold off;
+% Phase
+ax2 = subaxis(2,1,2);
+hold on;
+plot(freqs, 180/pi*angle(squeeze(freqresp(G_ms_flexible(['Dg' dir], ['Fm' dir]), freqs, 'Hz'))));
+plot(freqs, 180/pi*angle(squeeze(freqresp(G_ms_ty_rigid(['Dg' dir], ['Fm' dir]), freqs, 'Hz'))), '--');
+plot(freqs, 180/pi*angle(squeeze(freqresp(meas_sys(['Dm' dir], ['Fh' dir]), freqs, 'Hz'))), '.');
+set(gca,'xscale','log');
+ylim([-180, 180]);
+yticks([-180, -90, 0, 90, 180]);
+xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+hold off;
+linkaxes([ax1,ax2],'x');
diff --git a/identification/id_micro_station_comp_meas.m b/identification/id_micro_station_comp_meas.m
index 6a97556..6c02a84 100644
--- a/identification/id_micro_station_comp_meas.m
+++ b/identification/id_micro_station_comp_meas.m
@@ -16,85 +16,90 @@ load('2018_01_12.mat', 'm_object');
 
 %% Get Measurements Data
 opts = struct('freq_min', 10, 'est_backend', 'idfrd');
-meas_sys = getDynamicTFs(m_object, 'marble', 'hexa', {'tx', 'tx'}, opts);
+meas_sys = getDynamicTFs(m_object, 'marble', 'hexa', {{'tx', 'tx'},{'ty', 'ty'},{'tz', 'tz'}}, opts);
 
 %% Granite to Granite
-figure;
-% Amplitude
-ax1 = subaxis(2,1,1);
-hold on;
-plot(freqs, abs(squeeze(freqresp(G_ms('Dgx', 'Fgx'), freqs, 'Hz'))));
-plot(freqs, abs(squeeze(freqresp(meas_sys('Dmx', 'Fmx'), freqs, 'Hz'))), '.');
-set(gca,'xscale','log'); set(gca,'yscale','log');
-ylabel('Amplitude [m/N]');
-set(gca, 'XTickLabel',[]);
-legend({'Model', 'Meas.'});
-hold off;
-% Phase
-ax2 = subaxis(2,1,2);
-hold on;
-plot(freqs, 180/pi*angle(squeeze(freqresp(G_ms('Dgz', 'Fgz'), freqs, 'Hz'))));
-plot(freqs, 180/pi*angle(squeeze(freqresp(meas_sys('Dmx', 'Fmx'), freqs, 'Hz'))), '.');
-set(gca,'xscale','log');
-ylim([-180, 180]);
-yticks([-180, -90, 0, 90, 180]);
-xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
-hold off;
+for dir = 'xyz'
+    figure;
+    % Amplitude
+    ax1 = subaxis(2,1,1);
+    hold on;
+    plot(freqs, abs(squeeze(freqresp(G_ms(['Dg' dir], ['Fg' dir]), freqs, 'Hz'))));
+    plot(freqs, abs(squeeze(freqresp(meas_sys(['Dm' dir], ['Fm' dir]), freqs, 'Hz'))), '.');
+    set(gca,'xscale','log'); set(gca,'yscale','log');
+    ylabel('Amplitude [m/N]');
+    set(gca, 'XTickLabel',[]);
+    legend({'Model', 'Meas.'});
+    hold off;
+    % Phase
+    ax2 = subaxis(2,1,2);
+    hold on;
+    plot(freqs, 180/pi*angle(squeeze(freqresp(G_ms(['Dg' dir], ['Fg' dir]), freqs, 'Hz'))));
+    plot(freqs, 180/pi*angle(squeeze(freqresp(meas_sys(['Dm' dir], ['Fm' dir]), freqs, 'Hz'))), '.');
+    set(gca,'xscale','log');
+    ylim([-180, 180]);
+    yticks([-180, -90, 0, 90, 180]);
+    xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+    hold off;
+    linkaxes([ax1,ax2],'x');
 
-linkaxes([ax1,ax2],'x');
-
-if save_fig; exportFig('comp_meas_g_g', 'normal-normal', struct('path', 'identification')); end
+    if save_fig; exportFig(['comp_meas_g_g_' dir], 'normal-normal', struct('path', 'identification')); end
+end
 
 %% Hexapod to Hexapod
-figure;
-% Amplitude
-ax1 = subaxis(2,1,1);
-hold on;
-plot(freqs, abs(squeeze(freqresp(G_ms('Dmx', 'Fmx'), freqs, 'Hz'))));
-plot(freqs, abs(squeeze(freqresp(meas_sys('Dhx', 'Fhx'), freqs, 'Hz'))), '.');
-set(gca,'xscale','log'); set(gca,'yscale','log');
-ylabel('Amplitude [m/N]');
-set(gca, 'XTickLabel',[]);
-legend({'Model', 'Meas.'});
-hold off;
-% Phase
-ax2 = subaxis(2,1,2);
-hold on;
-plot(freqs, 180/pi*angle(squeeze(freqresp(G_ms('Dmx', 'Fmx'), freqs, 'Hz'))));
-plot(freqs, 180/pi*angle(squeeze(freqresp(meas_sys('Dhx', 'Fhx'), freqs, 'Hz'))), '.');
-set(gca,'xscale','log');
-ylim([-180, 180]);
-yticks([-180, -90, 0, 90, 180]);
-xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
-hold off;
+for dir = 'xyz'
+    figure;
+    % Amplitude
+    ax1 = subaxis(2,1,1);
+    hold on;
+    plot(freqs, abs(squeeze(freqresp(G_ms(['Dm' dir], ['Fm' dir]), freqs, 'Hz'))));
+    plot(freqs, abs(squeeze(freqresp(meas_sys(['Dh' dir], ['Fh' dir]), freqs, 'Hz'))), '.');
+    set(gca,'xscale','log'); set(gca,'yscale','log');
+    ylabel('Amplitude [m/N]');
+    set(gca, 'XTickLabel',[]);
+    legend({'Model', 'Meas.'});
+    hold off;
+    % Phase
+    ax2 = subaxis(2,1,2);
+    hold on;
+    plot(freqs, 180/pi*angle(squeeze(freqresp(G_ms(['Dm' dir], ['Fm' dir]), freqs, 'Hz'))));
+    plot(freqs, 180/pi*angle(squeeze(freqresp(meas_sys(['Dh' dir], ['Fh' dir]), freqs, 'Hz'))), '.');
+    set(gca,'xscale','log');
+    ylim([-180, 180]);
+    yticks([-180, -90, 0, 90, 180]);
+    xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+    hold off;
 
-linkaxes([ax1,ax2],'x');
-
-if save_fig; exportFig('comp_meas_m_m', 'normal-normal', struct('path', 'identification')); end
+    linkaxes([ax1,ax2],'x');
+    
+    if save_fig; exportFig(['comp_meas_m_m_' dir], 'normal-normal', struct('path', 'identification')); end
+end
 
 %% Hexapod to Granite
-figure;
-% Amplitude
-ax1 = subaxis(2,1,1);
-hold on;
-plot(freqs, abs(squeeze(freqresp(G_ms('Dmx', 'Fgx'), freqs, 'Hz'))));
-plot(freqs, abs(squeeze(freqresp(meas_sys('Dhx', 'Fmx'), freqs, 'Hz'))), '.');
-set(gca,'xscale','log'); set(gca,'yscale','log');
-ylabel('Amplitude [m/N]');
-set(gca, 'XTickLabel',[]);
-legend({'Model', 'Meas.'});
-hold off;
-% Phase
-ax2 = subaxis(2,1,2);
-hold on;
-plot(freqs, 180/pi*angle(squeeze(freqresp(G_ms('Dmx', 'Fgx'), freqs, 'Hz'))));
-plot(freqs, 180/pi*angle(squeeze(freqresp(meas_sys('Dhx', 'Fmx'), freqs, 'Hz'))), '.');
-set(gca,'xscale','log');
-ylim([-180, 180]);
-yticks([-180, -90, 0, 90, 180]);
-xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
-hold off;
+for dir = 'xyz'
+    figure;
+    % Amplitude
+    ax1 = subaxis(2,1,1);
+    hold on;
+    plot(freqs, abs(squeeze(freqresp(G_ms(['Dm' dir], ['Fg' dir]), freqs, 'Hz'))));
+    plot(freqs, abs(squeeze(freqresp(meas_sys(['Dh' dir], ['Fm' dir]), freqs, 'Hz'))), '.');
+    set(gca,'xscale','log'); set(gca,'yscale','log');
+    ylabel('Amplitude [m/N]');
+    set(gca, 'XTickLabel',[]);
+    legend({'Model', 'Meas.'});
+    hold off;
+    % Phase
+    ax2 = subaxis(2,1,2);
+    hold on;
+    plot(freqs, 180/pi*angle(squeeze(freqresp(G_ms(['Dm' dir], ['Fg' dir]), freqs, 'Hz'))));
+    plot(freqs, 180/pi*angle(squeeze(freqresp(meas_sys(['Dh' dir], ['Fm' dir]), freqs, 'Hz'))), '.');
+    set(gca,'xscale','log');
+    ylim([-180, 180]);
+    yticks([-180, -90, 0, 90, 180]);
+    xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+    hold off;
 
-linkaxes([ax1,ax2],'x');
+    linkaxes([ax1,ax2],'x');
 
-if save_fig; exportFig('comp_meas_m_g', 'normal-normal', struct('path', 'identification')); end
+    if save_fig; exportFig(['comp_meas_m_g_' dir], 'normal-normal', struct('path', 'identification')); end
+end
diff --git a/identification/id_nano_station.m b/identification/id_nano_station.m
index 2e10216..cd343e6 100644
--- a/identification/id_nano_station.m
+++ b/identification/id_nano_station.m
@@ -3,7 +3,7 @@ clear; close all; clc;
 
 %%
 K_iff = tf(zeros(6));
-save('./mat/K_iff.mat', 'K_iff');
+save('./mat/controllers.mat', 'K_iff', '-append');
 
 %% Light Sample
 initializeSample(struct('mass', 1));
diff --git a/identification/sim_micro_station.slx b/identification/sim_micro_station.slx
index b124c9b..a9cd0d8 100644
Binary files a/identification/sim_micro_station.slx and b/identification/sim_micro_station.slx differ
diff --git a/init_data.m b/init_data.m
index f4f8cbf..1994e97 100644
--- a/init_data.m
+++ b/init_data.m
@@ -52,6 +52,7 @@ initializeSample(struct('mass', 20));
 
 %% Controllers
 K = tf(zeros(6));
-save('./mat/controller.mat', 'K');
+
 K_iff = tf(zeros(6));
-save('./mat/K_iff.mat', 'K_iff');
\ No newline at end of file
+
+save('./mat/controllers.mat', 'K', 'K_iff');
diff --git a/init_perturbations.m b/init_perturbations.m
index f9e4498..65f8114 100644
--- a/init_perturbations.m
+++ b/init_perturbations.m
@@ -9,9 +9,8 @@ Wxg = 1e-5*(s/(2e2)^(1/3) + 2*pi*0.1)^3/(s + 2*pi*0.1)^3;
 Wxg = Wxg*(s/(0.5e6)^(1/3) + 2*pi*10)^3/(s + 2*pi*10)^3;
 Wxg = Wxg/(1+s/(2*pi*2000));
 
-save('./mat/weight_Wxg.mat', 'Wxg');
-
 %% Sensor Noise
 Wn = tf(1e-12);
 
-save('./mat/weight_Wn.mat', 'Wn');
+%% Save Weights
+save('./mat/perturbations.mat', 'Wxg', 'Wn');
\ No newline at end of file
diff --git a/init_simulation.m b/init_simulation.m
index cb4c7c2..ce08643 100644
--- a/init_simulation.m
+++ b/init_simulation.m
@@ -1,17 +1,17 @@
 %% Script that is run just before
 % the simulation is started
-% Load all the data used for the simulation
-load('./mat/sim_conf.mat', 'sim_conf');
+
+%% Load all the data used for the simulation
+load('./mat/sim_conf.mat');
 
 %% Load SolidWorks Data
-load('./mat/solids.mat', 'solids');
+load('./mat/solids.mat');
 
 %% Load data of each stage
-load('./mat/stages.mat', 'ground', 'granite', 'ty', 'ry', 'rz', 'micro_hexapod', 'axisc', 'nano_hexapod', 'mirror', 'sample');
+load('./mat/stages.mat');
 
 %% Load Signals Applied to the system
-load('./mat/inputs.mat', 'inputs');
+load('./mat/inputs.mat');
 
 %% Load Controller
-load('./mat/controller.mat', 'K');
-load('./mat/K_iff.mat', 'K_iff');
+load('./mat/controllers.mat');
diff --git a/initialize/initializeInputs.m b/initialize/initializeInputs.m
index a4aab4f..16d7eca 100644
--- a/initialize/initializeInputs.m
+++ b/initialize/initializeInputs.m
@@ -28,7 +28,7 @@ function [inputs] = initializeInputs(opts_param)
 
     %% Ground motion
     if islogical(opts.Dw) && opts.Dw == true
-        load('./mat/weight_Wxg.mat', 'Wxg');
+        load('./mat/perturbations.mat', 'Wxg');
         Dw = 1/sqrt(2)*100*random('norm', 0, 1, length(time_vector), 3);
         Dw(:, 1) = lsim(Wxg, Dw(:, 1), time_vector);
         Dw(:, 2) = lsim(Wxg, Dw(:, 2), time_vector);
@@ -64,17 +64,17 @@ function [inputs] = initializeInputs(opts_param)
     inputs.ry = timeseries(ry, time_vector);
 
     %% Spindle [rad]
-    if islogical(opts.Rz) && opts.Rz == true
-        Rz = 2*pi*0.5*time_vector;
-    elseif islogical(opts.Rz) && opts.Rz == false
-        Rz = zeros(length(time_vector), 1);
-    elseif isnumeric(opts.Rz) && length(opts.Rz) == 1
-        Rz = 2*pi*(opts.Rz/60)*time_vector;
+    if islogical(opts.rz) && opts.rz == true
+        rz = 2*pi*0.5*time_vector;
+    elseif islogical(opts.rz) && opts.rz == false
+        rz = zeros(length(time_vector), 1);
+    elseif isnumeric(opts.rz) && length(opts.rz) == 1
+        rz = 2*pi*(opts.rz/60)*time_vector;
     else
-        Rz = opts.Rz;
+        rz = opts.rz;
     end
 
-    inputs.Rz = timeseries(Rz, time_vector);
+    inputs.rz = timeseries(rz, time_vector);
 
     %% Micro Hexapod
     if islogical(opts.u_hexa) && opts.setpoint == true
@@ -85,19 +85,19 @@ function [inputs] = initializeInputs(opts_param)
         u_hexa = opts.u_hexa;
     end
 
-    inputs.micro_hexapod = timeseries(u_hexa, time_vector);
+    inputs.u_hexa = timeseries(u_hexa, time_vector);
 
     %% Center of gravity compensation
     if islogical(opts.mass) && opts.setpoint == true
-        Rm = zeros(length(time_vector), 2);
+        axisc = zeros(length(time_vector), 2);
     elseif islogical(opts.mass) && opts.setpoint == false
-        Rm = zeros(length(time_vector), 2);
-        Rm(:, 2) = pi*ones(length(time_vector), 1);
+        axisc = zeros(length(time_vector), 2);
+        axisc(:, 2) = pi*ones(length(time_vector), 1);
     else
-        Rm = opts.mass;
+        axisc = opts.mass;
     end
 
-    inputs.Rm = timeseries(Rm, time_vector);
+    inputs.axisc = timeseries(axisc, time_vector);
 
     %% Nano Hexapod
     if islogical(opts.n_hexa) && opts.setpoint == true
@@ -108,7 +108,7 @@ function [inputs] = initializeInputs(opts_param)
         n_hexa = opts.n_hexa;
     end
 
-    inputs.nano_hexapod = timeseries(n_hexa, time_vector);
+    inputs.n_hexa = timeseries(n_hexa, time_vector);
 
     %% Set point [m, rad]
     if islogical(opts.setpoint) && opts.setpoint == true
diff --git a/initialize/initializeMicroHexapod.m b/initialize/initializeMicroHexapod.m
index c389dc6..899d61f 100644
--- a/initialize/initializeMicroHexapod.m
+++ b/initialize/initializeMicroHexapod.m
@@ -12,9 +12,9 @@ function [] = initializeMicroHexapod(opts_param)
     %% Stewart Object
     micro_hexapod = struct();
     micro_hexapod.h        = 350; % Total height of the platform [mm]
-    micro_hexapod.jacobian = 265; % Point where the Jacobian is computed => Center of rotation [mm]
+    micro_hexapod.jacobian = 265; % Distance from the top platform to the Jacobian point [mm]
 
-    %% Bottom Plate
+    %% Bottom Plate - Mechanical Design
     BP = struct();
 
     BP.rad.int   = 110;   % Internal Radius [mm]
@@ -26,7 +26,7 @@ function [] = initializeMicroHexapod(opts_param)
     BP.color     = [0.6 0.6 0.6]; % Color [rgb]
     BP.shape     = [BP.rad.int BP.thickness; BP.rad.int 0; BP.rad.ext 0; BP.rad.ext BP.thickness];
 
-    %% Top Plate
+    %% Top Plate - Mechanical Design
     TP = struct();
 
     TP.rad.int   = 82;   % Internal Radius [mm]
@@ -38,7 +38,7 @@ function [] = initializeMicroHexapod(opts_param)
     TP.color     = [0.6 0.6 0.6]; % Color [rgb]
     TP.shape     = [TP.rad.int TP.thickness; TP.rad.int 0; TP.rad.ext 0; TP.rad.ext TP.thickness];
 
-    %% Leg
+    %% Struts
     Leg = struct();
 
     Leg.stroke     = 10e-3; % Maximum Stroke of each leg [m]
diff --git a/initialize/initializeRy.m b/initialize/initializeRy.m
index 10f5ba5..8fc923f 100644
--- a/initialize/initializeRy.m
+++ b/initialize/initializeRy.m
@@ -1,10 +1,24 @@
-function [ry] = initializeRy()
+function [ry] = initializeRy(opts_param)
+    %% Default values for opts
+    opts = struct('rigid', false);
+
+    %% Populate opts with input parameters
+    if exist('opts_param','var')
+        for opt = fieldnames(opts_param)'
+            opts.(opt{1}) = opts_param.(opt{1});
+        end
+    end
+
     %%
     ry = struct();
 
     ry.m = 200; % [kg]
 
-    ry.k.tilt = 1e4; % Rotation stiffness around y [N*m/deg]
+    if opts.rigid
+        ry.k.tilt = 1e10; % Rotation stiffness around y [N*m/deg]
+    else
+        ry.k.tilt = 1e4; % Rotation stiffness around y [N*m/deg]
+    end
 
     ry.k.h    = 357e6/4; % Stiffness in the direction of the guidance [N/m]
     ry.k.rad  = 555e6/4; % Stiffness in the top direction [N/m]
diff --git a/initialize/initializeRz.m b/initialize/initializeRz.m
index 76d81ef..f0ef2c9 100644
--- a/initialize/initializeRz.m
+++ b/initialize/initializeRz.m
@@ -1,18 +1,36 @@
-function [rz] = initializeRz()
+function [rz] = initializeRz(opts_param)
+    %% Default values for opts
+    opts = struct('rigid', false);
+
+    %% Populate opts with input parameters
+    if exist('opts_param','var')
+        for opt = fieldnames(opts_param)'
+            opts.(opt{1}) = opts_param.(opt{1});
+        end
+    end
+
     %%
     rz = struct();
 
+    % Estimated mass of the mooving part
     rz.m = 250; % [kg]
 
+    % Estimated stiffnesses
     rz.k.ax   = 2e9; % Axial Stiffness [N/m]
     rz.k.rad  = 7e8; % Radial Stiffness [N/m]
-    rz.k.rot  = 1e2; % Rotational Stiffness [N*m/deg]
-    rz.k.tilt = 1e2; % TODO [N*m/deg]
+    rz.k.rot  = 100e6*2*pi/360; % Rotational Stiffness [N*m/deg]
 
-    rz.c.ax   = 10*(1/5)*sqrt(rz.k.ax/rz.m);
-    rz.c.rad  = 10*(1/5)*sqrt(rz.k.rad/rz.m);
-    rz.c.tilt = 100*(1/1)*sqrt(rz.k.tilt/rz.m);
-    rz.c.rot  = 100*(1/1)*sqrt(rz.k.rot/rz.m);
+    if opts.rigid
+        rz.k.tilt = 1e10; % Vertical Rotational Stiffness [N*m/deg]
+    else
+        rz.k.tilt = 1e2; % TODO what value should I put? [N*m/deg]
+    end
+
+    % TODO
+    rz.c.ax   = 2*sqrt(rz.k.ax/rz.m);
+    rz.c.rad  = 2*sqrt(rz.k.rad/rz.m);
+    rz.c.tilt = 100*sqrt(rz.k.tilt/rz.m);
+    rz.c.rot  = 100*sqrt(rz.k.rot/rz.m);
 
     %% Save
     save('./mat/stages.mat', 'rz', '-append');
diff --git a/initialize/initializeSample.m b/initialize/initializeSample.m
index c2336b1..51e30d2 100644
--- a/initialize/initializeSample.m
+++ b/initialize/initializeSample.m
@@ -1,9 +1,9 @@
 function [] = initializeSample(opts_param)
     %% Default values for opts
-    sample = struct('radius', 100,...
-                    'height', 300,...
-                    'mass',   50,...
-                    'offset', 0,...
+    sample = struct('radius', 100, ...
+                    'height', 300, ...
+                    'mass',   50,  ...
+                    'offset', 0,   ...
                     'color',  [0.45, 0.45, 0.45] ...
     );
 
diff --git a/initialize/initializeTy.m b/initialize/initializeTy.m
index 7cef6e5..a3e7fab 100644
--- a/initialize/initializeTy.m
+++ b/initialize/initializeTy.m
@@ -1,10 +1,24 @@
-function [ty] = initializeTy()
+function [ty] = initializeTy(opts_param)
+    %% Default values for opts
+    opts = struct('rigid', false);
+
+    %% Populate opts with input parameters
+    if exist('opts_param','var')
+        for opt = fieldnames(opts_param)'
+            opts.(opt{1}) = opts_param.(opt{1});
+        end
+    end
+
     %%
     ty = struct();
 
     ty.m = 250; % [kg]
 
-    ty.k.ax  = 1e7/4; % Axial Stiffness for each of the 4 guidance (y) [N/m]
+    if opts.rigid
+        ty.k.ax  = 1e10; % Axial Stiffness for each of the 4 guidance (y) [N/m]
+    else
+        ty.k.ax  = 1e7/4; % Axial Stiffness for each of the 4 guidance (y) [N/m]
+    end
     ty.k.rad = 9e9/4; % Radial Stiffness for each of the 4 guidance (x-z) [N/m]
 
     ty.c.ax  = 100*(1/5)*sqrt(ty.k.ax/ty.m);
diff --git a/mat/config.mat b/mat/config.mat
index b2c9dc7..632d562 100644
Binary files a/mat/config.mat and b/mat/config.mat differ
diff --git a/mat/controllers.mat b/mat/controllers.mat
new file mode 100644
index 0000000..600a1b3
Binary files /dev/null and b/mat/controllers.mat differ
diff --git a/mat/id_micro_station.mat b/mat/id_micro_station.mat
index 8155ac3..6272084 100644
Binary files a/mat/id_micro_station.mat and b/mat/id_micro_station.mat differ
diff --git a/mat/id_micro_station_flexibility.mat b/mat/id_micro_station_flexibility.mat
new file mode 100644
index 0000000..e833e42
Binary files /dev/null and b/mat/id_micro_station_flexibility.mat differ
diff --git a/mat/inputs.mat b/mat/inputs.mat
index 6b7eb53..0b14007 100644
Binary files a/mat/inputs.mat and b/mat/inputs.mat differ
diff --git a/mat/perturbations.mat b/mat/perturbations.mat
new file mode 100644
index 0000000..6606250
Binary files /dev/null and b/mat/perturbations.mat differ
diff --git a/mat/sim_conf.mat b/mat/sim_conf.mat
index 8aaac97..d187942 100644
Binary files a/mat/sim_conf.mat and b/mat/sim_conf.mat differ
diff --git a/mat/solids.mat b/mat/solids.mat
index 334d828..d948fee 100644
Binary files a/mat/solids.mat and b/mat/solids.mat differ
diff --git a/mat/stages.mat b/mat/stages.mat
index 928fc2b..6c998a4 100644
Binary files a/mat/stages.mat and b/mat/stages.mat differ
diff --git a/mat/weight_Wxg.mat b/mat/weight_Wxg.mat
deleted file mode 100644
index a02769c..0000000
Binary files a/mat/weight_Wxg.mat and /dev/null differ
diff --git a/micro_station_test.slx b/micro_station_test.slx
deleted file mode 100644
index dbf6359..0000000
Binary files a/micro_station_test.slx and /dev/null differ
diff --git a/src/runSimulation.m b/src/runSimulation.m
index 40d7466..7959a28 100644
--- a/src/runSimulation.m
+++ b/src/runSimulation.m
@@ -19,10 +19,10 @@ function [] = runSimulation(sys_name, sys_mass, ctrl_type, act_damp)
     if strcmp(act_damp, 'iff')
         K_iff_crit = load('./mat/K_iff_crit.mat');
         K_iff = K_iff_crit.(sprintf('K_iff_%s_%s', sys_mass, sys_name)); %#ok
-        save('./mat/K_iff.mat', 'K_iff');
+        save('./mat/controllers.mat', 'K_iff', '-append');
     elseif strcmp(act_damp, 'none')
         K_iff = tf(zeros(6)); %#ok
-        save('./mat/K_iff.mat', 'K_iff');
+        save('./mat/controllers.mat', 'K_iff', '-append');
     end
 
     %%
@@ -33,7 +33,7 @@ function [] = runSimulation(sys_name, sys_mass, ctrl_type, act_damp)
     else
         error('sys_name should be pz or vc');
     end
-    
+
     if strcmp(sys_mass, 'light')
     	initializeSample(struct('mass', 1));
     elseif strcmp(sys_mass, 'heavy')