diff --git a/index.org b/index.org
index 50380ea..7763fb4 100644
--- a/index.org
+++ b/index.org
@@ -1856,7 +1856,17 @@ The file =mat/plant.mat= is accessible [[./mat/plant.mat][here]].
   figure; bode(sys_cl({'Vph', 'Vpv'}, {'Uch', 'Ucv'}), sys({'Vph', 'Vpv'}, {'Uch', 'Ucv'}))
 #+end_src
 
-* TODO Huddle Test
+* Huddle Test
+** Matlab Init                                              :noexport:ignore:
+#+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
+  <<matlab-dir>>
+#+end_src
+
+#+begin_src matlab :exports none :results silent :noweb yes
+  <<matlab-init>>
+#+end_src
+
+** Load data
 We load the data taken during the Huddle Test.
 #+begin_src matlab
   load('mat/data_huddle_test.mat', ...
@@ -1868,6 +1878,7 @@ We load the data taken during the Huddle Test.
        'Va');
 #+end_src
 
+** Pre-processing
 We remove the first second of data where everything is settling down.
 #+begin_src matlab
   t0 = 1;
@@ -1888,11 +1899,12 @@ We remove the first second of data where everything is settling down.
   t = t - t(1); % We start at t=0
 #+end_src
 
+** Time Domain Data
 #+begin_src matlab :exports none
   figure;
   hold on;
-  plot(t, Vph, 'DisplayName', '$Vp_h$');
-  plot(t, Vpv, 'DisplayName', '$Vp_v$');
+  plot(t, Uch, 'DisplayName', '$Vp_h$');
+  plot(t, Ucv, 'DisplayName', '$Vp_v$');
   hold off;
   xlabel('Time [s]');
   ylabel('Amplitude [V]');
@@ -1906,7 +1918,7 @@ We compute the Power Spectral Density of the horizontal and vertical positions o
   [psd_Vpv, ~] = pwelch(Vpv, hanning(ceil(1*fs)), [], [], fs);
 #+end_src
 
-#+begin_src matlab :results none
+#+begin_src matlab :exports none
   figure;
   hold on;
   plot(f, sqrt(psd_Vph), 'DisplayName', '$\Gamma_{Vp_h}$');
@@ -1936,7 +1948,7 @@ We compute the Power Spectral Density of the voltage across the inductance used
   [psd_Vcv, ~] = pwelch(Vcv, hanning(ceil(1*fs)), [], [], fs);
 #+end_src
 
-#+begin_src matlab :results none
+#+begin_src matlab :exports none
   figure;
   hold on;
   plot(f, sqrt(psd_Vch), 'DisplayName', '$\Gamma_{Vc_h}$');
@@ -2237,7 +2249,12 @@ The controllers can be downloaded [[./mat/K_newport.mat][here]].
   save('mat/K_newport.mat', 'Kn', 'Knd');
 #+end_src
 
-* Measurement of the non-repeatability
+* Measuement of the non-repeatability
+** Introduction                                                     :ignore:
+- Explanation of the procedure
+- List all sources of error and their effects on the Attocube measurement
+- Think about how to determine the value of the individual sources of error
+
 ** Matlab Init                                              :noexport:ignore:
 #+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
   <<matlab-dir>>
@@ -2247,37 +2264,305 @@ The controllers can be downloaded [[./mat/K_newport.mat][here]].
   <<matlab-init>>
 #+end_src
 
+#+begin_src matlab
+  fs = 1e4;
+#+end_src
+
 ** Data Load
 #+begin_src matlab
