diff --git a/matlab/src/initializeDisturbances.m b/matlab/src/initializeDisturbances.m
index 3e7cbc1..37a4add 100644
--- a/matlab/src/initializeDisturbances.m
+++ b/matlab/src/initializeDisturbances.m
@@ -1,183 +1,229 @@
-  function [] = initializeDisturbances(args)
-  % initializeDisturbances - Initialize the disturbances
-  %
-  % Syntax: [] = initializeDisturbances(args)
-  %
-  % Inputs:
-  %    - args -
+function [] = initializeDisturbances(args)
+% initializeDisturbances - Initialize the disturbances
+%
+% Syntax: [] = initializeDisturbances(args)
+%
+% Inputs:
+%    - args -
 
-  arguments
-      % Global parameter to enable or disable the disturbances
-      args.enable logical {mustBeNumericOrLogical} = true
-      % Ground Motion - X direction
-      args.Dwx logical {mustBeNumericOrLogical} = true
-      % Ground Motion - Y direction
-      args.Dwy logical {mustBeNumericOrLogical} = true
-      % Ground Motion - Z direction
-      args.Dwz logical {mustBeNumericOrLogical} = true
-      % Translation Stage - X direction
-      args.Fty_x logical {mustBeNumericOrLogical} = true
-      % Translation Stage - Z direction
-      args.Fty_z logical {mustBeNumericOrLogical} = true
-      % Spindle - X direction
-      args.Frz_x logical {mustBeNumericOrLogical} = true
-      % Spindle - Y direction
-      args.Frz_y logical {mustBeNumericOrLogical} = true
-      % Spindle - Z direction
-      args.Frz_z logical {mustBeNumericOrLogical} = true
-  end
+arguments
+    % Global parameter to enable or disable the disturbances
+    args.enable logical {mustBeNumericOrLogical} = true
+    % Ground Motion - X direction
+    args.Dwx logical {mustBeNumericOrLogical} = true
+    % Ground Motion - Y direction
+    args.Dwy logical {mustBeNumericOrLogical} = true
+    % Ground Motion - Z direction
+    args.Dwz logical {mustBeNumericOrLogical} = true
+    % Translation Stage - X direction
+    args.Fty_x logical {mustBeNumericOrLogical} = true
+    % Translation Stage - Z direction
+    args.Fty_z logical {mustBeNumericOrLogical} = true
+    % Spindle - X direction
+    args.Frz_x logical {mustBeNumericOrLogical} = true
+    % Spindle - Y direction
+    args.Frz_y logical {mustBeNumericOrLogical} = true
+    % Spindle - Z direction
+    args.Frz_z logical {mustBeNumericOrLogical} = true
+end
 
-  load('dist_psd.mat', 'dist_f');
+% Initialization of random numbers
+rng("shuffle");
 
-  dist_f.f = dist_f.f(2:end);
-  dist_f.psd_gm = dist_f.psd_gm(2:end);
-  dist_f.psd_ty = dist_f.psd_ty(2:end);
-  dist_f.psd_rz = dist_f.psd_rz(2:end);
+%% Ground Motion
+load('dist_psd.mat', 'dist_f');
 
-  Fs = 2*dist_f.f(end);      % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz]
-  N  = 2*length(dist_f.f);   % Number of Samples match the one of the wanted PSD
-  T0 = N/Fs;                 % Signal Duration [s]
-  df = 1/T0;                 % Frequency resolution of the DFT [Hz]
-                             % Also equal to (dist_f.f(2)-dist_f.f(1))
-  t = linspace(0, T0, N+1)'; % Time Vector [s]
-  Ts = 1/Fs;                 % Sampling Time [s]
+% Frequency Data
+Dw.f = dist_f.f(2:end);
+Dw.psd_x = dist_f.psd_gm(2:end);
+Dw.psd_y = dist_f.psd_gm(2:end);
+Dw.psd_z = dist_f.psd_gm(2:end);
 
-  phi = dist_f.psd_gm;
-  C = zeros(N/2,1);
-  for i = 1:N/2
-    C(i) = sqrt(phi(i)*df);
-  end
+% Time data
+Fs = 2*Dw.f(end);      % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz]
+N  = 2*length(Dw.f);   % Number of Samples match the one of the wanted PSD
+T0 = N/Fs;             % Signal Duration [s]
+Dw.t = linspace(0, T0, N+1)'; % Time Vector [s]
 
