2413 lines
69 KiB
Org Mode
2413 lines
69 KiB
Org Mode
#+TITLE: ESRF Double Crystal Monochromator - Feedback Controller
|
|
:DRAWER:
|
|
#+LANGUAGE: en
|
|
#+EMAIL: dehaeze.thomas@gmail.com
|
|
#+AUTHOR: Dehaeze Thomas
|
|
|
|
#+HTML_LINK_HOME: ../index.html
|
|
#+HTML_LINK_UP: ../index.html
|
|
|
|
#+HTML_HEAD: <link rel="stylesheet" type="text/css" href="https://research.tdehaeze.xyz/css/style.css"/>
|
|
#+HTML_HEAD: <script type="text/javascript" src="https://research.tdehaeze.xyz/js/script.js"></script>
|
|
|
|
#+BIND: org-latex-image-default-option "scale=1"
|
|
#+BIND: org-latex-image-default-width ""
|
|
|
|
#+LaTeX_CLASS: scrreprt
|
|
#+LaTeX_CLASS_OPTIONS: [a4paper, 10pt, DIV=12, parskip=full]
|
|
#+LaTeX_HEADER_EXTRA: \input{preamble.tex}
|
|
#+LATEX_HEADER_EXTRA: \bibliography{ref}
|
|
|
|
#+PROPERTY: header-args:matlab :session *MATLAB*
|
|
#+PROPERTY: header-args:matlab+ :comments org
|
|
#+PROPERTY: header-args:matlab+ :exports both
|
|
#+PROPERTY: header-args:matlab+ :results none
|
|
#+PROPERTY: header-args:matlab+ :tangle no
|
|
#+PROPERTY: header-args:matlab+ :eval no-export
|
|
#+PROPERTY: header-args:matlab+ :noweb yes
|
|
#+PROPERTY: header-args:matlab+ :mkdirp yes
|
|
#+PROPERTY: header-args:matlab+ :output-dir figs
|
|
|
|
#+PROPERTY: header-args:latex :headers '("\\usepackage{tikz}" "\\usepackage{import}" "\\import{$HOME/Cloud/tikz/org/}{config.tex}")
|
|
#+PROPERTY: header-args:latex+ :imagemagick t :fit yes
|
|
#+PROPERTY: header-args:latex+ :iminoptions -scale 100% -density 150
|
|
#+PROPERTY: header-args:latex+ :imoutoptions -quality 100
|
|
#+PROPERTY: header-args:latex+ :results file raw replace
|
|
#+PROPERTY: header-args:latex+ :buffer no
|
|
#+PROPERTY: header-args:latex+ :tangle no
|
|
#+PROPERTY: header-args:latex+ :eval no-export
|
|
#+PROPERTY: header-args:latex+ :exports results
|
|
#+PROPERTY: header-args:latex+ :mkdirp yes
|
|
#+PROPERTY: header-args:latex+ :output-dir figs
|
|
#+PROPERTY: header-args:latex+ :post pdf2svg(file=*this*, ext="png")
|
|
:END:
|
|
|
|
#+begin_export html
|
|
<hr>
|
|
<p>This report is also available as a <a href="./dcm-feedback-control.pdf">pdf</a>.</p>
|
|
<hr>
|
|
#+end_export
|
|
|
|
#+latex: \clearpage
|
|
|
|
* Introduction :ignore:
|
|
|
|
* Estimation of Sensitivity Function
|
|
** 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 :tangle no :noweb yes
|
|
<<m-init-path>>
|
|
#+end_src
|
|
|
|
#+begin_src matlab :eval no :noweb yes
|
|
<<m-init-path-tangle>>
|
|
#+end_src
|
|
|
|
#+begin_src matlab :noweb yes
|
|
<<m-init-other>>
|
|
#+end_src
|
|
|
|
** Load Data
|
|
Two scans are performed:
|
|
- =1.1= in mode B
|
|
- =3.1= in mode C
|
|
|
|
The difference between the two is that mode C adds the feedback controller.
|
|
|
|
#+begin_src matlab
|
|
%% Load Data of the new LUT method
|
|
Ts = 0.1;
|
|
|
|
ol_drx = 1e-9*double(h5read('xanes_0003.h5','/1.1/measurement/xtal_111_drx_filter')); % Rx [rad]
|
|
cl_drx = 1e-9*double(h5read('xanes_0003.h5','/3.1/measurement/xtal_111_drx_filter')); % Rx [rad]
|
|
|
|
ol_dry = 1e-9*double(h5read('xanes_0003.h5','/1.1/measurement/xtal_111_dry_filter')); % Ry [rad]
|
|
cl_dry = 1e-9*double(h5read('xanes_0003.h5','/3.1/measurement/xtal_111_dry_filter')); % Ry [rad]
|
|
|
|
t = linspace(Ts, Ts*length(ol_drx), length(ol_drx));
|
|
#+end_src
|
|
|
|
By comparison the frequency content of the crystal orientation errors between mode B and mode C, it is possible to estimate the Sensitivity transfer function (Figure [[fig:sensitivity_function_drx_est]]).
|
|
#+begin_src matlab
|
|
win = hanning(ceil(1/Ts));
|
|
|
|
[pxx_ol_drx, f] = pwelch(ol_drx, win, [], [], 1/Ts);
|
|
[pxx_cl_drx, ~] = pwelch(cl_drx, win, [], [], 1/Ts);
|
|
|
|
[pxx_ol_dry, ~] = pwelch(ol_dry, win, [], [], 1/Ts);
|
|
[pxx_cl_dry, ~] = pwelch(cl_dry, win, [], [], 1/Ts);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Estimation of Sensitivity function
|
|
figure;
|
|
hold on;
|
|
plot(f, sqrt(pxx_cl_drx)./sqrt(pxx_ol_drx), 'DisplayName', '$R_x$');
|
|
plot(f, sqrt(pxx_cl_dry)./sqrt(pxx_ol_dry), 'DisplayName', '$R_y$');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('Sensitivity Function');
|
|
legend('location', 'northwest');
|
|
#+end_src
|
|
|
|
#+begin_src matlab :tangle no :exports results :results file replace
|
|
exportFig('figs/sensitivity_function_drx_est.pdf', 'width', 'wide', 'height', 'normal');
|
|
#+end_src
|
|
|
|
#+name: fig:sensitivity_function_drx_est
|
|
#+caption: Estimation of the sensitivity transfer function magnitude
|
|
#+RESULTS:
|
|
[[file:figs/sensitivity_function_drx_est.png]]
|
|
|
|
|
|
** Controller
|
|
#+begin_src matlab
|
|
load('X_tal_cage_PID.mat', 'K');
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Bode plot for the controller
|
|
figure;
|
|
tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
|
|
|
|
ax1 = nexttile([2,1]);
|
|
hold on;
|
|
plot(freqs, abs(squeeze(freqresp(K(1,1), freqs, 'Hz'))));
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
ylabel('Amplitude [-]'); set(gca, 'XTickLabel',[]);
|
|
|
|
ax2 = nexttile;
|
|
hold on;
|
|
plot(freqs, 180/pi*angle(squeeze(freqresp(K(1,1), freqs, 'Hz'))));
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
hold off;
|
|
yticks(-360:90:360);
|
|
ylim([-180, 180]);
|
|
|
|
linkaxes([ax1,ax2],'x');
|
|
xlim([freqs(1), freqs(end)]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :tangle no :exports results :results file replace
|
|
exportFig('figs/bode_plot_cur_controller.pdf', 'width', 'wide', 'height', 'tall');
|
|
#+end_src
|
|
|
|
#+name: fig:bode_plot_cur_controller
|
|
#+caption: Bode Plot of the Controller
|
|
#+RESULTS:
|
|
[[file:figs/bode_plot_cur_controller.png]]
|
|
|
|
** Test
|
|
#+begin_src matlab
|
|
Ts = 5e-3;
|
|
cl_drx = 1e-9*double(h5read('xanes_0003.h5','/16.1/measurement/xtal_111_drx_filter')); % Rx [rad]
|
|
ol_drx = 1e-9*double(h5read('xanes_0003.h5','/18.1/measurement/xtal_111_drx_filter')); % Rx [rad]
|
|
|
|
t = linspace(Ts, Ts*length(ol_drx), length(ol_drx));
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
figure;
|
|
hold on;
|
|
plot(t, ol_drx)
|
|
plot(t, cl_drx)
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
win = hanning(ceil(10/Ts));
|
|
|
|
[pxx_ol_drx, f] = pwelch(ol_drx, win, [], [], 1/Ts);
|
|
[pxx_cl_drx, ~] = pwelch(cl_drx, win, [], [], 1/Ts);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Estimation of Sensitivity function
|
|
figure;
|
|
hold on;
|
|
plot(f, sqrt(pxx_cl_drx)./sqrt(pxx_ol_drx), 'DisplayName', '$R_x$');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('Sensitivity Function');
|
|
legend('location', 'northwest');
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Estimation of Sensitivity function
|
|
figure;
|
|
hold on;
|
|
plot(f, sqrt(pxx_ol_drx), 'DisplayName', '$R_x$ - OL');
|
|
plot(f, sqrt(pxx_cl_drx), 'DisplayName', '$R_x$ - CL');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('Amplitude Spectral Density [$m/\sqrt{Hz}$]');
|
|
legend('location', 'northwest');
|
|
#+end_src
|
|
|
|
* System Identification
|
|
** 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 :tangle no :noweb yes
|
|
<<m-init-path>>
|
|
#+end_src
|
|
|
|
#+begin_src matlab :eval no :noweb yes
|
|
<<m-init-path-tangle>>
|
|
#+end_src
|
|
|
|
#+begin_src matlab :noweb yes
|
|
<<m-init-other>>
|
|
#+end_src
|
|
|
|
** Identification ID24
|
|
#+begin_src matlab
|
|
load('test_id_id24_3.mat')
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
t = 1e-4*ones(size(fjpur, 1), 1);
|
|
|
|
ur.dz = fjpur(:,1) - mean(fjpur(:,1));
|
|
ur.dry = fjpur(:,2) - mean(fjpur(:,2));
|
|
ur.drx = fjpur(:,3) - mean(fjpur(:,3));
|
|
ur.u = fjpur(:,7) - mean(fjpur(:,7));
|
|
|
|
uh.dz = fjpuh(:,1) - mean(fjpuh(:,1));
|
|
uh.dry = fjpuh(:,2) - mean(fjpuh(:,2));
|
|
uh.drx = fjpuh(:,3) - mean(fjpuh(:,3));
|
|
uh.u = fjpuh(:,8) - mean(fjpuh(:,8));
|
|
|
|
d.dz = fjpd(:,1) - mean(fjpd(:,1));
|
|
d.dry = fjpd(:,2) - mean(fjpd(:,2));
|
|
d.drx = fjpd(:,3) - mean(fjpd(:,3));
|
|
d.u = fjpd(:,9) - mean(fjpd(:,9));
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
J_a_311 = [1, 0.14, -0.0675
|
|
1, 0.14, 0.1525
|
|
1, -0.14, 0.0425];
|
|
|
|
J_a_111 = [1, 0.14, -0.1525
|
|
1, 0.14, 0.0675
|
|
1, -0.14, -0.0425];
|
|
|
|
ur.y = [J_a_311 * [-ur.dz, ur.dry,-ur.drx]']';
|
|
uh.y = [J_a_311 * [-uh.dz, uh.dry,-uh.drx]']';
|
|
d.y = [J_a_311 * [-d.dz, d.dry, -d.drx]']';
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Sampling Time and Frequency
|
|
Ts = 1e-4; % [s]
|
|
Fs = 1/Ts; % [Hz]
|
|
|
|
% Hannning Windows
|
|
win = hanning(ceil(1*Fs));
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% And we get the frequency vector
|
|
[G_ur, f] = tfestimate(ur.u, ur.y, win, [], [], 1/Ts);
|
|
[G_uh, ~] = tfestimate(uh.u, uh.y, win, [], [], 1/Ts);
|
|
[G_d, ~] = tfestimate(d.u, d.y, win, [], [], 1/Ts);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
[coh_ur, ~] = mscohere(ur.u, ur.y, win, [], [], 1/Ts);
|
|
[coh_uh, ~] = mscohere(uh.u, uh.y, win, [], [], 1/Ts);
|
|
[coh_d, ~] = mscohere(d.u, d.y, win, [], [], 1/Ts);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Coherence
|
|
figure;
|
|
tiledlayout(1, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
|
|
hold on;
|
|
plot(f, coh_ur(:,1), 'DisplayName', '$u_r$');
|
|
plot(f, coh_uh(:,2), 'DisplayName', '$u_h$');
|
|
plot(f, coh_d( :,3), 'DisplayName', '$u_d$');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Coherence');
|
|
legend('location', 'southeast');
|
|
xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Bode plot of the DCM dynamics in the frame of the Fast Jacks
|
|
figure;
|
|
tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
|
|
|
|
ax1 = nexttile([2,1]);
|
|
hold on;
|
|
plot(f, abs(G_ur(:,1)), 'DisplayName', '$u_r$');
|
|
plot(f, abs(G_uh(:,2)), 'DisplayName', '$u_h$');
|
|
plot(f, abs(G_d( :,3)), 'DisplayName', '$u_d$');
|
|
plot(f, abs(G_ur(:,2)), 'color', [0, 0, 0, 0.2], ...
|
|
'DisplayName', 'Off Diagonal');
|
|
plot(f, abs(G_ur(:,3)), 'color', [0, 0, 0, 0.2], ...
|
|
'HandleVisibility', 'off');
|
|
plot(f, abs(G_uh(:,1)), 'color', [0, 0, 0, 0.2], ...
|
|
'HandleVisibility', 'off');
|
|
plot(f, abs(G_uh(:,3)), 'color', [0, 0, 0, 0.2], ...
|
|
'HandleVisibility', 'off');
|
|
plot(f, abs(G_d(:,1)), 'color', [0, 0, 0, 0.2], ...
|
|
'HandleVisibility', 'off');
|
|
plot(f, abs(G_d(:,2)), 'color', [0, 0, 0, 0.2], ...
|
|
'HandleVisibility', 'off');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
|
|
ylim([1e-2, 1e2]);
|
|
legend('location', 'northwest');
|
|
|
|
ax2 = nexttile;
|
|
hold on;
|
|
plot(f, 180/pi*angle(G_ur(:,1)));
|
|
plot(f, 180/pi*angle(G_uh(:,2)));
|
|
plot(f, 180/pi*angle(G_d(:,3)));
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
hold off;
|
|
yticks(-360:90:360);
|
|
% ylim([-45, 45]);
|
|
|
|
linkaxes([ax1,ax2],'x');
|
|
xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
** Identification
|
|
#+begin_src matlab
|
|
ur = load('FJPUR_step.mat');
|
|
uh = load('FJPUH_step.mat');
|
|
d = load('FJPD_step.mat');
|
|
#+end_src
|
|
|
|
| 1 | dz311 |
|
|
| 2 | dry311 |
|
|
| 3 | drx311 |
|
|
| 4 | dz111 |
|
|
| 5 | dry111 |
|
|
| 6 | drx111 |
|
|
| 7 | fjpur |
|
|
| 8 | fjpuh |
|
|
| 9 | fjpd |
|
|
| 10 | bragg |
|
|
|
|
#+begin_src matlab
|
|
ur.time = ur.time - ur.time(1);
|
|
ur.allValues(:, 1) = ur.allValues(:, 1) - mean(ur.allValues(ur.time<1, 1));
|
|
ur.allValues(:, 2) = ur.allValues(:, 2) - mean(ur.allValues(ur.time<1, 2));
|
|
ur.allValues(:, 3) = ur.allValues(:, 3) - mean(ur.allValues(ur.time<1, 3));
|
|
|
|
t_filt = ur.time > 48 & ur.time < 60;
|
|
|
|
ur.u = ur.allValues(t_filt, 7);
|
|
ur.y_111 = [-ur.allValues(t_filt, 1), ur.allValues(t_filt, 2), ur.allValues(t_filt, 3)];
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
uh.time = uh.time - uh.time(1);
|
|
uh.allValues(:, 1) = uh.allValues(:, 1) - mean(uh.allValues(uh.time<1, 1));
|
|
uh.allValues(:, 2) = uh.allValues(:, 2) - mean(uh.allValues(uh.time<1, 2));
|
|
uh.allValues(:, 3) = uh.allValues(:, 3) - mean(uh.allValues(uh.time<1, 3));
|
|
|
|
uh.u = uh.allValues(t_filt, 8);
|
|
uh.y_111 = [-uh.allValues(t_filt, 1), uh.allValues(t_filt, 2), uh.allValues(t_filt, 3)];
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
d.time = d.time - d.time(1);
|
|
d.allValues(:, 1) = d.allValues(:, 1) - mean(d.allValues(d.time<1, 1));
|
|
d.allValues(:, 2) = d.allValues(:, 2) - mean(d.allValues(d.time<1, 2));
|
|
d.allValues(:, 3) = d.allValues(:, 3) - mean(d.allValues(d.time<1, 3));
|
|
|
|
d.u = d.allValues(t_filt, 9);
|
|
d.y_111 = [-d.allValues(t_filt, 1), d.allValues(t_filt, 2), d.allValues(t_filt, 3)];
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
J_a_111 = [1, 0.14, -0.0675
|
|
1, 0.14, 0.1525
|
|
1, -0.14, 0.0425];
|
|
|
|
J_a_311 = [1, 0.14, -0.1525
|
|
1, 0.14, 0.0675
|
|
1, -0.14, -0.0425];
|
|
|
|
ur.y = [J_a_311 * ur.y_111']';
|
|
uh.y = [J_a_311 * uh.y_111']';
|
|
d.y = [J_a_311 * d.y_111']';
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Sampling Time and Frequency
|
|
Ts = 1e-4; % [s]
|
|
Fs = 1/Ts; % [Hz]
|
|
|
|
% Hannning Windows
|
|
win = hanning(ceil(5*Fs));
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% And we get the frequency vector
|
|
[G_ur, f] = tfestimate(ur.u, ur.y, win, [], [], 1/Ts);
|
|
[G_uh, ~] = tfestimate(uh.u, uh.y, win, [], [], 1/Ts);
|
|
[G_d, ~] = tfestimate(d.u, d.y, win, [], [], 1/Ts);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
[coh_ur, ~] = mscohere(ur.u, ur.y, win, [], [], 1/Ts);
|
|
[coh_uh, ~] = mscohere(uh.u, uh.y, win, [], [], 1/Ts);
|
|
[coh_d, ~] = mscohere(d.u, d.y, win, [], [], 1/Ts);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Coherence
|
|
figure;
|
|
tiledlayout(1, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
|
|
hold on;
|
|
plot(f, coh_ur(:,1), 'DisplayName', '$u_r$');
|
|
plot(f, coh_uh(:,2), 'DisplayName', '$u_h$');
|
|
plot(f, coh_d( :,3), 'DisplayName', '$u_d$');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Coherence');
|
|
legend('location', 'southeast');
|
|
xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :tangle no :exports results :results file replace
|
|
exportFig('figs/coherence_id_dcm_dyn.pdf', 'width', 'wide', 'height', 'normal');
|
|
#+end_src
|
|
|
|
#+name: fig:coherence_id_dcm_dyn
|
|
#+caption: Coherence
|
|
#+RESULTS:
|
|
[[file:figs/coherence_id_dcm_dyn.png]]
|
|
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Bode plot of the DCM dynamics in the frame of the Fast Jacks
|
|
figure;
|
|
tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
|
|
|
|
ax1 = nexttile([2,1]);
|
|
hold on;
|
|
plot(f, abs(G_ur(:,1)), 'DisplayName', '$u_r$');
|
|
plot(f, abs(G_uh(:,2)), 'DisplayName', '$u_h$');
|
|
plot(f, abs(G_d( :,3)), 'DisplayName', '$u_d$');
|
|
plot(f, abs(G_ur(:,2)), 'color', [0, 0, 0, 0.2], ...
|
|
'DisplayName', 'Off Diagonal');
|
|
plot(f, abs(G_ur(:,3)), 'color', [0, 0, 0, 0.2], ...
|
|
'HandleVisibility', 'off');
|
|
plot(f, abs(G_uh(:,1)), 'color', [0, 0, 0, 0.2], ...
|
|
'HandleVisibility', 'off');
|
|
plot(f, abs(G_uh(:,3)), 'color', [0, 0, 0, 0.2], ...
|
|
'HandleVisibility', 'off');
|
|
plot(f, abs(G_d(:,1)), 'color', [0, 0, 0, 0.2], ...
|
|
'HandleVisibility', 'off');
|
|
plot(f, abs(G_d(:,2)), 'color', [0, 0, 0, 0.2], ...
|
|
'HandleVisibility', 'off');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
|
|
ylim([1e-2, 1e2]);
|
|
legend('location', 'southeast');
|
|
|
|
ax2 = nexttile;
|
|
hold on;
|
|
plot(f, 180/pi*angle(G_ur(:,1)));
|
|
plot(f, 180/pi*angle(G_uh(:,2)));
|
|
plot(f, 180/pi*angle(G_d(:,3)));
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
hold off;
|
|
yticks(-360:90:360);
|
|
% ylim([-45, 45]);
|
|
|
|
linkaxes([ax1,ax2],'x');
|
|
xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :tangle no :exports results :results file replace
|
|
exportFig('figs/bode_plot_dcm_dynamics.pdf', 'width', 'wide', 'height', 'tall');
|
|
#+end_src
|
|
|
|
#+name: fig:bode_plot_dcm_dynamics
|
|
#+caption: Bode Plot of the DCM dynamics in the frame of the fast jack.
|
|
#+RESULTS:
|
|
[[file:figs/bode_plot_dcm_dynamics.png]]
|
|
|
|
|
|
#+begin_src matlab
|
|
%% Previously used controller
|
|
load('X_tal_cage_PID.mat', 'K');
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Controller design
|
|
s = tf('s');
|
|
|
|
% Lead
|
|
a = 4; % Amount of phase lead / width of the phase lead / high frequency gain
|
|
wc = 2*pi*20; % Frequency with the maximum phase lead [rad/s]
|
|
|
|
% Low Pass Filter
|
|
w0 = 2*pi*100; % Cut-off frequency [rad/s]
|
|
xi = 0.4; % Damping Ratio
|
|
|
|
Kb = eye(3)*(2*pi*20)^2/(s^2) *1/(sqrt(a))* (1 + s/(wc/sqrt(a)))/(1 + s/(wc*sqrt(a))) * 1/(1 + 2*xi/w0*s + s^2/w0^2);;
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Loop Gain
|
|
figure;
|
|
tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
|
|
|
|
ax1 = nexttile([2,1]);
|
|
hold on;
|
|
plot(f, abs(G_ur(:,1).*squeeze(freqresp(Kb(1,1), f, 'Hz'))), 'DisplayName', '$u_r$');
|
|
plot(f, abs(G_uh(:,2).*squeeze(freqresp(Kb(2,2), f, 'Hz'))), 'DisplayName', '$u_h$');
|
|
plot(f, abs(G_d( :,3).*squeeze(freqresp(Kb(3,3), f, 'Hz'))), 'DisplayName', '$d$');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
|
|
ylim([1e-4, 1e2]);
|
|
legend('location', 'northeast');
|
|
|
|
ax2 = nexttile;
|
|
hold on;
|
|
plot(f, 180/pi*angle(G_ur(:,1).*squeeze(freqresp(K(1,1), f, 'Hz'))));
|
|
plot(f, 180/pi*angle(G_uh(:,2).*squeeze(freqresp(K(2,2), f, 'Hz'))));
|
|
plot(f, 180/pi*angle(G_d(:,3).*squeeze(freqresp(K(3,3), f, 'Hz'))));
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
hold off;
|
|
yticks(-360:90:360);
|
|
|
|
linkaxes([ax1,ax2],'x');
|
|
xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :tangle no :exports results :results file replace
|
|
exportFig('figs/loop_gain_dcm_contr_simple.pdf', 'width', 'wide', 'height', 'tall');
|
|
#+end_src
|
|
|
|
#+name: fig:loop_gain_dcm_contr_simple
|
|
#+caption: Loop gain
|
|
#+RESULTS:
|
|
[[file:figs/loop_gain_dcm_contr_simple.png]]
|
|
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Loop Gain
|
|
figure;
|
|
tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
|
|
|
|
ax1 = nexttile([2,1]);
|
|
hold on;
|
|
plot(f, abs(G_d(:,3).*squeeze(freqresp(K(1,1), f, 'Hz'))));
|
|
plot(f, abs(G_d(:,3).*squeeze(freqresp(Kb(1,1), f, 'Hz'))));
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
|
|
ylim([1e-4, 1e2]);
|
|
|
|
ax2 = nexttile;
|
|
hold on;
|
|
plot(f, 180/pi*angle(G_d(:,3).*squeeze(freqresp(K(1,1), f, 'Hz'))));
|
|
plot(f, 180/pi*angle(G_d(:,3).*squeeze(freqresp(Kb(1,1), f, 'Hz'))));
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
hold off;
|
|
yticks(-360:90:360);
|
|
|
|
linkaxes([ax1,ax2],'x');
|
|
xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :tangle no :exports results :results file replace
|
|
exportFig('figs/loop_gain_diag_old_new_contr.pdf', 'width', 'wide', 'height', 'tall');
|
|
#+end_src
|
|
|
|
#+name: fig:loop_gain_diag_old_new_contr
|
|
#+caption: Loop gain
|
|
#+RESULTS:
|
|
[[file:figs/loop_gain_diag_old_new_contr.png]]
|
|
|
|
Compare Sensitivity functions
|
|
#+begin_src matlab :exports none
|
|
%% Sensitivity function
|
|
figure;
|
|
tiledlayout(1, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
|
|
|
|
hold on;
|
|
plot(f, 1./abs(1 + G_d(:,3).*squeeze(freqresp(K(1,1), f, 'Hz'))));
|
|
plot(f, 1./abs(1 + G_d(:,3).*squeeze(freqresp(Kb(1,1), f, 'Hz'))));
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('Amplitude');
|
|
ylim([1e-2, 1e1]); xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :tangle no :exports results :results file replace
|
|
exportFig('figs/sensitivity_comp.pdf', 'width', 'wide', 'height', 'normal');
|
|
#+end_src
|
|
|
|
#+name: fig:sensitivity_comp
|
|
#+caption: Comparison of sensitivity functions
|
|
#+RESULTS:
|
|
[[file:figs/sensitivity_comp.png]]
|
|
|
|
|
|
#+begin_src matlab
|
|
L = zeros(3, 3, length(f));
|
|
Lb = zeros(3, 3, length(f));
|
|
|
|
for i_f = 1:length(f)
|
|
L(:,:,i_f) = [G_ur(i_f,:); G_uh(i_f,:); G_d(i_f,:)]*freqresp(K , f(i_f), 'Hz');
|
|
Lb(:,:,i_f) = [G_ur(i_f,:); G_uh(i_f,:); G_d(i_f,:)]*freqresp(Kb, f(i_f), 'Hz');
|
|
end
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Compute the Eigenvalues of the loop gain
|
|
Ldet = zeros(3, length(f));
|
|
Ldetb = zeros(3, length(f));
|
|
|
|
% Loop gain
|
|
for i_f = 1:length(f)
|
|
Ldet( :, i_f) = eig(squeeze(L( :,:,i_f)));
|
|
Ldetb(:, i_f) = eig(squeeze(Lb(:,:,i_f)));
|
|
end
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Plot of the eigenvalues of L in the complex plane
|
|
figure;
|
|
hold on;
|
|
for i = 1:size(Ldet,1)
|
|
plot(real(Ldet(i,:)), imag(Ldet(i,:)), '.', 'color', colors(1, :));
|
|
plot(real(Ldet(i,:)), -imag(Ldet(i,:)), '.', 'color', colors(1, :));
|
|
|
|
plot(real(Ldetb(i,:)), imag(Ldetb(i,:)), '.', 'color', colors(2, :));
|
|
plot(real(Ldetb(i,:)), -imag(Ldetb(i,:)), '.', 'color', colors(2, :));
|
|
end
|
|
plot(-1, 0, 'kx', 'HandleVisibility', 'off');
|
|
hold off;
|
|
set(gca, 'XScale', 'lin'); set(gca, 'YScale', 'lin');
|
|
xlabel('Real'); ylabel('Imag');
|
|
xlim([-3, 1]); ylim([-2, 2]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :tangle no :exports results :results file replace
|
|
exportFig('figs/loci_loop_gain_comp_controllers.pdf', 'width', 'wide', 'height', 'tall');
|
|
#+end_src
|
|
|
|
#+caption: Root Locus
|
|
#+RESULTS:
|
|
[[file:figs/loci_loop_gain_comp_controllers.png]]
|
|
|
|
|
|
** Identification - New
|
|
#+begin_src matlab
|
|
ur = load('FJPUR_step_new.mat');
|
|
uh = load('FJPUH_step_new.mat');
|
|
d = load('FJPD_step_new.mat');
|
|
#+end_src
|
|
|
|
| 1 | dz311 |
|
|
| 2 | dry311 |
|
|
| 3 | drx311 |
|
|
| 4 | dz111 |
|
|
| 5 | dry111 |
|
|
| 6 | drx111 |
|
|
| 7 | fjpur |
|
|
| 8 | fjpuh |
|
|
| 9 | fjpd |
|
|
| 10 | bragg |
|
|
|
|
#+begin_src matlab
|
|
ur.time = ur.time - ur.time(1);
|
|
ur.allValues(:, 1) = ur.allValues(:, 1) - mean(ur.allValues(ur.time<0.1, 1));
|
|
ur.allValues(:, 2) = ur.allValues(:, 2) - mean(ur.allValues(ur.time<0.1, 2));
|
|
ur.allValues(:, 3) = ur.allValues(:, 3) - mean(ur.allValues(ur.time<0.1, 3));
|
|
|
|
t_filt = ur.time < 5;
|
|
|
|
ur.u = ur.allValues(t_filt, 7);
|
|
ur.y_111 = [-ur.allValues(t_filt, 1), ur.allValues(t_filt, 2), ur.allValues(t_filt, 3)];
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
uh.time = uh.time - uh.time(1);
|
|
uh.allValues(:, 1) = uh.allValues(:, 1) - mean(uh.allValues(uh.time<0.1, 1));
|
|
uh.allValues(:, 2) = uh.allValues(:, 2) - mean(uh.allValues(uh.time<0.1, 2));
|
|
uh.allValues(:, 3) = uh.allValues(:, 3) - mean(uh.allValues(uh.time<0.1, 3));
|
|
|
|
uh.u = uh.allValues(t_filt, 8);
|
|
uh.y_111 = [-uh.allValues(t_filt, 1), uh.allValues(t_filt, 2), uh.allValues(t_filt, 3)];
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
d.time = d.time - d.time(1);
|
|
d.allValues(:, 1) = d.allValues(:, 1) - mean(d.allValues(d.time<0.1, 1));
|
|
d.allValues(:, 2) = d.allValues(:, 2) - mean(d.allValues(d.time<0.1, 2));
|
|
d.allValues(:, 3) = d.allValues(:, 3) - mean(d.allValues(d.time<0.1, 3));
|
|
|
|
d.u = d.allValues(t_filt, 9);
|
|
d.y_111 = [-d.allValues(t_filt, 1), d.allValues(t_filt, 2), d.allValues(t_filt, 3)];
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
J_a_111 = [1, 0.14, -0.0675
|
|
1, 0.14, 0.1525
|
|
1, -0.14, 0.0425];
|
|
|
|
J_a_311 = [1, 0.14, -0.1525
|
|
1, 0.14, 0.0675
|
|
1, -0.14, -0.0425];
|
|
|
|
ur.y = [J_a_311 * ur.y_111']';
|
|
uh.y = [J_a_311 * uh.y_111']';
|
|
d.y = [J_a_311 * d.y_111']';
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Sampling Time and Frequency
|
|
Ts = 1e-4; % [s]
|
|
Fs = 1/Ts; % [Hz]
|
|
|
|
% Hannning Windows
|
|
win = hanning(ceil(5*Fs));
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% And we get the frequency vector
|
|
[G_ur, f] = tfestimate(ur.u, ur.y, win, [], [], 1/Ts);
|
|
[G_uh, ~] = tfestimate(uh.u, uh.y, win, [], [], 1/Ts);
|
|
[G_d, ~] = tfestimate(d.u, d.y, win, [], [], 1/Ts);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
[coh_ur, ~] = mscohere(ur.u, ur.y, win, [], [], 1/Ts);
|
|
[coh_uh, ~] = mscohere(uh.u, uh.y, win, [], [], 1/Ts);
|
|
[coh_d, ~] = mscohere(d.u, d.y, win, [], [], 1/Ts);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Coherence
|
|
figure;
|
|
tiledlayout(1, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
|
|
hold on;
|
|
plot(f, coh_ur(:,1), 'DisplayName', '$u_r$');
|
|
plot(f, coh_uh(:,2), 'DisplayName', '$u_h$');
|
|
plot(f, coh_d( :,3), 'DisplayName', '$u_d$');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Coherence');
|
|
legend('location', 'southeast');
|
|
xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Bode plot of the DCM dynamics in the frame of the Fast Jacks
|
|
figure;
|
|
tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
|
|
|
|
ax1 = nexttile([2,1]);
|
|
hold on;
|
|
plot(f, abs(G_ur(:,1)), 'DisplayName', '$u_r$');
|
|
plot(f, abs(G_uh(:,2)), 'DisplayName', '$u_h$');
|
|
plot(f, abs(G_d( :,3)), 'DisplayName', '$u_d$');
|
|
plot(f, abs(G_ur(:,2)), 'color', [0, 0, 0, 0.2], ...
|
|
'DisplayName', 'Off Diagonal');
|
|
plot(f, abs(G_ur(:,3)), 'color', [0, 0, 0, 0.2], ...
|
|
'HandleVisibility', 'off');
|
|
plot(f, abs(G_uh(:,1)), 'color', [0, 0, 0, 0.2], ...
|
|
'HandleVisibility', 'off');
|
|
plot(f, abs(G_uh(:,3)), 'color', [0, 0, 0, 0.2], ...
|
|
'HandleVisibility', 'off');
|
|
plot(f, abs(G_d(:,1)), 'color', [0, 0, 0, 0.2], ...
|
|
'HandleVisibility', 'off');
|
|
plot(f, abs(G_d(:,2)), 'color', [0, 0, 0, 0.2], ...
|
|
'HandleVisibility', 'off');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
|
|
ylim([1e-2, 1e2]);
|
|
legend('location', 'southeast');
|
|
|
|
ax2 = nexttile;
|
|
hold on;
|
|
plot(f, 180/pi*angle(G_ur(:,1)));
|
|
plot(f, 180/pi*angle(G_uh(:,2)));
|
|
plot(f, 180/pi*angle(G_d(:,3)));
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
hold off;
|
|
yticks(-360:90:360);
|
|
% ylim([-45, 45]);
|
|
|
|
linkaxes([ax1,ax2],'x');
|
|
xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Previously used controller
|
|
load('X_tal_cage_PID.mat', 'K');
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Controller design
|
|
s = tf('s');
|
|
|
|
% Lead
|
|
a = 4; % Amount of phase lead / width of the phase lead / high frequency gain
|
|
wc = 2*pi*20; % Frequency with the maximum phase lead [rad/s]
|
|
|
|
% Low Pass Filter
|
|
w0 = 2*pi*100; % Cut-off frequency [rad/s]
|
|
xi = 0.4; % Damping Ratio
|
|
|
|
Kb = eye(3)*(2*pi*20)^2/(s^2) *1/(sqrt(a))* (1 + s/(wc/sqrt(a)))/(1 + s/(wc*sqrt(a))) * 1/(1 + 2*xi/w0*s + s^2/w0^2);;
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Loop Gain
|
|
figure;
|
|
tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
|
|
|
|
ax1 = nexttile([2,1]);
|
|
hold on;
|
|
plot(f, abs(G_ur(:,1).*squeeze(freqresp(Kb(1,1), f, 'Hz'))), 'DisplayName', '$u_r$');
|
|
plot(f, abs(G_uh(:,2).*squeeze(freqresp(Kb(2,2), f, 'Hz'))), 'DisplayName', '$u_h$');
|
|
plot(f, abs(G_d( :,3).*squeeze(freqresp(Kb(3,3), f, 'Hz'))), 'DisplayName', '$d$');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
|
|
ylim([1e-4, 1e2]);
|
|
legend('location', 'northeast');
|
|
|
|
ax2 = nexttile;
|
|
hold on;
|
|
plot(f, 180/pi*angle(G_ur(:,1).*squeeze(freqresp(K(1,1), f, 'Hz'))));
|
|
plot(f, 180/pi*angle(G_uh(:,2).*squeeze(freqresp(K(2,2), f, 'Hz'))));
|
|
plot(f, 180/pi*angle(G_d(:,3).*squeeze(freqresp(K(3,3), f, 'Hz'))));
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
hold off;
|
|
yticks(-360:90:360);
|
|
|
|
linkaxes([ax1,ax2],'x');
|
|
xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
Compare Sensitivity functions
|
|
#+begin_src matlab :exports none
|
|
%% Sensitivity function
|
|
figure;
|
|
tiledlayout(1, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
|
|
|
|
|
|
hold on;
|
|
plot(f, 1./abs(1 + G_d(:,3).*squeeze(freqresp(K(1,1), f, 'Hz'))));
|
|
plot(f, 1./abs(1 + G_d(:,3).*squeeze(freqresp(Kb(1,1), f, 'Hz'))));
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('Amplitude');
|
|
ylim([1e-2, 1e1]); xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :tangle no :exports results :results file replace
|
|
exportFig('figs/sensitivity_comp.pdf', 'width', 'wide', 'height', 'normal');
|
|
#+end_src
|
|
|
|
#+name: fig:sensitivity_comp
|
|
#+caption: Comparison of sensitivity functions
|
|
#+RESULTS:
|
|
[[file:figs/sensitivity_comp.png]]
|
|
|
|
|
|
#+begin_src matlab
|
|
L = zeros(3, 3, length(f));
|
|
Lb = zeros(3, 3, length(f));
|
|
|
|
for i_f = 1:length(f)
|
|
L(:,:,i_f) = [G_ur(i_f,:); G_uh(i_f,:); G_d(i_f,:)]*freqresp(K , f(i_f), 'Hz');
|
|
Lb(:,:,i_f) = [G_ur(i_f,:); G_uh(i_f,:); G_d(i_f,:)]*freqresp(Kb, f(i_f), 'Hz');
|
|
end
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Compute the Eigenvalues of the loop gain
|
|
Ldet = zeros(3, length(f));
|
|
Ldetb = zeros(3, length(f));
|
|
|
|
% Loop gain
|
|
for i_f = 1:length(f)
|
|
Ldet( :, i_f) = eig(squeeze(L( :,:,i_f)));
|
|
Ldetb(:, i_f) = eig(squeeze(Lb(:,:,i_f)));
|
|
end
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Plot of the eigenvalues of L in the complex plane
|
|
figure;
|
|
hold on;
|
|
for i = 1:size(Ldet,1)
|
|
plot(real(Ldet(i,:)), imag(Ldet(i,:)), '.', 'color', colors(1, :));
|
|
plot(real(Ldet(i,:)), -imag(Ldet(i,:)), '.', 'color', colors(1, :));
|
|
|
|
plot(real(Ldetb(i,:)), imag(Ldetb(i,:)), '.', 'color', colors(2, :));
|
|
plot(real(Ldetb(i,:)), -imag(Ldetb(i,:)), '.', 'color', colors(2, :));
|
|
end
|
|
plot(-1, 0, 'kx', 'HandleVisibility', 'off');
|
|
hold off;
|
|
set(gca, 'XScale', 'lin'); set(gca, 'YScale', 'lin');
|
|
xlabel('Real'); ylabel('Imag');
|
|
xlim([-3, 1]); ylim([-2, 2]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :tangle no :exports results :results file replace
|
|
exportFig('figs/loci_loop_gain_comp_controllers.pdf', 'width', 'wide', 'height', 'tall');
|
|
#+end_src
|
|
|
|
#+caption: Root Locus
|
|
#+RESULTS:
|
|
[[file:figs/loci_loop_gain_comp_controllers.png]]
|
|
|
|
|
|
** Identification - White noise
|
|
#+begin_src matlab
|
|
ur = load('fjpur_white_noise.mat');
|
|
uh = load('fjpuh_white_noise.mat');
|
|
d = load('fjpd_white_noise.mat');
|
|
#+end_src
|
|
|
|
| 1 | dz111 |
|
|
| 2 | dry111 |
|
|
| 3 | drx111 |
|
|
| 4 | fjpur |
|
|
| 5 | fjpuh |
|
|
| 6 | fjpd |
|
|
| 7 | bragg |
|
|
|
|
#+begin_src matlab
|
|
ur.time = ur.time - ur.time(1);
|
|
|
|
ur.drx = ur.drx - mean(ur.drx);
|
|
ur.dry = ur.dry - mean(ur.dry);
|
|
ur.dz = ur.dz - mean(ur.dz);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
uh.time = uh.time - uh.time(1);
|
|
|
|
uh.drx = uh.drx - mean(uh.drx);
|
|
uh.dry = uh.dry - mean(uh.dry);
|
|
uh.dz = uh.dz - mean(uh.dz);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
d.time = d.time - d.time(1);
|
|
|
|
d.drx = d.drx - mean(d.drx);
|
|
d.dry = d.dry - mean(d.dry);
|
|
d.dz = d.dz - mean(d.dz);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
J_a_111 = [1, 0.14, -0.0675
|
|
1, 0.14, 0.1525
|
|
1, -0.14, 0.0425];
|
|
|
|
ur.y = [J_a_111 * [-ur.dz, ur.dry, ur.drx]']';
|
|
uh.y = [J_a_111 * [-uh.dz, uh.dry, uh.drx]']';
|
|
d.y = [J_a_111 * [-d.dz, d.dry, d.drx]']';
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Sampling Time and Frequency
|
|
Ts = 1e-4; % [s]
|
|
Fs = 1/Ts; % [Hz]
|
|
|
|
% Hannning Windows
|
|
win = hanning(ceil(0.5*Fs));
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% And we get the frequency vector
|
|
[G_ur, f] = tfestimate(ur.fjpur, ur.y, win, [], [], 1/Ts);
|
|
[G_uh, ~] = tfestimate(uh.fjpuh, uh.y, win, [], [], 1/Ts);
|
|
[G_d, ~] = tfestimate(d.fjpd, d.y, win, [], [], 1/Ts);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
[coh_ur, ~] = mscohere(ur.fjpur, ur.y, win, [], [], 1/Ts);
|
|
[coh_uh, ~] = mscohere(uh.fjpuh, uh.y, win, [], [], 1/Ts);
|
|
[coh_d, ~] = mscohere(d.fjpd, d.y, win, [], [], 1/Ts);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Coherence
|
|
figure;
|
|
tiledlayout(1, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
|
|
hold on;
|
|
plot(f, coh_ur(:,1), 'DisplayName', '$u_r$');
|
|
plot(f, coh_uh(:,2), 'DisplayName', '$u_h$');
|
|
plot(f, coh_d(:,3), 'DisplayName', '$d$');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Coherence');
|
|
legend('location', 'southeast');
|
|
xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :tangle no :exports results :results file replace
|
|
exportFig('figs/coherence_ident_noise.pdf', 'width', 'wide', 'height', 'normal');
|
|
#+end_src
|
|
|
|
#+name: fig:coherence_ident_noise
|
|
#+caption: description
|
|
#+RESULTS:
|
|
[[file:figs/coherence_ident_noise.png]]
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Bode plot of the DCM dynamics in the frame of the Fast Jacks
|
|
figure;
|
|
tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
|
|
|
|
ax1 = nexttile([2,1]);
|
|
hold on;
|
|
plot(f, abs(G_ur(:,1)), 'DisplayName', '$u_r$');
|
|
plot(f, abs(G_uh(:,2)), 'DisplayName', '$u_h$');
|
|
plot(f, abs(G_d( :,3)), 'DisplayName', '$u_d$');
|
|
plot(f, abs(G_ur(:,2)), 'color', [0, 0, 0, 0.2], ...
|
|
'DisplayName', 'Off Diagonal');
|
|
plot(f, abs(G_ur(:,3)), 'color', [0, 0, 0, 0.2], ...
|
|
'HandleVisibility', 'off');
|
|
plot(f, abs(G_uh(:,1)), 'color', [0, 0, 0, 0.2], ...
|
|
'HandleVisibility', 'off');
|
|
plot(f, abs(G_uh(:,3)), 'color', [0, 0, 0, 0.2], ...
|
|
'HandleVisibility', 'off');
|
|
plot(f, abs(G_d(:,1)), 'color', [0, 0, 0, 0.2], ...
|
|
'HandleVisibility', 'off');
|
|
plot(f, abs(G_d(:,2)), 'color', [0, 0, 0, 0.2], ...
|
|
'HandleVisibility', 'off');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
|
|
ylim([1e-2, 1e2]);
|
|
legend('location', 'northwest');
|
|
|
|
ax2 = nexttile;
|
|
hold on;
|
|
plot(f, 180/pi*angle(G_ur(:,1)));
|
|
plot(f, 180/pi*angle(G_uh(:,2)));
|
|
plot(f, 180/pi*angle(G_d(:,3)));
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
hold off;
|
|
yticks(-360:90:360);
|
|
% ylim([-45, 45]);
|
|
|
|
linkaxes([ax1,ax2],'x');
|
|
xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :tangle no :exports results :results file replace
|
|
exportFig('figs/bode_plot_ident_noise.pdf', 'width', 'full', 'height', 'tall');
|
|
#+end_src
|
|
|
|
#+name: fig:bode_plot_ident_noise
|
|
#+caption: Bode Plot of the DCM dynamics in the frame of the fast jack.
|
|
#+RESULTS:
|
|
[[file:figs/bode_plot_ident_noise.png]]
|
|
|
|
|
|
#+begin_src matlab
|
|
%% Previously used controller
|
|
load('X_tal_cage_PID.mat', 'K');
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Controller design
|
|
s = tf('s');
|
|
|
|
% Lead
|
|
a = 8; % Amount of phase lead / width of the phase lead / high frequency gain
|
|
wc = 2*pi*20; % Frequency with the maximum phase lead [rad/s]
|
|
|
|
% Low Pass Filter
|
|
w0 = 2*pi*80; % Cut-off frequency [rad/s]
|
|
xi = 0.4; % Damping Ratio
|
|
|
|
Kb = eye(3)*(2*pi*20)^2/(s^2) *1/(sqrt(a))* (1 + s/(wc/sqrt(a)))/(1 + s/(wc*sqrt(a))) * 1/(1 + 2*xi/w0*s + s^2/w0^2);;
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Loop Gain
|
|
figure;
|
|
tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
|
|
|
|
ax1 = nexttile([2,1]);
|
|
hold on;
|
|
plot(f, abs(G_ur(:,1).*squeeze(freqresp(Kb(1,1), f, 'Hz'))), 'DisplayName', '$u_r$');
|
|
plot(f, abs(G_uh(:,2).*squeeze(freqresp(Kb(2,2), f, 'Hz'))), 'DisplayName', '$u_h$');
|
|
plot(f, abs(G_d( :,3).*squeeze(freqresp(Kb(3,3), f, 'Hz'))), 'DisplayName', '$d$');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
|
|
ylim([1e-4, 1e2]);
|
|
legend('location', 'northeast');
|
|
|
|
ax2 = nexttile;
|
|
hold on;
|
|
plot(f, 180/pi*angle(G_ur(:,1).*squeeze(freqresp(Kb(1,1), f, 'Hz'))));
|
|
plot(f, 180/pi*angle(G_uh(:,2).*squeeze(freqresp(Kb(2,2), f, 'Hz'))));
|
|
plot(f, 180/pi*angle(G_d(:,3).*squeeze(freqresp(Kb(3,3), f, 'Hz'))));
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
yticks(-360:90:360);
|
|
|
|
linkaxes([ax1,ax2],'x');
|
|
xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :tangle no :exports results :results file replace
|
|
exportFig('figs/loop_gain_dcm_contr_simple.pdf', 'width', 'full', 'height', 'tall');
|
|
#+end_src
|
|
|
|
#+name: fig:loop_gain_dcm_contr_simple
|
|
#+caption: Loop gain
|
|
#+RESULTS:
|
|
[[file:figs/loop_gain_dcm_contr_simple.png]]
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Loop Gain
|
|
figure;
|
|
tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
|
|
|
|
ax1 = nexttile([2,1]);
|
|
hold on;
|
|
plot(f, abs(G_ur(:,1).*squeeze(freqresp(K(1,1), f, 'Hz'))), 'color', colors(1,:), 'DisplayName', 'Old K');
|
|
plot(f, abs(G_uh(:,2).*squeeze(freqresp(K(2,2), f, 'Hz'))), 'color', colors(1,:), 'HandleVisibility', 'off');
|
|
plot(f, abs(G_d(:,3).*squeeze(freqresp(K(3,3), f, 'Hz'))), 'color', colors(1,:), 'HandleVisibility', 'off');
|
|
plot(f, abs(G_ur(:,1).*squeeze(freqresp(Kb(1,1), f, 'Hz'))), 'color', colors(2,:), 'DisplayName', 'New K');
|
|
plot(f, abs(G_uh(:,2).*squeeze(freqresp(Kb(2,2), f, 'Hz'))), 'color', colors(2,:), 'HandleVisibility', 'off');
|
|
plot(f, abs(G_d(:,3).*squeeze(freqresp(Kb(3,3), f, 'Hz'))), 'color', colors(2,:), 'HandleVisibility', 'off');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
|
|
ylim([1e-4, 1e2]);
|
|
legend('location', 'northeast');
|
|
|
|
ax2 = nexttile;
|
|
hold on;
|
|
plot(f, 180/pi*angle(G_ur(:,1).*squeeze(freqresp(K(1,1), f, 'Hz'))), 'color', colors(1,:));
|
|
plot(f, 180/pi*angle(G_uh(:,2).*squeeze(freqresp(K(2,2), f, 'Hz'))), 'color', colors(1,:));
|
|
plot(f, 180/pi*angle(G_d(:,3).*squeeze(freqresp(K(3,3), f, 'Hz'))), 'color', colors(1,:));
|
|
plot(f, 180/pi*angle(G_ur(:,1).*squeeze(freqresp(Kb(1,1), f, 'Hz'))), 'color', colors(2,:));
|
|
plot(f, 180/pi*angle(G_uh(:,2).*squeeze(freqresp(Kb(2,2), f, 'Hz'))), 'color', colors(2,:));
|
|
plot(f, 180/pi*angle(G_d(:,3).*squeeze(freqresp(Kb(3,3), f, 'Hz'))), 'color', colors(2,:));
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
yticks(-360:90:360);
|
|
|
|
linkaxes([ax1,ax2],'x');
|
|
xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :tangle no :exports results :results file replace
|
|
exportFig('figs/loop_gain_diag_old_new_contr.pdf', 'width', 'full', 'height', 'tall');
|
|
#+end_src
|
|
|
|
#+name: fig:loop_gain_diag_old_new_contr
|
|
#+caption: Loop gain
|
|
#+RESULTS:
|
|
[[file:figs/loop_gain_diag_old_new_contr.png]]
|
|
|
|
Compare Sensitivity functions
|
|
#+begin_src matlab :exports none
|
|
%% Sensitivity function
|
|
figure;
|
|
tiledlayout(1, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
|
|
|
|
hold on;
|
|
plot(f, 1./abs(1 + G_d(:,3).*squeeze(freqresp(K(1,1), f, 'Hz'))), 'DisplayName', 'Old K');
|
|
plot(f, 1./abs(1 + G_d(:,3).*squeeze(freqresp(Kb(1,1), f, 'Hz'))), 'DisplayName', 'New K');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('Amplitude');
|
|
ylim([1e-2, 1e1]); xlim([1, 1e3]);
|
|
legend('location', 'southeast');
|
|
#+end_src
|
|
|
|
#+begin_src matlab :tangle no :exports results :results file replace
|
|
exportFig('figs/sensitivity_comp.pdf', 'width', 'full', 'height', 'tall');
|
|
#+end_src
|
|
|
|
#+name: fig:sensitivity_comp
|
|
#+caption: Comparison of sensitivity functions
|
|
#+RESULTS:
|
|
[[file:figs/sensitivity_comp.png]]
|
|
|
|
|
|
#+begin_src matlab
|
|
L = zeros(3, 3, length(f));
|
|
Lb = zeros(3, 3, length(f));
|
|
|
|
for i_f = 1:length(f)
|
|
L(:,:,i_f) = [G_ur(i_f,:); G_uh(i_f,:); G_d(i_f,:)]*freqresp(K , f(i_f), 'Hz');
|
|
Lb(:,:,i_f) = [G_ur(i_f,:); G_uh(i_f,:); G_d(i_f,:)]*freqresp(Kb, f(i_f), 'Hz');
|
|
end
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Compute the Eigenvalues of the loop gain
|
|
Ldet = zeros(3, length(f));
|
|
Ldetb = zeros(3, length(f));
|
|
|
|
% Loop gain
|
|
for i_f = 1:length(f)
|
|
Ldet( :, i_f) = eig(squeeze(L( :,:,i_f)));
|
|
Ldetb(:, i_f) = eig(squeeze(Lb(:,:,i_f)));
|
|
end
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Plot of the eigenvalues of L in the complex plane
|
|
figure;
|
|
hold on;
|
|
for i = 1:size(Ldet,1)
|
|
plot(real(Ldet(i,:)), imag(Ldet(i,:)), '.', 'color', colors(1, :));
|
|
plot(real(Ldet(i,:)), -imag(Ldet(i,:)), '.', 'color', colors(1, :));
|
|
|
|
plot(real(Ldetb(i,:)), imag(Ldetb(i,:)), '.', 'color', colors(2, :));
|
|
plot(real(Ldetb(i,:)), -imag(Ldetb(i,:)), '.', 'color', colors(2, :));
|
|
end
|
|
plot(-1, 0, 'kx', 'HandleVisibility', 'off');
|
|
hold off;
|
|
set(gca, 'XScale', 'lin'); set(gca, 'YScale', 'lin');
|
|
xlabel('Real'); ylabel('Imag');
|
|
xlim([-3, 1]); ylim([-2, 2]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :tangle no :exports results :results file replace
|
|
exportFig('figs/loci_loop_gain_comp_controllers.pdf', 'width', 'wide', 'height', 'tall');
|
|
#+end_src
|
|
|
|
#+caption: Root Locus
|
|
#+RESULTS:
|
|
[[file:figs/loci_loop_gain_comp_controllers.png]]
|
|
|
|
|
|
** Test
|
|
#+begin_src matlab
|
|
%% Notch
|
|
gm = 0.015;
|
|
xi = 0.1;
|
|
wn = 2*pi*208;
|
|
|
|
K_notch = (s^2 + 2*gm*xi*wn*s + wn^2)/(s^2 + 2*xi*wn*s + wn^2);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Double integrator
|
|
w0 = 2*pi*40;
|
|
K_int = (w0^2)/(s^2);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Lead
|
|
a = 3; % Amount of phase lead / width of the phase lead / high frequency gain
|
|
K_lead = 1/(sqrt(a))*(1 + s/(w0/sqrt(a)))/(1 + s/(w0*sqrt(a)));
|
|
K_lead = K_lead*K_lead;
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Low Pass Filter
|
|
w0 = 2*pi*120; % Cut-off frequency [rad/s]
|
|
xi = 0.3; % Damping Ratio
|
|
|
|
K_lpf = 1/(1 + 2*xi/w0*s + s^2/w0^2);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Diagonal controller
|
|
Kb = 0.8*eye(3)*K_notch*K_int*K_lead*K_lpf;
|
|
#+end_src
|
|
** New controller - Higher bandwidth
|
|
|
|
#+begin_src matlab
|
|
%% Previously used controller
|
|
load('X_tal_cage_PID.mat', 'K');
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Current Controller design
|
|
% Lead
|
|
a = 8; % Amount of phase lead / width of the phase lead / high frequency gain
|
|
wc = 2*pi*20; % Frequency with the maximum phase lead [rad/s]
|
|
|
|
% Low Pass Filter
|
|
w0 = 2*pi*80; % Cut-off frequency [rad/s]
|
|
xi = 0.4; % Damping Ratio
|
|
|
|
Kb_old = eye(3)*(2*pi*20)^2/(s^2) *1/(sqrt(a))* (1 + s/(wc/sqrt(a)))/(1 + s/(wc*sqrt(a))) * 1/(1 + 2*xi/w0*s + s^2/w0^2);;
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Notch
|
|
gm = 0.015;
|
|
xi = 0.2;
|
|
wn = 2*pi*208;
|
|
|
|
K_notch = (s^2 + 2*gm*xi*wn*s + wn^2)/(s^2 + 2*xi*wn*s + wn^2);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Double integrator
|
|
w0 = 2*pi*40;
|
|
K_int = (w0^2)/(s^2);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Lead
|
|
a = 3; % Amount of phase lead / width of the phase lead / high frequency gain
|
|
w0 = 2*pi*40;
|
|
K_lead = 1/(sqrt(a))*(1 + s/(w0/sqrt(a)))/(1 + s/(w0*sqrt(a)));
|
|
K_lead = K_lead*K_lead;
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Low Pass Filter
|
|
w0 = 2*pi*120; % Cut-off frequency [rad/s]
|
|
xi = 0.3; % Damping Ratio
|
|
|
|
K_lpf = 1/(1 + 2*xi/w0*s + s^2/w0^2);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Diagonal controller
|
|
Kb = 0.9*eye(3)*K_notch*K_int*K_lead*K_lpf;
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Loop Gain
|
|
figure;
|
|
tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
|
|
|
|
ax1 = nexttile([2,1]);
|
|
hold on;
|
|
plot(f, abs(G_ur(:,1).*squeeze(freqresp(K(1,1), f, 'Hz'))), 'DisplayName', 'PID');
|
|
plot(f, abs(G_ur(:,1).*squeeze(freqresp(Kb_old(1,1), f, 'Hz'))), 'DisplayName', 'Double Integrator + Lead');
|
|
plot(f, abs(G_ur(:,1).*squeeze(freqresp(Kb(1,1), f, 'Hz'))), 'DisplayName', 'Double Integrator + Lead + Notch');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
|
|
ylim([1e-3, 1e3]);
|
|
legend('location', 'southwest', 'FontSize', 8, 'NumColumns', 1);
|
|
|
|
ax2 = nexttile;
|
|
hold on;
|
|
plot(f, 180/pi*angle(G_ur(:,1).*squeeze(freqresp(K(1,1), f, 'Hz'))));
|
|
plot(f, 180/pi*angle(G_ur(:,1).*squeeze(freqresp(Kb_old(1,1), f, 'Hz'))));
|
|
plot(f, 180/pi*angle(G_ur(:,1).*squeeze(freqresp(Kb(1,1), f, 'Hz'))));
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
yticks(-360:90:360);
|
|
|
|
linkaxes([ax1,ax2],'x');
|
|
xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :tangle no :exports results :results file replace
|
|
exportFig('figs/loop_gain_compare.pdf', 'width', 'full', 'height', 'tall');
|
|
#+end_src
|
|
|
|
#+name: fig:loop_gain_compare
|
|
#+caption: description
|
|
#+RESULTS:
|
|
[[file:figs/loop_gain_compare.png]]
|
|
|
|
|
|
#+begin_src matlab
|
|
L = zeros(3, 3, length(f));
|
|
Lb = zeros(3, 3, length(f));
|
|
Lb_new = zeros(3, 3, length(f));
|
|
|
|
for i_f = 1:length(f)
|
|
L(:,:,i_f) = [G_ur(i_f,:); G_uh(i_f,:); G_d(i_f,:)]*freqresp(K , f(i_f), 'Hz');
|
|
Lb(:,:,i_f) = [G_ur(i_f,:); G_uh(i_f,:); G_d(i_f,:)]*freqresp(Kb_old, f(i_f), 'Hz');
|
|
Lb_new(:,:,i_f) = [G_ur(i_f,:); G_uh(i_f,:); G_d(i_f,:)]*freqresp(Kb, f(i_f), 'Hz');
|
|
end
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Compute the Eigenvalues of the loop gain
|
|
Ldet = zeros(3, length(f));
|
|
Ldetb = zeros(3, length(f));
|
|
Ldetb_new = zeros(3, length(f));
|
|
|
|
% Loop gain
|
|
for i_f = 1:length(f)
|
|
Ldet( :, i_f) = eig(squeeze(L( :,:,i_f)));
|
|
Ldetb(:, i_f) = eig(squeeze(Lb(:,:,i_f)));
|
|
Ldetb_new(:, i_f) = eig(squeeze(Lb_new(:,:,i_f)));
|
|
end
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Plot of the eigenvalues of L in the complex plane
|
|
figure;
|
|
hold on;
|
|
plot(real(Ldet(1,:)), imag(Ldet(1,:)), '.', 'color', colors(1, :), 'DisplayName', 'PID');
|
|
plot(real(Ldet(1,:)), -imag(Ldet(1,:)), '.', 'color', colors(1, :), 'HandleVisibility', 'off');
|
|
|
|
plot(real(Ldetb(1,:)), imag(Ldetb(1,:)), '.', 'color', colors(2, :), 'DisplayName', 'Double Integrator + Lead');
|
|
plot(real(Ldetb(1,:)), -imag(Ldetb(1,:)), '.', 'color', colors(2, :), 'HandleVisibility', 'off');
|
|
|
|
plot(real(Ldetb_new(1,:)), imag(Ldetb_new(1,:)), '.', 'color', colors(3, :), 'DisplayName', 'Double Integrator + Lead + Notch');
|
|
plot(real(Ldetb_new(1,:)), -imag(Ldetb_new(1,:)), '.', 'color', colors(3, :), 'HandleVisibility', 'off');
|
|
for i = 2:size(Ldet,1)
|
|
plot(real(Ldet(i,:)), imag(Ldet(i,:)), '.', 'color', colors(1, :), 'HandleVisibility', 'off');
|
|
plot(real(Ldet(i,:)), -imag(Ldet(i,:)), '.', 'color', colors(1, :), 'HandleVisibility', 'off');
|
|
|
|
plot(real(Ldetb(i,:)), imag(Ldetb(i,:)), '.', 'color', colors(2, :), 'HandleVisibility', 'off');
|
|
plot(real(Ldetb(i,:)), -imag(Ldetb(i,:)), '.', 'color', colors(2, :), 'HandleVisibility', 'off');
|
|
|
|
plot(real(Ldetb_new(i,:)), imag(Ldetb_new(i,:)), '.', 'color', colors(3, :), 'HandleVisibility', 'off');
|
|
plot(real(Ldetb_new(i,:)), -imag(Ldetb_new(i,:)), '.', 'color', colors(3, :), 'HandleVisibility', 'off');
|
|
end
|
|
plot(-1, 0, 'kx', 'HandleVisibility', 'off');
|
|
hold off;
|
|
set(gca, 'XScale', 'lin'); set(gca, 'YScale', 'lin');
|
|
xlabel('Real'); ylabel('Imag');
|
|
xlim([-3, 1]); ylim([-2, 2]);
|
|
axis square;
|
|
legend('location', 'southeast', 'FontSize', 8, 'NumColumns', 1);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :tangle no :exports results :results file replace
|
|
exportFig('figs/nyquist_compare.pdf', 'width', 'wide', 'height', 'tall');
|
|
#+end_src
|
|
|
|
#+name: fig:nyquist_compare
|
|
#+caption: Nyquist Plot
|
|
#+RESULTS:
|
|
[[file:figs/nyquist_compare.png]]
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Sensitivity function
|
|
figure;
|
|
tiledlayout(1, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
|
|
|
|
hold on;
|
|
plot(f, 1./abs(1 + G_d(:,3).*squeeze(freqresp(K(1,1), f, 'Hz'))), 'DisplayName', 'PID');
|
|
plot(f, 1./abs(1 + G_d(:,3).*squeeze(freqresp(Kb_old(1,1), f, 'Hz'))), 'DisplayName', '$K_{10Hz}$');
|
|
plot(f, 1./abs(1 + G_d(:,3).*squeeze(freqresp(Kb(1,1), f, 'Hz'))), 'DisplayName', '$K_{20Hz}$');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('Sensitivity Magnitude');
|
|
ylim([1e-3, 5]); xlim([1, 5e2]);
|
|
legend('location', 'southeast');
|
|
#+end_src
|
|
|
|
#+begin_src matlab :tangle no :exports results :results file replace
|
|
exportFig('figs/sensitivity_function_compare.pdf', 'width', 'normal', 'height', 'tall');
|
|
#+end_src
|
|
|
|
#+name: fig:sensitivity_function_compare
|
|
#+caption: description
|
|
#+RESULTS:
|
|
[[file:figs/sensitivity_function_compare.png]]
|
|
|
|
|
|
|
|
** Added gain
|
|
#+begin_src matlab
|
|
%% Notch
|
|
gm = 0.015;
|
|
xi = 0.2;
|
|
wn = 2*pi*208;
|
|
|
|
K_notch = (s^2 + 2*gm*xi*wn*s + wn^2)/(s^2 + 2*xi*wn*s + wn^2);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Double integrator
|
|
w0 = 2*pi*40;
|
|
K_int = (w0^2)/(s^2);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Lead
|
|
a = 3; % Amount of phase lead / width of the phase lead / high frequency gain
|
|
w0 = 2*pi*40;
|
|
K_lead = 1/(sqrt(a))*(1 + s/(w0/sqrt(a)))/(1 + s/(w0*sqrt(a)));
|
|
K_lead = K_lead*K_lead;
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Low Pass Filter
|
|
w0 = 2*pi*120; % Cut-off frequency [rad/s]
|
|
xi = 0.3; % Damping Ratio
|
|
|
|
K_lpf = 1/(1 + 2*xi/w0*s + s^2/w0^2);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
gm = 10;
|
|
xi = 0.02;
|
|
wn = 2*pi*15;
|
|
|
|
H = (s^2 + 2*gm*xi*wn*s + wn^2)/(s^2 + 2*xi*wn*s + wn^2);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Diagonal controller
|
|
Kb_gain = 0.9*eye(3)*H*K_notch*K_int*K_lead*K_lpf;
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
Lb_gain = zeros(3, 3, length(f));
|
|
|
|
for i_f = 1:length(f)
|
|
Lb_gain(:,:,i_f) = [G_ur(i_f,:); G_uh(i_f,:); G_d(i_f,:)]*freqresp(Kb_gain, f(i_f), 'Hz');
|
|
end
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Compute the Eigenvalues of the loop gain
|
|
Ldetb_gain = zeros(3, length(f));
|
|
|
|
% Loop gain
|
|
for i_f = 1:length(f)
|
|
Ldetb_gain(:, i_f) = eig(squeeze(Lb_gain(:,:,i_f)));
|
|
end
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Loop Gain
|
|
figure;
|
|
tiledlayout(3, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
|
|
|
|
ax1 = nexttile([2,1]);
|
|
hold on;
|
|
plot(f, abs(G_ur(:,1).*squeeze(freqresp(Kb(1,1), f, 'Hz'))), 'color', colors(3, :));
|
|
plot(f, abs(G_ur(:,1).*squeeze(freqresp(Kb_gain(1,1), f, 'Hz'))), 'color', colors(4, :));
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
ylabel('Amplitude'); set(gca, 'XTickLabel',[]);
|
|
ylim([1e-4, 1e2]);
|
|
|
|
ax2 = nexttile;
|
|
hold on;
|
|
plot(f, 180/pi*angle(G_ur(:,1).*squeeze(freqresp(Kb(1,1), f, 'Hz'))), 'color', colors(3, :));
|
|
plot(f, 180/pi*angle(G_ur(:,1).*squeeze(freqresp(Kb_gain(1,1), f, 'Hz'))), 'color', colors(4, :));
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'lin');
|
|
xlabel('Frequency [Hz]'); ylabel('Phase [deg]');
|
|
yticks(-360:90:360);
|
|
|
|
linkaxes([ax1,ax2],'x');
|
|
xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :tangle no :exports results :results file replace
|
|
exportFig('figs/loop_gain_compare_added_gain.pdf', 'width', 'wide', 'height', 'tall');
|
|
#+end_src
|
|
|
|
#+name: fig:loop_gain_compare_added_gain
|
|
#+caption: description
|
|
#+RESULTS:
|
|
[[file:figs/loop_gain_compare_added_gain.png]]
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Sensitivity function
|
|
figure;
|
|
tiledlayout(1, 1, 'TileSpacing', 'Compact', 'Padding', 'None');
|
|
|
|
hold on;
|
|
plot(f, 1./abs(1 + G_d(:,3).*squeeze(freqresp(K(1,1), f, 'Hz'))), 'DisplayName', 'Old K');
|
|
plot(f, 1./abs(1 + G_d(:,3).*squeeze(freqresp(Kb_old(1,1), f, 'Hz'))), 'DisplayName', 'Act K');
|
|
plot(f, 1./abs(1 + G_d(:,3).*squeeze(freqresp(Kb(1,1), f, 'Hz'))), 'DisplayName', 'New K');
|
|
plot(f, 1./abs(1 + G_d(:,3).*squeeze(freqresp(Kb_gain(1,1), f, 'Hz'))), 'DisplayName', 'New K');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('Amplitude');
|
|
ylim([1e-2, 1e1]); xlim([1, 1e3]);
|
|
legend('location', 'southeast');
|
|
#+end_src
|
|
|
|
#+begin_src matlab :tangle no :exports results :results file replace
|
|
exportFig('figs/sensitivity_new_gain_compare.pdf', 'width', 'wide', 'height', 'tall');
|
|
#+end_src
|
|
|
|
#+name: fig:sensitivity_new_gain_compare
|
|
#+caption: description
|
|
#+RESULTS:
|
|
[[file:figs/sensitivity_new_gain_compare.png]]
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Plot of the eigenvalues of L in the complex plane
|
|
figure;
|
|
hold on;
|
|
for i = 1:size(Ldet,1)
|
|
plot(real(Ldetb_new(i,:)), imag(Ldetb_new(i,:)), '.', 'color', colors(3, :));
|
|
plot(real(Ldetb_new(i,:)), -imag(Ldetb_new(i,:)), '.', 'color', colors(3, :));
|
|
plot(real(Ldetb_gain(i,:)), imag(Ldetb_gain(i,:)), '.', 'color', colors(4, :));
|
|
plot(real(Ldetb_gain(i,:)), -imag(Ldetb_gain(i,:)), '.', 'color', colors(4, :));
|
|
end
|
|
plot(-1, 0, 'kx', 'HandleVisibility', 'off');
|
|
hold off;
|
|
set(gca, 'XScale', 'lin'); set(gca, 'YScale', 'lin');
|
|
xlabel('Real'); ylabel('Imag');
|
|
xlim([-3, 1]); ylim([-2, 2]);
|
|
axis square;
|
|
#+end_src
|
|
|
|
#+begin_src matlab :tangle no :exports results :results file replace
|
|
exportFig('figs/nyquist_after_gain_frequency.pdf', 'width', 'wide', 'height', 'tall');
|
|
#+end_src
|
|
|
|
#+name: fig:nyquist_after_gain_frequency
|
|
#+caption: nyquist plot
|
|
#+RESULTS:
|
|
[[file:figs/nyquist_after_gain_frequency.png]]
|
|
|
|
|
|
* Noise Budgeting
|
|
** 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 :tangle no :noweb yes
|
|
<<m-init-path>>
|
|
#+end_src
|
|
|
|
#+begin_src matlab :eval no :noweb yes
|
|
<<m-init-path-tangle>>
|
|
#+end_src
|
|
|
|
#+begin_src matlab :noweb yes
|
|
<<m-init-other>>
|
|
#+end_src
|
|
|
|
** No Displacement
|
|
|
|
| 1 | dz311 |
|
|
| 2 | dry311 |
|
|
| 3 | drx311 |
|
|
| 4 | dz111 |
|
|
| 5 | dry111 |
|
|
| 6 | drx111 |
|
|
| 7 | fjpur |
|
|
| 8 | fjpuh |
|
|
| 9 | fjpd |
|
|
| 10 | bragg |
|
|
|
|
#+begin_src matlab
|
|
data_10_deg = load('no_mov_10.mat');
|
|
data_70_deg = load('no_mov_70.mat');
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
data_10_deg = extractDatData('no_mov_10.mat', ...
|
|
{"dz311", "dry311", "drx311", "dz", "dry", "drx", "fjpur", "fjpuh", "fjpd", "bragg"}, ...
|
|
[1e-9, 1e-9, 1e-9, 1e-9, 1e-9, 1e-9, 1e-8, 1e-8, 1e-8, pi/180]);
|
|
data_10_deg = processMeasData(data_10_deg);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
Ts = 1e-4;
|
|
t = Ts*[1:length(data_10_deg.bragg)];
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
data_10_deg.dz = data_10_deg.allValues(:,4) - mean(data_10_deg.allValues(:,4));
|
|
data_10_deg.dry = data_10_deg.allValues(:,5) - mean(data_10_deg.allValues(:,5));
|
|
data_10_deg.drx = data_10_deg.allValues(:,6) - mean(data_10_deg.allValues(:,6));
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Compute motion error in the frame of the fast jack
|
|
J_a_111 = [1, 0.14, -0.1525
|
|
1, 0.14, 0.0675
|
|
1, -0.14, 0.0425];
|
|
|
|
de_111 = [data_10_deg.dz'; data_10_deg.dry'; data_10_deg.drx'];
|
|
de_fj = J_a_111*de_111;
|
|
data_10_deg.fj_ur = de_fj(1,:)';
|
|
data_10_deg.fj_uh = de_fj(2,:)';
|
|
data_10_deg.fj_d = de_fj(3,:)';
|
|
|
|
de_111 = [data_70_deg.dz'; data_70_deg.dry'; data_70_deg.drx'];
|
|
de_fj = J_a_111*de_111;
|
|
data_70_deg.fj_ur = de_fj(1,:)';
|
|
data_70_deg.fj_uh = de_fj(2,:)';
|
|
data_70_deg.fj_d = de_fj(3,:)';
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
win = hanning(ceil(1/Ts));
|
|
|
|
[pxx_10_ur, f] = pwelch(data_10_deg.fj_ur, win, [], [], 1/Ts);
|
|
[pxx_70_ur, ~] = pwelch(data_70_deg.fj_ur, win, [], [], 1/Ts);
|
|
|
|
[pxx_10_uh, ~] = pwelch(data_10_deg.fj_uh, win, [], [], 1/Ts);
|
|
[pxx_70_uh, ~] = pwelch(data_70_deg.fj_uh, win, [], [], 1/Ts);
|
|
|
|
[pxx_10_d, ~] = pwelch(data_10_deg.fj_d, win, [], [], 1/Ts);
|
|
[pxx_70_d, ~] = pwelch(data_70_deg.fj_d, win, [], [], 1/Ts);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
CPS_10_ur = flip(-cumtrapz(flip(f), flip(pxx_10_ur)));
|
|
CPS_10_uh = flip(-cumtrapz(flip(f), flip(pxx_10_uh)));
|
|
CPS_10_d = flip(-cumtrapz(flip(f), flip(pxx_10_d)));
|
|
|
|
CPS_70_ur = flip(-cumtrapz(flip(f), flip(pxx_70_ur)));
|
|
CPS_70_uh = flip(-cumtrapz(flip(f), flip(pxx_70_uh)));
|
|
CPS_70_d = flip(-cumtrapz(flip(f), flip(pxx_70_d)));
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
figure;
|
|
hold on;
|
|
plot(f, sqrt(CPS_10_ur), '-' , 'color', colors(1, :), 'DisplayName', '10 deg - $u_r$')
|
|
plot(f, sqrt(CPS_70_ur), '--', 'color', colors(1, :), 'DisplayName', '70 deg - $u_r$')
|
|
plot(f, sqrt(CPS_10_uh), '-' , 'color', colors(2, :), 'DisplayName', '10 deg - $u_h$')
|
|
plot(f, sqrt(CPS_70_uh), '--', 'color', colors(2, :), 'DisplayName', '70 deg - $u_h$')
|
|
plot(f, sqrt(CPS_10_d), '-' , 'color', colors(3, :), 'DisplayName', '10 deg - $d$')
|
|
plot(f, sqrt(CPS_70_d), '--', 'color', colors(3, :), 'DisplayName', '70 deg - $d$')
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('ASD [$\frac{nrad}{\sqrt{Hz}}$]');
|
|
legend('location', 'northwest');
|
|
xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
figure;
|
|
hold on;
|
|
plot(f, sqrt(pxx_10_ur))
|
|
plot(f, sqrt(pxx_70_ur))
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('ASD [$\frac{nrad}{\sqrt{Hz}}$]');
|
|
legend('location', 'northwest');
|
|
xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
data_70_deg.dz = data_70_deg.allValues(:,4) - mean(data_70_deg.allValues(:,4));
|
|
data_70_deg.dry = data_70_deg.allValues(:,5) - mean(data_70_deg.allValues(:,5));
|
|
data_70_deg.drx = data_70_deg.allValues(:,6) - mean(data_70_deg.allValues(:,6));
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
win = hanning(ceil(1/Ts));
|
|
|
|
[pxx_10_drx, f] = pwelch(data_10_deg.drx, win, [], [], 1/Ts);
|
|
[pxx_70_drx, ~] = pwelch(data_70_deg.drx, win, [], [], 1/Ts);
|
|
|
|
[pxx_10_dry, ~] = pwelch(data_10_deg.dry, win, [], [], 1/Ts);
|
|
[pxx_70_dry, ~] = pwelch(data_70_deg.dry, win, [], [], 1/Ts);
|
|
|
|
[pxx_10_dz, ~] = pwelch(data_10_deg.dz, win, [], [], 1/Ts);
|
|
[pxx_70_dz, ~] = pwelch(data_70_deg.dz, win, [], [], 1/Ts);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
figure;
|
|
hold on;
|
|
plot(f, sqrt(pxx_10_drx))
|
|
plot(f, sqrt(pxx_70_drx))
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('ASD [$\frac{nrad}{\sqrt{Hz}}$]');
|
|
legend('location', 'northwest');
|
|
xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
figure;
|
|
hold on;
|
|
plot(f, sqrt(pxx_10_dry))
|
|
plot(f, sqrt(pxx_70_dry))
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('ASD [$\frac{nrad}{\sqrt{Hz}}$]');
|
|
legend('location', 'northwest');
|
|
xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
figure;
|
|
hold on;
|
|
plot(f, sqrt(pxx_10_dz))
|
|
plot(f, sqrt(pxx_70_dz))
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('ASD [$\frac{nm}{\sqrt{Hz}}$]');
|
|
legend('location', 'northwest');
|
|
xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
|
|
|
|
** Scans
|
|
|
|
| 1 | dz311 |
|
|
| 2 | dry311 |
|
|
| 3 | drx311 |
|
|
| 4 | dz111 |
|
|
| 5 | dry111 |
|
|
| 6 | drx111 |
|
|
| 7 | fjpur |
|
|
| 8 | fjpuh |
|
|
| 9 | fjpd |
|
|
| 10 | bragg |
|
|
|
|
#+begin_src matlab
|
|
%% Load Data
|
|
data_10_70_deg = extractDatData('thtraj_10_70.mat', ...
|
|
{"dz311", "dry311", "drx311", "dz", "dry", "drx", "fjur", "fjuh", "fjd", "bragg"}, ...
|
|
[1e-9, 1e-9, 1e-9, 1e-9, 1e-9, 1e-9, 1e-8, 1e-8, 1e-8, pi/180]);
|
|
|
|
Ts = 1e-4;
|
|
t = Ts*[1:length(data_10_70_deg.bragg)];
|
|
|
|
%% Actuator Jacobian
|
|
J_a_111 = [1, 0.14, -0.0675
|
|
1, 0.14, 0.1525
|
|
1, -0.14, 0.0425];
|
|
|
|
data_10_70_deg.ddz = 10.5e-3./(2*cos(data_10_70_deg.bragg)) - data_10_70_deg.dz;
|
|
|
|
%% Computation of the position of the FJ as measured by the interferometers
|
|
error = J_a_111 * [data_10_70_deg.ddz, data_10_70_deg.dry, data_10_70_deg.drx]';
|
|
|
|
data_10_70_deg.fjur_e = error(1,:)'; % [m]
|
|
data_10_70_deg.fjuh_e = error(2,:)'; % [m]
|
|
data_10_70_deg.fjd_e = error(3,:)'; % [m]
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Load Data
|
|
data_70_10_deg = extractDatData('thtraj_70_10.mat', ...
|
|
{"dz311", "dry311", "drx311", "dz", "dry", "drx", "fjur", "fjuh", "fjd", "bragg"}, ...
|
|
[1e-9, 1e-9, 1e-9, 1e-9, 1e-9, 1e-9, 1e-8, 1e-8, 1e-8, pi/180]);
|
|
%% Actuator Jacobian
|
|
J_a_111 = [1, 0.14, -0.0675
|
|
1, 0.14, 0.1525
|
|
1, -0.14, 0.0425];
|
|
|
|
data_70_10_deg.ddz = 10.5e-3./(2*cos(data_70_10_deg.bragg)) - data_70_10_deg.dz;
|
|
|
|
%% Computation of the position of the FJ as measured by the interferometers
|
|
error = J_a_111 * [data_70_10_deg.ddz, data_70_10_deg.dry, data_70_10_deg.drx]';
|
|
|
|
data_70_10_deg.fjur_e = error(1,:)'; % [m]
|
|
data_70_10_deg.fjuh_e = error(2,:)'; % [m]
|
|
data_70_10_deg.fjd_e = error(3,:)'; % [m]
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
win = hanning(ceil(1/Ts));
|
|
|
|
[pxx_10_70_ur, f] = pwelch(data_10_70_deg.fjur_e(t<100), win, [], [], 1/Ts);
|
|
[pxx_70_10_ur, ~] = pwelch(data_70_10_deg.fjur_e(t<100), win, [], [], 1/Ts);
|
|
|
|
[pxx_10_70_uh, ~] = pwelch(data_10_70_deg.fjuh_e(t<100), win, [], [], 1/Ts);
|
|
[pxx_70_10_uh, ~] = pwelch(data_70_10_deg.fjuh_e(t<100), win, [], [], 1/Ts);
|
|
|
|
[pxx_10_70_d, ~] = pwelch(data_10_70_deg.fjd_e(t<100), win, [], [], 1/Ts);
|
|
[pxx_70_10_d, ~] = pwelch(data_70_10_deg.fjd_e(t<100), win, [], [], 1/Ts);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
CPS_10_70_ur = flip(-cumtrapz(flip(f), flip(pxx_10_70_ur)));
|
|
CPS_10_70_uh = flip(-cumtrapz(flip(f), flip(pxx_10_70_uh)));
|
|
CPS_10_70_d = flip(-cumtrapz(flip(f), flip(pxx_10_70_d)));
|
|
|
|
CPS_70_10_ur = flip(-cumtrapz(flip(f), flip(pxx_70_10_ur)));
|
|
CPS_70_10_uh = flip(-cumtrapz(flip(f), flip(pxx_70_10_uh)));
|
|
CPS_70_10_d = flip(-cumtrapz(flip(f), flip(pxx_70_10_d)));
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Spectrogram
|
|
figure;
|
|
hold on;
|
|
pspectrum(data_10_70_deg.fjuh_e(t<100), 1e4, 'spectrogram', ...
|
|
'FrequencyResolution', 1e0, ...
|
|
'OverlapPercent', 99, ...
|
|
'FrequencyLimits', [1, 400]);
|
|
hold off;
|
|
% xlim([0.03, 1.14]); ylim([1, 400]);
|
|
% caxis([-160, -130])
|
|
title('');
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
figure;
|
|
hold on;
|
|
plot(f, 1e9*sqrt(CPS_10_70_ur))
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('ASD [$\frac{nrad}{\sqrt{Hz}}$]');
|
|
legend('location', 'northwest');
|
|
xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
figure;
|
|
hold on;
|
|
plot(f, sqrt(pxx_10_70_ur))
|
|
plot(f, sqrt(pxx_70_10_ur))
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('ASD [$\frac{nrad}{\sqrt{Hz}}$]');
|
|
legend('location', 'northwest');
|
|
xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
** Noise budgeting - No rotation
|
|
First, we look at the position errors when the bragg axis is not moving
|
|
|
|
#+begin_src matlab
|
|
%% Load Measurement Data
|
|
ol_data = load('FJPUR_step.mat');
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Pre-processing
|
|
ol_time = ol_data.time - ol_data.time(1);
|
|
|
|
ol_drx = ol_data.allValues(ol_time < 45, 6);
|
|
ol_dry = ol_data.allValues(ol_time < 45, 5);
|
|
ol_dz = ol_data.allValues(ol_time < 45, 4);
|
|
|
|
ol_drx = ol_drx - mean(ol_drx);
|
|
ol_dry = ol_dry - mean(ol_dry);
|
|
ol_dz = ol_dz - mean(ol_dz);
|
|
|
|
ol_time = ol_time(ol_time < 45);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
figure;
|
|
plot(ol_time, ol_drx)
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Parameters for Spectral Analysis
|
|
Ts = 1e-4;
|
|
win = hanning(ceil(1/Ts));
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Computation of the PSD
|
|
[pxx_ol_drx, f] = pwelch(ol_drx, win, [], [], 1/Ts);
|
|
[pxx_ol_dry, ~] = pwelch(ol_dry, win, [], [], 1/Ts);
|
|
[pxx_ol_dz, ~] = pwelch(ol_dz, win, [], [], 1/Ts);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Amplitude Spectral Density
|
|
figure;
|
|
hold on;
|
|
plot(f, sqrt(pxx_ol_drx), 'DisplayName', '$R_x$');
|
|
plot(f, sqrt(pxx_ol_dry), 'DisplayName', '$R_y$');
|
|
plot(f, sqrt(pxx_ol_dz), 'DisplayName', '$D_z$');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('ASD [$\frac{nm,nrad}{\sqrt{Hz}}$]');
|
|
legend('location', 'northwest');
|
|
xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :tangle no :exports results :results file replace
|
|
exportFig('figs/noise_budget_no_mov_asd.pdf', 'width', 'wide', 'height', 'normal');
|
|
#+end_src
|
|
|
|
#+name: fig:noise_budget_no_mov_asd
|
|
#+caption: Amplitude Spectral Density
|
|
#+RESULTS:
|
|
[[file:figs/noise_budget_no_mov_asd.png]]
|
|
|
|
#+begin_src matlab
|
|
%% Cumulative Spectral Density
|
|
CPS_drx = flip(-cumtrapz(flip(f), flip(pxx_ol_drx)));
|
|
CPS_dry = flip(-cumtrapz(flip(f), flip(pxx_ol_dry)));
|
|
CPS_dz = flip(-cumtrapz(flip(f), flip(pxx_ol_dz)));
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Cumulative Spectral Density
|
|
CPS_drx = cumtrapz(f, pxx_ol_drx);
|
|
CPS_dry = cumtrapz(f, pxx_ol_dry);
|
|
CPS_dz = cumtrapz(f, pxx_ol_dz);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Cumulative Spectral Density
|
|
figure;
|
|
hold on;
|
|
plot(f, sqrt(CPS_drx), 'DisplayName', '$R_x$');
|
|
plot(f, sqrt(CPS_dry), 'DisplayName', '$R_y$');
|
|
plot(f, sqrt(CPS_dz), 'DisplayName', '$D_z$');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('CAS [nm, nrad rms]');
|
|
legend('location', 'northwest');
|
|
xlim([1, 1e3]);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :tangle no :exports results :results file replace
|
|
exportFig('figs/noise_budget_no_mov_cas.pdf', 'width', 'wide', 'height', 'normal');
|
|
#+end_src
|
|
|
|
#+name: fig:noise_budget_no_mov_cas
|
|
#+caption: Cumulative Amplitude Spectrum
|
|
#+RESULTS:
|
|
[[file:figs/noise_budget_no_mov_cas.png]]
|
|
|
|
|
|
** TODO Noise budgeting - Bragg rotation
|
|
|
|
* Test Mode C
|
|
** 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 :tangle no :noweb yes
|
|
<<m-init-path>>
|
|
#+end_src
|
|
|
|
#+begin_src matlab :eval no :noweb yes
|
|
<<m-init-path-tangle>>
|
|
#+end_src
|
|
|
|
#+begin_src matlab :noweb yes
|
|
<<m-init-other>>
|
|
#+end_src
|
|
|
|
** Mode B and Mode C
|
|
#+begin_src matlab
|
|
data_B = extractDatData(sprintf("/home/thomas/mnt/data_id21/22Jan/blc13491/id21/test_regul_220119/%s","lut_const_fj_vel_19012022_1450.dat"), ...
|
|
{"bragg", "dz", "dry", "drx", "fjur", "fjuh", "fjd"}, ...
|
|
[pi/180, 1e-9, 1e-9, 1e-9, 1e-8, 1e-8, 1e-8]);
|
|
data_B = processMeasData(data_B);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
data_C = extractDatData(sprintf("/home/thomas/mnt/data_id21/22Jan/blc13491/id21/test_regul_220119/%s","lut_const_fj_vel_19012022_1454.dat"), ...
|
|
{"bragg", "dz", "dry", "drx", "fjur", "fjuh", "fjd"}, ...
|
|
[pi/180, 1e-9, 1e-9, 1e-9, 1e-8, 1e-8, 1e-8]);
|
|
data_C = processMeasData(data_C);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
figure;
|
|
hold on;
|
|
plot(180/pi*data_B.bragg, 1e9*data_B.drx)
|
|
hold off;
|
|
xlabel('Bragg Angle [deg]'); ylabel('DRX [nrad]');
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
figure;
|
|
hold on;
|
|
plot(180/pi*data_B.bragg, 1e9*data_B.fjur_e_filt)
|
|
plot(180/pi*data_C.bragg, 1e9*data_C.fjur_e_filt)
|
|
hold off;
|
|
xlabel('Bragg Angle [deg]'); ylabel('FJUR Error [nm]');
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
figure;
|
|
hold on;
|
|
plot(180/pi*data_B.bragg, 1e9*data_B.fjur_e)
|
|
plot(180/pi*data_C.bragg, 1e9*data_C.fjur_e)
|
|
hold off;
|
|
xlabel('Bragg Angle [deg]'); ylabel('FJUR Error [nm]');
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% FIR Filter
|
|
Fs = 1e4; % Sampling Frequency [Hz]
|
|
fir_order = 5000; % Filter's order
|
|
delay = fir_order/2; % Delay induced by the filter
|
|
B_fir = firls(5000, ... % Filter's order
|
|
[0 5/(Fs/2) 10/(Fs/2) 1], ... % Frequencies [Hz]
|
|
[1 1 0 0]); % Wanted Magnitudes
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Filtering all measured Fast Jack Position using the FIR filter
|
|
data_B.drx_filter = filter(B_fir, 1, data_B.drx);
|
|
data_B.drx_filter(1:end-delay) = data_B.drx_filter(delay+1:end);
|
|
|
|
data_C.drx_filter = filter(B_fir, 1, data_C.drx);
|
|
data_C.drx_filter(1:end-delay) = data_C.drx_filter(delay+1:end);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
figure;
|
|
hold on;
|
|
plot(180/pi*data_B.bragg, 1e9*data_B.drx_filter)
|
|
plot(180/pi*data_C.bragg, 1e9*data_C.drx_filter)
|
|
hold off;
|
|
xlabel('Bragg Angle [deg]'); ylabel('DRX [nrad]');
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
Ts = 1e-4;
|
|
win = hanning(ceil(1/Ts));
|
|
|
|
[pxx_B_drx, f] = pwelch(data_B.drx, win, [], [], 1/Ts);
|
|
[pxx_C_drx, ~] = pwelch(data_C.drx, win, [], [], 1/Ts);
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Estimation of Sensitivity function
|
|
figure;
|
|
hold on;
|
|
plot(f, sqrt(pxx_B_drx), 'DisplayName', 'B');
|
|
plot(f, sqrt(pxx_C_drx), 'DisplayName', 'C');
|
|
hold off;
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('ASD');
|
|
legend('location', 'northwest');
|
|
#+end_src
|
|
|
|
#+begin_src matlab :exports none
|
|
%% Estimation of Sensitivity function
|
|
figure;
|
|
plot(f, sqrt(pxx_C_drx)./sqrt(pxx_B_drx), 'k-');
|
|
set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
|
|
xlabel('Frequency [Hz]'); ylabel('Sensitivity Function');
|
|
legend('location', 'northwest');
|
|
#+end_src
|
|
|
|
* Export numerator and denominator
|
|
** 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 :tangle no :noweb yes
|
|
<<m-init-path>>
|
|
#+end_src
|
|
|
|
#+begin_src matlab :eval no :noweb yes
|
|
<<m-init-path-tangle>>
|
|
#+end_src
|
|
|
|
#+begin_src matlab :noweb yes
|
|
<<m-init-other>>
|
|
#+end_src
|
|
|
|
** Export
|
|
#+begin_src matlab
|
|
K_order = 10;
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
load('X_tal_cage_PID.mat', 'K');
|
|
K_order = order(K(1,1));
|
|
Kz = c2d(K(1,1)*(1 + s/2/pi/2e3)^(9-K_order)/(1 + s/2/pi/2e3)^(9-K_order), 1e-4);
|
|
[num, den] = tfdata(Kz, 'v');
|
|
|
|
formatSpec = '%.18e %.18e %.18e %.18e %.18e %.18e %.18e %.18e %.18e %.18e\n';
|
|
fileID = fopen('X_tal_cage_PID.dat', 'w');
|
|
fprintf(fileID, formatSpec, [num; den]');
|
|
fclose(fileID);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
load('X_tal_cage_PID_20Hz.mat', 'K');
|
|
K_order = order(K(1,1));
|
|
Kz = c2d(K(1,1)*(1 + s/2/pi/2e3)^(9-K_order)/(1 + s/2/pi/2e3)^(9-K_order), 1e-4);
|
|
[num, den] = tfdata(Kz, 'v');
|
|
|
|
formatSpec = '%.18e %.18e %.18e %.18e %.18e %.18e %.18e %.18e %.18e %.18e\n';
|
|
fileID = fopen('X_tal_cage_PID_20Hz.dat', 'w');
|
|
fprintf(fileID, formatSpec, [num; den]');
|
|
fclose(fileID);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
load('X_tal_cage_PID_40Hz.mat', 'K');
|
|
K_order = order(K(1,1));
|
|
Kz = c2d(K(1,1)*(1 + s/2/pi/2e3)^(9-K_order)/(1 + s/2/pi/2e3)^(9-K_order), 1e-4);
|
|
[num, den] = tfdata(Kz, 'v');
|
|
|
|
formatSpec = '%.18e %.18e %.18e %.18e %.18e %.18e %.18e %.18e %.18e %.18e\n';
|
|
fileID = fopen('X_tal_cage_PID_40Hz.dat', 'w');
|
|
fprintf(fileID, formatSpec, [num; den]');
|
|
fclose(fileID);
|
|
#+end_src
|
|
|
|
** Verify
|
|
#+begin_src matlab
|
|
K_data = importdata('X_tal_cage_PID_20Hz.dat');
|
|
K = tf(K_data(1,:), K_data(2,:), 1e-4)
|
|
#+end_src
|
|
|
|
* Helping Functions :noexport:
|
|
** Initialize Path
|
|
#+NAME: m-init-path
|
|
#+BEGIN_SRC matlab
|
|
%% Path for functions, data and scripts
|
|
addpath('./matlab/mat/'); % Path for data
|
|
addpath('./matlab/src/'); % Path for functions
|
|
addpath('./matlab/'); % Path for scripts
|
|
#+END_SRC
|
|
|
|
#+NAME: m-init-path-tangle
|
|
#+BEGIN_SRC matlab
|
|
%% Path for functions, data and scripts
|
|
addpath('./mat/'); % Path for data
|
|
addpath('./src/'); % Path for functions
|
|
#+END_SRC
|
|
|
|
** Initialize Simscape Model
|
|
#+NAME: m-init-simscape
|
|
#+begin_src matlab
|
|
#+end_src
|
|
|
|
** Initialize other elements
|
|
#+NAME: m-init-other
|
|
#+BEGIN_SRC matlab
|
|
%% Colors for the figures
|
|
colors = colororder;
|
|
|
|
%% Frequency Vector
|
|
freqs = logspace(0, 3, 1000);
|
|
#+END_SRC
|
|
|
|
** =extractDatData=
|
|
:PROPERTIES:
|
|
:header-args:matlab+: :tangle matlab/src/extractDatData.m
|
|
:header-args:matlab+: :comments none :mkdirp yes :eval no
|
|
:END:
|
|
<<sec:extractDatData>>
|
|
|
|
*** Function Description
|
|
#+begin_src matlab
|
|
function [data_struct] = extractDatData(dat_file, names, scale)
|
|
% extractDatData -
|
|
%
|
|
% Syntax: extractDatData(data_file, lut_file, args)
|
|
%
|
|
% Inputs:
|
|
% - data_file - Where to load the .mat file
|
|
% - lut_file - Where to save the .dat file
|
|
#+end_src
|
|
|
|
*** Load Data
|
|
#+begin_src matlab
|
|
%% Load Data
|
|
data_array = importdata(dat_file);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Initialize Struct
|
|
data_struct = struct();
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Populate Struct
|
|
for i = 1:length(names)
|
|
data_struct.(names{i}) = scale(i)*data_array(:,i);
|
|
end
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Add Time
|
|
data_struct.time = 1e-4*[1:1:length(data_struct.(names{1}))];
|
|
#+end_src
|
|
|
|
** =processMeasData=
|
|
:PROPERTIES:
|
|
:header-args:matlab+: :tangle matlab/src/processMeasData.m
|
|
:header-args:matlab+: :comments none :mkdirp yes :eval no
|
|
:END:
|
|
<<sec:processMeasData>>
|
|
|
|
*** Function Description
|
|
#+begin_src matlab
|
|
function [data] = processMeasData(data, args)
|
|
% processMeasData -
|
|
%
|
|
% Syntax: processMeasData(data_file, lut_file, args)
|
|
%
|
|
% Inputs:
|
|
% - data
|
|
#+end_src
|
|
|
|
*** Optional Parameters
|
|
#+begin_src matlab
|
|
arguments
|
|
data
|
|
args.fir_order (1,1) double {mustBeNumericOrLogical} = 5000
|
|
end
|
|
#+end_src
|
|
|
|
*** Process Data
|
|
#+begin_src matlab
|
|
%% Actuator Jacobian
|
|
J_a_111 = [1, 0.14, -0.0675
|
|
1, 0.14, 0.1525
|
|
1, -0.14, 0.0425];
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% FIR Filter
|
|
Fs = 1e4; % Sampling Frequency [Hz]
|
|
fir_order = args.fir_order; % Filter's order
|
|
delay = fir_order/2; % Delay induced by the filter
|
|
B_fir = firls(args.fir_order, ... % Filter's order
|
|
[0 25/(Fs/2) 34/(Fs/2) 1], ... % Frequencies [Hz]
|
|
[1 1 0 0]); % Wanted Magnitudes
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
data.ddz = 10.5e-3./(2*cos(data.bragg)) - data.dz;
|
|
|
|
%% Computation of the position of the FJ as measured by the interferometers
|
|
error = J_a_111 * [data.ddz, data.dry, data.drx]';
|
|
|
|
data.fjur_e = error(1,:)'; % [m]
|
|
data.fjuh_e = error(2,:)'; % [m]
|
|
data.fjd_e = error(3,:)'; % [m]
|
|
|
|
%% Filtering all measured Fast Jack Position using the FIR filter
|
|
data.fjur_e_filt = filter(B_fir, 1, data.fjur_e);
|
|
data.fjuh_e_filt = filter(B_fir, 1, data.fjuh_e);
|
|
data.fjd_e_filt = filter(B_fir, 1, data.fjd_e);
|
|
|
|
%% Compensation of the delay introduced by the FIR filter
|
|
data.fjur_e_filt(1:end-delay) = data.fjur_e_filt(delay+1:end);
|
|
data.fjuh_e_filt(1:end-delay) = data.fjuh_e_filt(delay+1:end);
|
|
data.fjd_e_filt( 1:end-delay) = data.fjd_e_filt( delay+1:end);
|
|
#+end_src
|
|
|
|
#+begin_src matlab
|
|
%% Re-sample data to have same data points in FJUR
|
|
[data.fjur_e_resampl, data.fjur_resampl] = resample(data.fjur_e_filt, data.fjur, 1/100e-9);
|
|
[data.fjuh_e_resampl, data.fjuh_resampl] = resample(data.fjuh_e_filt, data.fjuh, 1/100e-9);
|
|
[data.fjd_e_resampl, data.fjd_resampl] = resample(data.fjd_e_filt, data.fjd, 1/100e-9);
|
|
#+end_src
|
|
|
|
* Bibliography :ignore:
|
|
#+latex: \printbibliography
|