diff --git a/data/APA300ML_new/APA300ML_test_bench.slx b/data/APA300ML_new/APA300ML_test_bench.slx
deleted file mode 100644
index 9a74b75..0000000
Binary files a/data/APA300ML_new/APA300ML_test_bench.slx and /dev/null differ
diff --git a/data/strut_encoder/strut_encoder.slx.r2020a b/data/strut_encoder/strut_encoder.slx.r2020a
deleted file mode 100644
index 9afd9a6..0000000
Binary files a/data/strut_encoder/strut_encoder.slx.r2020a and /dev/null differ
diff --git a/doc/APA300ML.pdf b/doc/APA300ML.pdf
new file mode 100644
index 0000000..b7d2f01
Binary files /dev/null and b/doc/APA300ML.pdf differ
diff --git a/figs/apa300ml_comp_simpler_model.pdf b/figs/apa300ml_comp_simpler_model.pdf
index 65f5566..e126e57 100644
Binary files a/figs/apa300ml_comp_simpler_model.pdf and b/figs/apa300ml_comp_simpler_model.pdf differ
diff --git a/figs/apa300ml_comp_simpler_model.png b/figs/apa300ml_comp_simpler_model.png
index 725e098..aaff69e 100644
Binary files a/figs/apa300ml_comp_simpler_model.png and b/figs/apa300ml_comp_simpler_model.png differ
diff --git a/figs/apa300ml_comp_simpler_simscape.pdf b/figs/apa300ml_comp_simpler_simscape.pdf
index 5ed65b9..5c811d3 100644
Binary files a/figs/apa300ml_comp_simpler_simscape.pdf and b/figs/apa300ml_comp_simpler_simscape.pdf differ
diff --git a/figs/apa300ml_comp_simpler_simscape.png b/figs/apa300ml_comp_simpler_simscape.png
index f90a647..3b34db6 100644
Binary files a/figs/apa300ml_comp_simpler_simscape.png and b/figs/apa300ml_comp_simpler_simscape.png differ
diff --git a/figs/apa300ml_dvf_root_locus.pdf b/figs/apa300ml_dvf_root_locus.pdf
index 7f59327..6b5de5c 100644
Binary files a/figs/apa300ml_dvf_root_locus.pdf and b/figs/apa300ml_dvf_root_locus.pdf differ
diff --git a/figs/apa300ml_dvf_root_locus.png b/figs/apa300ml_dvf_root_locus.png
index 2cd3d7c..e8ee069 100644
Binary files a/figs/apa300ml_dvf_root_locus.png and b/figs/apa300ml_dvf_root_locus.png differ
diff --git a/figs/apa300ml_iff_plant.pdf b/figs/apa300ml_iff_plant.pdf
index f2c981f..c9426eb 100644
Binary files a/figs/apa300ml_iff_plant.pdf and b/figs/apa300ml_iff_plant.pdf differ
diff --git a/figs/apa300ml_iff_plant.png b/figs/apa300ml_iff_plant.png
index 2671c25..06eeaae 100644
Binary files a/figs/apa300ml_iff_plant.png and b/figs/apa300ml_iff_plant.png differ
diff --git a/figs/apa300ml_iff_root_locus.pdf b/figs/apa300ml_iff_root_locus.pdf
index 45f6b25..6576f07 100644
Binary files a/figs/apa300ml_iff_root_locus.pdf and b/figs/apa300ml_iff_root_locus.pdf differ
diff --git a/figs/apa300ml_iff_root_locus.png b/figs/apa300ml_iff_root_locus.png
index 31065b7..b71cc11 100644
Binary files a/figs/apa300ml_iff_root_locus.png and b/figs/apa300ml_iff_root_locus.png differ
diff --git a/figs/apa300ml_plant_dynamics.pdf b/figs/apa300ml_plant_dynamics.pdf
index eb09fea..5f0a5a9 100644
Binary files a/figs/apa300ml_plant_dynamics.pdf and b/figs/apa300ml_plant_dynamics.pdf differ
diff --git a/figs/apa300ml_plant_dynamics.png b/figs/apa300ml_plant_dynamics.png
index 31ef786..e854a81 100644
Binary files a/figs/apa300ml_plant_dynamics.png and b/figs/apa300ml_plant_dynamics.png differ
diff --git a/figs/apa300ml_resonance.pdf b/figs/apa300ml_resonance.pdf
index a6f736c..6e57d78 100644
Binary files a/figs/apa300ml_resonance.pdf and b/figs/apa300ml_resonance.pdf differ
diff --git a/figs/apa300ml_resonance.png b/figs/apa300ml_resonance.png
index 6e356a0..6758e7f 100644
Binary files a/figs/apa300ml_resonance.png and b/figs/apa300ml_resonance.png differ
diff --git a/index.org b/index.org
index ab32741..53e0370 100644
--- a/index.org
+++ b/index.org
@@ -34,796 +34,34 @@
 #+PROPERTY: header-args:latex+ :output-dir figs
 :END:
 