-  if args.Dwx && args.enable
-    rng(111);
+C = zeros(N/2,1);
+for i = 1:N/2
+    C(i) = sqrt(Dw.psd_x(i)/T0);
+end
+
+if args.Dwx && args.enable
     theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
     Cx = [0 ; C.*complex(cos(theta),sin(theta))];
     Cx = [Cx; flipud(conj(Cx(2:end)))];;
-    Dwx = N/sqrt(2)*ifft(Cx); % Ground Motion - x direction [m]
-  else
-    Dwx = zeros(length(t), 1);
-  end
+    Dw.x = N/sqrt(2)*ifft(Cx); % Ground Motion - x direction [m]
+else
+    Dw.x = zeros(length(Dw.t), 1);
+end
 
-  if args.Dwy && args.enable
-    rng(112);
+if args.Dwy && args.enable
     theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
     Cx = [0 ; C.*complex(cos(theta),sin(theta))];
     Cx = [Cx; flipud(conj(Cx(2:end)))];;
-    Dwy = N/sqrt(2)*ifft(Cx); % Ground Motion - y direction [m]
-  else
-    Dwy = zeros(length(t), 1);
-  end
+    Dw.y = N/sqrt(2)*ifft(Cx); % Ground Motion - y direction [m]
+else
+    Dw.y = zeros(length(Dw.t), 1);
+end
 
-  if args.Dwy && args.enable
-    rng(113);
+if args.Dwy && args.enable
     theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
     Cx = [0 ; C.*complex(cos(theta),sin(theta))];
     Cx = [Cx; flipud(conj(Cx(2:end)))];;
-    Dwz = N/sqrt(2)*ifft(Cx); % Ground Motion - z direction [m]
-  else
-    Dwz = zeros(length(t), 1);
-  end
+    Dw.z = N/sqrt(2)*ifft(Cx); % Ground Motion - z direction [m]
+else
+    Dw.z = zeros(length(Dw.t), 1);
+end
 
-  if args.Fty_x && args.enable
-    phi = dist_f.psd_ty; % TODO - we take here the vertical direction which is wrong but approximate
+load('dist_psd.mat', 'dist_f');
+dist_f.f = dist_f.f(2:end);
+dist_f.psd_gm = dist_f.psd_gm(2:end);
+dist_f.psd_ty = dist_f.psd_ty(2:end);
+dist_f.psd_rz = dist_f.psd_rz(2:end);
+
+%% Translation Stage
+load('dist_psd.mat', 'dist_f');
+
+% Frequency Data
+Ty.f = dist_f.f(2:end);
+Ty.psd_x = dist_f.psd_ty(2:end); % TODO - we take here the vertical direction which is wrong but approximate
+Ty.psd_z = dist_f.psd_ty(2:end);
+
+% Time data
+Fs = 2*Ty.f(end);      % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz]
+N  = 2*length(Ty.f);   % Number of Samples match the one of the wanted PSD
+T0 = N/Fs;             % Signal Duration [s]
+Ty.t = linspace(0, T0, N+1)'; % Time Vector [s]
+
+C = zeros(N/2,1);
+for i = 1:N/2
+    C(i) = sqrt(Ty.psd_x(i)/T0);
+end
+
+% Translation Stage - X
+if args.Fty_x && args.enable
+    phi = Ty.psd_x;
     C = zeros(N/2,1);
     for i = 1:N/2
-      C(i) = sqrt(phi(i)*df);
+        C(i) = sqrt(phi(i)/T0);
     end
-    rng(121);
     theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
     Cx = [0 ; C.*complex(cos(theta),sin(theta))];
     Cx = [Cx; flipud(conj(Cx(2:end)))];;
     u = N/sqrt(2)*ifft(Cx); % Disturbance Force Ty x [N]
-    Fty_x = u;
-  else
-    Fty_x = zeros(length(t), 1);
-  end
+    Ty.x = u;
+else
+    Ty.x = zeros(length(Ty.t), 1);
+end
 
-  if args.Fty_z && args.enable
-    phi = dist_f.psd_ty;
+% Translation Stage - Z
+if args.Fty_z && args.enable
+    phi = Ty.psd_z;
     C = zeros(N/2,1);
     for i = 1:N/2
