Add analysis of the damped plant with IFF

active_damping/index.org

@@ -66,8 +66,12 @@ The disturbances are:
 <<sec:undamped_system>>
 
 ** Introduction                                                      :ignore:
-We first look at the undamped system.
-The performance of this undamped system will be compared with the damped system using various techniques.
+In this section, we identify the dynamic of the system from forces applied in the nano-hexapod legs to the various sensors included in the nano-hexapod that could be use for Active Damping, namely:
+- A relative motion sensor, measuring the relative displacement of each of the leg
+- A force sensor measuring the total force transmitted to the top part of the leg in the direction of the leg
+- A absolute velocity sensor measuring the absolute velocity of the top part of the leg in the direction of the leg
+
+After that, a tomography experiment is simulation without any active damping techniques.
 
 ** Matlab Init                                              :noexport:ignore:
 #+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
@@ -115,6 +119,7 @@ We set the references to zero.
   initializeReferences();
 #+end_src
 
+No disturbance is included in the system.
 #+begin_src matlab
   initializeDisturbances('enable', false);
 #+end_src
@@ -143,17 +148,19 @@ First, we identify the dynamics of the system using the =linearize= function.
 
   %% Input/Output definition
   clear io; io_i = 1;
-  io(io_i) = linio([mdl, '/Fnl'],           1, 'openinput');              io_i = io_i + 1;
-  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Dnlm'); io_i = io_i + 1;
-  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Fnlm'); io_i = io_i + 1;
-  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Vlm');  io_i = io_i + 1;
+  io(io_i) = linio([mdl, '/Fnl'],           1, 'openinput');               io_i = io_i + 1; % Actuator Inputs
+  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Dnlm');  io_i = io_i + 1; % Relative Motion Outputs
+  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Fnlm');  io_i = io_i + 1; % Force Sensors
+  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Vlm');   io_i = io_i + 1; % Absolute Velocity Outputs
+  io(io_i) = linio([mdl, '/Compute Error in NASS base'], 2, 'openoutput'); io_i = io_i + 1; % Metrology Outputs
 
   %% Run the linearization
   G = linearize(mdl, io, 0.5, options);
   G.InputName  = {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'};
   G.OutputName = {'Dnlm1', 'Dnlm2', 'Dnlm3', 'Dnlm4', 'Dnlm5', 'Dnlm6', ...
                   'Fnlm1', 'Fnlm2', 'Fnlm3', 'Fnlm4', 'Fnlm5', 'Fnlm6', ...
-                  'Vnlm1', 'Vnlm2', 'Vnlm3', 'Vnlm4', 'Vnlm5', 'Vnlm6'};
+                  'Vnlm1', 'Vnlm2', 'Vnlm3', 'Vnlm4', 'Vnlm5', 'Vnlm6', ...
+                  'Dxn', 'Dyn', 'Dzn', 'Rxn', 'Ryn', 'Rzn'};
 #+end_src
 
 We then create transfer functions corresponding to the active damping plants.