-* Amplified Piezoelectric Actuator - 3D elements
-** Introduction                                                      :ignore:
-The idea here is to:
-- export a FEM of an amplified piezoelectric actuator from Ansys to Matlab
-- import it into a Simscape model
-- compare the obtained dynamics
-- add 10kg mass on top of the actuator and identify the dynamics
-- compare with results from Ansys where 10kg are directly added to the FEM
-
-** 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
-
-#+begin_src matlab
-  addpath('./data/piezo_amplified_3d/');
-#+end_src
-
-#+begin_src matlab :exports none
-  open('piezo_amplified_3d');
-#+end_src
-
-** Import Mass Matrix, Stiffness Matrix, and Interface Nodes Coordinates
-We first extract the stiffness and mass matrices.
-#+begin_src matlab
-  K = extractMatrix('piezo_amplified_3d_K.txt');
-  M = extractMatrix('piezo_amplified_3d_M.txt');
-#+end_src
-
-Then, we extract the coordinates of the interface nodes.
-#+begin_src matlab
-  [int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('piezo_amplified_3d.txt');
-#+end_src
-
-#+begin_src matlab
-  save('./mat/piezo_amplified_3d.mat', 'int_xyz', 'int_i', 'n_xyz', 'n_i', 'nodes', 'M', 'K');
-#+end_src
-
-** Output parameters
-#+begin_src matlab
-  load('./mat/piezo_amplified_3d.mat', 'int_xyz', 'int_i', 'n_xyz', 'n_i', 'nodes', 'M', 'K');
-#+end_src
-
-#+begin_src matlab :exports results :results value table replace :tangle no
-  data2orgtable([length(n_i); length(int_i); size(M,1) - 6*length(int_i); size(M,1)], {'Total number of Nodes', 'Number of interface Nodes', 'Number of Modes', 'Size of M and K matrices'}, {}, ' %.0f ');
-#+end_src
-
-#+RESULTS:
-| Total number of Nodes     | 168959 |
-| Number of interface Nodes |     13 |
-| Number of Modes           |     30 |
-| Size of M and K matrices  |    108 |
-
-#+name: fig:amplified_piezo_interface_nodes
-#+caption: Interface Nodes for the Amplified Piezo Actuator
-[[file:figs/amplified_piezo_interface_nodes.png]]
-
-#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
-  data2orgtable([[1:length(int_i)]', int_i, int_xyz], {}, {'Node i', 'Node Number', 'x [m]', 'y [m]', 'z [m]'}, ' %f ');
-#+end_src
-
-#+caption: Coordinates of the interface nodes
-#+RESULTS:
-| Node i | Node Number |  x [m] | y [m] | z [m] |
-|--------+-------------+--------+-------+-------|
-|    1.0 |    168947.0 |    0.0 |  0.03 |   0.0 |
-|    2.0 |    168949.0 |    0.0 | -0.03 |   0.0 |
-|    3.0 |    168950.0 | -0.035 |   0.0 |   0.0 |
-|    4.0 |    168951.0 | -0.028 |   0.0 |   0.0 |
-|    5.0 |    168952.0 | -0.021 |   0.0 |   0.0 |
-|    6.0 |    168953.0 | -0.014 |   0.0 |   0.0 |
-|    7.0 |    168954.0 | -0.007 |   0.0 |   0.0 |
-|    8.0 |    168955.0 |    0.0 |   0.0 |   0.0 |
-|    9.0 |    168956.0 |  0.007 |   0.0 |   0.0 |
-|   10.0 |    168957.0 |  0.014 |   0.0 |   0.0 |
-|   11.0 |    168958.0 |  0.021 |   0.0 |   0.0 |
-|   12.0 |    168959.0 |  0.035 |   0.0 |   0.0 |
-|   13.0 |    168960.0 |  0.028 |   0.0 |   0.0 |
-
-#+begin_src matlab :exports results :results value table replace :tangle no
-  data2orgtable(K(1:10, 1:10), {}, {}, ' %.1g ');
-#+end_src
-
-#+caption: First 10x10 elements of the Stiffness matrix
-#+RESULTS:
-| 300000000.0 |    -30000.0 |     8000.0 |    -200.0 |   -30.0 | -60000.0 |  20000000.0 |     -4000.0 |      500.0 |        8 |
-|    -30000.0 | 100000000.0 |      400.0 |      30.0 |   200.0 |       -1 |      4000.0 |  -8000000.0 |      800.0 |        7 |
-|      8000.0 |       400.0 | 50000000.0 | -800000.0 |  -300.0 |    -40.0 |       300.0 |       100.0 |  5000000.0 |  40000.0 |
-|      -200.0 |        30.0 |  -800000.0 |   20000.0 |       5 |        1 |       -10.0 |          -2 |   -40000.0 |   -300.0 |
-|       -30.0 |       200.0 |     -300.0 |         5 | 40000.0 |      0.3 |          -4 |       -10.0 |       40.0 |      0.4 |
-|    -60000.0 |          -1 |      -40.0 |         1 |     0.3 |   3000.0 |      7000.0 |         0.8 |         -1 |   0.0003 |
-|  20000000.0 |      4000.0 |      300.0 |     -10.0 |      -4 |   7000.0 | 300000000.0 |     20000.0 |     3000.0 |     80.0 |
-|     -4000.0 |  -8000000.0 |      100.0 |        -2 |   -10.0 |      0.8 |     20000.0 | 100000000.0 |    -4000.0 |   -100.0 |
-|       500.0 |       800.0 |  5000000.0 |  -40000.0 |    40.0 |       -1 |      3000.0 |     -4000.0 | 50000000.0 | 800000.0 |
-|           8 |           7 |    40000.0 |    -300.0 |     0.4 |   0.0003 |        80.0 |      -100.0 |   800000.0 |  20000.0 |
-
-
-#+begin_src matlab :exports results :results value table replace :tangle no
-  data2orgtable(M(1:10, 1:10), {}, {}, ' %.1g ');
-#+end_src
-
-#+caption: First 10x10 elements of the Mass matrix
-#+RESULTS:
-|   0.03 |  2e-06 |  -2e-07 |  1e-08 |  2e-08 | 0.0002 | -0.001 |  2e-07 |  -8e-08 | -9e-10 |
-|  2e-06 |   0.02 |  -5e-07 |  7e-09 |  3e-08 |  2e-08 | -3e-07 | 0.0003 |  -1e-08 |  1e-10 |
-| -2e-07 | -5e-07 |    0.02 | -9e-05 |  4e-09 | -1e-08 |  2e-07 | -2e-08 | -0.0006 | -5e-06 |
-|  1e-08 |  7e-09 |  -9e-05 |  1e-06 |  6e-11 |  4e-10 | -1e-09 |  3e-11 |   5e-06 |  3e-08 |
-|  2e-08 |  3e-08 |   4e-09 |  6e-11 |  1e-06 |  2e-10 | -2e-09 |  2e-10 |  -7e-09 | -4e-11 |
-| 0.0002 |  2e-08 |  -1e-08 |  4e-10 |  2e-10 |  2e-06 | -2e-06 | -1e-09 |  -7e-10 | -9e-12 |
-| -0.001 | -3e-07 |   2e-07 | -1e-09 | -2e-09 | -2e-06 |   0.03 | -2e-06 |  -1e-07 | -5e-09 |
-|  2e-07 | 0.0003 |  -2e-08 |  3e-11 |  2e-10 | -1e-09 | -2e-06 |   0.02 |  -8e-07 | -1e-08 |
-| -8e-08 | -1e-08 | -0.0006 |  5e-06 | -7e-09 | -7e-10 | -1e-07 | -8e-07 |    0.02 |  9e-05 |
-| -9e-10 |  1e-10 |  -5e-06 |  3e-08 | -4e-11 | -9e-12 | -5e-09 | -1e-08 |   9e-05 |  1e-06 |
-
-Using =K=, =M= and =int_xyz=, we can use the =Reduced Order Flexible Solid= simscape block.
-
-
-** Piezoelectric parameters
-Parameters for the APA95ML:
-
-#+begin_src matlab
-  d33 = 3e-10; % Strain constant [m/V]
-  n = 80; % Number of layers per stack
-  eT = 1.6e-7; % Permittivity under constant stress [F/m]
-  sD = 2e-11; % Elastic compliance under constant electric displacement [m2/N]
-  ka = 235e6; % Stack stiffness [N/m]
-  C = 5e-6; % Stack capactiance [F]
-#+end_src
-
-#+begin_src matlab
-  na = 2; % Number of stacks used as actuator
-  ns = 1; % Number of stacks used as force sensor
-#+end_src
-
-The ratio of the developed force to applied voltage is $d_{33} n k_a$ in [N/V].
-We denote this constant by $g_a$ and:
-\[ F_a = g_a V_a, \quad g_a = d_{33} n k_a \]
- #+begin_src matlab :results replace value
-  d33*(na*n)*(ka/(na + ns)) % [N/V]
-#+end_src
-
-#+RESULTS:
-: 3.76
-
-From cite:fleming14_desig_model_contr_nanop_system (page 123), the relation between relative displacement and generated voltage is:
-\[ V_s = \frac{d_{33}}{\epsilon^T s^D n} \Delta h \]
-where:
-- $V_s$: measured voltage [V]
-- $d_{33}$: strain constant [m/V]
-- $\epsilon^T$: permittivity under constant stress [F/m]
-- $s^D$: elastic compliance under constant electric displacement [m^2/N]
-- $n$: number of layers
-- $\Delta h$: relative displacement [m]
-#+begin_src matlab :results replace value
-  1e-6*d33/(eT*sD*ns*n) % [V/um]
-#+end_src
-
-#+RESULTS:
-: 1.1719
-
-** Identification of the Dynamics
-The flexible element is imported using the =Reduced Order Flexible Solid= simscape block.
-
-To model the actuator, an =Internal Force= block is added between the nodes 3 and 12.
-A =Relative Motion Sensor= block is added between the nodes 1 and 2 to measure the displacement and the amplified piezo.
-
-One mass is fixed at one end of the piezo-electric stack actuator, the other end is fixed to the world frame.
-
-We first set the mass to be zero.
-#+begin_src matlab
-  m = 0.01;
-#+end_src
-
-The dynamics is identified from the applied force to the measured relative displacement.
-#+begin_src matlab
-  %% Name of the Simulink File
-  mdl = 'piezo_amplified_3d';
-
-  %% Input/Output definition
-  clear io; io_i = 1;
-  io(io_i) = linio([mdl, '/F'], 1, 'openinput');  io_i = io_i + 1;
-  io(io_i) = linio([mdl, '/y'], 1, 'openoutput'); io_i = io_i + 1;
-
-  Gh = -linearize(mdl, io);
-#+end_src
-
-Then, we add 10Kg of mass:
-#+begin_src matlab
-  m = 5;
-#+end_src
-
-And the dynamics is identified.
-
-The two identified dynamics are compared in Figure [[fig:dynamics_act_disp_comp_mass]].
-#+begin_src matlab
-  %% Name of the Simulink File
-  mdl = 'piezo_amplified_3d';
-
-  %% Input/Output definition
-  clear io; io_i = 1;
-  io(io_i) = linio([mdl, '/F'], 1, 'openinput');  io_i = io_i + 1;
-  io(io_i) = linio([mdl, '/y'], 1, 'openoutput'); io_i = io_i + 1;
-
-  Ghm = -linearize(mdl, io);
-#+end_src
-
-#+begin_src matlab :exports none
-  freqs = logspace(0, 4, 5000);
-
-  figure;
-
-  ax1 = subplot(2,1,1);
-  hold on;
-  plot(freqs, abs(squeeze(freqresp(Gh, freqs, 'Hz'))), '-');
-  plot(freqs, abs(squeeze(freqresp(Ghm, freqs, 'Hz'))), '-');
-  hold off;
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
-  hold off;
-
-  ax2 = subplot(2,1,2);
-  hold on;
-  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gh, freqs, 'Hz')))), '-', ...
-       'DisplayName', '$m = 0kg$');
-  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Ghm, freqs, 'Hz')))), '-', ...
-       'DisplayName', '$m = 10kg$');
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
-  yticks(-360:90:360);
-  ylim([-360 0]);
-  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
-  hold off;
-  linkaxes([ax1,ax2],'x');
-  xlim([freqs(1), freqs(end)]);
-  legend('location', 'southwest');
-#+end_src
-
-#+begin_src matlab :tangle no :exports results :results file replace
-  exportFig('figs/dynamics_act_disp_comp_mass.pdf', 'width', 'full', 'height', 'full');
-#+end_src
-
-#+name: fig:dynamics_act_disp_comp_mass
-#+caption: Dynamics from $F$ to $d$ without a payload and with a 10kg payload
-#+RESULTS:
-[[file:figs/dynamics_act_disp_comp_mass.png]]
-
-** Comparison with Ansys
-Let's import the results from an Harmonic response analysis in Ansys.
-#+begin_src matlab
-  Gresp0 = readtable('FEA_HarmResponse_00kg.txt');
-  Gresp10 = readtable('FEA_HarmResponse_10kg.txt');
-#+end_src
-
-The obtained dynamics from the Simscape model and from the Ansys analysis are compare in Figure [[fig:dynamics_force_disp_comp_anasys]].
-
-#+begin_src matlab :exports none
-  freqs = logspace(1, 5, 1000);
-
-  figure;
-
-  ax1 = subplot(2,1,1);
-  hold on;
-  set(gca,'ColorOrderIndex',1)
-  plot(freqs, abs(squeeze(freqresp(Gh, freqs, 'Hz'))), '-');
-  set(gca,'ColorOrderIndex',1)
-  plot(Gresp0{:, 2}, 1e-3*Gresp0{:, 3}, '--');
-
-  set(gca,'ColorOrderIndex',2)
-  plot(freqs, abs(squeeze(freqresp(Ghm, freqs, 'Hz'))), '-');
-  set(gca,'ColorOrderIndex',2)
-  plot(Gresp0{:, 2}, 1e-3*Gresp10{:, 3}, '--');
-  hold off;
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
-  hold off;
-
-  ax2 = subplot(2,1,2);
-  hold on;
-  set(gca,'ColorOrderIndex',1)
-  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gh, freqs, 'Hz')))), '-', ...
-       'DisplayName', '$m = 0kg$, Simscape');
-  set(gca,'ColorOrderIndex',1)
-  plot(Gresp0{:, 2}, 180/pi*unwrap(pi/180*Gresp0{:, 4}), '--', ...
-       'DisplayName', '$m = 0kg$, Ansys');
-  set(gca,'ColorOrderIndex',2)
-  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Ghm, freqs, 'Hz')))), '-', ...
-       'DisplayName', '$m = 10kg$, Simscape');
-  set(gca,'ColorOrderIndex',2)
-  plot(Gresp0{:, 2}, 180/pi*unwrap(pi/180*Gresp10{:, 4}), '--', ...
-       'DisplayName', '$m = 10kg$, Ansys');
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
-  yticks(-360:90:360);
-  ylim([-390 30]);
-  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
-  hold off;
-  legend('location', 'southwest');
-
-  linkaxes([ax1,ax2],'x');
-  xlim([freqs(1), freqs(end)]);
-#+end_src
-
-#+begin_src matlab :tangle no :exports results :results file replace
-  exportFig('figs/dynamics_force_disp_comp_anasys.pdf', 'width', 'full', 'height', 'full');
-#+end_src
-
-#+name: fig:dynamics_force_disp_comp_anasys
-#+caption: Comparison of the obtained dynamics using Simscape with the harmonic response analysis using Ansys
-#+RESULTS:
-[[file:figs/dynamics_force_disp_comp_anasys.png]]
-
-** Force Sensor
-The dynamics is identified from internal forces applied between nodes 3 and 11 to the relative displacement of nodes 11 and 13.
-
-The obtained dynamics is shown in Figure [[fig:dynamics_force_force_sensor_comp_mass]].
-
-#+begin_src matlab
-  m = 0;
-#+end_src
-
-#+begin_src matlab
-  %% Name of the Simulink File
-  mdl = 'piezo_amplified_3d';
-
-  %% Input/Output definition
-  clear io; io_i = 1;
-  io(io_i) = linio([mdl, '/Fa'], 1, 'openinput');  io_i = io_i + 1;
-  io(io_i) = linio([mdl, '/Fs'], 1, 'openoutput'); io_i = io_i + 1;
-
-  Gf = linearize(mdl, io);
-#+end_src
-
-#+begin_src matlab
-  m = 10;
-#+end_src
-
-#+begin_src matlab
-  %% Name of the Simulink File
-  mdl = 'piezo_amplified_3d';
-
-  %% Input/Output definition
-  clear io; io_i = 1;
-  io(io_i) = linio([mdl, '/Fa'], 1, 'openinput');  io_i = io_i + 1;
-  io(io_i) = linio([mdl, '/Fs'], 1, 'openoutput'); io_i = io_i + 1;
-
-  Gfm = linearize(mdl, io);
-#+end_src
-
-#+begin_src matlab :exports none
-  freqs = logspace(1, 5, 1000);
-
-  figure;
-
-  ax1 = subplot(2,1,1);
-  hold on;
-  plot(freqs, abs(squeeze(freqresp(Gf, freqs, 'Hz'))), '-');
-  plot(freqs, abs(squeeze(freqresp(Gfm, freqs, 'Hz'))), '-');
-  hold off;
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
-  hold off;
-
-  ax2 = subplot(2,1,2);
-  hold on;
-  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gf, freqs, 'Hz')))), '-', ...
-       'DisplayName', '$m = 0kg$');
-  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gfm, freqs, 'Hz')))), '-', ...
-       'DisplayName', '$m = 10kg$');
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
-  yticks(-360:90:360);
-  ylim([-390 30]);
-  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
-  hold off;
-  linkaxes([ax1,ax2],'x');
-  xlim([freqs(1), freqs(end)]);
-  legend('location', 'southwest');
-#+end_src
-
-#+begin_src matlab :tangle no :exports results :results file replace
-  exportFig('figs/dynamics_force_force_sensor_comp_mass.pdf', 'width', 'full', 'height', 'full');
-#+end_src
-
-#+name: fig:dynamics_force_force_sensor_comp_mass
-#+caption: Dynamics from $F$ to $F_m$ for $m=0$ and $m = 10kg$
-#+RESULTS:
-[[file:figs/dynamics_force_force_sensor_comp_mass.png]]
-
-** Distributed Actuator
-#+begin_src matlab
-  m = 0;
-#+end_src
-
-The dynamics is identified from the applied force to the measured relative displacement.
-#+begin_src matlab
-  %% Name of the Simulink File
-  mdl = 'piezo_amplified_3d_distri';
-
-  %% Input/Output definition
-  clear io; io_i = 1;
-  io(io_i) = linio([mdl, '/F'], 1, 'openinput');  io_i = io_i + 1;
-  io(io_i) = linio([mdl, '/y'], 1, 'openoutput'); io_i = io_i + 1;
-
-  Gd = linearize(mdl, io);
-#+end_src
-
-Then, we add 10Kg of mass:
-#+begin_src matlab
-  m = 10;
-#+end_src
-
-And the dynamics is identified.
-#+begin_src matlab
-  %% Name of the Simulink File
-  mdl = 'piezo_amplified_3d_distri';
-
-  %% Input/Output definition
-  clear io; io_i = 1;
-  io(io_i) = linio([mdl, '/F'], 1, 'openinput');  io_i = io_i + 1;
-  io(io_i) = linio([mdl, '/y'], 1, 'openoutput'); io_i = io_i + 1;
-
-  Gdm = linearize(mdl, io);
-#+end_src
-
-#+begin_src matlab :exports none
-  freqs = logspace(1, 5, 5000);
-
-  figure;
-
-  ax1 = subplot(2,1,1);
-  hold on;
-  plot(freqs, abs(squeeze(freqresp(Gh, freqs, 'Hz'))), '-');
-  plot(freqs, abs(squeeze(freqresp(Ghm, freqs, 'Hz'))), '-');
-  set(gca,'ColorOrderIndex',1)
-  plot(freqs, abs(squeeze(freqresp(Gd, freqs, 'Hz'))), '--');
-  plot(freqs, abs(squeeze(freqresp(Gdm, freqs, 'Hz'))), '--');
-  hold off;
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
-  hold off;
-
-  ax2 = subplot(2,1,2);
-  hold on;
-  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gh, freqs, 'Hz')))), '-', ...
-       'DisplayName', '$m = 0kg$');
-  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Ghm, freqs, 'Hz')))), '-', ...
-       'DisplayName', '$m = 10kg$');
-  set(gca,'ColorOrderIndex',1)
-  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gd, freqs, 'Hz')))), '--', ...
-       'DisplayName', '$m = 0kg$, distri');
-  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gdm, freqs, 'Hz')))), '--', ...
-       'DisplayName', '$m = 10kg$, distri');
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
-  yticks(-360:90:360);
-  ylim([-360 0]);
-  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
-  hold off;
-  linkaxes([ax1,ax2],'x');
-  xlim([freqs(1), freqs(end)]);
-  legend('location', 'southwest');
-#+end_src
-
-** Distributed Actuator and Force Sensor
-#+begin_src matlab
-  m = 0;
-#+end_src
-
-#+begin_src matlab
-  %% Name of the Simulink File
-  mdl = 'piezo_amplified_3d_distri_act_sens';
-
-  %% Input/Output definition
-  clear io; io_i = 1;
-  io(io_i) = linio([mdl, '/F'], 1, 'openinput');  io_i = io_i + 1;
-  io(io_i) = linio([mdl, '/Fm'], 1, 'openoutput'); io_i = io_i + 1;
-
-  Gfd = linearize(mdl, io);
-#+end_src
-
-#+begin_src matlab
-  m = 10;
-#+end_src
-
-#+begin_src matlab
-  %% Name of the Simulink File
-  mdl = 'piezo_amplified_3d_distri_act_sens';
-
-  %% Input/Output definition
-  clear io; io_i = 1;
-  io(io_i) = linio([mdl, '/F'], 1, 'openinput');  io_i = io_i + 1;
-  io(io_i) = linio([mdl, '/Fm'], 1, 'openoutput'); io_i = io_i + 1;
-
-  Gfdm = linearize(mdl, io);
-#+end_src
-
-#+begin_src matlab :exports none
-  freqs = logspace(1, 5, 1000);
-
-  figure;
-
-  ax1 = subplot(2,1,1);
-  hold on;
-  plot(freqs, abs(squeeze(freqresp(Gf, freqs, 'Hz'))), '-');
-  plot(freqs, abs(squeeze(freqresp(Gfm, freqs, 'Hz'))), '-');
-  set(gca,'ColorOrderIndex',1)
-  plot(freqs, abs(squeeze(freqresp(Gfd, freqs, 'Hz'))), '--');
-  plot(freqs, abs(squeeze(freqresp(Gfdm, freqs, 'Hz'))), '--');
-  hold off;
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
-  hold off;
-
-  ax2 = subplot(2,1,2);
-  hold on;
-  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gf, freqs, 'Hz')))), '-', ...
-       'DisplayName', '$m = 0kg$');
-  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gfm, freqs, 'Hz')))), '-', ...
-       'DisplayName', '$m = 10kg$');
-  set(gca,'ColorOrderIndex',1)
-  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gfd, freqs, 'Hz')))), '--', ...
-       'DisplayName', '$m = 0kg$, distri');
-  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gfdm, freqs, 'Hz')))), '--', ...
-       'DisplayName', '$m = 10kg$, distri');
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
-  yticks(-360:90:360);
-  ylim([-390 30]);
-  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
-  hold off;
-  linkaxes([ax1,ax2],'x');
-  xlim([freqs(1), freqs(end)]);
-  legend('location', 'southwest');
-#+end_src
-
-** Dynamics from input voltage to displacement
-#+begin_src matlab
-  m = 5;
-#+end_src
-
-And the dynamics is identified.
-
-The two identified dynamics are compared in Figure [[fig:dynamics_act_disp_comp_mass]].
-#+begin_src matlab
-  %% Name of the Simulink File
-  mdl = 'piezo_amplified_3d';
-
-  %% Input/Output definition
-  clear io; io_i = 1;
-  io(io_i) = linio([mdl, '/V'], 1, 'openinput');  io_i = io_i + 1;
-  io(io_i) = linio([mdl, '/y'], 1, 'openoutput'); io_i = io_i + 1;
-
-  G = -linearize(mdl, io);
-#+end_src
-
-#+begin_src matlab :exports none
-  freqs = logspace(1, 4, 5000);
-
-  figure;
-
-  ax1 = subplot(2,1,1);
-  hold on;
-  plot(freqs, abs(squeeze(freqresp(G, freqs, 'Hz'))));
-  hold off;
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
-  hold off;
-
-  ax2 = subplot(2,1,2);
-  hold on;
-  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G, freqs, 'Hz')))));
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
-  yticks(-360:90:360);
-  ylim([-360 0]);
-  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
-  hold off;
-  linkaxes([ax1,ax2],'x');
-  % xlim([10, 5e3]);
-#+end_src
-
-#+begin_src matlab
-  save('../test-bench-apa/mat/fem_model_5kg.mat', 'G')
-#+end_src
-
-** Dynamics from input voltage to output voltage
-#+begin_src matlab
-  m = 5;
-#+end_src
-
-#+begin_src matlab
-  %% Name of the Simulink File
-  mdl = 'piezo_amplified_3d';
-
-  %% Input/Output definition
-  clear io; io_i = 1;
-  io(io_i) = linio([mdl, '/Fa'], 1, 'openinput');  io_i = io_i + 1;
-  io(io_i) = linio([mdl, '/dL'], 1, 'openoutput'); io_i = io_i + 1;
-
-  G = -linearize(mdl, io);
-#+end_src
-
-#+begin_src matlab :exports none
-  freqs = logspace(1, 4, 5000);
-
-  figure;
-
-  ax1 = subplot(2,1,1);
-  hold on;
-  plot(freqs, abs(squeeze(freqresp(G(2,1), freqs, 'Hz'))));
-  plot(freqs, abs(squeeze(freqresp(G(3,1), freqs, 'Hz'))));
-  plot(freqs, abs(squeeze(freqresp(G(3,1) + G(3,2), freqs, 'Hz'))));
-  hold off;
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
-  hold off;
-
-  ax2 = subplot(2,1,2);
-  hold on;
-  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G(2,1), freqs, 'Hz')))));
-  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G(3,1), freqs, 'Hz')))));
-  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G(3,1) + G(3,2), freqs, 'Hz')))));
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
-  yticks(-360:90:360);
-  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
-  hold off;
-  linkaxes([ax1,ax2],'x');
-  xlim([10, 1e4]);
-#+end_src
-
-** Identification for a simpler model
-The goal in this section is to identify the parameters of a simple APA model from the FEM.
-This can be useful is a lower order model is to be used for simulations.
-
-The presented model is based on cite:souleille18_concep_activ_mount_space_applic.
-
-The model represents the Amplified Piezo Actuator (APA) from Cedrat-Technologies (Figure [[fig:souleille18_model_piezo]]).
-The parameters are shown in the table below.
-
-#+name: fig:souleille18_model_piezo
-#+caption: Picture of an APA100M from Cedrat Technologies. Simplified model of a one DoF payload mounted on such isolator
-[[file:./figs/souleille18_model_piezo.png]]
-
-#+caption: Parameters used for the model of the APA 100M
-|       | Meaning                                                        |
-|-------+----------------------------------------------------------------|
-| $k_e$ | Stiffness used to adjust the pole of the isolator              |
-| $k_1$ | Stiffness of the metallic suspension when the stack is removed |
-| $k_a$ | Stiffness of the actuator                                      |
-| $c_1$ | Added viscous damping                                          |
-
-The goal is to determine $k_e$, $k_a$ and $k_1$ so that the simplified model fits the FEM model.
-
-\[ \alpha = \frac{x_1}{f}(\omega=0) = \frac{\frac{k_e}{k_e + k_a}}{k_1 + \frac{k_e k_a}{k_e + k_a}} \]
-\[ \beta = \frac{x_1}{F}(\omega=0) = \frac{1}{k_1 + \frac{k_e k_a}{k_e + k_a}} \]
-
-If we can fix $k_a$, we can determine $k_e$ and $k_1$ with:
-\[ k_e = \frac{k_a}{\frac{\beta}{\alpha} - 1} \]
-\[ k_1 = \frac{1}{\beta} - \frac{k_e k_a}{k_e + k_a} \]
-
-#+begin_src matlab :exports none
-  m = 10;
-#+end_src
-
-#+begin_src matlab :exports none
-  %% Name of the Simulink File
-  mdl = 'piezo_amplified_3d_identification';
-
-  %% Input/Output definition
-  clear io; io_i = 1;
-  io(io_i) = linio([mdl, '/Fd'], 1, 'openinput');  io_i = io_i + 1; % External Vertical Force [N]
-  io(io_i) = linio([mdl, '/w'],  1, 'openinput');  io_i = io_i + 1; % Base Motion [m]
-  io(io_i) = linio([mdl, '/Fa'], 1, 'openinput');  io_i = io_i + 1; % Actuator Force [N]
-  io(io_i) = linio([mdl, '/z'],  1, 'openoutput'); io_i = io_i + 1; % Vertical Displacement [m]
-  io(io_i) = linio([mdl, '/Vs'], 1, 'openoutput'); io_i = io_i + 1; % Force Sensor [V]
-
-  G = linearize(mdl, io);
-
-  G.InputName = {'Fd', 'w', 'Fa'};
-  G.OutputName = {'y', 'Fs'};
-#+end_src
-
-From the identified dynamics, compute $\alpha$ and $\beta$
-#+begin_src matlab
-  alpha = abs(dcgain(G('y', 'Fa')));
-  beta  = abs(dcgain(G('y', 'Fd')));
-#+end_src
-
-$k_a$ is estimated using the following formula:
-#+begin_src matlab
-  ka = 0.9/abs(dcgain(G('y', 'Fa')));
-#+end_src
-The factor can be adjusted to better match the curves.
-
-Then $k_e$ and $k_1$ are computed.
-#+begin_src matlab
-  ke = ka/(beta/alpha - 1);
-  k1 = 1/beta - ke*ka/(ke + ka);
-#+end_src
-
-#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
-  data2orgtable(1e-6*[ka; ke; k1], {'ka', 'ke', 'k1'}, {'Value [N/um]'}, ' %.1f ');
-#+end_src
-
-#+RESULTS:
-|    | Value [N/um] |
-|----+--------------|
-| ka |         54.9 |
-| ke |         25.1 |
-| k1 |          4.3 |
-
-The damping in the system is adjusted to match the FEM model if necessary.
-#+begin_src matlab
-  c1 = 1e2;
-#+end_src
-
-Analytical model of the simpler system:
-#+begin_src matlab
-  Ga = 1/(m*s^2 + k1 + c1*s + ke*ka/(ke + ka)) * ...
-       [ 1 ,              k1 + c1*s + ke*ka/(ke + ka)  , ke/(ke + ka) ;
-        -ke*ka/(ke + ka), ke*ka/(ke + ka)*m*s^2 ,       -ke/(ke + ka)*(m*s^2 + c1*s + k1)];
-
-  Ga.InputName = {'Fd', 'w', 'Fa'};
-  Ga.OutputName = {'y', 'Fs'};
-#+end_src
-
-Adjust the DC gain for the force sensor:
-#+begin_src matlab
-  F_gain = dcgain(G('Fs', 'Fd'))/dcgain(Ga('Fs', 'Fd'));
-#+end_src
-
-#+begin_src matlab :exports none
-  freqs = logspace(0, 5, 1000);
-
-  figure;
-
-  ax1 = subplot(2, 2, 1);
-  title('$\displaystyle \frac{x_1}{f}$')
-  hold on;
-  plot(freqs, abs(squeeze(freqresp(G( 'y', 'Fa'), freqs, 'Hz'))));
-  plot(freqs, abs(squeeze(freqresp(Ga('y', 'Fa'), freqs, 'Hz'))));
-  hold off;
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude [m/N]');xlabel('Frequency [Hz]');
-
-  ax2 = subplot(2, 2, 2);
-  title('$\displaystyle \frac{x_1}{F}$')
-  hold on;
-  plot(freqs, abs(squeeze(freqresp(G( 'y', 'Fd'), freqs, 'Hz'))));
-  plot(freqs, abs(squeeze(freqresp(Ga('y', 'Fd'), freqs, 'Hz'))));
-  hold off;
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude [m/N]');xlabel('Frequency [Hz]');
-
-  ax3 = subplot(2, 2, 3);
-  title('$\displaystyle \frac{F_s}{f}$')
-  hold on;
-  plot(freqs, abs(squeeze(freqresp(G( 'Fs', 'Fa'), freqs, 'Hz'))));
-  plot(freqs, abs(squeeze(freqresp(F_gain*Ga('Fs', 'Fa'), freqs, 'Hz'))));
-  hold off;
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude [m/N]');xlabel('Frequency [Hz]');
-
-  ax4 = subplot(2, 2, 4);
-  title('$\displaystyle \frac{F_s}{F}$')
-  hold on;
-  plot(freqs, abs(squeeze(freqresp(G( 'Fs', 'Fd'), freqs, 'Hz'))));
-  plot(freqs, abs(squeeze(freqresp(F_gain*Ga('Fs', 'Fd'), freqs, 'Hz'))));
-  hold off;
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
-
-  linkaxes([ax1,ax2,ax3,ax4],'x');
-#+end_src
-
-#+begin_src matlab :tangle no :exports results :results file replace
-  exportFig('figs/apa95ml_comp_simpler_model.pdf', 'width', 'full', 'height', 'full');
-#+end_src
-
-#+name: fig:apa95ml_comp_simpler_model
-#+caption: Comparison of the Dynamics between the FEM model and the simplified one
-#+RESULTS:
-[[file:figs/apa95ml_comp_simpler_model.png]]
-
-We save the parameters of the simplified model for the APA95ML:
-#+begin_src matlab
-  save('./mat/APA95ML_simplified_model.mat', 'ka', 'ke', 'k1', 'c1', 'F_gain');
-#+end_src
+* Introduction                                                        :ignore:
+
+In this document, Finite Element Models (FEM) of parts of the Nano-Hexapod are developed and integrated into Simscape for dynamical analysis.
+
+- Section [[sec:APA300ML]]:
+  A super-element of the Amplified Piezoelectric Actuator APA300ML used for the NASS is exported using Ansys and imported in Simscape.
+  The static and dynamical properties of the APA300ML are then estimated using the Simscape model.
+- Section [[sec:first_flexible_joint]]:
+  A first geometry of a Flexible joint is modelled and its characteristics are identified from the Stiffness matrix as well as from the Simscape model.
+- Section [[sec:optimized_flexible_joint]]:
+  An optimized flexible joint is developed for the Nano-Hexapod and is then imported in a Simscape model.
+- Section [[sec:integral_force_feedback]]:
+- Section [[sec:strut_fem]]:
+  A super element of a complete strut is exported.
 
 * APA300ML