-      C(i) = sqrt(phi(i)*df);
+        C(i) = sqrt(phi(i)/T0);
     end
-    rng(122);
     theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
     Cx = [0 ; C.*complex(cos(theta),sin(theta))];
     Cx = [Cx; flipud(conj(Cx(2:end)))];;
     u = N/sqrt(2)*ifft(Cx); % Disturbance Force Ty z [N]
-    Fty_z = u;
-  else
-    Fty_z = zeros(length(t), 1);
-  end
+    Ty.z = u;
+else
+    Ty.z = zeros(length(Ty.t), 1);
+end
 
-  % if args.Frz_x && args.enable
-  %   phi = dist_f.psd_rz;
-  %   C = zeros(N/2,1);
-  %   for i = 1:N/2
-  %     C(i) = sqrt(phi(i)*df);
-  %   end
-  %   rng(131);
-  %   theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
-  %   Cx = [0 ; C.*complex(cos(theta),sin(theta))];
-  %   Cx = [Cx; flipud(conj(Cx(2:end)))];;
-  %   u = N/sqrt(2)*ifft(Cx); % Disturbance Force Rz z [N]
-  %   Frz_x = u;
-  % else
-    Frz_x = zeros(length(t), 1);
-  % end
+%% Translation Stage
+load('dist_psd.mat', 'dist_f');
 
-  % if args.Frz_y && args.enable
-  %   phi = dist_f.psd_rz;
-  %   C = zeros(N/2,1);
-  %   for i = 1:N/2
-  %     C(i) = sqrt(phi(i)*df);
-  %   end
-  %   rng(131);
-  %   theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
-  %   Cx = [0 ; C.*complex(cos(theta),sin(theta))];
-  %   Cx = [Cx; flipud(conj(Cx(2:end)))];;
-  %   u = N/sqrt(2)*ifft(Cx); % Disturbance Force Rz z [N]
-  %   Frz_z = u;
-  % else
-    Frz_y = zeros(length(t), 1);
-  % end
+% Frequency Data
+Rz.f = dist_f.f(2:end);
+Rz.psd_x = dist_f.psd_rz(2:end); % TODO - we take here the vertical direction which is wrong but approximate
+Rz.psd_y = dist_f.psd_rz(2:end); % TODO - we take here the vertical direction which is wrong but approximate
+Rz.psd_z = dist_f.psd_rz(2:end);
 
-  if args.Frz_z && args.enable
-    phi = dist_f.psd_rz;
+% Time data
+Fs = 2*Rz.f(end);      % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz]
+N  = 2*length(Rz.f);   % Number of Samples match the one of the wanted PSD
+T0 = N/Fs;             % Signal Duration [s]
+Rz.t = linspace(0, T0, N+1)'; % Time Vector [s]
+
+C = zeros(N/2,1);
+for i = 1:N/2
+    C(i) = sqrt(Rz.psd_x(i)/T0);
+end
+
+% Translation Stage - X
+if args.Frz_x && args.enable
+    phi = Rz.psd_x;
     C = zeros(N/2,1);
     for i = 1:N/2
-      C(i) = sqrt(phi(i)*df);
+        C(i) = sqrt(phi(i)/T0);
+    end
+    theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
+    Cx = [0 ; C.*complex(cos(theta),sin(theta))];
+    Cx = [Cx; flipud(conj(Cx(2:end)))];;
+    u = N/sqrt(2)*ifft(Cx); % Disturbance Force Rz x [N]
+    Rz.x = u;
+else
+    Rz.x = zeros(length(Rz.t), 1);
+end
+
+% Translation Stage - Y
+if args.Frz_y && args.enable
+    phi = Rz.psd_y;
+    C = zeros(N/2,1);
+    for i = 1:N/2
+        C(i) = sqrt(phi(i)/T0);
+    end
+    theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
+    Cx = [0 ; C.*complex(cos(theta),sin(theta))];
+    Cx = [Cx; flipud(conj(Cx(2:end)))];;
+    u = N/sqrt(2)*ifft(Cx); % Disturbance Force Rz y [N]
+    Rz.y = u;
+else
+    Rz.y = zeros(length(Rz.t), 1);
+end
+
+% Translation Stage - Z
+if args.Frz_z && args.enable
+    phi = Rz.psd_z;
+    C = zeros(N/2,1);
+    for i = 1:N/2
+        C(i) = sqrt(phi(i)/T0);
     end