@@ -163,9 +170,16 @@ We then create transfer functions corresponding to the active damping plants.
   G_ine = minreal(G({'Vnlm1', 'Vnlm2', 'Vnlm3', 'Vnlm4', 'Vnlm5', 'Vnlm6'}, {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'}));
 #+end_src
 
+#+begin_src matlab
+  load('mat/stages.mat', 'nano_hexapod');
+  G_cart = minreal(G({'Dxn', 'Dyn', 'Dzn', 'Rxn', 'Ryn', 'Rzn'}, {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'}))*inv(nano_hexapod.J');
+  G_cart.InputName = {'Fnx', 'Fny', 'Fnz', 'Mnx', 'Mny', 'Mnz'};
+#+end_src
+
 And we save them for further analysis.
 #+begin_src matlab
   save('./active_damping/mat/undamped_plants.mat', 'G_iff', 'G_dvf', 'G_ine');
+  save('./active_damping/mat/cart_plants.mat', 'G_cart');
 #+end_src
 
 *** Obtained Plants for Active Damping
@@ -207,7 +221,7 @@ And we save them for further analysis.
 #+end_src
 
 #+NAME: fig:nass_active_damping_iff_plant
-#+CAPTION: =G_iff=: IFF Plant ([[./figs/nass_active_damping_iff_plant.png][png]], [[./figs/nass_active_damping_iff_plant.pdf][pdf]])
+#+CAPTION: =G_iff=: Transfer functions from forces applied in the actuators to the force sensor in each actuator ([[./figs/nass_active_damping_iff_plant.png][png]], [[./figs/nass_active_damping_iff_plant.pdf][pdf]])
 [[file:figs/nass_active_damping_iff_plant.png]]
 
 #+begin_src matlab :exports none
@@ -244,8 +258,8 @@ And we save them for further analysis.
 #+end_src
 
 #+NAME: fig:nass_active_damping_dvf_plant
-#+CAPTION: =G_dvf=: Plant for Direct Velocity Feedback ([[./figs/nass_active_damping_dvf_plant.png][png]], [[./figs/nass_active_damping_dvf_plant.pdf][pdf]])
-[[file:figs/nass_active_damping_ine_plant.png]]
+#+CAPTION: =G_dvf=: Transfer functions from forces applied in the actuators to the relative motion sensor in each actuator ([[./figs/nass_active_damping_dvf_plant.png][png]], [[./figs/nass_active_damping_dvf_plant.pdf][pdf]])
+[[file:figs/nass_active_damping_dvf_plant.png]]
 
 #+begin_src matlab :exports none
   freqs = logspace(0, 3, 1000);
@@ -281,7 +295,7 @@ And we save them for further analysis.
 #+end_src
 
 #+NAME: fig:nass_active_damping_inertial_plant
-#+CAPTION: =G_ine=: Inertial Feedback Plant ([[./figs/nass_active_damping_inertial_plant.png][png]], [[./figs/nass_active_damping_inertial_plant.pdf][pdf]])
+#+CAPTION: =G_ine=: Transfer functions from forces applied in the actuators to the geophone located in each leg measuring the absolute velocity of the top part of the leg in the direction of the leg ([[./figs/nass_active_damping_inertial_plant.png][png]], [[./figs/nass_active_damping_inertial_plant.pdf][pdf]])
 [[file:figs/nass_active_damping_inertial_plant.png]]
 
 ** Tomography Experiment
@@ -508,7 +522,8 @@ We identify the dynamics for the following sample mass.
   ax2 = subplot(2, 1, 2);
   hold on;
   for i = 1:length(Gm_iff)
-    plot(freqs, 180/pi*angle(squeeze(freqresp(Gm_iff{i}('Fnlm1', 'Fnl1'), freqs, 'Hz'))), 'DisplayName', sprintf('$M = %.0f$ [kg]', masses(i)));
+    plot(freqs, 180/pi*angle(squeeze(freqresp(Gm_iff{i}('Fnlm1', 'Fnl1'), freqs, 'Hz'))), ...
+         'DisplayName', sprintf('$M = %.0f$ [kg]', masses(i)));
   end
   hold off;
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
@@ -1688,7 +1703,16 @@ Also, for the Inertial Sensor, a RHP zero may appear when the spindle is rotatin
 #+end_note
 
 ** Introduction                                                      :ignore:
-Integral Force Feedback is applied on the simscape model.
+Here, we study the use of *Integral Force Feedback* (IFF) to actively damp the resonances.
+
+The IFF control is applied in a decentralized way: there is on controller for each leg.
+
+The control architecture is represented in figure [[fig:iff_1dof]] where one of the 6 nano-hexapod legs is represented.
+
+#+name: fig:iff_1dof
+#+caption: Integral Force Feedback applied to a 1dof system
+#+RESULTS:
+[[file:figs/iff_1dof.png]]
 
 ** Matlab Init                                              :noexport:ignore:
 #+begin_src matlab :tangle no :exports none :results silent :noweb yes :var current_dir=(file-name-directory buffer-file-name)
@@ -1713,9 +1737,10 @@ Integral Force Feedback is applied on the simscape model.
 
 ** Control Design
 *** Plant
-Let's load the previously indentified undamped plant:
+Let's load the previously identified undamped plant:
 #+begin_src matlab
   load('./active_damping/mat/undamped_plants.mat', 'G_iff');
+  load('./active_damping/mat/plants_variable.mat', 'masses', 'Gm_iff');
 #+end_src
 
 Let's look at the transfer function from actuator forces in the nano-hexapod to the force sensor in the nano-hexapod legs for all 6 pairs of actuator/sensor (figure [[fig:iff_plant]]).
@@ -1727,8 +1752,8 @@ Let's look at the transfer function from actuator forces in the nano-hexapod to
 
   ax1 = subplot(2, 1, 1);
   hold on;
-  for i=1:6
-    plot(freqs, abs(squeeze(freqresp(G_iff(['Fnlm', num2str(i)], ['Fnl', num2str(i)]), freqs, 'Hz'))));
+  for i=1:length(masses)
+    plot(freqs, abs(squeeze(freqresp(-Gm_iff{i}('Fnlm1', 'Fnl1'), freqs, 'Hz'))));
   end
   hold off;
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
@@ -1736,14 +1761,16 @@ Let's look at the transfer function from actuator forces in the nano-hexapod to
 
   ax2 = subplot(2, 1, 2);
   hold on;
-  for i=1:6
-    plot(freqs, 180/pi*angle(squeeze(freqresp(G_iff(['Fnlm', num2str(i)], ['Fnl', num2str(i)]), freqs, 'Hz'))));
+  for i=1:length(masses)
+    plot(freqs, 180/pi*angle(squeeze(freqresp(-Gm_iff{i}('Fnlm1', 'Fnl1'), freqs, 'Hz'))), ...
+         'DisplayName', sprintf('$M = %.0f$ [kg]', masses(i)));
   end
   hold off;
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
   ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
   ylim([-180, 180]);
   yticks([-180, -90, 0, 90, 180]);
+  legend('location', 'southwest');
 
   linkaxes([ax1,ax2],'x');
 #+end_src
@@ -1760,7 +1787,8 @@ Let's look at the transfer function from actuator forces in the nano-hexapod to
 *** Control Design
 The controller for each pair of actuator/sensor is:
 #+begin_src matlab
-  K_iff = 1000/s;
+  w0 = 2*pi*50;
+  K_iff = -5000/s * (s/w0)/(1 + s/w0);
 #+end_src
 
 The corresponding loop gains are shown in figure [[fig:iff_open_loop]].
@@ -1772,8 +1800,8 @@ The corresponding loop gains are shown in figure [[fig:iff_open_loop]].
 
   ax1 = subplot(2, 1, 1);
   hold on;
-  for i=1:6
-    plot(freqs, abs(squeeze(freqresp(K_iff*G_iff(['Fnlm', num2str(i)], ['Fnl', num2str(i)]), freqs, 'Hz'))));
+  for i=1:length(masses)
+    plot(freqs, abs(squeeze(freqresp(K_iff*Gm_iff{i}('Fnlm1', 'Fnl1'), freqs, 'Hz'))));
   end
   hold off;
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
@@ -1781,14 +1809,16 @@ The corresponding loop gains are shown in figure [[fig:iff_open_loop]].
 
   ax2 = subplot(2, 1, 2);
   hold on;
-  for i=1:6
-    plot(freqs, 180/pi*angle(squeeze(freqresp(K_iff*G_iff(['Fnlm', num2str(i)], ['Fnl', num2str(i)]), freqs, 'Hz'))));
+  for i=1:length(masses)
+    plot(freqs, 180/pi*angle(squeeze(freqresp(K_iff*Gm_iff{i}('Fnlm1', 'Fnl1'), freqs, 'Hz'))), ...
+         'DisplayName', sprintf('$M = %.0f$ [kg]', masses(i)));
   end
   hold off;
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
   ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
   ylim([-180, 180]);
   yticks([-180, -90, 0, 90, 180]);
+  legend('location', 'southwest');
 
   linkaxes([ax1,ax2],'x');
 #+end_src
@@ -1814,94 +1844,38 @@ We save the controller for further analysis.
   save('./active_damping/mat/K_iff.mat', 'K_iff');
 #+end_src
 
-*** IFF with High Pass Filter
-#+begin_src matlab
-  w_hpf = 2*pi*10; % Cut-off frequency for the high pass filter [rad/s]
-
-  K_iff = 2*pi*200/s * (s/w_hpf)/(s/w_hpf + 1);
-#+end_src
-
-The corresponding loop gains are shown in figure [[fig:iff_hpf_open_loop]].
-#+begin_src matlab :exports none
-  freqs = logspace(0, 3, 1000);
-
-  figure;
-
-  ax1 = subplot(2, 1, 1);
-  hold on;
-  for i=1:6
-    plot(freqs, abs(squeeze(freqresp(K_iff*G_iff(['Fnlm', num2str(i)], ['Fnl', num2str(i)]), freqs, 'Hz'))));
-  end
-  hold off;
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude [N/N]'); set(gca, 'XTickLabel',[]);
-
-  ax2 = subplot(2, 1, 2);
-  hold on;
-  for i=1:6
-    plot(freqs, 180/pi*angle(squeeze(freqresp(K_iff*G_iff(['Fnlm', num2str(i)], ['Fnl', num2str(i)]), freqs, 'Hz'))));
-  end
-  hold off;
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
-  ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
-  ylim([-180, 180]);
-  yticks([-180, -90, 0, 90, 180]);
-
-  linkaxes([ax1,ax2],'x');
-#+end_src
-
-#+HEADER: :tangle no :exports results :results none :noweb yes
-#+begin_src matlab :var filepath="figs/iff_hpf_open_loop.pdf" :var figsize="full-tall" :post pdf2svg(file=*this*, ext="png")
-  <<plt-matlab>>
-#+end_src
-
-#+NAME: fig:iff_hpf_open_loop
-#+CAPTION: Loop Gain for the Integral Force Feedback with an High pass filter ([[./figs/iff_hpf_open_loop.png][png]], [[./figs/iff_hpf_open_loop.pdf][pdf]])
-[[file:figs/iff_hpf_open_loop.png]]
-
-We create the diagonal controller and we add a minus sign as we have a positive
-feedback architecture.
-#+begin_src matlab
-  K_iff = -K_iff*eye(6);
-#+end_src
-
-We save the controller for further analysis.
-#+begin_src matlab
-  save('./active_damping/mat/K_iff_hpf.mat', 'K_iff');
-#+end_src
-
 ** TODO Identification of the damped plant                         :noexport:
-*** Initialize the Simulation
+*** Initialize the Simulation                                     :noexport:
 We initialize all the stages with the default parameters.
 #+begin_src matlab
-    initializeGround();
-    initializeGranite();
-    initializeTy();
-    initializeRy();
-    initializeRz();
-    initializeMicroHexapod();
-    initializeAxisc();
-    initializeMirror();
+  initializeGround();
+  initializeGranite();
+  initializeTy();
+  initializeRy();
+  initializeRz();
+  initializeMicroHexapod();
+  initializeAxisc();
+  initializeMirror();
 #+end_src
 
-The nano-hexapod is a piezoelectric hexapod and the sample has a mass of 50kg.
+No disturbances.
+#+begin_src matlab
+  initializeDisturbances('enable', false);
+#+end_src
+
+The nano-hexapod is a piezoelectric hexapod.
 #+begin_src matlab
   initializeNanoHexapod('actuator', 'piezo');
   initializeSample('mass', 50);
 #+end_src
 
-We set the references to zero.
-#+begin_src matlab
-  initializeReferences();
-#+end_src
-
-And all the controllers are set to 0 except for the IFF.
+And all the controllers are set to 0.
 #+begin_src matlab
   K = tf(zeros(6));
   save('./mat/controllers.mat', 'K', '-append');
   K_ine = tf(zeros(6));
   save('./mat/controllers.mat', 'K_ine', '-append');
-  K_iff = K_iff;
+  load('./active_damping/mat/K_iff.mat', 'K_iff');
   save('./mat/controllers.mat', 'K_iff', '-append');
   K_dvf = tf(zeros(6));
   save('./mat/controllers.mat', 'K_dvf', '-append');
@@ -1919,26 +1893,24 @@ First, we identify the dynamics of the system using the =linearize= function.
 
   %% Input/Output definition
   clear io; io_i = 1;
-  io(io_i) = linio([mdl, '/Fnl'],           1, 'openinput');              io_i = io_i + 1;
-  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Dnlm'); io_i = io_i + 1;
-  io(io_i) = linio([mdl, '/Micro-Station'], 3, 'openoutput', [], 'Fnlm'); io_i = io_i + 1;
+  io(io_i) = linio([mdl, '/Fnl'],           1, 'openinput');               io_i = io_i + 1;
+  io(io_i) = linio([mdl, '/Compute Error in NASS base'], 2, 'openoutput'); io_i = io_i + 1;
 
   %% Run the linearization
-  G = linearize(mdl, io, options);
+  G = linearize(mdl, io, 0.5, options);
   G.InputName  = {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'};
-  G.OutputName = {'Dnlm1', 'Dnlm2', 'Dnlm3', 'Dnlm4', 'Dnlm5', 'Dnlm6', ...
-                  'Fnlm1', 'Fnlm2', 'Fnlm3', 'Fnlm4', 'Fnlm5', 'Fnlm6'};
+  G.OutputName = {'Dxn', 'Dyn', 'Dzn', 'Rxn', 'Ryn', 'Rzn'};
 #+end_src
 
-We then create transfer functions corresponding to the active damping plants.
 #+begin_src matlab
-  G_iff = minreal(G({'Fnlm1', 'Fnlm2', 'Fnlm3', 'Fnlm4', 'Fnlm5', 'Fnlm6'}, {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'}));
-  % G_rmc = minreal(G({'Dnlm1', 'Dnlm2', 'Dnlm3', 'Dnlm4', 'Dnlm5', 'Dnlm6'}, {'Fnl1', 'Fnl2', 'Fnl3', 'Fnl4', 'Fnl5', 'Fnl6'}));
+  load('mat/stages.mat', 'nano_hexapod');
+  G_cart_iff = G*inv(nano_hexapod.J');
+  G_cart_iff.InputName = {'Fnx', 'Fny', 'Fnz', 'Mnx', 'Mny', 'Mnz'};
 #+end_src
 
 And we save them for further analysis.
 #+begin_src matlab
-  save('./active_damping/mat/plants.mat', 'G_iff', '-append');
+  save('./active_damping/mat/cart_plants.mat', 'G_cart_iff', '-append');
 #+end_src
 
 *** TODO Sensitivity to disturbances
@@ -2022,8 +1994,11 @@ As shown on figure [[fig:sensitivity_dist_iff]]:
 #+CAPTION: Sensitivity to force disturbances in various stages when IFF is applied ([[./figs/sensitivity_dist_stages_iff.png][png]], [[./figs/sensitivity_dist_stages_iff.pdf][pdf]])
 [[file:figs/sensitivity_dist_stages_iff.png]]
 
-*** TODO Damped Plant
+*** Damped Plant
 Now, look at the new damped plant to control.
+#+begin_src matlab
+  load('./active_damping/mat/cart_plants.mat', 'G_cart', 'G_cart_iff');
+#+end_src
 
 It damps the plant (resonance of the nano hexapod as well as other resonances) as shown in figure [[fig:plant_iff_damped]].
 
@@ -2034,59 +2009,57 @@ It damps the plant (resonance of the nano hexapod as well as other resonances) a
 
   ax1 = subplot(2, 2, 1);
   hold on;
-  plot(freqs, abs(squeeze(freqresp(G.G_cart('Dx', 'Fnx'), freqs, 'Hz'))));
-  plot(freqs, abs(squeeze(freqresp(G.G_cart('Dy', 'Fny'), freqs, 'Hz'))));
-  plot(freqs, abs(squeeze(freqresp(G.G_cart('Dz', 'Fnz'), freqs, 'Hz'))));
+  plot(freqs, abs(squeeze(freqresp(G_cart('Dxn', 'Fnx'), freqs, 'Hz'))));
+  plot(freqs, abs(squeeze(freqresp(G_cart('Dyn', 'Fny'), freqs, 'Hz'))));
+  plot(freqs, abs(squeeze(freqresp(G_cart('Dzn', 'Fnz'), freqs, 'Hz'))));
   set(gca,'ColorOrderIndex',1);
-  plot(freqs, abs(squeeze(freqresp(G_iff.G_cart('Dx', 'Fnx'), freqs, 'Hz'))), '--');
-  plot(freqs, abs(squeeze(freqresp(G_iff.G_cart('Dy', 'Fny'), freqs, 'Hz'))), '--');
-  plot(freqs, abs(squeeze(freqresp(G_iff.G_cart('Dz', 'Fnz'), freqs, 'Hz'))), '--');
+  plot(freqs, abs(squeeze(freqresp(G_cart_iff('Dxn', 'Fnx'), freqs, 'Hz'))), '--');
+  plot(freqs, abs(squeeze(freqresp(G_cart_iff('Dyn', 'Fny'), freqs, 'Hz'))), '--');
+  plot(freqs, abs(squeeze(freqresp(G_cart_iff('Dzn', 'Fnz'), freqs, 'Hz'))), '--');
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
   ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
 
   ax2 = subplot(2, 2, 2);
   hold on;
-  plot(freqs, abs(squeeze(freqresp(G.G_cart('Rx', 'Mnx'), freqs, 'Hz'))));
-  plot(freqs, abs(squeeze(freqresp(G.G_cart('Ry', 'Mny'), freqs, 'Hz'))));
-  plot(freqs, abs(squeeze(freqresp(G.G_cart('Rz', 'Mnz'), freqs, 'Hz'))));
+  plot(freqs, abs(squeeze(freqresp(G_cart('Rxn', 'Mnx'), freqs, 'Hz'))));
+  plot(freqs, abs(squeeze(freqresp(G_cart('Ryn', 'Mny'), freqs, 'Hz'))));
+  plot(freqs, abs(squeeze(freqresp(G_cart('Rzn', 'Mnz'), freqs, 'Hz'))));
   set(gca,'ColorOrderIndex',1);
-  plot(freqs, abs(squeeze(freqresp(G_iff.G_cart('Rx', 'Mnx'), freqs, 'Hz'))), '--');
-  plot(freqs, abs(squeeze(freqresp(G_iff.G_cart('Ry', 'Mny'), freqs, 'Hz'))), '--');
-  plot(freqs, abs(squeeze(freqresp(G_iff.G_cart('Rz', 'Mnz'), freqs, 'Hz'))), '--');
+  plot(freqs, abs(squeeze(freqresp(G_cart_iff('Rxn', 'Mnx'), freqs, 'Hz'))), '--');
+  plot(freqs, abs(squeeze(freqresp(G_cart_iff('Ryn', 'Mny'), freqs, 'Hz'))), '--');
+  plot(freqs, abs(squeeze(freqresp(G_cart_iff('Rzn', 'Mnz'), freqs, 'Hz'))), '--');
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
   ylabel('Amplitude [rad/(Nm)]'); xlabel('Frequency [Hz]');
 
   ax3 = subplot(2, 2, 3);
   hold on;
-  plot(freqs, 180/pi*angle(squeeze(freqresp(G.G_cart('Dx', 'Fnx'), freqs, 'Hz'))), 'DisplayName', '$\left|D_x / F_{n,x}\right|$');
-  plot(freqs, 180/pi*angle(squeeze(freqresp(G.G_cart('Dy', 'Fny'), freqs, 'Hz'))), 'DisplayName', '$\left|D_y / F_{n,y}\right|$');
-  plot(freqs, 180/pi*angle(squeeze(freqresp(G.G_cart('Dz', 'Fnz'), freqs, 'Hz'))), 'DisplayName', '$\left|D_z / F_{n,z}\right|$');
+  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_cart('Dxn', 'Fnx'), freqs, 'Hz')))), 'DisplayName', '$\left|D_x / F_{n,x}\right|$');
+  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_cart('Dyn', 'Fny'), freqs, 'Hz')))), 'DisplayName', '$\left|D_y / F_{n,y}\right|$');
+  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_cart('Dzn', 'Fnz'), freqs, 'Hz')))), 'DisplayName', '$\left|D_z / F_{n,z}\right|$');
   set(gca,'ColorOrderIndex',1);