+:PROPERTIES:
+:header-args:matlab+: :tangle matlab/APA300ML.m
+:END:
+<<sec:APA300ML>>
 ** Introduction                                                      :ignore:
 
+In this section, the Amplified Piezoelectric Actuator APA300ML ([[file:doc/APA300ML.pdf][doc]]) is modeled using a Finite Element Software.
+Then a /super element/ is exported and imported in Simscape where its dynamic is studied.
+
+A 3D view of the Amplified Piezoelectric Actuator (APA300ML) is shown in Figure [[fig:apa300ml_ansys]].
+The remote point used are also shown in this figure.
+
 #+name: fig:apa300ml_ansys
 #+caption: Ansys FEM of the APA300ML
 [[file:figs/apa300ml_ansys.jpg]]
@@ -837,12 +75,17 @@ We save the parameters of the simplified model for the APA95ML:
   <<matlab-init>>
 #+end_src
 
-#+begin_src matlab
-  addpath('./data/APA300ML_new/');
+#+begin_src matlab :tangle no
+  addpath('matlab/');
+  addpath('matlab/APA300ML/');
 #+end_src
 
-#+begin_src matlab :exports none
-  open('APA300ML_test_bench');
+#+begin_src matlab :eval no
+  addpath('APA300ML/');
+#+end_src
+
+#+begin_src matlab
+  open('APA300ML.slx');
 #+end_src
 
 ** Import Mass Matrix, Stiffness Matrix, and Interface Nodes Coordinates
@@ -852,46 +95,6 @@ We first extract the stiffness and mass matrices.
   M = readmatrix('mat_M.CSV');
 #+end_src
 
-Then, we extract the coordinates of the interface nodes.
-#+begin_src matlab
-  [int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('out_nodes_3D.txt');
-#+end_src
-
-#+begin_src matlab
-  save('./mat/APA300ML.mat', 'int_xyz', 'int_i', 'n_xyz', 'n_i', 'nodes', 'M', 'K');
-#+end_src
-
-** Output parameters
-#+begin_src matlab
-  load('./mat/APA300ML.mat', 'int_xyz', 'int_i', 'n_xyz', 'n_i', 'nodes', 'M', 'K');
-#+end_src
-
-#+begin_src matlab :exports results :results value table replace :tangle no
-  data2orgtable([length(n_i); length(int_i); size(M,1) - 6*length(int_i); size(M,1)], {'Total number of Nodes', 'Number of interface Nodes', 'Number of Modes', 'Size of M and K matrices'}, {}, ' %.0f ');
-#+end_src
-
-#+RESULTS:
-| Total number of Nodes     |   7 |
-| Number of interface Nodes |   7 |
-| Number of Modes           | 120 |
-| Size of M and K matrices  | 162 |
-
-#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
-  data2orgtable([[1:length(int_i)]', int_i, int_xyz], {}, {'Node i', 'Node Number', 'x [m]', 'y [m]', 'z [m]'}, ' %f ');
-#+end_src
-
-#+caption: Coordinates of the interface nodes
-#+RESULTS:
-| Node i | Node Number |   x [m] | y [m] |  z [m] |
-|--------+-------------+---------+-------+--------|
-|    1.0 |    697783.0 |     0.0 |   0.0 | -0.015 |
-|    2.0 |    697784.0 |     0.0 |   0.0 |  0.015 |
-|    3.0 |    697785.0 | -0.0325 |   0.0 |    0.0 |
-|    4.0 |    697786.0 | -0.0125 |   0.0 |    0.0 |
-|    5.0 |    697787.0 | -0.0075 |   0.0 |    0.0 |
-|    6.0 |    697788.0 |  0.0125 |   0.0 |    0.0 |
-|    7.0 |    697789.0 |  0.0325 |   0.0 |    0.0 |
-
 #+begin_src matlab :exports results :results value table replace :tangle no
   data2orgtable(K(1:10, 1:10), {}, {}, ' %.1g ');
 #+end_src
@@ -927,50 +130,103 @@ Then, we extract the coordinates of the interface nodes.
 |   6e-07 |   1e-08 |  6e-05 | -4e-10 |  1e-09 | -4e-10 |  -3e-07 |  -4e-07 |  0.009 | -2e-10 |
 |  -3e-09 |  -9e-07 |  1e-10 |  3e-09 | -8e-12 | -5e-12 |  -2e-08 |   2e-05 | -2e-10 |  3e-07 |
 
-Using =K=, =M= and =int_xyz=, we can use the =Reduced Order Flexible Solid= simscape block.
+
+Then, we extract the coordinates of the interface nodes.
+#+begin_src matlab
+  [int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('out_nodes_3D.txt');
+#+end_src
+
+#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
+  data2orgtable([[1:length(int_i)]', int_i, int_xyz], {}, {'Node i', 'Node Number', 'x [m]', 'y [m]', 'z [m]'}, ' %f ');
+#+end_src
+
+#+caption: Coordinates of the interface nodes
+#+RESULTS:
+| Node i | Node Number |   x [m] | y [m] |  z [m] |
+|--------+-------------+---------+-------+--------|
+|    1.0 |    697783.0 |     0.0 |   0.0 | -0.015 |
+|    2.0 |    697784.0 |     0.0 |   0.0 |  0.015 |
+|    3.0 |    697785.0 | -0.0325 |   0.0 |    0.0 |
+|    4.0 |    697786.0 | -0.0125 |   0.0 |    0.0 |
+|    5.0 |    697787.0 | -0.0075 |   0.0 |    0.0 |
+|    6.0 |    697788.0 |  0.0125 |   0.0 |    0.0 |
+|    7.0 |    697789.0 |  0.0325 |   0.0 |    0.0 |
+
+#+begin_src matlab :exports results :results value table replace :tangle no
+  data2orgtable([length(n_i); length(int_i); size(M,1) - 6*length(int_i); size(M,1)], {'Total number of Nodes', 'Number of interface Nodes', 'Number of Modes', 'Size of M and K matrices'}, {}, ' %.0f ');
+#+end_src
+
+#+caption: Some extracted parameters of the FEM
+#+RESULTS:
+| Total number of Nodes     |   7 |
+| Number of interface Nodes |   7 |
+| Number of Modes           | 120 |
+| Size of M and K matrices  | 162 |
+
+Using =K=, =M= and =int_xyz=, we can now use the =Reduced Order Flexible Solid= simscape block.
 
 ** Piezoelectric parameters
-Parameters for the APA300ML:
-
+In order to make the conversion from applied voltage to generated force or from the strain to the generated voltage, we need to defined some parameters corresponding to the piezoelectric material:
 #+begin_src matlab
-  d33 = 3e-10; % Strain constant [m/V]
-  n = 80; % Number of layers per stack
-  eT = 1.6e-8; % Permittivity under constant stress [F/m]
-  sD = 2e-11; % Elastic compliance under constant electric displacement [m2/N]
-  ka = 235e6; % Stack stiffness [N/m]
-  C = 5e-6; % Stack capactiance [F]
+  d33 = 300e-12; % Strain constant [m/V]
+  n   = 80;      % Number of layers per stack
+  eT  = 1.6e-8;  % Permittivity under constant stress [F/m]
+  sD  = 1e-11;   % Compliance under constant electric displacement [m2/N]
+  ka  = 235e6;   % Stack stiffness [N/m]
+  C   = 5e-6;    % Stack capactiance [F]
 #+end_src
 
-#+begin_src matlab
-  na = 2; % Number of stacks used as actuator
-  ns = 1; % Number of stacks used as force sensor
-#+end_src
-
-The ratio of the developed force to applied voltage is $d_{33} n k_a$ in [N/V].
-We denote this constant by $g_a$ and:
-\[ F_a = g_a V_a, \quad g_a = d_{33} n k_a \]
-#+begin_src matlab :results replace value
-  d33*(na*n)*(ka/(na + ns)) % [N/V]
-#+end_src
-
-#+RESULTS:
-: 3.76
-
-From cite:fleming14_desig_model_contr_nanop_system (page 123), the relation between relative displacement and generated voltage is:
-\[ V_s = \frac{d_{33}}{\epsilon^T s^D n} \Delta h \]
+The ratio of the developed force to applied voltage is:
+#+name: eq:piezo_voltage_to_force
+\begin{equation}
+  F_a = g_a V_a, \quad g_a = d_{33} n k_a
+\end{equation}
 where:
-- $V_s$: measured voltage [V]
-- $d_{33}$: strain constant [m/V]
-- $\epsilon^T$: permittivity under constant stress [F/m]
-- $s^D$: elastic compliance under constant electric displacement [m^2/N]
-- $n$: number of layers
-- $\Delta h$: relative displacement [m]
+- $F_a$: developed force in [N]
+- $n$: number of layers of the actuator stack
+- $d_{33}$: strain constant in [m/V]
+- $k_a$: actuator stack stiffness in [N/m]
+- $V_a$: applied voltage in [V]
+
+If we take the numerical values, we obtain:
 #+begin_src matlab :results replace value
-  1e-6*d33/(eT*sD*ns*n) % [V/um]
+  d33*n*ka % [N/V]
 #+end_src
 
 #+RESULTS:
-: 11.719
+: 5.64
+
+From cite:fleming14_desig_model_contr_nanop_system (page 123), the relation between relative displacement of the sensor stack and generated voltage is:
+#+name: eq:piezo_strain_to_voltage
+\begin{equation}
+  V_s = \frac{d_{33}}{\epsilon^T s^D n} \Delta h
+\end{equation}
+where:
+- $V_s$: measured voltage in [V]
+- $d_{33}$: strain constant in [m/V]
+- $\epsilon^T$: permittivity under constant stress in [F/m]
+- $s^D$: elastic compliance under constant electric displacement in [m^2/N]
+- $n$: number of layers of the sensor stack
+- $\Delta h$: relative displacement in [m]
+
+If we take the numerical values, we obtain:
+#+begin_src matlab :results replace value
+  1e-6*d33/(eT*sD*n) % [V/um]
+#+end_src
+
+#+RESULTS:
+: 23.438
+
+** Simscape Model
+The flexible element is imported using the =Reduced Order Flexible Solid= simscape block.
+
+Let's say we use two stacks as a force sensor and one stack as an actuator:
+- A =Relative Motion Sensor= block is added between the nodes A and C
+- An =Internal Force= block is added between the remote points E and B
+
+The interface nodes are shown in Figure [[fig:apa300ml_ansys]].
+
+One mass is fixed at one end of the piezo-electric stack actuator (remove point F), the other end is fixed to the world frame (remote point G).
 
 ** Identification of the APA Characteristics
 *** Stiffness
@@ -982,12 +238,12 @@ The transfer function from vertical external force to the relative vertical disp
 
 #+begin_src matlab :exports none
   %% Name of the Simulink File
-  mdl = 'APA300ML_test_bench';
+  mdl = 'APA300ML';
 
   %% Input/Output definition
   clear io; io_i = 1;
   io(io_i) = linio([mdl, '/Fd'], 1, 'openinput');  io_i = io_i + 1;
-  io(io_i) = linio([mdl, '/y'], 1, 'openoutput'); io_i = io_i + 1;
+  io(io_i) = linio([mdl, '/z'],  1, 'openoutput'); io_i = io_i + 1;
 
   G = linearize(mdl, io);
 #+end_src
@@ -1028,11 +284,11 @@ This is also the case for the FEM model (Figure [[fig:apa300ml_resonance]]).
 [[file:figs/apa300ml_resonance.png]]
 
 *** Amplification factor
-The amplification factor is the ratio of the axial displacement to the stack displacement.
+The amplification factor is the ratio of the vertical displacement to the stack displacement.
 
 #+begin_src matlab :exports none
   %% Name of the Simulink File
-  mdl = 'APA300ML_test_bench';
+  mdl = 'APA300ML';
 
   %% Input/Output definition
   clear io; io_i = 1;
@@ -1051,7 +307,7 @@ The ratio of the two displacement is computed from the FEM model.
 #+RESULTS:
 : 5.0749
 
-If we take the ratio of the piezo height and length (approximation of the amplification factor):
+This is actually correct and approximately corresponds to the ratio of the piezo height and length:
 #+begin_src matlab :results replace value
   75/15
 #+end_src
@@ -1078,15 +334,8 @@ with:
 
 This is exactly the specified stroke in the data-sheet.
 
-** Identification of the Dynamics
-The flexible element is imported using the =Reduced Order Flexible Solid= simscape block.
-
-To model the actuator, an =Internal Force= block is added between the nodes 3 and 12.
-A =Relative Motion Sensor= block is added between the nodes 1 and 2 to measure the displacement and the amplified piezo.
-
-One mass is fixed at one end of the piezo-electric stack actuator, the other end is fixed to the world frame.
-
-We first set the mass to be zero.
+** Identification of the Dynamics from actuator to replace displacement
+We first set the mass to be approximately zero.
 #+begin_src matlab :exports none
   m = 0.01;
 #+end_src
@@ -1094,7 +343,7 @@ We first set the mass to be zero.
 The dynamics is identified from the applied force to the measured relative displacement.
 #+begin_src matlab :exports none
   %% Name of the Simulink File
-  mdl = 'APA300ML_test_bench';
+  mdl = 'APA300ML';
 
   %% Input/Output definition
   clear io; io_i = 1;
@@ -1111,7 +360,7 @@ The same dynamics is identified for a payload mass of 10Kg.
 
 #+begin_src matlab :exports none
   %% Name of the Simulink File
-  mdl = 'APA300ML_test_bench';
+  mdl = 'APA300ML';
 
   %% Input/Output definition
   clear io; io_i = 1;
@@ -1125,8 +374,9 @@ The same dynamics is identified for a payload mass of 10Kg.
   freqs = logspace(0, 4, 5000);
 
   figure;
+  tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
 
-  ax1 = subplot(2,1,1);
+  ax1 = nexttile([2,1]);
   hold on;
   plot(freqs, abs(squeeze(freqresp(Gh, freqs, 'Hz'))), '-');
   plot(freqs, abs(squeeze(freqresp(Ghm, freqs, 'Hz'))), '-');
@@ -1135,7 +385,7 @@ The same dynamics is identified for a payload mass of 10Kg.
   ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
   hold off;
 
-  ax2 = subplot(2,1,2);
+  ax2 = nexttile;
   hold on;
   plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gh, freqs, 'Hz')))), '-', ...
        'DisplayName', '$m = 0kg$');
@@ -1152,7 +402,7 @@ The same dynamics is identified for a payload mass of 10Kg.
 #+end_src
 
 #+begin_src matlab :tangle no :exports results :results file replace
-  exportFig('figs/apa300ml_plant_dynamics.pdf', 'width', 'full', 'height', 'full');
+  exportFig('figs/apa300ml_plant_dynamics.pdf', 'width', 'wide', 'height', 'tall');
 #+end_src
 
 #+name: fig:apa300ml_plant_dynamics
@@ -1160,17 +410,46 @@ The same dynamics is identified for a payload mass of 10Kg.
 #+RESULTS:
 [[file:figs/apa300ml_plant_dynamics.png]]
 
-** IFF
-Let's use 2 stacks as actuators and 1 stack as force sensor.
+The root locus corresponding to Direct Velocity Feedback with a mass of 10kg is shown in Figure [[fig:apa300ml_dvf_root_locus]].
+#+begin_src matlab :exports none
+  figure;
 
-The transfer function from actuator to sensors is identified and shown in Figure [[fig:apa300ml_iff_plant]].
+  gains = logspace(0, 5, 500);
+
+  hold on;
+  plot(real(pole(Ghm)),  imag(pole(G)), 'kx');
+  plot(real(tzero(Ghm)),  imag(tzero(G)), 'ko');
+  for k = 1:length(gains)
+      cl_poles = pole(feedback(Ghm, gains(k)*s));
+      plot(real(cl_poles), imag(cl_poles), 'k.');
+  end
+  hold off;
+  axis square;
+  xlim([-500, 10]); ylim([0, 510]);
+
+  xlabel('Real Part'); ylabel('Imaginary Part');
+#+end_src
+
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/apa300ml_dvf_root_locus.pdf', 'width', 'wide', 'height', 'tall');
+#+end_src
+
+#+name: fig:apa300ml_dvf_root_locus
+#+caption: Root Locus for Direct Velocity Feedback
+#+RESULTS:
+[[file:figs/apa300ml_dvf_root_locus.png]]
+
+** Identification of the Dynamics from actuator to force sensor
+Let's use 2 stacks as a force sensor and 1 stack as force actuator.
+
+The transfer function from actuator voltage to sensor voltage is identified and shown in Figure [[fig:apa300ml_iff_plant]].
 #+begin_src matlab :exports none
   m = 10;
 #+end_src
 
 #+begin_src matlab :exports none
   %% Name of the Simulink File
-  mdl = 'APA300ML_test_bench';
+  mdl = 'APA300ML';
 
   %% Input/Output definition
   clear io; io_i = 1;
@@ -1184,8 +463,9 @@ The transfer function from actuator to sensors is identified and shown in Figure
   freqs = logspace(0, 4, 5000);
 
   figure;
+  tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
 
-  ax1 = subplot(2,1,1);
+  ax1 = nexttile([2,1]);
   hold on;
   plot(freqs, abs(squeeze(freqresp(Giff, freqs, 'Hz'))), '-');
   hold off;
@@ -1193,7 +473,7 @@ The transfer function from actuator to sensors is identified and shown in Figure
   ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
   hold off;
 
-  ax2 = subplot(2,1,2);
+  ax2 = nexttile;
   hold on;
   plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Giff, freqs, 'Hz')))), '-');
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
@@ -1206,7 +486,7 @@ The transfer function from actuator to sensors is identified and shown in Figure
 #+end_src
 
 #+begin_src matlab :tangle no :exports results :results file replace