-    rng(131);
     theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
     Cx = [0 ; C.*complex(cos(theta),sin(theta))];
     Cx = [Cx; flipud(conj(Cx(2:end)))];;
     u = N/sqrt(2)*ifft(Cx); % Disturbance Force Rz z [N]
-    Frz_z = u;
-  else
-    Frz_z = zeros(length(t), 1);
-  end
+    Rz.z = u;
+else
+    Rz.z = zeros(length(Rz.t), 1);
+end
 
-  u = zeros(length(t), 6);
-  Fd = u;
+u = zeros(100, 6);
+Fd = u;
 
-  Dwx    = Dwx   - Dwx(1);
-  Dwy    = Dwy   - Dwy(1);
-  Dwz    = Dwz   - Dwz(1);
-  Fty_x  = Fty_x - Fty_x(1);
-  Fty_z  = Fty_z - Fty_z(1);
-  Frz_z  = Frz_z - Frz_z(1);
+Dw.x   = Dw.x   - Dw.x(1);
+Dw.y   = Dw.y   - Dw.y(1);
+Dw.z   = Dw.z   - Dw.z(1);
+Ty.x   = Ty.x   - Ty.x(1);
+Ty.z   = Ty.z   - Ty.z(1);
+Rz.x   = Rz.x   - Rz.x(1);
+Rz.y   = Rz.y   - Rz.y(1);
+Rz.z   = Rz.z   - Rz.z(1);
 
 if exist('./mat', 'dir')
     if exist('./mat/nass_disturbances.mat', 'file')
-        save('mat/nass_disturbances.mat', 'Dwx', 'Dwy', 'Dwz', 'Fty_x', 'Fty_z', 'Frz_x', 'Frz_y', 'Frz_z', 'Fd', 'Ts', 't', 'args', '-append');
+        save('mat/nass_disturbances.mat', 'Dw', 'Ty', 'Rz', 'Fd', 'args', '-append');
     else
-        save('mat/nass_disturbances.mat', 'Dwx', 'Dwy', 'Dwz', 'Fty_x', 'Fty_z', 'Frz_x', 'Frz_y', 'Frz_z', 'Fd', 'Ts', 't', 'args');
+        save('mat/nass_disturbances.mat', 'Dw', 'Ty', 'Rz', 'Fd', 'args');
     end
 elseif exist('./matlab', 'dir')
     if exist('./matlab/mat/nass_disturbances.mat', 'file')
-        save('matlab/mat/nass_disturbances.mat', 'Dwx', 'Dwy', 'Dwz', 'Fty_x', 'Fty_z', 'Frz_x', 'Frz_y', 'Frz_z', 'Fd', 'Ts', 't', 'args', '-append');
+        save('matlab/mat/nass_disturbances.mat', 'Dw', 'Ty', 'Rz', 'Fd', 'args', '-append');
     else
-        save('matlab/mat/nass_disturbances.mat', 'Dwx', 'Dwy', 'Dwz', 'Fty_x', 'Fty_z', 'Frz_x', 'Frz_y', 'Frz_z', 'Fd', 'Ts', 't', 'args');
+        save('matlab/mat/nass_disturbances.mat', 'Dw', 'Ty', 'Rz', 'Fd', 'args');
     end
 end
diff --git a/matlab/ustation_simscape.slx b/matlab/ustation_simscape.slx
index bbabff0..eb5a532 100644
Binary files a/matlab/ustation_simscape.slx and b/matlab/ustation_simscape.slx differ
diff --git a/simscape-micro-station.org b/simscape-micro-station.org
index 02bc89f..670c7c3 100644
--- a/simscape-micro-station.org
+++ b/simscape-micro-station.org
@@ -299,8 +299,11 @@ Procedure:
 
 [[*=initializeDisturbances=: Initialize Disturbances][=initializeDisturbances=: Initialize Disturbances]]
 
-- [ ] It is suppose in this script that all disturbances have the same frequency vectors, and therefore the same time vector...
-- [ ] See how to deal with that
+- [X] It is suppose in this script that all disturbances have the same frequency vectors, and therefore the same time vector...
+  It does not anymore
+- [X] See how to deal with that
+
+Be able to pass custom =.mat= files (one mat file per disturbance)?
 
 - [ ] Ground motion, X, Y and Z
 - [ ] Ty stage, X and Z