-  plot(freqs, 180/pi*angle(squeeze(freqresp(G_iff.G_cart('Dx', 'Fnx'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
-  plot(freqs, 180/pi*angle(squeeze(freqresp(G_iff.G_cart('Dy', 'Fny'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
-  plot(freqs, 180/pi*angle(squeeze(freqresp(G_iff.G_cart('Dz', 'Fnz'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
+  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_cart_iff('Dxn', 'Fnx'), freqs, 'Hz')))), '--', 'HandleVisibility', 'off');
+  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_cart_iff('Dyn', 'Fny'), freqs, 'Hz')))), '--', 'HandleVisibility', 'off');
+  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_cart_iff('Dzn', 'Fnz'), freqs, 'Hz')))), '--', 'HandleVisibility', 'off');
   hold off;
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
   ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
-  ylim([-180, 180]);
-  yticks([-180, -90, 0, 90, 180]);
-  legend('location', 'northwest');
+  yticks([-540:180:540]);
+  legend('location', 'southwest');
 
   ax4 = subplot(2, 2, 4);
   hold on;
-  plot(freqs, 180/pi*angle(squeeze(freqresp(G.G_cart('Rx', 'Mnx'), freqs, 'Hz'))), 'DisplayName', '$\left|R_x / M_{n,x}\right|$');
-  plot(freqs, 180/pi*angle(squeeze(freqresp(G.G_cart('Ry', 'Mny'), freqs, 'Hz'))), 'DisplayName', '$\left|R_y / M_{n,y}\right|$');
-  plot(freqs, 180/pi*angle(squeeze(freqresp(G.G_cart('Rz', 'Mnz'), freqs, 'Hz'))), 'DisplayName', '$\left|R_z / M_{n,z}\right|$');
+  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_cart('Rxn', 'Mnx'), freqs, 'Hz')))), 'DisplayName', '$\left|R_x / M_{n,x}\right|$');
+  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_cart('Ryn', 'Mny'), freqs, 'Hz')))), 'DisplayName', '$\left|R_y / M_{n,y}\right|$');
+  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_cart('Rzn', 'Mnz'), freqs, 'Hz')))), 'DisplayName', '$\left|R_z / M_{n,z}\right|$');
   set(gca,'ColorOrderIndex',1);