-  load('mat/data_rep_1.mat', ...
+  uh = load('mat/data_rep_h.mat', ...
+       't', 'Uch', 'Ucv', ...
+       'Unh', 'Unv', ...
+       'Vph', 'Vpv', ...
+       'Vnh', 'Vnv', ...
+       'Va');
+
+  uv = load('mat/data_rep_v.mat', ...
        't', 'Uch', 'Ucv', ...
        'Unh', 'Unv', ...
        'Vph', 'Vpv', ...
-       'Vch', 'Vcv', ...
        'Vnh', 'Vnv', ...
        'Va');
 #+end_src
 
 #+begin_src matlab
-  t0 = 5;
+  % Let's start one second after the first command in the system
+  i0 = find(uh.Unh ~= 0, 1) + fs;
+  iend = i0+fs*floor((length(uh.t)-i0)/fs);
 
-  Uch(t<t0) = [];
-  Ucv(t<t0) = [];
-  Unh(t<t0) = [];
-  Unv(t<t0) = [];
-  Vph(t<t0) = [];
-  Vpv(t<t0) = [];
-  Vch(t<t0) = [];
-  Vcv(t<t0) = [];
-  Vnh(t<t0) = [];
-  Vnv(t<t0) = [];
-  Va(t<t0)  = [];
-  t(t<t0)   = [];
+  uh.Uch([1:i0-1, iend:end]) = [];
+  uh.Ucv([1:i0-1, iend:end]) = [];
+  uh.Unh([1:i0-1, iend:end]) = [];
+  uh.Unv([1:i0-1, iend:end]) = [];
+  uh.Vph([1:i0-1, iend:end]) = [];
+  uh.Vpv([1:i0-1, iend:end]) = [];
+  uh.Vnh([1:i0-1, iend:end]) = [];
+  uh.Vnv([1:i0-1, iend:end]) = [];
+  uh.Va ([1:i0-1, iend:end]) = [];
+  uh.t  ([1:i0-1, iend:end]) = [];
 
-  t = t - t(1); % We start at t=0
+  % We reset the time t
+  uh.t = uh.t - uh.t(1);
 #+end_src
 
-** TODO Some Plots
+#+begin_src matlab
+  % Let's start one second after the first command in the system
+  i0 = find(uv.Unv ~= 0, 1) + fs;
+  iend = i0+fs*floor((length(uv.t)-i0)/fs);
+
+  uv.Uch([1:i0-1, iend:end]) = [];
+  uv.Ucv([1:i0-1, iend:end]) = [];
+  uv.Unh([1:i0-1, iend:end]) = [];
+  uv.Unv([1:i0-1, iend:end]) = [];
+  uv.Vph([1:i0-1, iend:end]) = [];
+  uv.Vpv([1:i0-1, iend:end]) = [];
+  uv.Vnh([1:i0-1, iend:end]) = [];
+  uv.Vnv([1:i0-1, iend:end]) = [];
+  uv.Va ([1:i0-1, iend:end]) = [];
+  uv.t  ([1:i0-1, iend:end]) = [];
+
+  % We reset the time t
+  uv.t = uv.t - uv.t(1);
+#+end_src
+
+** Verify Tracking Angle Error
+Let's verify that the positioning error of the beam is small and what could be the effect on the distance measured by the intereferometer.
+
+#+begin_src matlab
+  load('./mat/plant.mat', 'Gd');
+#+end_src
+
+#+begin_src matlab
+  figure;
+  ax1 = subplot(1, 2, 1);
+  hold on;
+  plot(uh.t(1:2*fs), 1e6*uh.Vph(1:2*fs)/freqresp(Gd(1,1), 0), 'DisplayName', '$\theta_{h}$');
+  plot(uh.t(1:2*fs), 1e6*uh.Vpv(1:2*fs)/freqresp(Gd(2,2), 0), 'DisplayName', '$\theta_{v}$');
+  hold off;
+  xlabel('Time [s]'); ylabel('$\theta$ [$\mu$ rad]');
+  title('Newport Tilt - Horizontal Direction');
+  legend();
+
+  ax2 = subplot(1, 2, 2);
+  hold on;
+  plot(uv.t(1:2*fs), 1e6*uv.Vph(1:2*fs)/freqresp(Gd(1,1), 0), 'DisplayName', '$\theta_{h}$');
+  plot(uv.t(1:2*fs), 1e6*uv.Vpv(1:2*fs)/freqresp(Gd(2,2), 0), 'DisplayName', '$\theta_{v}$');
+  hold off;
+  xlabel('Time [s]'); ylabel('$\theta$ [$\mu$ rad]');
+  title('Newport Tilt - Vertical Direction');
+  legend();
+
+  linkaxes([ax1,ax2],'xy');
+#+end_src
+
+Let's compute the PSD of the error to see the frequency content.
+
+#+begin_src matlab
+  [psd_UhRh, f] = pwelch(uh.Vph/freqresp(Gd(1,1), 0), hanning(ceil(1*fs)), [], [], fs);
+  [psd_UhRv, ~] = pwelch(uh.Vpv/freqresp(Gd(2,2), 0), hanning(ceil(1*fs)), [], [], fs);
+  [psd_UvRh, ~] = pwelch(uv.Vph/freqresp(Gd(1,1), 0), hanning(ceil(1*fs)), [], [], fs);
+  [psd_UvRv, ~] = pwelch(uv.Vpv/freqresp(Gd(2,2), 0), hanning(ceil(1*fs)), [], [], fs);
+#+end_src
+
+#+begin_src matlab :exports none
+  figure;
+  ax1 = subplot(1, 2, 1);
+  hold on;
+  plot(f, sqrt(psd_UhRh), 'DisplayName', '$\Gamma_{\theta_h}$');
+  plot(f, sqrt(psd_UhRv), 'DisplayName', '$\Gamma_{\theta_v}$');
+  hold off;
+  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('ASD $\left[\frac{rad}{\sqrt{Hz}}\right]$')
+  legend('Location', 'southwest');
+  title('Newport Tilt - Horizontal Direction');
+
+  ax2 = subplot(1, 2, 2);
+  hold on;
+  plot(f, sqrt(psd_UvRh), 'DisplayName', '$\Gamma_{\theta_h}$');
+  plot(f, sqrt(psd_UvRv), 'DisplayName', '$\Gamma_{\theta_v}$');
+  hold off;
+  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('ASD $\left[\frac{rad}{\sqrt{Hz}}\right]$')
+  legend('Location', 'southwest');
+  title('Newport Tilt - Vertical Direction');
+
+  linkaxes([ax1,ax2],'xy');
+  xlim([1, 1000]);
+#+end_src
+
+Let's convert that to errors in distance
+
+\[ \Delta L = L^\prime - L = \frac{L}{\cos(\alpha)} - L \approx \frac{L \alpha^2}{2} \]
+with
+- $L$ is the nominal distance traveled by the beam
+- $L^\prime$ is the distance traveled by the beam with an angle error
+- $\alpha$ is the angle error
+
+#+begin_src matlab
+  L = 0.1; % [m]
+#+end_src
+
+#+begin_src matlab
+  [psd_UhLh, f] = pwelch(0.5*L*(uh.Vph/freqresp(Gd(1,1), 0)).^2, hanning(ceil(1*fs)), [], [], fs);
+  [psd_UhLv, ~] = pwelch(0.5*L*(uh.Vpv/freqresp(Gd(2,2), 0)).^2, hanning(ceil(1*fs)), [], [], fs);
+  [psd_UvLh, ~] = pwelch(0.5*L*(uv.Vph/freqresp(Gd(1,1), 0)).^2, hanning(ceil(1*fs)), [], [], fs);
+  [psd_UvLv, ~] = pwelch(0.5*L*(uv.Vpv/freqresp(Gd(2,2), 0)).^2, hanning(ceil(1*fs)), [], [], fs);
+#+end_src
+
+#+begin_src matlab :exports none
+  figure;
+  ax1 = subplot(1, 2, 1);
+  hold on;
+  plot(f, sqrt(psd_UhLh), 'DisplayName', '$\Gamma_{L_h}$');
+  plot(f, sqrt(psd_UhLv), 'DisplayName', '$\Gamma_{L_v}$');
+  hold off;
+  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('ASD $\left[\frac{m}{\sqrt{Hz}}\right]$')
+  legend('Location', 'southwest');
+  title('Newport Tilt - Horizontal Direction');
+
+  ax2 = subplot(1, 2, 2);
+  hold on;
+  plot(f, sqrt(psd_UvLh), 'DisplayName', '$\Gamma_{L_h}$');
+  plot(f, sqrt(psd_UvLv), 'DisplayName', '$\Gamma_{L_v}$');
+  hold off;
+  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('ASD $\left[\frac{m}{\sqrt{Hz}}\right]$')
+  legend('Location', 'southwest');
+  title('Newport Tilt - Vertical Direction');
+
+  linkaxes([ax1,ax2],'xy');
+  xlim([1, 1000]);
+#+end_src
+
+Now, compare that with the PSD of the measured distance by the interferometer.
+#+begin_src matlab
+  [psd_Lh, f] = pwelch(uh.Va, hanning(ceil(1*fs)), [], [], fs);
+  [psd_Lv, ~] = pwelch(uv.Va, hanning(ceil(1*fs)), [], [], fs);
+#+end_src
+
+#+begin_src matlab :exports none
+  figure;
+  ax1 = subplot(1, 2, 1);
+  hold on;
+  plot(f, sqrt(psd_UhLh), 'DisplayName', '$\Gamma_{L_h}$');
+  plot(f, sqrt(psd_UhLv), 'DisplayName', '$\Gamma_{L_v}$');
+  plot(f, sqrt(psd_Lh), '--k', 'DisplayName', '$\Gamma_{L_h}$');
+  hold off;
+  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('ASD $\left[\frac{m}{\sqrt{Hz}}\right]$')
+  legend('Location', 'southwest');
+  title('Newport Tilt - Horizontal Direction');
+
+  ax2 = subplot(1, 2, 2);
+  hold on;
+  plot(f, sqrt(psd_UvLh), 'DisplayName', '$\Gamma_{L_h}$');
+  plot(f, sqrt(psd_UvLv), 'DisplayName', '$\Gamma_{L_v}$');
+  plot(f, sqrt(psd_Lv), '--k', 'DisplayName', '$\Gamma_{L_h}$');
+  hold off;
+  set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('ASD $\left[\frac{m}{\sqrt{Hz}}\right]$')
+  legend('Location', 'southwest');
+  title('Newport Tilt - Vertical Direction');
+
+  linkaxes([ax1,ax2],'xy');
+  xlim([1, 1000]);
+#+end_src
+
+** Processing
+First, we get the mean value as measured by the interferometer for each value of the Newport angle.
+#+begin_src matlab
+  Vahm = mean(reshape(uh.Va, [fs floor(length(uh.t)/fs)]),2);
+  Unhm = mean(reshape(uh.Unh, [fs floor(length(uh.t)/fs)]),2);
+
+  Vavm = mean(reshape(uv.Va, [fs floor(length(uv.t)/fs)]),2);
+  Unvm = mean(reshape(uv.Unv, [fs floor(length(uv.t)/fs)]),2);
+#+end_src
+
+#+begin_src matlab :exports none
+  figure;
+  ax1 = subplot(1, 2, 1);
+  hold on;
+  plot(uh.Unh, uh.Va);
+  plot(Unhm, Vahm)
+  hold off;
+  xlabel('$V_{n,h}$ [V]'); ylabel('$V_a$ [m]');
+
+  ax2 = subplot(1, 2, 2);
+  hold on;
+  plot(uv.Unv, uv.Va);
+  plot(Unvm, Vavm)
+  hold off;
+  xlabel('$V_{n,v}$ [V]'); ylabel('$V_a$ [m]');
+#+end_src
+
+#+begin_src matlab :exports none
+  figure;
+  ax1 = subplot(1, 2, 1);
+  hold on;
+  plot(uh.Unh, 1e9*(uh.Va - repmat(Vahm, length(uh.t)/length(Vahm),1)));
+  hold off;
+  xlabel('$V_{n,h}$ [V]'); ylabel('$V_a$ [nm]');
+
+  ax2 = subplot(1, 2, 2);
+  hold on;
+  plot(uv.Unv, 1e9*(uv.Va - repmat(Vavm, length(uv.t)/length(Vavm),1)));
+  hold off;
+  xlabel('$V_{n,v}$ [V]'); ylabel('$V_a$ [nm]');
+
+  linkaxes([ax1,ax2],'xy');
+  ylim([-100 100]);
+#+end_src
+
+We then subtract
+#+begin_src matlab
+  figure;
+  hold on;
+  plot(uh.Unh, 1e9*(uh.Va - repmat(Vam, length(uh.t)/length(Vam),1)))
+  hold off;
+  xlabel('$V_{n,h}$ [V]'); ylabel('$V_a$ [nm]');
+#+end_src
+
+** Some Plots
+#+begin_src matlab
+  figure;
+  hold on;
+  plot(uh.Unh, uh.Va);
+  plot(Unhm, Vam)
+  hold off;
+  xlabel('$V_{n,h}$ [V]'); ylabel('$V_a$ [m]');
+#+end_src
+
+#+begin_src matlab :exports none
+  figure;
+  ax1 = subplot(1, 2, 1);
+  hold on;
+  plot(uh.Vnh(1:fs/2), uh.Va(1:fs/2));
+  hold off;
+  xlabel('$V_{n,h}$ [V]'); ylabel('$V_a$ [m]');
+
+  ax2 = subplot(1, 2, 2);
+  hold on;
+  plot(uv.Vnv, uv.Va);
+  hold off;
+  xlabel('$V_{n,v}$ [V]'); ylabel('$V_a$ [m]');
+#+end_src
+
+
+#+begin_src matlab :exports none
+  figure;
+  ax1 = subplot(1, 2, 1);
+  hold on;
+  plot(uh.Vnh, uh.Va);
+  hold off;
+  xlabel('$V_{n,h}$ [V]'); ylabel('$V_a$ [m]');
+
+  ax2 = subplot(1, 2, 2);
+  hold on;
+  plot(uv.Vnv, uv.Va);
+  hold off;
+  xlabel('$V_{n,v}$ [V]'); ylabel('$V_a$ [m]');
+#+end_src
 
 ** Repeatability
 #+begin_src matlab :exports none
diff --git a/matlab/test_cercalo.slx b/matlab/test_cercalo.slx
index 9745183..d8c36b3 100644
Binary files a/matlab/test_cercalo.slx and b/matlab/test_cercalo.slx differ