@@ -2486,8 +2489,8 @@ initializeDisturbances(...
     'Dwz',   true, ... % Ground Motion - Z direction
     'Fty_x', false, ... % Translation Stage - X direction
     'Fty_z', false, ... % Translation Stage - Z direction
-    'Frz_x', false, ... % Spindle - X direction
-    'Frz_y', false, ... % Spindle - Y direction
+    'Frz_x', true, ... % Spindle - X direction
+    'Frz_y', true, ... % Spindle - Y direction
     'Frz_z', true);    % Spindle - Z direction
 
 initializeReferences(...
@@ -2501,9 +2504,9 @@ tomo_align_dist = simout;
 #+begin_src matlab :exports none
 figure;
 hold on;
-plot(tomo_align_dist.y.x.Time, tomo_align_dist.y.x.Data)
-plot(tomo_align_dist.y.y.Time, tomo_align_dist.y.y.Data)
-plot(tomo_align_dist.y.z.Time, tomo_align_dist.y.z.Data)
+plot(tomo_align_dist.y.x.Time, 1e6*tomo_align_dist.y.x.Data)
+plot(tomo_align_dist.y.y.Time, 1e6*tomo_align_dist.y.y.Data)
+plot(tomo_align_dist.y.z.Time, 1e6*tomo_align_dist.y.z.Data)
 hold off;
 #+end_src
 
@@ -2540,18 +2543,16 @@ hold off;
 #+end_src
 
 ** Raster Scans with the translation stage
-<<sec:ty_scans>>
 #+begin_src matlab
-  initializeReferences('Dy_type', 'triangular', 'Dy_amplitude', 10e-3, 'Dy_period', 1);
-  sim(mdl);
-  ty_scan_triangle = simout;
-  % save('./mat/experiment_tomography.mat', 'ty_scan_triangle', '-append');
+initializeReferences('Dy_type', 'triangular', 'Dy_amplitude', 10e-3, 'Dy_period', 1);
+sim(mdl);
+ty_scan_triangle = simout;
+% save('./mat/experiment_tomography.mat', 'ty_scan_triangle', '-append');
 #+end_src
 
 * Conclusion
 <<sec:uniaxial_conclusion>>
 
-
 * Bibliography                                                        :ignore:
 #+latex: \printbibliography[heading=bibintoc,title={Bibliography}]
 
@@ -2979,200 +2980,228 @@ arguments
 end
 #+end_src
 
-*** Load Data
 #+begin_src matlab
-load('dist_psd.mat', 'dist_f');
+% Initialization of random numbers
+rng("shuffle");
 #+end_src
 
+*** Ground Motion
+#+begin_src matlab
+%% Ground Motion
+load('dist_psd.mat', 'dist_f');
+
+% Frequency Data
+Dw.f = dist_f.f(2:end);
+Dw.psd_x = dist_f.psd_gm(2:end);
+Dw.psd_y = dist_f.psd_gm(2:end);
+Dw.psd_z = dist_f.psd_gm(2:end);
+
+% Time data
+Fs = 2*Dw.f(end);      % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz]
+N  = 2*length(Dw.f);   % Number of Samples match the one of the wanted PSD
+T0 = N/Fs;             % Signal Duration [s]
+Dw.t = linspace(0, T0, N+1)'; % Time Vector [s]
+
+C = zeros(N/2,1);
+for i = 1:N/2
+    C(i) = sqrt(Dw.psd_x(i)/T0);
+end
+
+if args.Dwx && args.enable
+    theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
+    Cx = [0 ; C.*complex(cos(theta),sin(theta))];
+    Cx = [Cx; flipud(conj(Cx(2:end)))];;
+    Dw.x = N/sqrt(2)*ifft(Cx); % Ground Motion - x direction [m]
+else
+    Dw.x = zeros(length(Dw.t), 1);
+end
+
+if args.Dwy && args.enable
+    theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
+    Cx = [0 ; C.*complex(cos(theta),sin(theta))];
+    Cx = [Cx; flipud(conj(Cx(2:end)))];;
+    Dw.y = N/sqrt(2)*ifft(Cx); % Ground Motion - y direction [m]
+else
+    Dw.y = zeros(length(Dw.t), 1);
+end
+
+if args.Dwy && args.enable
+    theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
+    Cx = [0 ; C.*complex(cos(theta),sin(theta))];
+    Cx = [Cx; flipud(conj(Cx(2:end)))];;
+    Dw.z = N/sqrt(2)*ifft(Cx); % Ground Motion - z direction [m]
+else
+    Dw.z = zeros(length(Dw.t), 1);
+end
+#+end_src
+
+*** Translation Stage
+
+
 We remove the first frequency point that usually is very large.
 #+begin_src matlab :exports none