-  plot(freqs, 180/pi*angle(squeeze(freqresp(G_iff.G_cart('Rx', 'Mnx'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
-  plot(freqs, 180/pi*angle(squeeze(freqresp(G_iff.G_cart('Ry', 'Mny'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
-  plot(freqs, 180/pi*angle(squeeze(freqresp(G_iff.G_cart('Rz', 'Mnz'), freqs, 'Hz'))), '--', 'HandleVisibility', 'off');
+  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_cart_iff('Rxn', 'Mnx'), freqs, 'Hz')))), '--', 'HandleVisibility', 'off');
+  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_cart_iff('Ryn', 'Mny'), freqs, 'Hz')))), '--', 'HandleVisibility', 'off');
+  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G_cart_iff('Rzn', 'Mnz'), freqs, 'Hz')))), '--', 'HandleVisibility', 'off');
   hold off;
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
   ylabel('Phase [deg]'); xlabel('Frequency [Hz]');
-  ylim([-180, 180]);
-  yticks([-180, -90, 0, 90, 180]);
-  legend('location', 'northwest');
+  yticks([-540:180:540]);
+  legend('location', 'southwest');
 
   linkaxes([ax1,ax2,ax3,ax4],'x');
 #+end_src
@@ -2110,10 +2083,12 @@ However, it increases coupling at low frequency (figure [[fig:plant_iff_coupling
     for iy = 1:6
       subplot(6, 6, (ix-1)*6 + iy);
       hold on;
-      plot(freqs, abs(squeeze(freqresp(G.G_cart(ix, iy), freqs, 'Hz'))), 'k-');
-      plot(freqs, abs(squeeze(freqresp(G_iff.G_cart(ix, iy), freqs, 'Hz'))), 'k--');
+      plot(freqs, abs(squeeze(freqresp(G_cart(ix, iy), freqs, 'Hz'))), 'k-');
+      plot(freqs, abs(squeeze(freqresp(G_cart_iff(ix, iy), freqs, 'Hz'))), 'k--');
       set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-      ylim([1e-12, 1e-5]);
+      ylim([1e-13, 1e-4]);
+      xticks([])
+      yticks([])
     end
   end
 #+end_src
@@ -3673,8 +3648,8 @@ This Matlab function is accessible [[file:src/prepareTomographyExperiment.m][her
 #+begin_src matlab
   arguments
       args.nass_actuator       char   {mustBeMember(args.nass_actuator,{'piezo', 'lorentz'})} = 'piezo'
-      args.sample_mass   (1,1) double {mustBeNumeric, mustBePositive} = 50
-      args.Ry_period     (1,1) double {mustBeNumeric, mustBePositive} = 1
+      args.sample_mass   (1,1) double {mustBeNumeric, mustBePositive} = 50 % [kg]
+      args.Rz_period     (1,1) double {mustBeNumeric, mustBePositive} = 1 % [s]
   end
 #+end_src
 
@@ -3702,7 +3677,7 @@ The nano-hexapod is a piezoelectric hexapod and the sample has a mass of 50kg.
 
 We set the references to zero.
 #+begin_src matlab
-  initializeReferences('Rz_type', 'rotating', 'Rz_period', args.Ry_period);
+  initializeReferences('Rz_type', 'rotating', 'Rz_period', args.Rz_period);
 #+end_src
 
 And all the controllers are set to 0.
@@ -3762,7 +3737,7 @@ And all the controllers are set to 0.
 #+end_src
 
 #+begin_src matlab
-  lzoad('mat/stages.mat', 'nano_hexapod');
+  load('mat/stages.mat', 'nano_hexapod');
   G_cart = G*inv(nano_hexapod.J');
   G_cart.InputName = {'Fnx', 'Fny', 'Fnz', 'Mnx', 'Mny', 'Mnz'};
 #+end_src