-  exportFig('figs/apa300ml_iff_plant.pdf', 'width', 'full', 'height', 'full');
+  exportFig('figs/apa300ml_iff_plant.pdf', 'width', 'wide', 'height', 'tall');
 #+end_src
 
 #+name: fig:apa300ml_iff_plant
@@ -1243,271 +523,6 @@ For root locus corresponding to IFF is shown in Figure [[fig:apa300ml_iff_root_l
 #+RESULTS:
 [[file:figs/apa300ml_iff_root_locus.png]]
 
-** DVF
-Now the dynamics from the stack actuator to the relative motion sensor is identified and shown in Figure [[fig:apa300ml_dvf_plant]].
-#+begin_src matlab :exports none
-  m = 10;
-#+end_src
-
-#+begin_src matlab :exports none
-  %% Name of the Simulink File
-  mdl = 'APA300ML_test_bench';
-
-  %% Input/Output definition
-  clear io; io_i = 1;
-  io(io_i) = linio([mdl, '/F'], 1, 'openinput');  io_i = io_i + 1;
-  io(io_i) = linio([mdl, '/y'], 1, 'openoutput'); io_i = io_i + 1;
-
-  G = -linearize(mdl, io);
-#+end_src
-
-#+begin_src matlab :exports none
-  freqs = logspace(0, 4, 5000);
-
-  figure;
-
-  ax1 = subplot(2,1,1);
-  hold on;
-  plot(freqs, abs(squeeze(freqresp(G, freqs, 'Hz'))), '-');
-  hold off;
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
-  hold off;
-
-  ax2 = subplot(2,1,2);
-  hold on;
-  plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(G, freqs, 'Hz')))), '-');
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
-  yticks(-360:90:360);
-  ylim([-360 0]);
-  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
-  hold off;
-  linkaxes([ax1,ax2],'x');
-  xlim([freqs(1), freqs(end)]);
-#+end_src
-
-#+begin_src matlab :tangle no :exports results :results file replace
-  exportFig('figs/apa300ml_dvf_plant.pdf', 'width', 'full', 'height', 'full');
-#+end_src
-
-#+name: fig:apa300ml_dvf_plant
-#+caption: Transfer function from stack actuator to relative motion sensor
-#+RESULTS:
-[[file:figs/apa300ml_dvf_plant.png]]
-
-The root locus for DVF is shown in Figure [[fig:apa300ml_dvf_root_locus]].
-
-#+begin_src matlab :exports none
-  figure;
-
-  gains = logspace(0, 5, 500);
-
-  hold on;
-  plot(real(pole(G)),  imag(pole(G)), 'kx');
-  plot(real(tzero(G)),  imag(tzero(G)), 'ko');
-  for k = 1:length(gains)
-      cl_poles = pole(feedback(G, gains(k)*s));
-      plot(real(cl_poles), imag(cl_poles), 'k.');
-  end
-  hold off;
-  axis square;
-  xlim([-500, 10]); ylim([0, 510]);
-
-  xlabel('Real Part'); ylabel('Imaginary Part');
-#+end_src
-
-#+begin_src matlab :tangle no :exports results :results file replace
-  exportFig('figs/apa300ml_dvf_root_locus.pdf', 'width', 'wide', 'height', 'tall');
-#+end_src
-
-#+name: fig:apa300ml_dvf_root_locus
-#+caption: Root Locus for Direct Velocity Feedback
-#+RESULTS:
-[[file:figs/apa300ml_dvf_root_locus.png]]
-
-** TODO Sensor Fusion                                              :noexport:
-- [ ] What is the goal of that? Special control properties, lower the sensor noise?
-
-Use the relative motion sensor at low frequency and the force sensor at high frequency.
-#+begin_src matlab
-  m = 10;
-#+end_src
-
-#+begin_src matlab
-  %% Name of the Simulink File
-  mdl = 'APA300ML_test_bench';
-
-  %% Input/Output definition
-  clear io; io_i = 1;
-  io(io_i) = linio([mdl, '/Va'], 1, 'openinput');  io_i = io_i + 1;
-  io(io_i) = linio([mdl, '/Vs'], 1, 'openoutput'); io_i = io_i + 1;
-  io(io_i) = linio([mdl, '/y'],  1, 'openoutput'); io_i = io_i + 1;
-
-  G = linearize(mdl, io);
-#+end_src
-
-#+begin_src matlab :exports none
-  freqs = logspace(1, 5, 1000);
-
-  Gresp = squeeze(freqresp(G, freqs, 'Hz'));
-
-  figure;
-
-  ax1 = subplot(2,1,1);
-  hold on;
-  plot(freqs, abs(Gresp(1, :)), '-');
-  plot(freqs, abs(Gresp(2, :)), '-');
-  hold off;
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
-  hold off;
-
-  ax2 = subplot(2,1,2);
-  hold on;
-  plot(freqs, 180/pi*unwrap(angle(Gresp(1, :))), '-', ...
-       'DisplayName', 'Force');
-  plot(freqs, 180/pi*unwrap(angle(Gresp(2, :))), '-', ...
-       'DisplayName', 'Displacement');
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
-  yticks(-360:90:360);
-  ylim([-390 30]);
-  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
-  hold off;
-  linkaxes([ax1,ax2],'x');
-  xlim([freqs(1), freqs(end)]);
-  legend('location', 'southwest');
-#+end_src
-
-Merge around 1kHz
-#+begin_src matlab :exports none
-  freqs = logspace(1, 5, 1000);
-
-  figure;
-
-  ax1 = subplot(2,1,1);
-  hold on;
-  plot(freqs, abs(Gresp(1, :).*squeeze(freqresp(0.0065/s^2, freqs, 'Hz'))'), '-');
-  plot(freqs, abs(Gresp(2, :)), '-');
-  hold off;
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
-  hold off;
-
-  ax2 = subplot(2,1,2);
-  hold on;
-  plot(freqs, 180/pi*unwrap(angle(Gresp(1, :).*squeeze(freqresp(1/s^2, freqs, 'Hz'))')), '-', ...
-       'DisplayName', 'Force');
-  plot(freqs, 180/pi*unwrap(angle(Gresp(2, :))), '-', ...
-       'DisplayName', 'Displacement');
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
-  yticks(-360:90:360);
-  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
-  hold off;
-  linkaxes([ax1,ax2],'x');
-  xlim([freqs(1), freqs(end)]);
-  legend('location', 'southwest');
-#+end_src
-
-LPF and HPF
-#+begin_src matlab
-  n = 3; w0 = 2*pi*1e3; G0 = 1/10; G1 = 1e5; Gc = 1/2;
-  W1 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n;
-
-  n = 3; w0 = 2*pi*1e3; G0 = 1e5; G1 = 0.1; Gc = 1/2;
-  W2 = (((1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (G0/Gc)^(1/n))/((1/G1)^(1/n)*(1/w0)*sqrt((1-(G0/Gc)^(2/n))/(1-(Gc/G1)^(2/n)))*s + (1/Gc)^(1/n)))^n;
-#+end_src
-
-#+begin_src matlab
-  P = [W1 -W1;
-       0   W2;
-       1   0];
-      
-  [H2, ~, gamma, ~] = hinfsyn(P, 1, 1,'TOLGAM', 0.001, 'METHOD', 'ric', 'DISPLAY', 'on');
-  H1 = 1 - H2;
-#+end_src
-
-#+begin_src matlab :exports none
-  freqs = logspace(1, 5, 1000);
-
-  Gresp = squeeze(freqresp(G, freqs, 'Hz'));
-
-  figure;
-
-  ax1 = subplot(2,1,1);
-  hold on;
-  plot(freqs, abs(Gresp(1, :).*squeeze(freqresp(0.0065/s^2*H2, freqs, 'Hz'))'), '-');
-  plot(freqs, abs(Gresp(2, :).*squeeze(freqresp(H1, freqs, 'Hz'))'), '-');
-  plot(freqs, abs(Gresp(2, :).*squeeze(freqresp(H1, freqs, 'Hz'))' + Gresp(1, :).*squeeze(freqresp(0.0065/s^2*H2, freqs, 'Hz'))'), 'k--');
-  hold off;
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
-  hold off;
-
-  ax2 = subplot(2,1,2);
-  hold on;
-  plot(freqs, 180/pi*unwrap(angle(Gresp(1, :).*squeeze(freqresp(0.0065/s^2*H2, freqs, 'Hz'))')), '-', ...
-       'DisplayName', 'Force');
-  plot(freqs, 180/pi*unwrap(angle(Gresp(2, :).*squeeze(freqresp(H1, freqs, 'Hz'))')), '-', ...
-       'DisplayName', 'Displacement');
-  plot(freqs, 180/pi*unwrap(angle(Gresp(2, :).*squeeze(freqresp(H1, freqs, 'Hz'))' + Gresp(1, :).*squeeze(freqresp(0.0065/s^2*H2, freqs, 'Hz'))')), 'k--', ...
-       'DisplayName', 'Super Sensor');
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
-  yticks(-360:90:360);
-  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
-  hold off;
-  linkaxes([ax1,ax2],'x');
-  xlim([freqs(1), freqs(end)]);
-  legend('location', 'southwest');
-#+end_src
-
-#+begin_src matlab :exports none
-  Gss = -zpk(G(1)*0.0065/s^2*H2 + G(2)*H1);
-
-  freqs = logspace(1, 5, 1000);
-
-  figure;
-
-  ax1 = subplot(2,1,1);
-  hold on;
-  plot(freqs, abs(squeeze(freqresp(Gss, freqs, 'Hz'))'), '-');
-  hold off;
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
-  hold off;
-
-  ax2 = subplot(2,1,2);
-  hold on;
-  plot(freqs, 180/pi*(angle(squeeze(freqresp(Gss, freqs, 'Hz'))')), '-', ...
-       'DisplayName', 'Force');
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
-  yticks(-360:90:360);
-  xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
-  hold off;
-  linkaxes([ax1,ax2],'x');
-  xlim([freqs(1), freqs(end)]);
-  legend('location', 'southwest');
-#+end_src
-
-Root locus
-#+begin_src matlab
-  figure;
-
-  gains = logspace(4, 7, 100);
-
-  hold on;
-  plot(real(pole(Gss)),  imag(pole(Gss)), 'kx');
-  plot(real(tzero(Gss)),  imag(tzero(Gss)), 'ko');
-  for k = 1:length(gains)
-      cl_poles = pole(feedback(Gss, gains(k)*s));
-      plot(real(cl_poles), imag(cl_poles), 'k.');
-  end
-  hold off;
-  axis square;
-  xlim([-500, 10]); ylim([0, 510]);
-
-  xlabel('Real Part'); ylabel('Imaginary Part');
-#+end_src
-
 ** Identification for a simpler model
 The goal in this section is to identify the parameters of a simple APA model from the FEM.
 This can be useful is a lower order model is to be used for simulations.
@@ -1521,7 +536,7 @@ The parameters are shown in the table below.
 #+caption: Picture of an APA100M from Cedrat Technologies. Simplified model of a one DoF payload mounted on such isolator
 [[file:./figs/souleille18_model_piezo.png]]
 
-#+caption: Parameters used for the model of the APA 100M
+#+caption:Parameters used for the model of the APA 100M
 |       | Meaning                                                        |
 |-------+----------------------------------------------------------------|
 | $k_e$ | Stiffness used to adjust the pole of the isolator              |
@@ -1544,7 +559,7 @@ If we can fix $k_a$, we can determine $k_e$ and $k_1$ with:
 
 #+begin_src matlab :exports none
   %% Name of the Simulink File
-  mdl = 'APA300ML_test_bench';
+  mdl = 'APA300ML';
 
   %% Input/Output definition
   clear io; io_i = 1;
@@ -1595,7 +610,7 @@ The damping in the system is adjusted to match the FEM model if necessary.
   c1 = 1e2;
 #+end_src
 
-Analytical model of the simpler system:
+The analytical model of the simpler system is defined below:
 #+begin_src matlab
   Ga = 1/(m*s^2 + k1 + c1*s + ke*ka/(ke + ka)) * ...
        [ 1 ,              k1 + c1*s + ke*ka/(ke + ka)  , ke/(ke + ka) ;
@@ -1605,69 +620,78 @@ Analytical model of the simpler system:
   Ga.OutputName = {'y', 'Fs'};
 #+end_src
 
-Adjust the DC gain for the force sensor:
+And the DC gain is adjusted for the force sensor:
 #+begin_src matlab
   F_gain = dcgain(G('Fs', 'Fd'))/dcgain(Ga('Fs', 'Fd'));
 #+end_src
 
+The dynamics of the FEM model and the simpler model are compared in Figure [[fig:apa300ml_comp_simpler_model]].
+
 #+begin_src matlab :exports none
   freqs = logspace(0, 5, 1000);
 
   figure;
+  tiledlayout(2, 3, 'TileSpacing', 'None', 'Padding', 'None');
 
-  ax1 = subplot(2, 3, 1);
-  title('$\displaystyle \frac{x_1}{w}$')
+  ax1 = nexttile;
   hold on;
   plot(freqs, abs(squeeze(freqresp(G( 'y', 'w'), freqs, 'Hz'))));
   plot(freqs, abs(squeeze(freqresp(Ga('y', 'w'), freqs, 'Hz'))));
   hold off;
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude [m/m]');xlabel('Frequency [Hz]');
+  set(gca, 'XTickLabel',[]);
+  ylabel('$x_1/w$ [m/m]');
+  ylim([1e-6, 1e2]);
 
-  ax2 = subplot(2, 3, 2);
-  title('$\displaystyle \frac{x_1}{f}$')
+  ax2 = nexttile;
   hold on;
   plot(freqs, abs(squeeze(freqresp(G( 'y', 'Fa'), freqs, 'Hz'))));
   plot(freqs, abs(squeeze(freqresp(Ga('y', 'Fa'), freqs, 'Hz'))));
   hold off;
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude [m/N]');xlabel('Frequency [Hz]');
+  set(gca, 'XTickLabel',[]);
+  ylabel('$x_1/f$ [m/N]');
+  ylim([1e-14, 1e-6]);
 
-  ax3 = subplot(2, 3, 3);
-  title('$\displaystyle \frac{x_1}{F}$')
+  ax3 = nexttile;
   hold on;
   plot(freqs, abs(squeeze(freqresp(G( 'y', 'Fd'), freqs, 'Hz'))));
   plot(freqs, abs(squeeze(freqresp(Ga('y', 'Fd'), freqs, 'Hz'))));
   hold off;
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude [m/N]');xlabel('Frequency [Hz]');
+  set(gca, 'XTickLabel',[]);
+  ylabel('$x_1/F$ [m/N]');
+  ylim([1e-14, 1e-4]);
 
-  ax4 = subplot(2, 3, 4);
-  title('$\displaystyle \frac{F_s}{w}$')
+  ax4 = nexttile;
   hold on;
   plot(freqs, abs(squeeze(freqresp(G( 'Fs', 'w'), freqs, 'Hz'))));
   plot(freqs, abs(squeeze(freqresp(F_gain*Ga('Fs', 'w'), freqs, 'Hz'))));
   hold off;
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude [m/m]');xlabel('Frequency [Hz]');
+  xlabel('Frequency [Hz]');
+  ylabel('$F_s/w$ [m/m]');
+  ylim([1e2, 1e8]);
 
-  ax5 = subplot(2, 3, 5);
-  title('$\displaystyle \frac{F_s}{f}$')
+  ax5 = nexttile;
   hold on;
   plot(freqs, abs(squeeze(freqresp(G( 'Fs', 'Fa'), freqs, 'Hz'))));
   plot(freqs, abs(squeeze(freqresp(F_gain*Ga('Fs', 'Fa'), freqs, 'Hz'))));
   hold off;
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude [m/N]');xlabel('Frequency [Hz]');
+  xlabel('Frequency [Hz]');
+  ylabel('$F_s/f$ [m/N]');
+  ylim([1e-4, 1e1]);
 
-  ax6 = subplot(2, 3, 6);
-  title('$\displaystyle \frac{F_s}{F}$')
+  ax6 = nexttile;
   hold on;
   plot(freqs, abs(squeeze(freqresp(G( 'Fs', 'Fd'), freqs, 'Hz'))));
   plot(freqs, abs(squeeze(freqresp(F_gain*Ga('Fs', 'Fd'), freqs, 'Hz'))));
   hold off;
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
+  xlabel('Frequency [Hz]');
+  ylabel('$F_s/F$ [m/N]');
+  ylim([1e-7, 1e2]);
 
   linkaxes([ax1,ax2,ax3,ax4,ax5,ax6],'x');
 #+end_src
@@ -1681,10 +705,12 @@ Adjust the DC gain for the force sensor:
 #+RESULTS:
 [[file:figs/apa300ml_comp_simpler_model.png]]
 
-We now compare the FEM model with the simplified simscape model.
+The simplified model has also been implemented in Simscape.
+
+The dynamics of the Simscape simplified model is identified and compared with the FEM one in Figure [[fig:apa300ml_comp_simpler_simscape]].
 #+begin_src matlab :exports none
   %% Name of the Simulink File
-  mdl = 'APA300ML_test_bench_simplified';
+  mdl = 'APA300ML_simplified';
 
   %% Input/Output definition
   clear io; io_i = 1;
@@ -1704,60 +730,67 @@ We now compare the FEM model with the simplified simscape model.
   freqs = logspace(0, 5, 1000);
 
   figure;
+  tiledlayout(2, 3, 'TileSpacing', 'None', 'Padding', 'None');
 
-  ax1 = subplot(2, 3, 1);
-  title('$\displaystyle \frac{x_1}{w}$')
+  ax1 = nexttile;
   hold on;
   plot(freqs, abs(squeeze(freqresp(G( 'y', 'w'), freqs, 'Hz'))));
   plot(freqs, abs(squeeze(freqresp(Gs('y', 'w'), freqs, 'Hz'))));
   hold off;
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude [m/m]');xlabel('Frequency [Hz]');
+  set(gca, 'XTickLabel',[]);
+  ylabel('$x_1/w$ [m/m]');
+  ylim([1e-6, 1e2]);
 
-  ax2 = subplot(2, 3, 2);
-  title('$\displaystyle \frac{x_1}{f}$')
+  ax2 = nexttile;
   hold on;
   plot(freqs, abs(squeeze(freqresp(G( 'y', 'Fa'), freqs, 'Hz'))));
   plot(freqs, abs(squeeze(freqresp(Gs('y', 'Fa'), freqs, 'Hz'))));
   hold off;
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude [m/N]');xlabel('Frequency [Hz]');
+  set(gca, 'XTickLabel',[]);
+  ylabel('$x_1/f$ [m/N]');
+  ylim([1e-14, 1e-6]);
 
-  ax3 = subplot(2, 3, 3);
-  title('$\displaystyle \frac{x_1}{F}$')
+  ax3 = nexttile;
   hold on;
   plot(freqs, abs(squeeze(freqresp(G( 'y', 'Fd'), freqs, 'Hz'))));
   plot(freqs, abs(squeeze(freqresp(Gs('y', 'Fd'), freqs, 'Hz'))));
   hold off;
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude [m/N]');xlabel('Frequency [Hz]');
+  set(gca, 'XTickLabel',[]);
+  ylabel('$x_1/F$ [m/N]');
+  ylim([1e-14, 1e-4]);
 
-  ax4 = subplot(2, 3, 4);
-  title('$\displaystyle \frac{F_s}{w}$')
+  ax4 = nexttile;
   hold on;
   plot(freqs, abs(squeeze(freqresp(G( 'Fs', 'w'), freqs, 'Hz'))));
-  plot(freqs, abs(squeeze(freqresp(Gs('Fs', 'w'), freqs, 'Hz'))));
+  plot(freqs, abs(squeeze(freqresp(F_gain*Gs('Fs', 'w'), freqs, 'Hz'))));
   hold off;
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude [m/m]');xlabel('Frequency [Hz]');
+  xlabel('Frequency [Hz]');
+  ylabel('$F_s/w$ [m/m]');
+  ylim([1e2, 1e8]);
 
-  ax5 = subplot(2, 3, 5);
-  title('$\displaystyle \frac{F_s}{f}$')
+  ax5 = nexttile;
   hold on;
   plot(freqs, abs(squeeze(freqresp(G( 'Fs', 'Fa'), freqs, 'Hz'))));
-  plot(freqs, abs(squeeze(freqresp(Gs('Fs', 'Fa'), freqs, 'Hz'))));
+  plot(freqs, abs(squeeze(freqresp(F_gain*Gs('Fs', 'Fa'), freqs, 'Hz'))));
   hold off;
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude [m/N]');xlabel('Frequency [Hz]');
+  xlabel('Frequency [Hz]');
+  ylabel('$F_s/f$ [m/N]');
+  ylim([1e-4, 1e1]);
 
-  ax6 = subplot(2, 3, 6);
-  title('$\displaystyle \frac{F_s}{F}$')
+  ax6 = nexttile;
   hold on;
   plot(freqs, abs(squeeze(freqresp(G( 'Fs', 'Fd'), freqs, 'Hz'))));
-  plot(freqs, abs(squeeze(freqresp(Gs('Fs', 'Fd'), freqs, 'Hz'))));
+  plot(freqs, abs(squeeze(freqresp(F_gain*Gs('Fs', 'Fd'), freqs, 'Hz'))));
   hold off;
   set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  ylabel('Amplitude [m/N]'); xlabel('Frequency [Hz]');
+  xlabel('Frequency [Hz]');
+  ylabel('$F_s/F$ [m/N]');
+  ylim([1e-7, 1e2]);
 
   linkaxes([ax1,ax2,ax3,ax4,ax5,ax6],'x');
 #+end_src
@@ -1771,142 +804,8 @@ We now compare the FEM model with the simplified simscape model.
 #+RESULTS:
 [[file:figs/apa300ml_comp_simpler_simscape.png]]
 
-We save the parameters of the simplified model for the APA300ML:
-#+begin_src matlab
-  save('./mat/APA300ML_simplified_model.mat', 'ka', 'ke', 'k1', 'c1', 'F_gain');
-#+end_src
-
-** Identification of the stiffness properties
-*** APA Alone
-#+begin_src matlab :exports none
-  m = 10;
-#+end_src
-
-#+begin_src matlab :exports none
-  %% Name of the Simulink File
-  mdl = 'APA300ML_characterisation';
-
-  %% Input/Output definition
-  clear io; io_i = 1;
-  io(io_i) = linio([mdl, '/Fd'], 1, 'openinput');  io_i = io_i + 1; % External Vertical Force [N]
-  io(io_i) = linio([mdl, '/D'],  1, 'openoutput'); io_i = io_i + 1; % Displacement/Rotation [m]
-
-  G = linearize(mdl, io);
-#+end_src
-
-#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
-  data2orgtable([1e-6./dcgain(G(1,1)), 1e-6./dcgain(G(2,2)), 1e-6./dcgain(G(3,3)), 1./dcgain(G(4,4)), 1./dcgain(G(5,5)), 1./dcgain(G(6,6))]', {'Kx [N/um]', 'Ky [N/um]', 'Kz [N/um]', 'Rx [Nm/rad]', 'Ry [Nm/rad]', 'Rz [Nm/rad]'}, {'*Caracteristics*', '*Value*'}, ' %.1f ');
-#+end_src
-
-#+RESULTS:
-| *Caracteristics* | *Value* |
-|------------------+---------|
-| Kx [N/um]        |     0.8 |
-| Ky [N/um]        |     1.6 |
-| Kz [N/um]        |     1.8 |
-| Rx [Nm/rad]      |    71.4 |
-| Ry [Nm/rad]      |   148.2 |
-| Rz [Nm/rad]      |  4241.8 |
-
-*** See how the global stiffness is changing with the flexible joints
-#+begin_src matlab
-  flex = load('./mat/flexor_ID16.mat', 'int_xyz', 'int_i', 'n_xyz', 'n_i', 'nodes', 'M', 'K');
-#+end_src
-
-#+begin_src matlab :exports none
-  %% Name of the Simulink File
-  mdl = 'APA300ML_characterisation_with_joints';
-
-  %% Input/Output definition
-  clear io; io_i = 1;
-  io(io_i) = linio([mdl, '/Fd'], 1, 'openinput');  io_i = io_i + 1; % External Vertical Force [N]
-  io(io_i) = linio([mdl, '/D'],  1, 'openoutput'); io_i = io_i + 1; % Displacement/Rotation [m]
-
-  Gf = linearize(mdl, io);
-#+end_src
-
-#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
-  data2orgtable([1e-6./dcgain(Gf(1,1)), 1e-6./dcgain(Gf(2,2)), 1e-6./dcgain(Gf(3,3)), 1./dcgain(Gf(4,4)), 1./dcgain(Gf(5,5)), 1./dcgain(Gf(6,6))]', {'Kx [N/um]', 'Ky [N/um]', 'Kz [N/um]', 'Rx [Nm/rad]', 'Ry [Nm/rad]', 'Rz [Nm/rad]'}, {'*Caracteristic*', '*Value*'}, ' %.1f ');
-#+end_src
-
-#+RESULTS:
-| *Caracteristic* | *Value* |
-|-----------------+---------|
-| Kx [N/um]       |     0.0 |
-| Ky [N/um]       |     0.0 |
-| Kz [N/um]       |     1.8 |
-| Rx [Nm/rad]     |   722.9 |
-| Ry [Nm/rad]     |   129.6 |
-| Rz [Nm/rad]     |   115.3 |
-
-
-#+begin_src matlab
-  freqs = logspace(-2, 5, 1000);
-
-  figure;
-  hold on;
-  plot(freqs, abs(squeeze(freqresp(G(2,2), freqs, 'Hz'))), '-', 'DisplayName', 'APA');
-  plot(freqs, abs(squeeze(freqresp(Gf(2,2), freqs, 'Hz'))), '-', 'DisplayName', 'Flex');
-  hold off;
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  xlabel('Frequency [Hz]'); ylabel('Amplitude [m/N]');
-  hold off;
-  legend('location', 'northeast');
-#+end_src
-
-#+begin_src matlab
-  freqs = logspace(-2, 5, 1000);
-
-  figure;
-  hold on;
-  plot(freqs, abs(squeeze(freqresp(G(3,3), freqs, 'Hz'))), '-', 'DisplayName', 'APA');
-  plot(freqs, abs(squeeze(freqresp(Gf(3,3), freqs, 'Hz'))), '-', 'DisplayName', 'Flex');
-  hold off;
-  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
-  xlabel('Frequency [Hz]'); ylabel('Amplitude [m/N]');
-  hold off;
-  legend('location', 'northeast');
-#+end_src
-
-** Effect of APA300ML in the flexibility of the leg
-#+begin_src matlab :exports none
-  m = 10;
-#+end_src
-
-#+begin_src matlab :exports none
-  %% Name of the Simulink File
-  mdl = 'APA300ML_flex_joints';
-
-  %% Input/Output definition
-  clear io; io_i = 1;
-  io(io_i) = linio([mdl, '/Fd'], 1, 'openinput');  io_i = io_i + 1; % External Vertical Force [N]
-  io(io_i) = linio([mdl, '/D'],  1, 'openoutput'); io_i = io_i + 1; % Displacement/Rotation [m]
-#+end_src
-
-#+begin_src matlab :exports none
-  G = linearize(mdl, io);
-#+end_src
-
-#+begin_src matlab :exports none
-  Gf = linearize(mdl, io);
-#+end_src
-
-#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
-  data2orgtable([[1e-6./dcgain(G(1,1)), 1e-6./dcgain(G(2,2)), 1e-6./dcgain(G(3,3)), 1./dcgain(G(4,4)), 1./dcgain(G(5,5)), 1./dcgain(G(6,6))]', [1e-6./dcgain(Gf(1,1)), 1e-6./dcgain(Gf(2,2)), 1e-6./dcgain(Gf(3,3)), 1./dcgain(Gf(4,4)), 1./dcgain(Gf(5,5)), 1./dcgain(Gf(6,6))]'], {'Kx [N/um]', 'Ky [N/um]', 'Kz [N/um]', 'Rx [Nm/rad]', 'Ry [Nm/rad]', 'Rz [Nm/rad]'}, {'*Caracteristic*', '*Rigid APA*', '*Flexible APA*'}, ' %.3f ');
-#+end_src
-
-
-#+RESULTS:
-| *Caracteristic* | *Rigid APA* | *Flexible APA* |
-|-----------------+-------------+----------------|
-| Kx [N/um]       |       0.018 |          0.019 |
-| Ky [N/um]       |       0.018 |          0.018 |
-| Kz [N/um]       |        60.0 |          2.647 |
-| Rx [Nm/rad]     |      16.705 |        557.682 |
-| Ry [Nm/rad]     |      16.535 |        185.939 |
-| Rz [Nm/rad]     |       118.0 |        114.803 |
-
-* Flexible Joint
+* First Flexible Joint Geometry
+<<sec:first_flexible_joint>>
 ** Introduction                                                      :ignore:
 The studied flexor is shown in Figure [[fig:flexor_id16_screenshot]].
 
@@ -2066,7 +965,6 @@ And we find the same parameters as the one estimated from the Stiffness matrix.
 | Torsion Stiffness Rz [Nm/rad] |     236 |              238 |
 
 ** Simpler Model
-
 Let's now model the flexible joint with a "perfect" Bushing joint as shown in Figure [[fig:flexible_joint_simscape]].
 
 #+name: fig:flexible_joint_simscape
@@ -2153,8 +1051,10 @@ The two obtained dynamics are compared in Figure
 #+RESULTS:
 [[file:figs/flexor_ID16_compare_bushing_joint.png]]
 
-* Optimal Flexible Joint
+* Optimized Flexible Joint
+<<sec:optimized_flexible_joint>>
 ** Introduction                                                      :ignore:
+
 #+name: fig:optimal_flexor
 #+caption: Flexor studied
 [[file:data/flexor_circ_025/CS.jpg]]
@@ -2393,7 +1293,8 @@ The two obtained dynamics are compared in Figure
 #+RESULTS:
 [[file:figs/flexor_ID16_compare_bushing_joint.png]]
 
-* Integral Force Feedback with Amplified Piezo
+* Force Sensor and Integral Force Feedback with an Amplified Piezoelectric Actuator
+<<sec:integral_force_feedback>>
 ** Introduction                                                      :ignore:
 In this section, we try to replicate the results obtained in cite:souleille18_concep_activ_mount_space_applic.
 
@@ -2654,6 +1555,7 @@ The controller is defined and the loop gain is shown in Figure [[fig:piezo_ampli
 [[file:figs/souleille18_results.png]]
 
 * Complete Strut with Encoder
+<<sec:strut_fem>>
 ** Introduction
 
 #+name: fig:strut_encoder_points3
@@ -2847,6 +1749,137 @@ The same dynamics is identified for a payload mass of 10Kg.
   % legend('location', 'southwest');
 #+end_src
 
+* To Order                                                          :noexport:
+** TODO Identification of the stiffness properties of the APA300ML
+*** APA Alone
+#+begin_src matlab :exports none
+  m = 10;
+#+end_src
+
+#+begin_src matlab :exports none
+  %% Name of the Simulink File
+  mdl = 'APA300ML_characterisation';
+
+  %% Input/Output definition
+  clear io; io_i = 1;
+  io(io_i) = linio([mdl, '/Fd'], 1, 'openinput');  io_i = io_i + 1; % External Vertical Force [N]
+  io(io_i) = linio([mdl, '/D'],  1, 'openoutput'); io_i = io_i + 1; % Displacement/Rotation [m]
+
+  G = linearize(mdl, io);
+#+end_src
+
+#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
+  data2orgtable([1e-6./dcgain(G(1,1)), 1e-6./dcgain(G(2,2)), 1e-6./dcgain(G(3,3)), 1./dcgain(G(4,4)), 1./dcgain(G(5,5)), 1./dcgain(G(6,6))]', {'Kx [N/um]', 'Ky [N/um]', 'Kz [N/um]', 'Rx [Nm/rad]', 'Ry [Nm/rad]', 'Rz [Nm/rad]'}, {'*Caracteristics*', '*Value*'}, ' %.1f ');
+#+end_src
+
+#+RESULTS:
+| *Caracteristics* | *Value* |
+|------------------+---------|
+| Kx [N/um]        |     0.8 |
+| Ky [N/um]        |     1.6 |
+| Kz [N/um]        |     1.8 |
+| Rx [Nm/rad]      |    71.4 |
+| Ry [Nm/rad]      |   148.2 |
+| Rz [Nm/rad]      |  4241.8 |
+
+*** See how the global stiffness is changing with the flexible joints
+#+begin_src matlab
+  flex = load('./mat/flexor_ID16.mat', 'int_xyz', 'int_i', 'n_xyz', 'n_i', 'nodes', 'M', 'K');
+#+end_src
+
+#+begin_src matlab :exports none
+  %% Name of the Simulink File
+  mdl = 'APA300ML_characterisation_with_joints';
+
+  %% Input/Output definition
+  clear io; io_i = 1;
+  io(io_i) = linio([mdl, '/Fd'], 1, 'openinput');  io_i = io_i + 1; % External Vertical Force [N]
+  io(io_i) = linio([mdl, '/D'],  1, 'openoutput'); io_i = io_i + 1; % Displacement/Rotation [m]
+
+  Gf = linearize(mdl, io);
+#+end_src
+
+#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
+  data2orgtable([1e-6./dcgain(Gf(1,1)), 1e-6./dcgain(Gf(2,2)), 1e-6./dcgain(Gf(3,3)), 1./dcgain(Gf(4,4)), 1./dcgain(Gf(5,5)), 1./dcgain(Gf(6,6))]', {'Kx [N/um]', 'Ky [N/um]', 'Kz [N/um]', 'Rx [Nm/rad]', 'Ry [Nm/rad]', 'Rz [Nm/rad]'}, {'*Caracteristic*', '*Value*'}, ' %.1f ');
+#+end_src
+
+#+RESULTS:
+| *Caracteristic* | *Value* |
+|-----------------+---------|
+| Kx [N/um]       |     0.0 |
+| Ky [N/um]       |     0.0 |
+| Kz [N/um]       |     1.8 |
+| Rx [Nm/rad]     |   722.9 |
+| Ry [Nm/rad]     |   129.6 |
+| Rz [Nm/rad]     |   115.3 |
+
+
+#+begin_src matlab
+  freqs = logspace(-2, 5, 1000);
+
+  figure;
+  hold on;
+  plot(freqs, abs(squeeze(freqresp(G(2,2), freqs, 'Hz'))), '-', 'DisplayName', 'APA');
+  plot(freqs, abs(squeeze(freqresp(Gf(2,2), freqs, 'Hz'))), '-', 'DisplayName', 'Flex');
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('Amplitude [m/N]');
+  hold off;
+  legend('location', 'northeast');
+#+end_src
+
+#+begin_src matlab
+  freqs = logspace(-2, 5, 1000);
+
+  figure;
+  hold on;
+  plot(freqs, abs(squeeze(freqresp(G(3,3), freqs, 'Hz'))), '-', 'DisplayName', 'APA');
+  plot(freqs, abs(squeeze(freqresp(Gf(3,3), freqs, 'Hz'))), '-', 'DisplayName', 'Flex');
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  xlabel('Frequency [Hz]'); ylabel('Amplitude [m/N]');
+  hold off;
+  legend('location', 'northeast');
+#+end_src
+
+** TODO Effect of APA300ML in the flexibility of the leg
+#+begin_src matlab :exports none
+  m = 10;
+#+end_src
+
+#+begin_src matlab :exports none
+  %% Name of the Simulink File
+  mdl = 'APA300ML_flex_joints';
+
+  %% Input/Output definition
+  clear io; io_i = 1;
+  io(io_i) = linio([mdl, '/Fd'], 1, 'openinput');  io_i = io_i + 1; % External Vertical Force [N]
+  io(io_i) = linio([mdl, '/D'],  1, 'openoutput'); io_i = io_i + 1; % Displacement/Rotation [m]
+#+end_src
+
+#+begin_src matlab :exports none
+  G = linearize(mdl, io);
+#+end_src
+
+#+begin_src matlab :exports none
+  Gf = linearize(mdl, io);
+#+end_src
+
+#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
+  data2orgtable([[1e-6./dcgain(G(1,1)), 1e-6./dcgain(G(2,2)), 1e-6./dcgain(G(3,3)), 1./dcgain(G(4,4)), 1./dcgain(G(5,5)), 1./dcgain(G(6,6))]', [1e-6./dcgain(Gf(1,1)), 1e-6./dcgain(Gf(2,2)), 1e-6./dcgain(Gf(3,3)), 1./dcgain(Gf(4,4)), 1./dcgain(Gf(5,5)), 1./dcgain(Gf(6,6))]'], {'Kx [N/um]', 'Ky [N/um]', 'Kz [N/um]', 'Rx [Nm/rad]', 'Ry [Nm/rad]', 'Rz [Nm/rad]'}, {'*Caracteristic*', '*Rigid APA*', '*Flexible APA*'}, ' %.3f ');
+#+end_src
+
+
+#+RESULTS:
+| *Caracteristic* | *Rigid APA* | *Flexible APA* |
+|-----------------+-------------+----------------|
+| Kx [N/um]       |       0.018 |          0.019 |
+| Ky [N/um]       |       0.018 |          0.018 |
+| Kz [N/um]       |        60.0 |          2.647 |
+| Rx [Nm/rad]     |      16.705 |        557.682 |
+| Ry [Nm/rad]     |      16.535 |        185.939 |
+| Rz [Nm/rad]     |       118.0 |        114.803 |
+
 * Bibliography                                                        :ignore:
 bibliographystyle:unsrt
 bibliography:ref.bib
diff --git a/mat/APA300ML.mat b/mat/APA300ML.mat
index 560064a..9bf5907 100644
Binary files a/mat/APA300ML.mat and b/mat/APA300ML.mat differ
diff --git a/matlab/APA300ML.m b/matlab/APA300ML.m
new file mode 100644
index 0000000..e8f4984
--- /dev/null
+++ b/matlab/APA300ML.m
@@ -0,0 +1,592 @@
+%% Clear Workspace and Close figures
+clear; close all; clc;
+
+%% Intialize Laplace variable
+s = zpk('s');
+
+addpath('APA300ML/');
+
+open('APA300ML.slx');
+
+% Import Mass Matrix, Stiffness Matrix, and Interface Nodes Coordinates
+% We first extract the stiffness and mass matrices.
+
+K = readmatrix('mat_K.CSV');
+M = readmatrix('mat_M.CSV');
+
+
+
+% #+caption: First 10x10 elements of the Mass matrix
+% #+RESULTS:
+% |    0.01 |  -2e-06 |  1e-06 |  6e-09 |  5e-05 | -5e-09 | -0.0005 |  -7e-07 |  6e-07 | -3e-09 |
+% |  -2e-06 |    0.01 |  8e-07 | -2e-05 | -8e-09 |  2e-09 |  -9e-07 | -0.0002 |  1e-08 | -9e-07 |
+% |   1e-06 |   8e-07 |  0.009 |  5e-10 |  1e-09 | -1e-09 |  -5e-07 |   3e-08 |  6e-05 |  1e-10 |
+% |   6e-09 |  -2e-05 |  5e-10 |  3e-07 |  2e-11 | -3e-12 |   3e-09 |   9e-07 | -4e-10 |  3e-09 |
+% |   5e-05 |  -8e-09 |  1e-09 |  2e-11 |  6e-07 | -4e-11 |  -1e-06 |  -2e-09 |  1e-09 | -8e-12 |
+% |  -5e-09 |   2e-09 | -1e-09 | -3e-12 | -4e-11 |  1e-07 |  -2e-09 |  -1e-09 | -4e-10 | -5e-12 |
+% | -0.0005 |  -9e-07 | -5e-07 |  3e-09 | -1e-06 | -2e-09 |    0.01 |   1e-07 | -3e-07 | -2e-08 |
+% |  -7e-07 | -0.0002 |  3e-08 |  9e-07 | -2e-09 | -1e-09 |   1e-07 |    0.01 | -4e-07 |  2e-05 |
+% |   6e-07 |   1e-08 |  6e-05 | -4e-10 |  1e-09 | -4e-10 |  -3e-07 |  -4e-07 |  0.009 | -2e-10 |
+% |  -3e-09 |  -9e-07 |  1e-10 |  3e-09 | -8e-12 | -5e-12 |  -2e-08 |   2e-05 | -2e-10 |  3e-07 |
+
+
+% Then, we extract the coordinates of the interface nodes.
+
+[int_xyz, int_i, n_xyz, n_i, nodes] = extractNodes('out_nodes_3D.txt');
+
+% Piezoelectric parameters
+% In order to make the conversion from applied voltage to generated force or from the strain to the generated voltage, we need to defined some parameters corresponding to the piezoelectric material:
+
+d33 = 300e-12; % Strain constant [m/V]
+n   = 80;      % Number of layers per stack
+eT  = 1.6e-8;  % Permittivity under constant stress [F/m]
+sD  = 1e-11;   % Compliance under constant electric displacement [m2/N]
+ka  = 235e6;   % Stack stiffness [N/m]
+C   = 5e-6;    % Stack capactiance [F]
+
+
+
+% The ratio of the developed force to applied voltage is:
+% #+name: eq:piezo_voltage_to_force
+% \begin{equation}
+%   F_a = g_a V_a, \quad g_a = d_{33} n k_a
+% \end{equation}
+% where:
+% - $F_a$: developed force in [N]
+% - $n$: number of layers of the actuator stack
+% - $d_{33}$: strain constant in [m/V]
+% - $k_a$: actuator stack stiffness in [N/m]
+% - $V_a$: applied voltage in [V]
+
+% If we take the numerical values, we obtain:
+
+d33*n*ka % [N/V]
+
+
+
+% #+RESULTS:
+% : 5.64
+
+% From cite:fleming14_desig_model_contr_nanop_system (page 123), the relation between relative displacement of the sensor stack and generated voltage is:
+% #+name: eq:piezo_strain_to_voltage
+% \begin{equation}
+%   V_s = \frac{d_{33}}{\epsilon^T s^D n} \Delta h
+% \end{equation}
+% where:
+% - $V_s$: measured voltage in [V]
+% - $d_{33}$: strain constant in [m/V]
+% - $\epsilon^T$: permittivity under constant stress in [F/m]
+% - $s^D$: elastic compliance under constant electric displacement in [m^2/N]
+% - $n$: number of layers of the sensor stack
+% - $\Delta h$: relative displacement in [m]
+
+% If we take the numerical values, we obtain:
+
+1e-6*d33/(eT*sD*n) % [V/um]
+
+% Stiffness
+
+m = 0.001;
+
+
+
+% The transfer function from vertical external force to the relative vertical displacement is identified.
+
+
+%% Name of the Simulink File
+mdl = 'APA300ML';
+
+%% Input/Output definition
+clear io; io_i = 1;
+io(io_i) = linio([mdl, '/Fd'], 1, 'openinput');  io_i = io_i + 1;
+io(io_i) = linio([mdl, '/z'],  1, 'openoutput'); io_i = io_i + 1;
+
+G = linearize(mdl, io);
+
+
+
+% The inverse of its DC gain is the axial stiffness of the APA:
+
+1e-6/dcgain(G) % [N/um]
+
+% Resonance Frequency
+% The resonance frequency is specified to be between 650Hz and 840Hz.
+% This is also the case for the FEM model (Figure [[fig:apa300ml_resonance]]).
+
+
+freqs = logspace(2, 4, 5000);
+
+figure;
+hold on;
+plot(freqs, abs(squeeze(freqresp(G, freqs, 'Hz'))));
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('Frequency [Hz]'); ylabel('Amplitude');
+hold off;
+
+% Amplification factor
+% The amplification factor is the ratio of the vertical displacement to the stack displacement.
+
+
+%% Name of the Simulink File
+mdl = 'APA300ML';
+
+%% Input/Output definition
+clear io; io_i = 1;
+io(io_i) = linio([mdl, '/F'], 1, 'openinput');  io_i = io_i + 1;
+io(io_i) = linio([mdl, '/z'], 1, 'openoutput'); io_i = io_i + 1;
+io(io_i) = linio([mdl, '/d'], 1, 'openoutput'); io_i = io_i + 1;
+
+G = linearize(mdl, io);
+
+
+
+% The ratio of the two displacement is computed from the FEM model.
+
+abs(dcgain(G(1,1))./dcgain(G(2,1)))
+
+
+
+% #+RESULTS:
+% : 5.0749
+
+% This is actually correct and approximately corresponds to the ratio of the piezo height and length:
+
+75/15
+
+% Stroke
+
+% Estimation of the actuator stroke:
+% \[ \Delta H = A n \Delta L \]
+% with:
+% - $\Delta H$ Axial Stroke of the APA
+% - $A$ Amplification factor (5 for the APA300ML)
+% - $n$ Number of stack used
+% - $\Delta L$ Stroke of the stack (0.1% of its length)
+
+
+1e6 * 5 * 3 * 20e-3 * 0.1e-2
+
+% Identification of the Dynamics from actuator to replace displacement
+% We first set the mass to be approximately zero.
+
+m = 0.01;
+
+
+
+% The dynamics is identified from the applied force to the measured relative displacement.
+
+%% Name of the Simulink File
+mdl = 'APA300ML';
+
+%% Input/Output definition
+clear io; io_i = 1;
+io(io_i) = linio([mdl, '/F'], 1, 'openinput');  io_i = io_i + 1;
+io(io_i) = linio([mdl, '/z'], 1, 'openoutput'); io_i = io_i + 1;
+
+Gh = -linearize(mdl, io);
+
+
+
+% The same dynamics is identified for a payload mass of 10Kg.
+
+m = 10;
+
+%% Name of the Simulink File
+mdl = 'APA300ML';
+
+%% Input/Output definition
+clear io; io_i = 1;
+io(io_i) = linio([mdl, '/F'], 1, 'openinput');  io_i = io_i + 1;
+io(io_i) = linio([mdl, '/z'], 1, 'openoutput'); io_i = io_i + 1;
+
+Ghm = -linearize(mdl, io);
+
+freqs = logspace(0, 4, 5000);
+
+figure;
+tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
+ax1 = nexttile([2,1]);
+hold on;
+plot(freqs, abs(squeeze(freqresp(Gh, freqs, 'Hz'))), '-');
+plot(freqs, abs(squeeze(freqresp(Ghm, freqs, 'Hz'))), '-');
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
+hold off;
+
+ax2 = nexttile;
+hold on;
+plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Gh, freqs, 'Hz')))), '-', ...
+     'DisplayName', '$m = 0kg$');
+plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Ghm, freqs, 'Hz')))), '-', ...
+     'DisplayName', '$m = 10kg$');
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
+yticks(-360:90:360);
+ylim([-360 0]);
+xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+hold off;
+linkaxes([ax1,ax2],'x');
+xlim([freqs(1), freqs(end)]);
+legend('location', 'southwest');
+
+
+
+% #+name: fig:apa300ml_plant_dynamics
+% #+caption: Transfer function from forces applied by the stack to the axial displacement of the APA
+% #+RESULTS:
+% [[file:figs/apa300ml_plant_dynamics.png]]
+
+% The root locus corresponding to Direct Velocity Feedback with a mass of 10kg is shown in Figure [[fig:apa300ml_dvf_root_locus]].
+
+figure;
+
+gains = logspace(0, 5, 500);
+
+hold on;
+plot(real(pole(Ghm)),  imag(pole(G)), 'kx');
+plot(real(tzero(Ghm)),  imag(tzero(G)), 'ko');
+for k = 1:length(gains)
+    cl_poles = pole(feedback(Ghm, gains(k)*s));
+    plot(real(cl_poles), imag(cl_poles), 'k.');
+end
+hold off;
+axis square;
+xlim([-500, 10]); ylim([0, 510]);
+
+xlabel('Real Part'); ylabel('Imaginary Part');
+
+% Identification of the Dynamics from actuator to force sensor
+% Let's use 2 stacks as a force sensor and 1 stack as force actuator.
+
+% The transfer function from actuator voltage to sensor voltage is identified and shown in Figure [[fig:apa300ml_iff_plant]].
+
+m = 10;
+
+%% Name of the Simulink File
+mdl = 'APA300ML';
+
+%% Input/Output definition
+clear io; io_i = 1;
+io(io_i) = linio([mdl, '/Va'], 1, 'openinput');  io_i = io_i + 1;
+io(io_i) = linio([mdl, '/Vs'], 1, 'openoutput'); io_i = io_i + 1;
+
+Giff = -linearize(mdl, io);
+
+freqs = logspace(0, 4, 5000);
+
+figure;
+tiledlayout(3, 1, 'TileSpacing', 'None', 'Padding', 'None');
+
+ax1 = nexttile([2,1]);
+hold on;
+plot(freqs, abs(squeeze(freqresp(Giff, freqs, 'Hz'))), '-');
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
+hold off;
+
+ax2 = nexttile;
+hold on;
+plot(freqs, 180/pi*unwrap(angle(squeeze(freqresp(Giff, freqs, 'Hz')))), '-');
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
+yticks(-360:90:360);
+ylim([-180 180]);
+xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
+hold off;
+linkaxes([ax1,ax2],'x');
+xlim([freqs(1), freqs(end)]);
+
+
+
+% #+name: fig:apa300ml_iff_plant
+% #+caption: Transfer function from actuator to force sensor
+% #+RESULTS:
+% [[file:figs/apa300ml_iff_plant.png]]
+
+% For root locus corresponding to IFF is shown in Figure [[fig:apa300ml_iff_root_locus]].
+
+figure;
+
+gains = logspace(0, 5, 500);
+
+hold on;
+plot(real(pole(Giff)),  imag(pole(Giff)), 'kx');
+plot(real(tzero(Giff)),  imag(tzero(Giff)), 'ko');
+for k = 1:length(gains)
+    cl_poles = pole(feedback(Giff, gains(k)/s));
+    plot(real(cl_poles), imag(cl_poles), 'k.');
+end
+hold off;
+axis square;
+xlim([-500, 10]); ylim([0, 510]);
+
+xlabel('Real Part'); ylabel('Imaginary Part');
+
+% Identification for a simpler model
+% The goal in this section is to identify the parameters of a simple APA model from the FEM.
+% This can be useful is a lower order model is to be used for simulations.
+
+% The presented model is based on cite:souleille18_concep_activ_mount_space_applic.
+
+% The model represents the Amplified Piezo Actuator (APA) from Cedrat-Technologies (Figure [[fig:souleille18_model_piezo]]).
+% The parameters are shown in the table below.
+
+% #+name: fig:souleille18_model_piezo
+% #+caption: Picture of an APA100M from Cedrat Technologies. Simplified model of a one DoF payload mounted on such isolator
+% [[file:./figs/souleille18_model_piezo.png]]
+
+% #+caption:Parameters used for the model of the APA 100M
+% |       | Meaning                                                        |
+% |-------+----------------------------------------------------------------|
+% | $k_e$ | Stiffness used to adjust the pole of the isolator              |
+% | $k_1$ | Stiffness of the metallic suspension when the stack is removed |
+% | $k_a$ | Stiffness of the actuator                                      |
+% | $c_1$ | Added viscous damping                                          |
+
+% The goal is to determine $k_e$, $k_a$ and $k_1$ so that the simplified model fits the FEM model.
+
+% \[ \alpha = \frac{x_1}{f}(\omega=0) = \frac{\frac{k_e}{k_e + k_a}}{k_1 + \frac{k_e k_a}{k_e + k_a}} \]
+% \[ \beta = \frac{x_1}{F}(\omega=0) = \frac{1}{k_1 + \frac{k_e k_a}{k_e + k_a}} \]
+
+% If we can fix $k_a$, we can determine $k_e$ and $k_1$ with:
+% \[ k_e = \frac{k_a}{\frac{\beta}{\alpha} - 1} \]
+% \[ k_1 = \frac{1}{\beta} - \frac{k_e k_a}{k_e + k_a} \]
+
+
+m = 10;
+
+%% Name of the Simulink File
+mdl = 'APA300ML';
+
+%% Input/Output definition
+clear io; io_i = 1;
+io(io_i) = linio([mdl, '/Fd'], 1, 'openinput');  io_i = io_i + 1; % External Vertical Force [N]
+io(io_i) = linio([mdl, '/w'],  1, 'openinput');  io_i = io_i + 1; % Base Motion [m]
+io(io_i) = linio([mdl, '/Fa'], 1, 'openinput');  io_i = io_i + 1; % Actuator Force [N]
+io(io_i) = linio([mdl, '/z'],  1, 'openoutput'); io_i = io_i + 1; % Vertical Displacement [m]
+io(io_i) = linio([mdl, '/Vs'], 1, 'openoutput'); io_i = io_i + 1; % Force Sensor [V]
+io(io_i) = linio([mdl, '/d'],  1, 'openoutput'); io_i = io_i + 1; % Stack Displacement [m]
+
+G = linearize(mdl, io);
+
+G.InputName = {'Fd', 'w', 'Fa'};
+G.OutputName = {'y', 'Fs', 'd'};
+
+
+
+% From the identified dynamics, compute $\alpha$ and $\beta$
+
+alpha = abs(dcgain(G('y', 'Fa')));
+beta  = abs(dcgain(G('y', 'Fd')));
+
+
+
+% $k_a$ is estimated using the following formula:
+
+ka = 0.8/abs(dcgain(G('y', 'Fa')));
+
+
+% The factor can be adjusted to better match the curves.
+
+% Then $k_e$ and $k_1$ are computed.
+
+ke = ka/(beta/alpha - 1);
+k1 = 1/beta - ke*ka/(ke + ka);
+
+
+
+% #+RESULTS:
+% |    | Value [N/um] |
+% |----+--------------|
+% | ka |         40.5 |
+% | ke |          1.5 |
+% | k1 |          0.4 |
+
+% The damping in the system is adjusted to match the FEM model if necessary.
+
+c1 = 1e2;
+
+
+
+% The analytical model of the simpler system is defined below:
+
+Ga = 1/(m*s^2 + k1 + c1*s + ke*ka/(ke + ka)) * ...
+     [ 1 ,              k1 + c1*s + ke*ka/(ke + ka)  , ke/(ke + ka) ;
+      -ke*ka/(ke + ka), ke*ka/(ke + ka)*m*s^2 ,       -ke/(ke + ka)*(m*s^2 + c1*s + k1)];
+
+Ga.InputName = {'Fd', 'w', 'Fa'};
+Ga.OutputName = {'y', 'Fs'};
+
+
+
+% And the DC gain is adjusted for the force sensor:
+
+F_gain = dcgain(G('Fs', 'Fd'))/dcgain(Ga('Fs', 'Fd'));
+
+
+
+% The dynamics of the FEM model and the simpler model are compared in Figure [[fig:apa300ml_comp_simpler_model]].
+
+
+freqs = logspace(0, 5, 1000);
+
+figure;
+tiledlayout(2, 3, 'TileSpacing', 'None', 'Padding', 'None');
+
+ax1 = nexttile;
+hold on;
+plot(freqs, abs(squeeze(freqresp(G( 'y', 'w'), freqs, 'Hz'))));
+plot(freqs, abs(squeeze(freqresp(Ga('y', 'w'), freqs, 'Hz'))));
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+set(gca, 'XTickLabel',[]);
+ylabel('$x_1/w$ [m/m]');
+ylim([1e-6, 1e2]);
+
+ax2 = nexttile;
+hold on;
+plot(freqs, abs(squeeze(freqresp(G( 'y', 'Fa'), freqs, 'Hz'))));
+plot(freqs, abs(squeeze(freqresp(Ga('y', 'Fa'), freqs, 'Hz'))));
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+set(gca, 'XTickLabel',[]);
+ylabel('$x_1/f$ [m/N]');
+ylim([1e-14, 1e-6]);
+
+ax3 = nexttile;
+hold on;
+plot(freqs, abs(squeeze(freqresp(G( 'y', 'Fd'), freqs, 'Hz'))));
+plot(freqs, abs(squeeze(freqresp(Ga('y', 'Fd'), freqs, 'Hz'))));
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+set(gca, 'XTickLabel',[]);
+ylabel('$x_1/F$ [m/N]');
+ylim([1e-14, 1e-4]);
+
+ax4 = nexttile;
+hold on;
+plot(freqs, abs(squeeze(freqresp(G( 'Fs', 'w'), freqs, 'Hz'))));
+plot(freqs, abs(squeeze(freqresp(F_gain*Ga('Fs', 'w'), freqs, 'Hz'))));
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('Frequency [Hz]');
+ylabel('$F_s/w$ [m/m]');
+ylim([1e2, 1e8]);
+
+ax5 = nexttile;
+hold on;
+plot(freqs, abs(squeeze(freqresp(G( 'Fs', 'Fa'), freqs, 'Hz'))));
+plot(freqs, abs(squeeze(freqresp(F_gain*Ga('Fs', 'Fa'), freqs, 'Hz'))));
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('Frequency [Hz]');
+ylabel('$F_s/f$ [m/N]');
+ylim([1e-4, 1e1]);
+
+ax6 = nexttile;
+hold on;
+plot(freqs, abs(squeeze(freqresp(G( 'Fs', 'Fd'), freqs, 'Hz'))));
+plot(freqs, abs(squeeze(freqresp(F_gain*Ga('Fs', 'Fd'), freqs, 'Hz'))));
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('Frequency [Hz]');
+ylabel('$F_s/F$ [m/N]');
+ylim([1e-7, 1e2]);
+
+linkaxes([ax1,ax2,ax3,ax4,ax5,ax6],'x');
+
+
+
+% #+name: fig:apa300ml_comp_simpler_model
+% #+caption: Comparison of the Dynamics between the FEM model and the simplified one
+% #+RESULTS:
+% [[file:figs/apa300ml_comp_simpler_model.png]]
+
+% The simplified model has also been implemented in Simscape.
+
+% The dynamics of the Simscape simplified model is identified and compared with the FEM one in Figure [[fig:apa300ml_comp_simpler_simscape]].
+
+%% Name of the Simulink File
+mdl = 'APA300ML_simplified';
+
+%% Input/Output definition
+clear io; io_i = 1;
+io(io_i) = linio([mdl, '/Fd'], 1, 'openinput');  io_i = io_i + 1; % External Vertical Force [N]
+io(io_i) = linio([mdl, '/w'],  1, 'openinput');  io_i = io_i + 1; % Base Motion [m]
+io(io_i) = linio([mdl, '/Fa'], 1, 'openinput');  io_i = io_i + 1; % Actuator Force [N]
+io(io_i) = linio([mdl, '/y'],  1, 'openoutput'); io_i = io_i + 1; % Vertical Displacement [m]
+io(io_i) = linio([mdl, '/Fs'], 1, 'openoutput'); io_i = io_i + 1; % Force Sensor [V]
+
+Gs = linearize(mdl, io);
+
+Gs.InputName = {'Fd', 'w', 'Fa'};
+Gs.OutputName = {'y', 'Fs'};
+
+freqs = logspace(0, 5, 1000);
+
+figure;
+tiledlayout(2, 3, 'TileSpacing', 'None', 'Padding', 'None');
+
+ax1 = nexttile;
+hold on;
+plot(freqs, abs(squeeze(freqresp(G( 'y', 'w'), freqs, 'Hz'))));
+plot(freqs, abs(squeeze(freqresp(Gs('y', 'w'), freqs, 'Hz'))));
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+set(gca, 'XTickLabel',[]);
+ylabel('$x_1/w$ [m/m]');
+ylim([1e-6, 1e2]);
+
+ax2 = nexttile;
+hold on;
+plot(freqs, abs(squeeze(freqresp(G( 'y', 'Fa'), freqs, 'Hz'))));
+plot(freqs, abs(squeeze(freqresp(Gs('y', 'Fa'), freqs, 'Hz'))));
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+set(gca, 'XTickLabel',[]);
+ylabel('$x_1/f$ [m/N]');
+ylim([1e-14, 1e-6]);
+
+ax3 = nexttile;
+hold on;
+plot(freqs, abs(squeeze(freqresp(G( 'y', 'Fd'), freqs, 'Hz'))));
+plot(freqs, abs(squeeze(freqresp(Gs('y', 'Fd'), freqs, 'Hz'))));
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+set(gca, 'XTickLabel',[]);
+ylabel('$x_1/F$ [m/N]');
+ylim([1e-14, 1e-4]);
+
+ax4 = nexttile;
+hold on;
+plot(freqs, abs(squeeze(freqresp(G( 'Fs', 'w'), freqs, 'Hz'))));
+plot(freqs, abs(squeeze(freqresp(F_gain*Gs('Fs', 'w'), freqs, 'Hz'))));
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('Frequency [Hz]');
+ylabel('$F_s/w$ [m/m]');
+ylim([1e2, 1e8]);
+
+ax5 = nexttile;
+hold on;
+plot(freqs, abs(squeeze(freqresp(G( 'Fs', 'Fa'), freqs, 'Hz'))));
+plot(freqs, abs(squeeze(freqresp(F_gain*Gs('Fs', 'Fa'), freqs, 'Hz'))));
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('Frequency [Hz]');
+ylabel('$F_s/f$ [m/N]');
+ylim([1e-4, 1e1]);
+
+ax6 = nexttile;
+hold on;
+plot(freqs, abs(squeeze(freqresp(G( 'Fs', 'Fd'), freqs, 'Hz'))));
+plot(freqs, abs(squeeze(freqresp(F_gain*Gs('Fs', 'Fd'), freqs, 'Hz'))));
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('Frequency [Hz]');
+ylabel('$F_s/F$ [m/N]');
+ylim([1e-7, 1e2]);
+
+linkaxes([ax1,ax2,ax3,ax4,ax5,ax6],'x');
diff --git a/matlab/APA300ML.slx b/matlab/APA300ML.slx
new file mode 100644
index 0000000..f119adc
Binary files /dev/null and b/matlab/APA300ML.slx differ
diff --git a/data/APA300ML_new/APA300ML.step b/matlab/APA300ML/APA300ML.step
similarity index 100%
rename from data/APA300ML_new/APA300ML.step
rename to matlab/APA300ML/APA300ML.step
diff --git a/data/APA300ML_new/CS.jpg b/matlab/APA300ML/CS.jpg
similarity index 100%
rename from data/APA300ML_new/CS.jpg
rename to matlab/APA300ML/CS.jpg
diff --git a/data/APA300ML_new/MDoF.jpg b/matlab/APA300ML/MDoF.jpg
similarity index 100%
rename from data/APA300ML_new/MDoF.jpg
rename to matlab/APA300ML/MDoF.jpg
diff --git a/data/APA300ML_new/cp112.sub b/matlab/APA300ML/cp112.sub
similarity index 100%
rename from data/APA300ML_new/cp112.sub
rename to matlab/APA300ML/cp112.sub
diff --git a/data/APA300ML_new/mat_K.CSV b/matlab/APA300ML/mat_K.CSV
similarity index 100%
rename from data/APA300ML_new/mat_K.CSV
rename to matlab/APA300ML/mat_K.CSV
diff --git a/data/APA300ML_new/mat_K.matrix b/matlab/APA300ML/mat_K.matrix
similarity index 100%
rename from data/APA300ML_new/mat_K.matrix
rename to matlab/APA300ML/mat_K.matrix
diff --git a/data/APA300ML_new/mat_M.CSV b/matlab/APA300ML/mat_M.CSV
similarity index 100%
rename from data/APA300ML_new/mat_M.CSV
rename to matlab/APA300ML/mat_M.CSV
diff --git a/data/APA300ML_new/mat_M.matrix b/matlab/APA300ML/mat_M.matrix
similarity index 100%
rename from data/APA300ML_new/mat_M.matrix
rename to matlab/APA300ML/mat_M.matrix
diff --git a/data/APA300ML_new/out_nodes_3D.txt b/matlab/APA300ML/out_nodes_3D.txt
similarity index 100%
rename from data/APA300ML_new/out_nodes_3D.txt
rename to matlab/APA300ML/out_nodes_3D.txt
diff --git a/data/APA300ML_new/APA300ML_test_bench_simplified.slx b/matlab/APA300ML_simplified.slx
similarity index 100%
rename from data/APA300ML_new/APA300ML_test_bench_simplified.slx
rename to matlab/APA300ML_simplified.slx
diff --git a/data/APA300ML_new/APA300ML_characterisation.slx b/matlab/to-order/APA300ML/to-order/APA300ML_characterisation.slx
similarity index 100%
rename from data/APA300ML_new/APA300ML_characterisation.slx
rename to matlab/to-order/APA300ML/to-order/APA300ML_characterisation.slx
diff --git a/data/APA300ML/APA300ML_characterisation_with_joints.slx b/matlab/to-order/APA300ML/to-order/APA300ML_characterisation_with_joints.slx
similarity index 100%
rename from data/APA300ML/APA300ML_characterisation_with_joints.slx
rename to matlab/to-order/APA300ML/to-order/APA300ML_characterisation_with_joints.slx
diff --git a/data/APA300ML_new/APA300ML_old.step b/matlab/to-order/APA300ML/to-order/APA300ML_old.step
similarity index 100%
rename from data/APA300ML_new/APA300ML_old.step
rename to matlab/to-order/APA300ML/to-order/APA300ML_old.step
diff --git a/data/APA300ML_new/APA300ML_old_old.step b/matlab/to-order/APA300ML/to-order/APA300ML_old_old.step
similarity index 100%
rename from data/APA300ML_new/APA300ML_old_old.step
rename to matlab/to-order/APA300ML/to-order/APA300ML_old_old.step
diff --git a/data/APA300ML/APA300ML_characterisation.slx b/matlab/to-order/APA300ML_old/APA300ML_characterisation.slx
similarity index 100%
rename from data/APA300ML/APA300ML_characterisation.slx
rename to matlab/to-order/APA300ML_old/APA300ML_characterisation.slx
diff --git a/data/APA300ML_new/APA300ML_characterisation_with_joints.slx b/matlab/to-order/APA300ML_old/APA300ML_characterisation_with_joints.slx
similarity index 100%
rename from data/APA300ML_new/APA300ML_characterisation_with_joints.slx
rename to matlab/to-order/APA300ML_old/APA300ML_characterisation_with_joints.slx
diff --git a/data/APA300ML/APA300ML_flex_joints.slx b/matlab/to-order/APA300ML_old/APA300ML_flex_joints.slx
similarity index 100%
rename from data/APA300ML/APA300ML_flex_joints.slx
rename to matlab/to-order/APA300ML_old/APA300ML_flex_joints.slx
diff --git a/data/APA300ML/Model-Mdof.jpg b/matlab/to-order/APA300ML_old/Model-Mdof.jpg
similarity index 100%
rename from data/APA300ML/Model-Mdof.jpg
rename to matlab/to-order/APA300ML_old/Model-Mdof.jpg
diff --git a/data/APA300ML/Model.jpg b/matlab/to-order/APA300ML_old/Model.jpg
similarity index 100%
rename from data/APA300ML/Model.jpg
rename to matlab/to-order/APA300ML_old/Model.jpg
diff --git a/data/APA300ML/Nodes_MDoF_NLIST_MLIST.txt b/matlab/to-order/APA300ML_old/Nodes_MDoF_NLIST_MLIST.txt
similarity index 100%
rename from data/APA300ML/Nodes_MDoF_NLIST_MLIST.txt
rename to matlab/to-order/APA300ML_old/Nodes_MDoF_NLIST_MLIST.txt
diff --git a/data/APA300ML/mat_K-48modes-7MDoF.matrix b/matlab/to-order/APA300ML_old/mat_K-48modes-7MDoF.matrix
similarity index 100%
rename from data/APA300ML/mat_K-48modes-7MDoF.matrix
rename to matlab/to-order/APA300ML_old/mat_K-48modes-7MDoF.matrix
diff --git a/data/APA300ML/mat_K-80modes-7MDoF.matrix b/matlab/to-order/APA300ML_old/mat_K-80modes-7MDoF.matrix
similarity index 100%
rename from data/APA300ML/mat_K-80modes-7MDoF.matrix
rename to matlab/to-order/APA300ML_old/mat_K-80modes-7MDoF.matrix
diff --git a/data/APA300ML/mat_M-48modes-7MDoF.matrix b/matlab/to-order/APA300ML_old/mat_M-48modes-7MDoF.matrix
similarity index 100%
rename from data/APA300ML/mat_M-48modes-7MDoF.matrix
rename to matlab/to-order/APA300ML_old/mat_M-48modes-7MDoF.matrix
diff --git a/data/APA300ML/mat_M-80modes-7MDoF.matrix b/matlab/to-order/APA300ML_old/mat_M-80modes-7MDoF.matrix
similarity index 100%
rename from data/APA300ML/mat_M-80modes-7MDoF.matrix
rename to matlab/to-order/APA300ML_old/mat_M-80modes-7MDoF.matrix
diff --git a/data/flexor_ID16/22921061.STEP b/matlab/to-order/flexor_ID16/22921061.STEP
similarity index 100%
rename from data/flexor_ID16/22921061.STEP
rename to matlab/to-order/flexor_ID16/22921061.STEP
diff --git a/data/flexor_ID16/cp91_2MDoF_6modes.sub b/matlab/to-order/flexor_ID16/data/cp91_2MDoF_6modes.sub
similarity index 100%
rename from data/flexor_ID16/cp91_2MDoF_6modes.sub
rename to matlab/to-order/flexor_ID16/data/cp91_2MDoF_6modes.sub
diff --git a/data/flexor_ID16/mat_K_6modes_2MDoF.matrix b/matlab/to-order/flexor_ID16/data/mat_K_6modes_2MDoF.matrix
similarity index 100%
rename from data/flexor_ID16/mat_K_6modes_2MDoF.matrix
rename to matlab/to-order/flexor_ID16/data/mat_K_6modes_2MDoF.matrix
diff --git a/data/flexor_ID16/mat_M_6modes_2MDoF.matrix b/matlab/to-order/flexor_ID16/data/mat_M_6modes_2MDoF.matrix
similarity index 100%
rename from data/flexor_ID16/mat_M_6modes_2MDoF.matrix
rename to matlab/to-order/flexor_ID16/data/mat_M_6modes_2MDoF.matrix
diff --git a/data/flexor_ID16/out_nodes_3D.txt b/matlab/to-order/flexor_ID16/data/out_nodes_3D.txt
similarity index 100%
rename from data/flexor_ID16/out_nodes_3D.txt
rename to matlab/to-order/flexor_ID16/data/out_nodes_3D.txt
diff --git a/data/flexor_ID16/flexor_ID16.STEP b/matlab/to-order/flexor_ID16/flexor_ID16.STEP
similarity index 100%
rename from data/flexor_ID16/flexor_ID16.STEP
rename to matlab/to-order/flexor_ID16/flexor_ID16.STEP
diff --git a/data/flexor_ID16/flexor_ID16.slx b/matlab/to-order/flexor_ID16/flexor_ID16.slx
similarity index 100%
rename from data/flexor_ID16/flexor_ID16.slx
rename to matlab/to-order/flexor_ID16/flexor_ID16.slx
diff --git a/data/flexor_ID16/flexor_ID16_simplified.slx b/matlab/to-order/flexor_ID16/flexor_ID16_simplified.slx
similarity index 100%
rename from data/flexor_ID16/flexor_ID16_simplified.slx
rename to matlab/to-order/flexor_ID16/flexor_ID16_simplified.slx
diff --git a/data/flexor_circ_025/Flex-DC-0.25.step b/matlab/to-order/flexor_circ_025/Flex-DC-0.25.step
similarity index 100%
rename from data/flexor_circ_025/Flex-DC-0.25.step
rename to matlab/to-order/flexor_circ_025/Flex-DC-0.25.step
diff --git a/data/flexor_circ_025/mat_K.CSV b/matlab/to-order/flexor_circ_025/data/mat_K.CSV
similarity index 100%
rename from data/flexor_circ_025/mat_K.CSV
rename to matlab/to-order/flexor_circ_025/data/mat_K.CSV
diff --git a/data/flexor_circ_025/mat_K.matrix b/matlab/to-order/flexor_circ_025/data/mat_K.matrix
similarity index 100%
rename from data/flexor_circ_025/mat_K.matrix
rename to matlab/to-order/flexor_circ_025/data/mat_K.matrix
diff --git a/data/flexor_circ_025/mat_M.CSV b/matlab/to-order/flexor_circ_025/data/mat_M.CSV
similarity index 100%
rename from data/flexor_circ_025/mat_M.CSV
rename to matlab/to-order/flexor_circ_025/data/mat_M.CSV
diff --git a/data/flexor_circ_025/mat_M.matrix b/matlab/to-order/flexor_circ_025/data/mat_M.matrix
similarity index 100%
rename from data/flexor_circ_025/mat_M.matrix
rename to matlab/to-order/flexor_circ_025/data/mat_M.matrix
diff --git a/data/flexor_circ_025/out_nodes_3D.txt b/matlab/to-order/flexor_circ_025/data/out_nodes_3D.txt
similarity index 100%
rename from data/flexor_circ_025/out_nodes_3D.txt
rename to matlab/to-order/flexor_circ_025/data/out_nodes_3D.txt
diff --git a/data/flexor_circ_025/CS.jpg b/matlab/to-order/flexor_circ_025/figs/CS.jpg
similarity index 100%
rename from data/flexor_circ_025/CS.jpg
rename to matlab/to-order/flexor_circ_025/figs/CS.jpg
diff --git a/data/flexor_circ_025/Geometry.jpg b/matlab/to-order/flexor_circ_025/figs/Geometry.jpg
similarity index 100%
rename from data/flexor_circ_025/Geometry.jpg
rename to matlab/to-order/flexor_circ_025/figs/Geometry.jpg
diff --git a/data/flexor_circ_025/MDoF.jpg b/matlab/to-order/flexor_circ_025/figs/MDoF.jpg
similarity index 100%
rename from data/flexor_circ_025/MDoF.jpg
rename to matlab/to-order/flexor_circ_025/figs/MDoF.jpg
diff --git a/data/flexor_circ_025/flexor_025.slx b/matlab/to-order/flexor_circ_025/flexor_025.slx
similarity index 100%
rename from data/flexor_circ_025/flexor_025.slx
rename to matlab/to-order/flexor_circ_025/flexor_025.slx
diff --git a/data/flexor_circ_025/flexor_025_simplified.slx b/matlab/to-order/flexor_circ_025/flexor_025_simplified.slx
similarity index 100%
rename from data/flexor_circ_025/flexor_025_simplified.slx
rename to matlab/to-order/flexor_circ_025/flexor_025_simplified.slx
diff --git a/data/flexor_circ_050/Flex-DC-0.50.step b/matlab/to-order/flexor_circ_050/Flex-DC-0.50.step
similarity index 100%
rename from data/flexor_circ_050/Flex-DC-0.50.step
rename to matlab/to-order/flexor_circ_050/Flex-DC-0.50.step
diff --git a/data/flexor_circ_050/mat_K.CSV b/matlab/to-order/flexor_circ_050/data/mat_K.CSV
similarity index 100%
rename from data/flexor_circ_050/mat_K.CSV
rename to matlab/to-order/flexor_circ_050/data/mat_K.CSV
diff --git a/data/flexor_circ_050/mat_K.matrix b/matlab/to-order/flexor_circ_050/data/mat_K.matrix
similarity index 100%
rename from data/flexor_circ_050/mat_K.matrix
rename to matlab/to-order/flexor_circ_050/data/mat_K.matrix
diff --git a/data/flexor_circ_050/mat_M.CSV b/matlab/to-order/flexor_circ_050/data/mat_M.CSV
similarity index 100%
rename from data/flexor_circ_050/mat_M.CSV
rename to matlab/to-order/flexor_circ_050/data/mat_M.CSV
diff --git a/data/flexor_circ_050/mat_M.matrix b/matlab/to-order/flexor_circ_050/data/mat_M.matrix
similarity index 100%
rename from data/flexor_circ_050/mat_M.matrix
rename to matlab/to-order/flexor_circ_050/data/mat_M.matrix
diff --git a/data/flexor_circ_050/out_nodes_3D.txt b/matlab/to-order/flexor_circ_050/data/out_nodes_3D.txt
similarity index 100%
rename from data/flexor_circ_050/out_nodes_3D.txt
rename to matlab/to-order/flexor_circ_050/data/out_nodes_3D.txt
diff --git a/data/flexor_circ_050/CS.jpg b/matlab/to-order/flexor_circ_050/figs/CS.jpg
similarity index 100%
rename from data/flexor_circ_050/CS.jpg
rename to matlab/to-order/flexor_circ_050/figs/CS.jpg
diff --git a/data/flexor_circ_050/Geometry.jpg b/matlab/to-order/flexor_circ_050/figs/Geometry.jpg
similarity index 100%
rename from data/flexor_circ_050/Geometry.jpg
rename to matlab/to-order/flexor_circ_050/figs/Geometry.jpg
diff --git a/data/flexor_circ_050/MDoF.jpg b/matlab/to-order/flexor_circ_050/figs/MDoF.jpg
similarity index 100%
rename from data/flexor_circ_050/MDoF.jpg
rename to matlab/to-order/flexor_circ_050/figs/MDoF.jpg
diff --git a/data/piezo_amplified_3d/FEA_HarmResponse_00kg.txt b/matlab/to-order/piezo_amplified_3d/data/FEA_HarmResponse_00kg.txt
similarity index 100%
rename from data/piezo_amplified_3d/FEA_HarmResponse_00kg.txt
rename to matlab/to-order/piezo_amplified_3d/data/FEA_HarmResponse_00kg.txt
diff --git a/data/piezo_amplified_3d/FEA_HarmResponse_10kg.txt b/matlab/to-order/piezo_amplified_3d/data/FEA_HarmResponse_10kg.txt
similarity index 100%
rename from data/piezo_amplified_3d/FEA_HarmResponse_10kg.txt
rename to matlab/to-order/piezo_amplified_3d/data/FEA_HarmResponse_10kg.txt
diff --git a/data/piezo_amplified_3d/mat_K_108dof.matrix b/matlab/to-order/piezo_amplified_3d/data/mat_K_108dof.matrix
similarity index 100%
rename from data/piezo_amplified_3d/mat_K_108dof.matrix
rename to matlab/to-order/piezo_amplified_3d/data/mat_K_108dof.matrix
diff --git a/data/piezo_amplified_3d/mat_M_108dof.matrix b/matlab/to-order/piezo_amplified_3d/data/mat_M_108dof.matrix
similarity index 100%
rename from data/piezo_amplified_3d/mat_M_108dof.matrix
rename to matlab/to-order/piezo_amplified_3d/data/mat_M_108dof.matrix
diff --git a/data/piezo_amplified_3d/out_nodes_3D.txt b/matlab/to-order/piezo_amplified_3d/data/out_nodes_3D.txt
similarity index 100%
rename from data/piezo_amplified_3d/out_nodes_3D.txt
rename to matlab/to-order/piezo_amplified_3d/data/out_nodes_3D.txt
diff --git a/data/piezo_amplified_3d/piezo_amplified_3d.txt b/matlab/to-order/piezo_amplified_3d/data/piezo_amplified_3d.txt
similarity index 100%
rename from data/piezo_amplified_3d/piezo_amplified_3d.txt
rename to matlab/to-order/piezo_amplified_3d/data/piezo_amplified_3d.txt
diff --git a/data/piezo_amplified_3d/piezo_amplified_3d_K.txt b/matlab/to-order/piezo_amplified_3d/data/piezo_amplified_3d_K.txt
similarity index 100%
rename from data/piezo_amplified_3d/piezo_amplified_3d_K.txt
rename to matlab/to-order/piezo_amplified_3d/data/piezo_amplified_3d_K.txt
diff --git a/data/piezo_amplified_3d/piezo_amplified_3d_M.txt b/matlab/to-order/piezo_amplified_3d/data/piezo_amplified_3d_M.txt
similarity index 100%
rename from data/piezo_amplified_3d/piezo_amplified_3d_M.txt
rename to matlab/to-order/piezo_amplified_3d/data/piezo_amplified_3d_M.txt
diff --git a/data/piezo_amplified_3d/piezo_amplified_3d.slx b/matlab/to-order/piezo_amplified_3d/piezo_amplified_3d.slx
similarity index 100%
rename from data/piezo_amplified_3d/piezo_amplified_3d.slx
rename to matlab/to-order/piezo_amplified_3d/piezo_amplified_3d.slx
diff --git a/data/piezo_amplified_3d/piezo_amplified_3d_distri.slx b/matlab/to-order/piezo_amplified_3d/piezo_amplified_3d_distri.slx
similarity index 100%
rename from data/piezo_amplified_3d/piezo_amplified_3d_distri.slx
rename to matlab/to-order/piezo_amplified_3d/piezo_amplified_3d_distri.slx
diff --git a/data/piezo_amplified_3d/piezo_amplified_3d_distri_act_sens.slx b/matlab/to-order/piezo_amplified_3d/piezo_amplified_3d_distri_act_sens.slx
similarity index 100%
rename from data/piezo_amplified_3d/piezo_amplified_3d_distri_act_sens.slx
rename to matlab/to-order/piezo_amplified_3d/piezo_amplified_3d_distri_act_sens.slx
diff --git a/data/piezo_amplified_3d/piezo_amplified_3d_identification.slx b/matlab/to-order/piezo_amplified_3d/piezo_amplified_3d_identification.slx
similarity index 100%
rename from data/piezo_amplified_3d/piezo_amplified_3d_identification.slx
rename to matlab/to-order/piezo_amplified_3d/piezo_amplified_3d_identification.slx
diff --git a/data/piezo_amplified_3d/solve_modal_00kg.out b/matlab/to-order/piezo_amplified_3d/solve_modal_00kg.out
similarity index 100%
rename from data/piezo_amplified_3d/solve_modal_00kg.out
rename to matlab/to-order/piezo_amplified_3d/solve_modal_00kg.out
diff --git a/data/piezo_amplified_3d/solve_modal_10kg.out b/matlab/to-order/piezo_amplified_3d/solve_modal_10kg.out
similarity index 100%
rename from data/piezo_amplified_3d/solve_modal_10kg.out
rename to matlab/to-order/piezo_amplified_3d/solve_modal_10kg.out
diff --git a/data/piezo_amplified_IFF/FEA_HarmResponse_00kg.txt b/matlab/to-order/piezo_amplified_IFF/data/FEA_HarmResponse_00kg.txt
similarity index 100%
rename from data/piezo_amplified_IFF/FEA_HarmResponse_00kg.txt
rename to matlab/to-order/piezo_amplified_IFF/data/FEA_HarmResponse_00kg.txt
diff --git a/data/piezo_amplified_IFF/FEA_HarmResponse_10kg.txt b/matlab/to-order/piezo_amplified_IFF/data/FEA_HarmResponse_10kg.txt
similarity index 100%
rename from data/piezo_amplified_IFF/FEA_HarmResponse_10kg.txt
rename to matlab/to-order/piezo_amplified_IFF/data/FEA_HarmResponse_10kg.txt
diff --git a/data/piezo_amplified_IFF/piezo_amplified_3d.txt b/matlab/to-order/piezo_amplified_IFF/data/piezo_amplified_3d.txt
similarity index 100%
rename from data/piezo_amplified_IFF/piezo_amplified_3d.txt
rename to matlab/to-order/piezo_amplified_IFF/data/piezo_amplified_3d.txt
diff --git a/data/piezo_amplified_IFF/piezo_amplified_3d_K.txt b/matlab/to-order/piezo_amplified_IFF/data/piezo_amplified_3d_K.txt
similarity index 100%
rename from data/piezo_amplified_IFF/piezo_amplified_3d_K.txt
rename to matlab/to-order/piezo_amplified_IFF/data/piezo_amplified_3d_K.txt
diff --git a/data/piezo_amplified_IFF/piezo_amplified_3d_M.txt b/matlab/to-order/piezo_amplified_IFF/data/piezo_amplified_3d_M.txt
similarity index 100%
rename from data/piezo_amplified_IFF/piezo_amplified_3d_M.txt
rename to matlab/to-order/piezo_amplified_IFF/data/piezo_amplified_3d_M.txt
diff --git a/data/piezo_amplified_IFF/piezo_amplified_IFF.txt b/matlab/to-order/piezo_amplified_IFF/data/piezo_amplified_IFF.txt
similarity index 100%
rename from data/piezo_amplified_IFF/piezo_amplified_IFF.txt
rename to matlab/to-order/piezo_amplified_IFF/data/piezo_amplified_IFF.txt
diff --git a/data/piezo_amplified_IFF/piezo_amplified_IFF_K.txt b/matlab/to-order/piezo_amplified_IFF/data/piezo_amplified_IFF_K.txt
similarity index 100%
rename from data/piezo_amplified_IFF/piezo_amplified_IFF_K.txt
rename to matlab/to-order/piezo_amplified_IFF/data/piezo_amplified_IFF_K.txt
diff --git a/data/piezo_amplified_IFF/piezo_amplified_IFF_M.txt b/matlab/to-order/piezo_amplified_IFF/data/piezo_amplified_IFF_M.txt
similarity index 100%
rename from data/piezo_amplified_IFF/piezo_amplified_IFF_M.txt
rename to matlab/to-order/piezo_amplified_IFF/data/piezo_amplified_IFF_M.txt
diff --git a/data/piezo_amplified_IFF/piezo_amplified_3d.slx b/matlab/to-order/piezo_amplified_IFF/piezo_amplified_3d.slx
similarity index 100%
rename from data/piezo_amplified_IFF/piezo_amplified_3d.slx
rename to matlab/to-order/piezo_amplified_IFF/piezo_amplified_3d.slx
diff --git a/data/piezo_amplified_IFF/piezo_amplified_3d_distri.slx b/matlab/to-order/piezo_amplified_IFF/piezo_amplified_3d_distri.slx
similarity index 100%
rename from data/piezo_amplified_IFF/piezo_amplified_3d_distri.slx
rename to matlab/to-order/piezo_amplified_IFF/piezo_amplified_3d_distri.slx
diff --git a/data/piezo_amplified_IFF/piezo_amplified_3d_distri_act_sens.slx b/matlab/to-order/piezo_amplified_IFF/piezo_amplified_3d_distri_act_sens.slx
similarity index 100%
rename from data/piezo_amplified_IFF/piezo_amplified_3d_distri_act_sens.slx
rename to matlab/to-order/piezo_amplified_IFF/piezo_amplified_3d_distri_act_sens.slx
diff --git a/data/piezo_amplified_IFF/piezo_amplified_IFF.slx b/matlab/to-order/piezo_amplified_IFF/piezo_amplified_IFF.slx
similarity index 100%
rename from data/piezo_amplified_IFF/piezo_amplified_IFF.slx
rename to matlab/to-order/piezo_amplified_IFF/piezo_amplified_IFF.slx
diff --git a/data/piezo_amplified_IFF/solve_modal_00kg.out b/matlab/to-order/piezo_amplified_IFF/solve_modal_00kg.out
similarity index 100%
rename from data/piezo_amplified_IFF/solve_modal_00kg.out
rename to matlab/to-order/piezo_amplified_IFF/solve_modal_00kg.out
diff --git a/data/piezo_amplified_IFF/solve_modal_10kg.out b/matlab/to-order/piezo_amplified_IFF/solve_modal_10kg.out
similarity index 100%
rename from data/piezo_amplified_IFF/solve_modal_10kg.out
rename to matlab/to-order/piezo_amplified_IFF/solve_modal_10kg.out
diff --git a/data/strut_encoder/cp143.sub b/matlab/to-order/strut_encoder/data/cp143.sub
similarity index 100%
rename from data/strut_encoder/cp143.sub
rename to matlab/to-order/strut_encoder/data/cp143.sub
diff --git a/data/strut_encoder/mat_K.CSV b/matlab/to-order/strut_encoder/data/mat_K.CSV
similarity index 100%
rename from data/strut_encoder/mat_K.CSV
rename to matlab/to-order/strut_encoder/data/mat_K.CSV
diff --git a/data/strut_encoder/mat_K.matrix b/matlab/to-order/strut_encoder/data/mat_K.matrix
similarity index 100%
rename from data/strut_encoder/mat_K.matrix
rename to matlab/to-order/strut_encoder/data/mat_K.matrix
diff --git a/data/strut_encoder/mat_M.CSV b/matlab/to-order/strut_encoder/data/mat_M.CSV
similarity index 100%
rename from data/strut_encoder/mat_M.CSV
rename to matlab/to-order/strut_encoder/data/mat_M.CSV
diff --git a/data/strut_encoder/mat_M.matrix b/matlab/to-order/strut_encoder/data/mat_M.matrix
similarity index 100%
rename from data/strut_encoder/mat_M.matrix
rename to matlab/to-order/strut_encoder/data/mat_M.matrix
diff --git a/data/strut_encoder/out_nodes_3D.txt b/matlab/to-order/strut_encoder/data/out_nodes_3D.txt
similarity index 100%
rename from data/strut_encoder/out_nodes_3D.txt
rename to matlab/to-order/strut_encoder/data/out_nodes_3D.txt
diff --git a/data/strut_encoder/global-cs1.jpg b/matlab/to-order/strut_encoder/figs/global-cs1.jpg
similarity index 100%
rename from data/strut_encoder/global-cs1.jpg
rename to matlab/to-order/strut_encoder/figs/global-cs1.jpg
diff --git a/data/strut_encoder/global-cs2.jpg b/matlab/to-order/strut_encoder/figs/global-cs2.jpg
similarity index 100%
rename from data/strut_encoder/global-cs2.jpg
rename to matlab/to-order/strut_encoder/figs/global-cs2.jpg
diff --git a/data/strut_encoder/points1.jpg b/matlab/to-order/strut_encoder/figs/points1.jpg
similarity index 100%
rename from data/strut_encoder/points1.jpg
rename to matlab/to-order/strut_encoder/figs/points1.jpg
diff --git a/data/strut_encoder/points2.jpg b/matlab/to-order/strut_encoder/figs/points2.jpg
similarity index 100%
rename from data/strut_encoder/points2.jpg
rename to matlab/to-order/strut_encoder/figs/points2.jpg
diff --git a/data/strut_encoder/points3.jpg b/matlab/to-order/strut_encoder/figs/points3.jpg
similarity index 100%
rename from data/strut_encoder/points3.jpg
rename to matlab/to-order/strut_encoder/figs/points3.jpg
diff --git a/data/strut_encoder/strut_encoder.slx b/matlab/to-order/strut_encoder/strut_encoder.slx
similarity index 100%
rename from data/strut_encoder/strut_encoder.slx
rename to matlab/to-order/strut_encoder/strut_encoder.slx
diff --git a/data/strut_encoder/strut_encoder.step b/matlab/to-order/strut_encoder/strut_encoder.step
similarity index 100%
rename from data/strut_encoder/strut_encoder.step
rename to matlab/to-order/strut_encoder/strut_encoder.step