+load('dist_psd.mat', 'dist_f');
 dist_f.f = dist_f.f(2:end);
 dist_f.psd_gm = dist_f.psd_gm(2:end);
 dist_f.psd_ty = dist_f.psd_ty(2:end);
 dist_f.psd_rz = dist_f.psd_rz(2:end);
 #+end_src
 
-*** Parameters
-We define some parameters that will be used in the algorithm.
 #+begin_src matlab
-Fs = 2*dist_f.f(end);      % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz]
-N  = 2*length(dist_f.f);   % Number of Samples match the one of the wanted PSD
-T0 = N/Fs;                 % Signal Duration [s]
-df = 1/T0;                 % Frequency resolution of the DFT [Hz]
-                           % Also equal to (dist_f.f(2)-dist_f.f(1))
-t = linspace(0, T0, N+1)'; % Time Vector [s]
-Ts = 1/Fs;                 % Sampling Time [s]
-#+end_src
+%% Translation Stage
+load('dist_psd.mat', 'dist_f');
+
+% Frequency Data
+Ty.f = dist_f.f(2:end);
+Ty.psd_x = dist_f.psd_ty(2:end); % TODO - we take here the vertical direction which is wrong but approximate
+Ty.psd_z = dist_f.psd_ty(2:end);
+
+% Time data
+Fs = 2*Ty.f(end);      % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz]
+N  = 2*length(Ty.f);   % Number of Samples match the one of the wanted PSD
+T0 = N/Fs;             % Signal Duration [s]
+Ty.t = linspace(0, T0, N+1)'; % Time Vector [s]
 
-*** Ground Motion
-#+begin_src matlab
-phi = dist_f.psd_gm;
 C = zeros(N/2,1);
 for i = 1:N/2
-    C(i) = sqrt(phi(i)*df);
+    C(i) = sqrt(Ty.psd_x(i)/T0);
 end
-#+end_src
 
-#+begin_src matlab
-if args.Dwx && args.enable
-    rng(111);
-    theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
-    Cx = [0 ; C.*complex(cos(theta),sin(theta))];
-    Cx = [Cx; flipud(conj(Cx(2:end)))];;
-    Dwx = N/sqrt(2)*ifft(Cx); % Ground Motion - x direction [m]
-else
-    Dwx = zeros(length(t), 1);
-end
-#+end_src
-
-#+begin_src matlab
-if args.Dwy && args.enable
-    rng(112);
-    theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
-    Cx = [0 ; C.*complex(cos(theta),sin(theta))];
-    Cx = [Cx; flipud(conj(Cx(2:end)))];;
-    Dwy = N/sqrt(2)*ifft(Cx); % Ground Motion - y direction [m]
-else
-    Dwy = zeros(length(t), 1);
-end
-#+end_src
-
-#+begin_src matlab
-if args.Dwy && args.enable
-    rng(113);
-    theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
-    Cx = [0 ; C.*complex(cos(theta),sin(theta))];
-    Cx = [Cx; flipud(conj(Cx(2:end)))];;
-    Dwz = N/sqrt(2)*ifft(Cx); % Ground Motion - z direction [m]
-else
-    Dwz = zeros(length(t), 1);
-end
-#+end_src
-
-*** Translation Stage - X direction
-#+begin_src matlab
+% Translation Stage - X
 if args.Fty_x && args.enable
-    phi = dist_f.psd_ty; % TODO - we take here the vertical direction which is wrong but approximate
+    phi = Ty.psd_x;
     C = zeros(N/2,1);
     for i = 1:N/2
-        C(i) = sqrt(phi(i)*df);
+        C(i) = sqrt(phi(i)/T0);
     end
