diff --git a/simscape-micro-station.org b/simscape-micro-station.org
index ccc0db2..1601ed1 100644
--- a/simscape-micro-station.org
+++ b/simscape-micro-station.org
@@ -164,9 +164,8 @@ initializeController(  'type', 'open-loop');
 
 #+end_src
 
-** TODO [#B] Make good "init" for the Simscape model
-
-** TODO [#B] Just keep smallest number of variants for each stage
+** DONE [#B] Just keep smallest number of variants for each stage
+CLOSED: [2024-10-30 Wed 16:15]
 
 - [ ] none
 - [ ] rigid
@@ -178,10 +177,41 @@ initializeController(  'type', 'open-loop');
 SCHEDULED: <2024-10-30 Wed>
 
 - [ ] Find the compliance measurements
+  - The one from modal analysis
+  - [[file:~/Cloud/work-projects/ID31-NASS/matlab/nass-measurements/2018-01-12 - Marc]] : Compliance measurement in X,Y,Z with geophone, marble not glued
+  - [[file:~/Cloud/work-projects/ID31-NASS/matlab/nass-measurements/2018-10-12 - Marc]]: same but in the hutch with glued marble
+  - *08/2020*: [[file:~/Cloud/work-projects/ID31-NASS/matlab/nass-measurements/micro-station-compliance/index.org::+TITLE: Compliance Measurement of the Micro Station]]
 - [ ] See if it matches somehow the current model
 - [ ] If not, see if model parameters can be tuned to have better match
       For instance from values here: file:/home/thomas/Cloud/meetings/esrf-meetings/2018-04-24-Simscape-Model/2018-04-24-Simscape-Model.pdf
 
+** TODO [#B] Make good "init" for the Simscape model
+
+** TODO [#C] Could see the effect of each stage on the compliance
+
+- [ ] Put =rigid= mode one by one to see the effect
+
+** TODO [#C] Make a comparison of the measured vibrations of the micro-station with the vibrations of the simscape model of the micro-station
+Do we have a correlation? At least in the frequency domain?
+
+Procedure:
+- [ ] Take the time domain measurement of the vibrations due to spindle or translation stage
+- [ ] Take the estimated PSD of the estimated disturbance force
+- [ ] Simulate the Simscape model without the nano-hexapod (same conditions as when measuring the disturbances) and add only the considered disturbance
+- [ ] Save the position errors due to the disturbance
+- [ ] Compare with the measured vibrations
+
+
+** TODO [#C] Add two perturbations: =Frz_x= and =Frz_y=
+Maybe I can estimate them from the measurements that was made on the spindle?
+
+** TODO [#C] Determine which Degree-Of-Freedom to keep and which to constrain
+
+For instance, for now the granite can not rotate, but in reality, the modes may be linked to the granite's rotation.
+By constraining more DoF, the simulation will be faster and the obtain state space will have a lower order.
+
+** TODO [#B] Compute eigenvalues of the model to see if we have similar frequencies than the modal analysis?
+
 * Introduction                                                        :ignore:
 
 
@@ -894,7 +924,385 @@ end
   end
 #+end_src
 
-** Obtained Compliance of the Micro-Station
+** Measured micro-station compliance
+*** Introduction                                                    :ignore:
+
+The most important dynamical characteristic of the micro-station that should be well modeled is its compliance.
+This is what can impact the nano-hexapod dynamics.
+
+- [ ] Add schematic of the experiment with Micro-Hexapod top platform, location of accelerometers, of impacts, etc...
+
+- 4 3-axis accelerometers
+- 10 hammer impacts on the micro-hexapod top plaftorm
+- *Was the rotation compensation axis present?* (I don't think so)
+
+*** Position of inertial sensors on top of the micro-hexapod
+Orientation is relative to the frame determined by the X-ray
+
+| *Num* | *Position* | *Orientation* | *Sensibility* | *Channels* |
+|-------+------------+---------------+---------------+------------|
+|     1 | [0, +A, 0] | [x, y, z]     | 1V/g          |        1-3 |
+|     2 | [-B, 0, 0] | [x, y, z]     | 1V/g          |        4-6 |
+|     3 | [0, -A, 0] | [x, y, z]     | 0.1V/g        |        7-9 |
+|     4 | [+B, 0, 0] | [x, y, z]     | 1V/g          |      10-12 |
+
+Instrumented Hammer:
+- Channel 13
+- Sensibility: 230 uV/N
+
+| Acc Number | Dir | Channel Number |
+|------------+-----+----------------|
+|          1 | x   |              1 |
+|          1 | y   |              2 |
+|          1 | z   |              3 |
+|          2 | x   |              4 |
+|          2 | y   |              5 |
+|          2 | z   |              6 |
+|          3 | x   |              7 |
+|          3 | y   |              8 |
+|          3 | z   |              9 |
+|          4 | x   |             10 |
+|          4 | y   |             11 |
+|          4 | z   |             12 |
+|     Hammer |     |             13 |
+
+From the acceleration measurement of the 4 accelerometers, we can compute the translations and rotations:
+|       | *Formula*                | *Formula (numbers)*    |
+|-------+--------------------------+-----------------------|
+| $D_x$ | (1x + 2x + 3x + 4x)/4    | (1 + 4 + 7 + 10)/4    |
+| $D_y$ | (1y + 2y + 3y + 4y)/4    | (2 + 5 + 8 + 11)/4    |
+| $D_z$ | (1z + 2z + 3z + 4z)/4    | (3 + 6 + 9 + 12)/4    |
+| $R_x$ | (1z - 3z)/A              | (1 - 9)/A             |
+| $R_y$ | (2z - 4z)/B              | (6 - 12)/B            |
+| $R_z$ | (3x - 1x)/A, (4y - 2y)/B | (7 - 1)/A, (11 - 5)/B |
+
+dL = J X
+#+begin_src matlab
+d = 0.14;
+J = [1  0  0  0  0 -d;
+     0  1  0  0  0  0;
+     0  0  1  d  0  0;
+     1  0  0  0  0  0;
+     0  1  0  0  0 -d;
+     0  0  1  0  d  0;
+     1  0  0  0  0  d;
+     0  1  0  0  0  0;
+     0  0  1 -d  0  0;
+     1  0  0  0  0  0;
+     0  1  0  0  0  d;
+     0  0  1  0 -d  0];
+
+J_inv = pinv(J);
+#+end_src
+
+|     | Dx | Dy | Dz | Rx | Ry | Rz |
+|-----+----+----+----+----+----+----|
+| a1x |  1 |  0 |  0 |  0 |  0 | -d |
+| a1y |  0 |  1 |  0 |  0 |  0 |  0 |
+| a1z |  0 |  0 |  1 |  d |  0 |  0 |
+| a2x |  1 |  0 |  0 |  0 |  0 |  0 |
+| a2y |  0 |  1 |  0 |  0 |  0 | -d |
+| a2z |  0 |  0 |  1 |  0 |  d |  0 |
+| a3x |  1 |  0 |  0 |  0 |  0 |  d |
+| a3y |  0 |  1 |  0 |  0 |  0 |  0 |
+| a3z |  0 |  0 |  1 | -d |  0 |  0 |
+| a4x |  1 |  0 |  0 |  0 |  0 |  0 |
+| a4y |  0 |  1 |  0 |  0 |  0 |  d |
+| a4z |  0 |  0 |  1 |  0 | -d |  0 |
+
+#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
+data2orgtable(J_inv, {'Dx', 'Dy', 'Dz', 'Rx', 'Ry', 'Rz'}, {'d1x', 'd1y', 'd1z', 'd2x', 'd2y', 'd2z', 'd3x', 'd3y', 'd3z', 'd4x', 'd4y', 'd4z'}, ' %.5f ');
+#+end_src
+
+#+RESULTS:
+|    |      d1x |  d1y |     d1z |  d2x |      d2y |     d2z |     d3x |  d3y |      d3z |  d4x |     d4y |      d4z |
+|----+----------+------+---------+------+----------+---------+---------+------+----------+------+---------+----------|
+| Dx |     0.25 |  0.0 |     0.0 | 0.25 |      0.0 |     0.0 |    0.25 |  0.0 |      0.0 | 0.25 |     0.0 |      0.0 |
+| Dy |      0.0 | 0.25 |     0.0 |  0.0 |     0.25 |     0.0 |     0.0 | 0.25 |      0.0 |  0.0 |    0.25 |      0.0 |
+| Dz |      0.0 |  0.0 |    0.25 |  0.0 |      0.0 |    0.25 |     0.0 |  0.0 |     0.25 |  0.0 |     0.0 |     0.25 |
+| Rx |      0.0 |  0.0 | 3.57143 |  0.0 |      0.0 |    -0.0 |     0.0 |  0.0 | -3.57143 |  0.0 |     0.0 |     -0.0 |
+| Ry |      0.0 |  0.0 |     0.0 |  0.0 |      0.0 | 3.57143 |     0.0 |  0.0 |      0.0 |  0.0 |     0.0 | -3.57143 |
+| Rz | -1.78571 |  0.0 |    -0.0 |  0.0 | -1.78571 |     0.0 | 1.78571 |  0.0 |     -0.0 |  0.0 | 1.78571 |     -0.0 |
+
+*** Hammer blow position/orientation
+
+| *Num* | *Direction* | *Position* | Accelerometer position | Jacobian number |
+|-------+-------------+------------+------------------------+-----------------|
+|     1 | -Y          | [0, +A, 0] |                      1 |              -2 |
+|     2 | -Z          | [0, +A, 0] |                      1 |              -3 |
+|     3 | X           | [-B, 0, 0] |                      2 |               4 |
+|     4 | -Z          | [-B, 0, 0] |                      2 |              -6 |
+|     5 | Y           | [0, -A, 0] |                      3 |               8 |
+|     6 | -Z          | [0, -A, 0] |                      3 |              -9 |
+|     7 | -X          | [+B, 0, 0] |                      4 |             -10 |
+|     8 | -Z          | [+B, 0, 0] |                      4 |             -12 |
+|     9 | -X          | [0, -A, 0] |                      3 |              -7 |
+|    10 | -X          | [0, +A, 0] |                      1 |              -1 |
+
+From hammer blows to pure forces / torques:
+|       | *Formula*    | Alternative |
+|-------+--------------+-------------|
+| $F_x$ | +3           |          -7 |
+| $F_y$ | -1           |          +5 |
+| $F_z$ | -(2 + 6)/2   |  -(4 + 8)/2 |
+| $M_x$ | A/2*(2 - 6)  |             |
+| $M_y$ | B/2*(8 - 4)  |             |
+| $M_z$ | A/2*(10 - 9) |             |
+
+F = J' tau
+#+begin_src matlab
+Jf = [0 -1  0  0  0  0;
+      0  0 -1 -d  0  0;
+      1  0  0  0  0  0;
+      0  0 -1  0 -d  0;
+      0  1  0  0  0  0;
+      0  0 -1  d  0  0;
+     -1  0  0  0  0  0;
+      0  0 -1  0  d  0;
+     -1  0  0  0  0 -d;
+     -1  0  0  0  0  d]';
+Jf_inv = pinv(Jf);
+#+end_src
+
+#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
+data2orgtable(Jf, {'Fx', 'Fy', 'Fz', 'Mx', 'My', 'Mz'}, {'-F1y', '-F1z', '+F2x', '-F2z', '+F3y', '-F3z', '-F4x', '-F4z', '-F3x', '-F1x'}, ' %.2f ');
+#+end_src
+
+#+RESULTS:
+|    | -F1y |  -F1z | +F2x |  -F2z | +F3y | -F3z | -F4x | -F4z |  -F3x | -F1x |
+|----+------+-------+------+-------+------+------+------+------+-------+------|
+| Fx |  0.0 |   0.0 |  1.0 |   0.0 |  0.0 |  0.0 | -1.0 |  0.0 |  -1.0 | -1.0 |
+| Fy | -1.0 |   0.0 |  0.0 |   0.0 |  1.0 |  0.0 |  0.0 |  0.0 |   0.0 |  0.0 |
+| Fz |  0.0 |  -1.0 |  0.0 |  -1.0 |  0.0 | -1.0 |  0.0 | -1.0 |   0.0 |  0.0 |
+| Mx |  0.0 | -0.14 |  0.0 |   0.0 |  0.0 | 0.14 |  0.0 |  0.0 |   0.0 |  0.0 |
+| My |  0.0 |   0.0 |  0.0 | -0.14 |  0.0 |  0.0 |  0.0 | 0.14 |   0.0 |  0.0 |
+| Mz |  0.0 |   0.0 |  0.0 |   0.0 |  0.0 |  0.0 |  0.0 |  0.0 | -0.14 | 0.14 |
+
+*** Compute FRF
+Raw measurements are in file:/home/thomas/Cloud/work-projects/ID31-NASS/matlab/nass-measurements/micro-station-compliance/data/record
+
+For each measurement (10):
+- Find the location of the 10 impacts based on "track13"
+- Then for the 12 accelerometer data, compute the FRF, and average them for the 10 impacts
+- Maybe have to take into account the sensitivity, etc...
+
+#+begin_src matlab
+raw_data_path = '/home/thomas/Cloud/work-projects/ID31-NASS/matlab/nass-measurements/micro-station-compliance/data/record/';
+data = [
+    load(sprintf('%s/Measurement1.mat', raw_data_path)), ...
+    load(sprintf('%s/Measurement2.mat', raw_data_path)), ...
+    load(sprintf('%s/Measurement3.mat', raw_data_path)), ...
+    load(sprintf('%s/Measurement4.mat', raw_data_path)), ...
+    load(sprintf('%s/Measurement5.mat', raw_data_path)), ...
+    load(sprintf('%s/Measurement6.mat', raw_data_path)), ...
+    load(sprintf('%s/Measurement7.mat', raw_data_path)), ...
+    load(sprintf('%s/Measurement8.mat', raw_data_path)), ...
+    load(sprintf('%s/Measurement9.mat', raw_data_path)), ...
+    load(sprintf('%s/Measurement10.mat', raw_data_path))];
+#+end_src
+
+#+begin_src matlab
+Ts = 1e-3; % Sampling Time [s]
+Nfft = floor(1/Ts); % Number of points for the FFT computation
+% win = hanning(Nfft); % Hanning window
+win = ones(Nfft, 1); % Rectangular window
+% Get the frequency vector
+[~, f] = tfestimate(data(1).Track13, data(1).Track1, win, [], Nfft, 1/Ts);
+#+end_src
+
+#+begin_src matlab
+pre_n = floor(0.1/Ts);
+post_n = Nfft - pre_n - 1;
+#+end_src
+
+#+begin_src matlab
+G_raw = zeros(12,10,length(f));
+for i = 1:10
+    % Find the impacts
+    [~, impacts_i] = find(diff(data(i).Track13 > 50)==1);
+    % Only keep the first 10 impacts if there are more than 10 impacts
+    if length(impacts_i)>10
+        impacts_i(11:end) = [];
+    end
+    % Initialize the FRF for the current experiment
+    G_impact = zeros(12,length(f));
+    for impact_i = impacts_i
+        i_start = impacts_i - pre_n;
+        i_end = impacts_i + post_n;
+        data_in = data(i).Track13(i_start:i_end); % [N]
+        % Remove hammer DC offset
+        data_in  = data_in  - mean(data_in(end-pre_n:end));
+        % Combine all outputs [m/s^2]
+        data_out = [data(i).Track1( i_start:i_end); ...
+                    data(i).Track2( i_start:i_end); ...
+                    data(i).Track3( i_start:i_end); ...
+                    data(i).Track4( i_start:i_end); ...
+                    data(i).Track5( i_start:i_end); ...
+                    data(i).Track6( i_start:i_end); ...
+                    data(i).Track7( i_start:i_end); ...
+                    data(i).Track8( i_start:i_end); ...
+                    data(i).Track9( i_start:i_end); ...
+                    data(i).Track10(i_start:i_end); ...
+                    data(i).Track11(i_start:i_end); ...
+                    data(i).Track12(i_start:i_end)];
+
+        [frf, ~] = tfestimate(data_in, data_out', win, [], Nfft, 1/Ts);
+        G_raw(:,i,:) = frf'./(-(2*pi*f').^2);
+    end
+end
+#+end_src
+
+#+begin_src matlab
+%% Compute transfer function in cartesian frame
+G_compl_b = pagemtimes(J_inv, pagemtimes(G_raw, Jf_inv));
+#+end_src
+
+#+begin_src matlab
+colors = colororder;
+figure;
+hold on;
+% plot(freqs, abs(squeeze(G_compl(1,1,:))), 'color', colors(1,:), 'DisplayName', '$C_x/F_x$')
+% plot(freqs, abs(squeeze(G_compl(2,2,:))), 'color', colors(2,:), 'DisplayName', '$C_y/F_y$')
+% plot(freqs, abs(squeeze(G_compl(3,3,:))), 'color', colors(3,:), 'DisplayName', '$C_z/F_z$')
+plot(f, abs(squeeze(G_compl_b(1,1,:))), '-', 'color', colors(1,:), 'DisplayName', '$C_x/F_x$')
+plot(f, abs(squeeze(G_compl_b(2,2,:))), '-', 'color', colors(2,:), 'DisplayName', '$C_y/F_y$')
+plot(f, abs(squeeze(G_compl_b(3,3,:))), '-', 'color', colors(3,:), 'DisplayName', '$C_z/F_z$')
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('Frequency [Hz]'); ylabel('Magnitude [m/N]');
+xlim([30, 300]); ylim([1e-9, 2e-6]);
+leg = legend('location', 'southwest', 'FontSize', 8, 'NumColumns', 1);
+leg.ItemTokenSize(1) = 15;
+#+end_src
+
+*** Load Data
+#+begin_src matlab
+% Load data
+data = [load('data/Measurement1.mat'), ...
+        load('data/Measurement2.mat'), ...
+        load('data/Measurement3.mat'), ...
+        load('data/Measurement4.mat'), ...
+        load('data/Measurement5.mat'), ...
+        load('data/Measurement6.mat'), ...
+        load('data/Measurement7.mat'), ...
+        load('data/Measurement8.mat'), ...
+        load('data/Measurement9.mat'), ...
+        load('data/Measurement10.mat')];
+#+end_src
+
+#+begin_src matlab
+% Frequency Vector
+freqs = m3.FFT1_H1_1_13_X_Val;
+w = 2*pi*freqs;
+
+% 12 outputs, 10 inputs
+G_raw = zeros(12, 10, length(freqs));
+
+for j = 1:10
+    for i = 1:12
+        G_raw(i,j,:) = data(j).(sprintf("FFT1_H1_%i_13_Y_ReIm", i))./(-w.^2);
+    end
+end
+#+end_src
+
+#+begin_src matlab
+%% Compute transfer function in cartesian frame
+G_compl = pagemtimes(J_inv, pagemtimes(G_raw, Jf_inv));
+#+end_src
+
+#+begin_src matlab
+colors = colororder;
+figure;
+hold on;
+plot(freqs, abs(squeeze(G_compl(1,1,:))), 'color', colors(1,:), 'DisplayName', '$C_x/F_x$')
+plot(freqs, abs(squeeze(G_compl(2,2,:))), 'color', colors(2,:), 'DisplayName', '$C_y/F_y$')
+plot(freqs, abs(squeeze(G_compl(3,3,:))), 'color', colors(3,:), 'DisplayName', '$C_z/F_z$')
+plot(f, abs(squeeze(G_compl_b(1,1,:))), '.', 'color', colors(1,:), 'DisplayName', '$C_x/F_x$')
+plot(f, abs(squeeze(G_compl_b(2,2,:))), '.', 'color', colors(2,:), 'DisplayName', '$C_y/F_y$')
+plot(f, abs(squeeze(G_compl_b(3,3,:))), '.', 'color', colors(3,:), 'DisplayName', '$C_z/F_z$')
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('Frequency [Hz]'); ylabel('Magnitude [m/N]');
+xlim([30, 300]); ylim([1e-9, 2e-6]);
+leg = legend('location', 'southwest', 'FontSize', 8, 'NumColumns', 1);
+leg.ItemTokenSize(1) = 15;
+#+end_src
+
+#+begin_src matlab
+colors = colororder;
+figure;
+hold on;
+plot(freqs, abs(squeeze(G_compl(4,4,:))), 'color', colors(1,:), 'DisplayName', '$R_x/M_x$')
+plot(freqs, abs(squeeze(G_compl(5,5,:))), 'color', colors(2,:), 'DisplayName', '$R_y/M_y$')
+plot(freqs, abs(squeeze(G_compl(6,6,:))), 'color', colors(3,:), 'DisplayName', '$R_z/M_z$')
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('Frequency [Hz]'); ylabel('Magnitude [m/N]');
+xlim([30, 300]); ylim([5e-7, 5e-5]);
+leg = legend('location', 'northeast', 'FontSize', 8, 'NumColumns', 1);
+leg.ItemTokenSize(1) = 15;
+#+end_src
+
+*** Diagonal Dynamics
+#+begin_src matlab
+figure;
+ax1 = subplot(2,1,1);
+hold on;
+plot(freqs, abs(squeeze(G(1,1,:)), '-', 'DisplayName', '$D_x/F_x$')
+plot(freqs, abs(squeeze(G(2,2,:)), '-', 'DisplayName', '$D_y/F_y$')
+plot(freqs, abs(squeeze(G(3,3,:)), '-', 'DisplayName', '$D_z/F_z$')
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('Frequency [Hz]'); ylabel('Magnitude [m/N]');
+ylim([1e-9, 2e-6]);
+legend('location', 'southwest');
+xlim([30, 300]);
+#+end_src
+
+#+begin_src matlab
+figure;
+ax1 = subplot(2,1,1);
+hold on;
+plot(freqs, abs(squeeze(G(4,4,:)), '-', 'DisplayName', '$R_x/M_x$')
+plot(freqs, abs(squeeze(G(5,5,:)), '-', 'DisplayName', '$R_y/M_y$')
+plot(freqs, abs(squeeze(G(6,6,:)), '-', 'DisplayName', '$R_z/M_z$')
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('Frequency [Hz]'); ylabel('Magnitude [rad/Nm]');
+ylim([1e-7, 2e-4]);
+legend('location', 'southwest');
+xlim([30, 300]);
+#+end_src
+
+*** Equivalent Stiffness and Mass Estimation
+
+#+begin_src matlab
+K = [1e7, 1e7, 2e8, 5e7, 3e7, 2e7];
+f_res = [125, 135, 390, 335, 335, 160];
+#+end_src
+
+#+begin_src matlab
+M = [20, 20, 20, 11, 7, 20];
+f_res_est = sqrt(K./M)./(2*pi);
+#+end_src
+
+Here is the inertia / stiffness to the granite that can represent the micro-station compliance dynamics:
+#+begin_src matlab :exports results :results value table replace :tangle no :post addhdr(*this*)
+data2orgtable([K'], {'x', 'y', 'z', 'Rx', 'Ry', 'Rz'}, {'Stiffness', 'Inertia'}, ' %.1g ');
+#+end_src
+
+
+#+RESULTS:
+| Stiffness |     Inertia |
+|-----------+-------------|
+| x         |  10000000.0 |
+| y         |  10000000.0 |
+| z         | 200000000.0 |
+| Rx        |  50000000.0 |
+| Ry        |  30000000.0 |
+| Rz        |  20000000.0 |
+
+** Compare with the Model
 #+begin_src matlab
 %% Initialize simulation with default parameters (flexible elements)
 initializeGround();
@@ -946,6 +1354,176 @@ legend('location', 'northwest');
 
 - [ ] Comparison with the estimated (or measured) compliance
 
+#+begin_src matlab
+  figure;
+  ax1 = subplot(2,1,1);
+  hold on;
+  plot(freqs, abs(squeeze(G(1,1,:))./(-w.^2)), '.')
+  plot(freqs, abs(squeeze(G(2,2,:))./(-w.^2)), '.')
+  plot(freqs, abs(squeeze(G(3,3,:))./(-w.^2)), '.')
+  set(gca,'ColorOrderIndex',1);
+  plot(freqs, abs(squeeze(freqresp(Gm(1,1,:), freqs, 'Hz'))), '-')
+  plot(freqs, abs(squeeze(freqresp(Gm(2,2,:), freqs, 'Hz'))), '-')
+  plot(freqs, abs(squeeze(freqresp(Gm(3,3,:), freqs, 'Hz'))), '-')
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('Magnitude [m/N]'); set(gca, 'XTickLabel',[]);
+  ylim([1e-9, 2e-6]);
+
+  ax2 = subplot(2,1,2);
+  hold on;
+  plot(freqs, 180/pi*angle(squeeze(G(1,1,:))./(-w.^2)), '.', 'DisplayName', '$D_x/F_x$')
+  plot(freqs, 180/pi*angle(squeeze(G(2,2,:))./(-w.^2)), '.', 'DisplayName', '$D_y/F_y$')
+  plot(freqs, 180/pi*angle(squeeze(G(3,3,:))./(-w.^2)), '.', 'DisplayName', '$D_z/F_z$')
+  set(gca,'ColorOrderIndex',1);
+  plot(freqs, 180/pi*angle(squeeze(freqresp(Gm(1,1,:), freqs, 'Hz'))), '-', 'HandleVisibility', 'off')
+  plot(freqs, 180/pi*angle(squeeze(freqresp(Gm(2,2,:), freqs, 'Hz'))), '-', 'HandleVisibility', 'off')
+  plot(freqs, 180/pi*angle(squeeze(freqresp(Gm(3,3,:), freqs, 'Hz'))), '-', 'HandleVisibility', 'off')
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
+  xlabel('Freqency [Hz]'); ylabel('Phase [deg]');
+  ylim([-180, 180]);
+  yticks([-180, -90, 0, 90, 180]);
+  legend('location', 'southwest');
+
+  linkaxes([ax1,ax2],'x');
+  xlim([30, 300]);
+#+end_src
+
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/compliance_diagonal_translations_comp_model.pdf', 'width', 'full', 'height', 'full');
+#+end_src
+
+#+name: fig:compliance_diagonal_translations_comp_model
+#+caption: Dynamics from Forces to Translations
+#+RESULTS:
+[[file:figs/compliance_diagonal_translations_comp_model.png]]
+
+#+begin_src matlab
+  figure;
+  ax1 = subplot(2,1,1);
+  hold on;
+  plot(freqs, abs(squeeze(G(4,4,:))./(-w.^2)), '.')
+  plot(freqs, abs(squeeze(G(5,5,:))./(-w.^2)), '.')
+  plot(freqs, abs(squeeze(G(6,6,:))./(-w.^2)), '.')
+  set(gca,'ColorOrderIndex',1);
+  plot(freqs, abs(squeeze(freqresp(Gm(4,4,:), freqs, 'Hz'))), '-')
+  plot(freqs, abs(squeeze(freqresp(Gm(5,5,:), freqs, 'Hz'))), '-')
+  plot(freqs, abs(squeeze(freqresp(Gm(6,6,:), freqs, 'Hz'))), '-')
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('Magnitude [rad/Nm]'); set(gca, 'XTickLabel',[]);
+  % ylim([1e-9, 2e-6]);
+
+  ax2 = subplot(2,1,2);
+  hold on;
+  plot(freqs, 180/pi*angle(squeeze(G(4,4,:))./(-w.^2)), '.', 'DisplayName', '$R_x/M_x$')
+  plot(freqs, 180/pi*angle(squeeze(G(5,5,:))./(-w.^2)), '.', 'DisplayName', '$R_y/M_y$')
+  plot(freqs, 180/pi*angle(squeeze(G(6,6,:))./(-w.^2)), '.', 'DisplayName', '$R_z/M_z$')
+  set(gca,'ColorOrderIndex',1);
+  plot(freqs, 180/pi*angle(squeeze(freqresp(Gm(4,4,:), freqs, 'Hz'))), '-', 'HandleVisibility', 'off')
+  plot(freqs, 180/pi*angle(squeeze(freqresp(Gm(5,5,:), freqs, 'Hz'))), '-', 'HandleVisibility', 'off')
+  plot(freqs, 180/pi*angle(squeeze(freqresp(Gm(6,6,:), freqs, 'Hz'))), '-', 'HandleVisibility', 'off')
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
+  xlabel('Freqency [Hz]'); ylabel('Phase [deg]');
+  ylim([-180, 180]);
+  yticks([-180, -90, 0, 90, 180]);
+  legend('location', 'southwest');
+
+  linkaxes([ax1,ax2],'x');
+  xlim([30, 300]);
+#+end_src
+
+#+begin_src matlab :tangle no :exports results :results file replace
+  exportFig('figs/compliance_diagonal_rotations_comp_model.pdf', 'width', 'full', 'height', 'full');
+#+end_src
+
+#+name: fig:compliance_diagonal_rotations_comp_model
+#+caption: Dynamics from Torques to Rotations
+#+RESULTS:
+[[file:figs/compliance_diagonal_rotations_comp_model.png]]
+
+|           | Stiffness | Unit     |
+|-----------+-----------+----------|
+| $K_x$     |       1e7 | [N/m]    |
+| $K_y$     |       1e7 | [N/m]    |
+| $K_z$     |       2e8 | [N/m]    |
+| $K_{R_x}$ |       5e7 | [Nm/rad] |
+| $K_{R_y}$ |       3e7 | [Nm/rad] |
+| $K_{R_z}$ |       2e7 | [Nm/rad] |
+
+*** Coupling Dynamics
+#+begin_src matlab
+  figure;
+  ax1 = subplot(2,1,1);
+  hold on;
+  plot(freqs, abs(squeeze(G(1,1,:))./(-w.^2)), '.')
+  plot(freqs, abs(squeeze(G(2,1,:))./(-w.^2)), '.')
+  plot(freqs, abs(squeeze(G(3,1,:))./(-w.^2)), '.')
+  set(gca,'ColorOrderIndex',1);
+  plot(freqs, abs(squeeze(freqresp(Gm(1,1,:), freqs, 'Hz'))), '-')
+  plot(freqs, abs(squeeze(freqresp(Gm(2,1,:), freqs, 'Hz'))), '-')
+  plot(freqs, abs(squeeze(freqresp(Gm(3,1,:), freqs, 'Hz'))), '-')
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('Magnitude [m/N]'); set(gca, 'XTickLabel',[]);
+  ylim([1e-9, 2e-6]);
+
+  ax2 = subplot(2,1,2);
+  hold on;
+  plot(freqs, 180/pi*angle(squeeze(G(1,1,:))./(-w.^2)), '.', 'DisplayName', '$D_x/F_x$')
+  plot(freqs, 180/pi*angle(squeeze(G(2,1,:))./(-w.^2)), '.', 'DisplayName', '$D_y/F_x$')
+  plot(freqs, 180/pi*angle(squeeze(G(3,1,:))./(-w.^2)), '.', 'DisplayName', '$D_z/F_x$')
+  set(gca,'ColorOrderIndex',1);
+  plot(freqs, 180/pi*angle(squeeze(freqresp(Gm(1,1,:), freqs, 'Hz'))), '-', 'HandleVisibility', 'off')
+  plot(freqs, 180/pi*angle(squeeze(freqresp(Gm(2,1,:), freqs, 'Hz'))), '-', 'HandleVisibility', 'off')
+  plot(freqs, 180/pi*angle(squeeze(freqresp(Gm(3,1,:), freqs, 'Hz'))), '-', 'HandleVisibility', 'off')
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
+  xlabel('Freqency [Hz]'); ylabel('Phase [deg]');
+  ylim([-180, 180]);
+  yticks([-180, -90, 0, 90, 180]);
+  legend('location', 'southwest');
+
+  linkaxes([ax1,ax2],'x');
+  xlim([30, 300]);
+#+end_src
+
+#+begin_src matlab
+  figure;
+  ax1 = subplot(2,1,1);
+  hold on;
+  plot(freqs, abs(squeeze(G(5,1,:))./(-w.^2)), '.')
+  plot(freqs, abs(squeeze(G(4,2,:))./(-w.^2)), '.')
+  set(gca,'ColorOrderIndex',1);
+  plot(freqs, abs(squeeze(freqresp(Gm(5,1,:), freqs, 'Hz'))), '-')
+  plot(freqs, abs(squeeze(freqresp(Gm(4,2,:), freqs, 'Hz'))), '-')
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+  ylabel('Magnitude [m/N]'); set(gca, 'XTickLabel',[]);
+  ylim([1e-9, 2e-6]);
+
+  ax2 = subplot(2,1,2);
+  hold on;
+  plot(freqs, 180/pi*angle(squeeze(G(5,1,:))./(-w.^2)), '.', 'DisplayName', '$R_y/F_x$')
+  plot(freqs, 180/pi*angle(squeeze(G(4,2,:))./(-w.^2)), '.', 'DisplayName', '$R_x/F_y$')
+  set(gca,'ColorOrderIndex',1);
+  plot(freqs, 180/pi*angle(squeeze(freqresp(Gm(5,1,:), freqs, 'Hz'))), '-', 'HandleVisibility', 'off')
+  plot(freqs, 180/pi*angle(squeeze(freqresp(Gm(4,2,:), freqs, 'Hz'))), '-', 'HandleVisibility', 'off')
+  hold off;
+  set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
+  xlabel('Freqency [Hz]'); ylabel('Phase [deg]');
+  ylim([-180, 180]);
+  yticks([-180, -90, 0, 90, 180]);
+  legend('location', 'southwest');
+
+  linkaxes([ax1,ax2],'x');
+  xlim([30, 300]);
+#+end_src
+
+
+
 ** Conclusion
 For such a complex system, we believe that the Simscape Model represents the dynamics of the system with enough fidelity.