Reworked the ground motion page

ground-motion/matlab/ground_meas_id31.m

@@ -0,0 +1,152 @@
+%% Clear Workspace and Close figures
+clear; close all; clc;
+
+%% Intialize Laplace variable
+s = zpk('s');
+
+% Load data
+
+data = load('mat/data_028.mat', 'data'); data = data.data;
+
+% Time domain plots
+
+figure;
+hold on;
+plot(data(:, 3), data(:, 1));
+hold off;
+xlabel('Time [s]'); ylabel('Voltage [V]');
+xlim([0, 100]);
+
+% Computation of the ASD of the measured voltage
+
+dt = data(2, 3) - data(1, 3);
+
+Fs = 1/dt;
+win = hanning(ceil(10*Fs));
+
+[px_dc, f] = pwelch(data(:, 1), win, [], [], Fs);
+
+figure;
+hold on;
+plot(f, sqrt(px_dc));
+hold off;
+set(gca, 'xscale', 'log');
+set(gca, 'yscale', 'log');
+xlabel('Frequency [Hz]'); ylabel('Amplitude Spectral Density $\left[\frac{V}{\sqrt{Hz}}\right]$')
+xlim([0.1, 500]);
+
+% Scaling to take into account the sensibility of the geophone and the voltage amplifier
+% The Geophone used are L22. Their sensibility is shown on figure [[fig:geophone_sensibility]].
+
+
+S0 = 88; % Sensitivity [V/(m/s)]
+f0 = 2; % Cut-off frequency [Hz]
+
+S = S0*(s/2/pi/f0)/(1+s/2/pi/f0);
+
+figure;
+bodeFig({S}, logspace(-1, 2, 1000));
+ylabel('Amplitude $\left[\frac{V}{m/s}\right]$')
+
+
+
+% #+NAME: fig:geophone_sensibility
+% #+CAPTION: Sensibility of the Geophone
+% #+RESULTS: fig:geophone_sensibility
+% [[file:figs/geophone_sensibility.png]]
+
+% We also take into account the gain of the electronics which is here set to be $60dB$.
+
+
+G0_db = 60; % [dB]
+
+G0 = 10^(G0_db/20); % [abs]
+
+
+
+% We divide the ASD measured (in $\text{V}/\sqrt{\text{Hz}}$) by the gain of the voltage amplifier to obtain the ASD of the voltage across the geophone.
+% We further divide the result by the sensibility of the Geophone to obtain the ASD of the velocity in $m/s/\sqrt{Hz}$.
+
+
+scaling = 1./squeeze(abs(freqresp(G0*S, f, 'Hz')));
+
+% Computation of the ASD of the velocity
+% The ASD of the measured velocity is shown on figure [[fig:ground_motion_id31_asd_velocity]].
+
+
+figure;
+hold on;
+plot(f, sqrt(px_dc).*scaling);
+hold off;
+set(gca, 'xscale', 'log');
+set(gca, 'yscale', 'log');
+xlabel('Frequency [Hz]'); ylabel('ASD of the measured Velocity $\left[\frac{m/s}{\sqrt{Hz}}\right]$')
+xlim([0.1, 500]);
+
+
+
+% #+NAME: fig:ground_motion_id31_asd_velocity
+% #+CAPTION: Amplitude Spectral Density of the Velocity
+% #+RESULTS: fig:ground_motion_id31_asd_velocity
+% [[file:figs/ground_motion_id31_asd_velocity.png]]
+
+% We also plot the ASD in displacement (figure [[fig:ground_motion_id31_asd_displacement]]);
+
+
+figure;
+hold on;
+plot(f, (sqrt(px_dc).*scaling)./(2*pi*f));
+hold off;
+set(gca, 'xscale', 'log'); set(gca, 'yscale', 'log');
+xlabel('Frequency [Hz]'); ylabel('ASD of the displacement $\left[\frac{m}{\sqrt{Hz}}\right]$')
+xlim([0.1, 500]);
+
+
+
+% #+NAME: fig:ground_motion_id31_asd_displacement
+% #+CAPTION: Amplitude Spectral Density of the Displacement
+% #+RESULTS: fig:ground_motion_id31_asd_displacement
+% [[file:figs/ground_motion_id31_asd_displacement.png]]
+% And also in $\frac{{\mu u}^2}{Hz}$ (figure [[fig:ground_motion_id31_psd_displacement]]).
+
+figure;
+hold on;
+plot(f, ((sqrt(px_dc).*scaling)./(2*pi*f).*1e6).^2);
+hold off;
+set(gca, 'xscale', 'log');
+set(gca, 'yscale', 'log');
+xlabel('Frequency [Hz]'); ylabel('PSD of the measured displacement $\left[\frac{{ \mu m }^2}{Hz}\right]$')
+xlim([0.1, 500]);
+
+% Load the measurement data
+% First we load the measurement data.
+% Here we have one measurement of the floor motion made at the ESRF in 2018, and one measurement made at CERN.
+
+id09 = load('./mat/id09_floor_september2018.mat');
+cern = load('./mat/ground_motion_dist.mat');
+
+% Compute PSD of the measurements
+% We compute the Power Spectral Densities of the measurements.
+
+Fs_id09 = 1/(id09.time3(2)-id09.time3(1));
+win_id09 = hanning(ceil(10*Fs_id09));
+[id09_pxx, id09_f] = pwelch(1e-6*id09.x_y_z(:, 3), win_id09, [], [], Fs_id09);
+
+Fs_cern = 1/(cern.gm.time(2)-cern.gm.time(1));
+win_cern = hanning(ceil(10*Fs_cern));
+[cern_pxx, cern_f] = pwelch(cern.gm.signal, win_cern, [], [], Fs_cern);
+
+% Compare PSD of Cern, ID09 and ID31
+% And we compare all the measurements (figure [[fig:ground_motion_compare]]).
+
+
+figure;
+hold on;
+plot(id09_f, id09_pxx, 'DisplayName', 'ID09');
+plot(cern_f, cern_pxx, 'DisplayName', 'CERN');
+plot(f, ((sqrt(px_dc).*scaling)./(2*pi*f)).^2, 'k', 'DisplayName', 'ID31');
+hold off;
+set(gca, 'XScale', 'log'); set(gca, 'YScale', 'log');
+xlabel('Frequency [Hz]'); ylabel('PSD [$m^2/Hz$]');
+legend('Location', 'northeast');
+xlim([0.1, 500]);