-    rng(121);
     theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
     Cx = [0 ; C.*complex(cos(theta),sin(theta))];
     Cx = [Cx; flipud(conj(Cx(2:end)))];;
     u = N/sqrt(2)*ifft(Cx); % Disturbance Force Ty x [N]
-    Fty_x = u;
+    Ty.x = u;
 else
-    Fty_x = zeros(length(t), 1);
+    Ty.x = zeros(length(Ty.t), 1);
 end
-#+end_src
 
-*** Translation Stage - Z direction
-#+begin_src matlab
+% Translation Stage - Z
 if args.Fty_z && args.enable
-    phi = dist_f.psd_ty;
+    phi = Ty.psd_z;
     C = zeros(N/2,1);
     for i = 1:N/2
-        C(i) = sqrt(phi(i)*df);
+        C(i) = sqrt(phi(i)/T0);
     end
-    rng(122);
     theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
     Cx = [0 ; C.*complex(cos(theta),sin(theta))];
     Cx = [Cx; flipud(conj(Cx(2:end)))];;
     u = N/sqrt(2)*ifft(Cx); % Disturbance Force Ty z [N]
-    Fty_z = u;
+    Ty.z = u;
 else
-    Fty_z = zeros(length(t), 1);
+    Ty.z = zeros(length(Ty.t), 1);
 end
 #+end_src
 
-*** Spindle - X direction
+*** Spindle
 #+begin_src matlab
-% if args.Frz_x && args.enable
-%   phi = dist_f.psd_rz;
-%   C = zeros(N/2,1);
-%   for i = 1:N/2
-%     C(i) = sqrt(phi(i)*df);
-%   end
-%   rng(131);
-%   theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
-%   Cx = [0 ; C.*complex(cos(theta),sin(theta))];
-%   Cx = [Cx; flipud(conj(Cx(2:end)))];;
-%   u = N/sqrt(2)*ifft(Cx); % Disturbance Force Rz z [N]
-%   Frz_x = u;
-% else
-Frz_x = zeros(length(t), 1);
-% end
-#+end_src
+%% Translation Stage
+load('dist_psd.mat', 'dist_f');
 
-*** Spindle - Y direction
-#+begin_src matlab
-% if args.Frz_y && args.enable
-%   phi = dist_f.psd_rz;
-%   C = zeros(N/2,1);
-%   for i = 1:N/2
-%     C(i) = sqrt(phi(i)*df);
-%   end
-%   rng(131);
-%   theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
-%   Cx = [0 ; C.*complex(cos(theta),sin(theta))];
-%   Cx = [Cx; flipud(conj(Cx(2:end)))];;
-%   u = N/sqrt(2)*ifft(Cx); % Disturbance Force Rz z [N]
-%   Frz_z = u;
-% else
-Frz_y = zeros(length(t), 1);
-% end
-#+end_src
+% Frequency Data
+Rz.f = dist_f.f(2:end);
+Rz.psd_x = dist_f.psd_rz(2:end); % TODO - we take here the vertical direction which is wrong but approximate
+Rz.psd_y = dist_f.psd_rz(2:end); % TODO - we take here the vertical direction which is wrong but approximate
+Rz.psd_z = dist_f.psd_rz(2:end);
 
-*** Spindle - Z direction
-#+begin_src matlab
-if args.Frz_z && args.enable
-    phi = dist_f.psd_rz;
+% Time data
+Fs = 2*Rz.f(end);      % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz]
+N  = 2*length(Rz.f);   % Number of Samples match the one of the wanted PSD
+T0 = N/Fs;             % Signal Duration [s]
+Rz.t = linspace(0, T0, N+1)'; % Time Vector [s]
+
+C = zeros(N/2,1);
+for i = 1:N/2
+    C(i) = sqrt(Rz.psd_x(i)/T0);
+end
+
+% Translation Stage - X
+if args.Frz_x && args.enable
+    phi = Rz.psd_x;
     C = zeros(N/2,1);
     for i = 1:N/2
-        C(i) = sqrt(phi(i)*df);
+        C(i) = sqrt(phi(i)/T0);
+    end
+    theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
+    Cx = [0 ; C.*complex(cos(theta),sin(theta))];
+    Cx = [Cx; flipud(conj(Cx(2:end)))];;
+    u = N/sqrt(2)*ifft(Cx); % Disturbance Force Rz x [N]
+    Rz.x = u;
+else
+    Rz.x = zeros(length(Rz.t), 1);
+end
+
+% Translation Stage - Y
+if args.Frz_y && args.enable
+    phi = Rz.psd_y;
+    C = zeros(N/2,1);
+    for i = 1:N/2
+        C(i) = sqrt(phi(i)/T0);
+    end
+    theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
+    Cx = [0 ; C.*complex(cos(theta),sin(theta))];
+    Cx = [Cx; flipud(conj(Cx(2:end)))];;
+    u = N/sqrt(2)*ifft(Cx); % Disturbance Force Rz y [N]
+    Rz.y = u;
+else
+    Rz.y = zeros(length(Rz.t), 1);
+end
+
+% Translation Stage - Z
+if args.Frz_z && args.enable
+    phi = Rz.psd_z;
+    C = zeros(N/2,1);
+    for i = 1:N/2
+        C(i) = sqrt(phi(i)/T0);
     end
-    rng(131);
     theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
     Cx = [0 ; C.*complex(cos(theta),sin(theta))];
     Cx = [Cx; flipud(conj(Cx(2:end)))];;
     u = N/sqrt(2)*ifft(Cx); % Disturbance Force Rz z [N]
-    Frz_z = u;
+    Rz.z = u;
 else
-    Frz_z = zeros(length(t), 1);
+    Rz.z = zeros(length(Rz.t), 1);
 end
 #+end_src
 
 *** Direct Forces
 #+begin_src matlab
-u = zeros(length(t), 6);
+u = zeros(100, 6);
 Fd = u;
 #+end_src
 
 *** Set initial value to zero
 #+begin_src matlab
-Dwx    = Dwx   - Dwx(1);
-Dwy    = Dwy   - Dwy(1);
-Dwz    = Dwz   - Dwz(1);
-Fty_x  = Fty_x - Fty_x(1);
-Fty_z  = Fty_z - Fty_z(1);
-Frz_z  = Frz_z - Frz_z(1);
+Dw.x   = Dw.x   - Dw.x(1);
+Dw.y   = Dw.y   - Dw.y(1);
+Dw.z   = Dw.z   - Dw.z(1);
+Ty.x   = Ty.x   - Ty.x(1);
+Ty.z   = Ty.z   - Ty.z(1);
+Rz.x   = Rz.x   - Rz.x(1);
+Rz.y   = Rz.y   - Rz.y(1);
+Rz.z   = Rz.z   - Rz.z(1);
 #+end_src
 
 *** Save the Structure
 #+begin_src matlab
 if exist('./mat', 'dir')
     if exist('./mat/nass_disturbances.mat', 'file')
-        save('mat/nass_disturbances.mat', 'Dwx', 'Dwy', 'Dwz', 'Fty_x', 'Fty_z', 'Frz_x', 'Frz_y', 'Frz_z', 'Fd', 'Ts', 't', 'args', '-append');
+        save('mat/nass_disturbances.mat', 'Dw', 'Ty', 'Rz', 'Fd', 'args', '-append');
     else
-        save('mat/nass_disturbances.mat', 'Dwx', 'Dwy', 'Dwz', 'Fty_x', 'Fty_z', 'Frz_x', 'Frz_y', 'Frz_z', 'Fd', 'Ts', 't', 'args');
+        save('mat/nass_disturbances.mat', 'Dw', 'Ty', 'Rz', 'Fd', 'args');
     end
 elseif exist('./matlab', 'dir')
     if exist('./matlab/mat/nass_disturbances.mat', 'file')
-        save('matlab/mat/nass_disturbances.mat', 'Dwx', 'Dwy', 'Dwz', 'Fty_x', 'Fty_z', 'Frz_x', 'Frz_y', 'Frz_z', 'Fd', 'Ts', 't', 'args', '-append');
+        save('matlab/mat/nass_disturbances.mat', 'Dw', 'Ty', 'Rz', 'Fd', 'args', '-append');
     else
-        save('matlab/mat/nass_disturbances.mat', 'Dwx', 'Dwy', 'Dwz', 'Fty_x', 'Fty_z', 'Frz_x', 'Frz_y', 'Frz_z', 'Fd', 'Ts', 't', 'args');
+        save('matlab/mat/nass_disturbances.mat', 'Dw', 'Ty', 'Rz', 'Fd', 'args');
     end
 end
 #+end_src