Add all data files

matlab/src/circlefit.m

@@ -0,0 +1,21 @@
+function   [xc,yc,R,a] = circlefit(x,y)
+%
+%   [xc yx R] = circfit(x,y)
+%
+%   fits a circle  in x,y plane in a more accurate
+%   (less prone to ill condition )
+%  procedure than circfit2 but using more memory
+%  x,y are column vector where (x(i),y(i)) is a measured point
+%
+%  result is center point (yc,xc) and radius R
+%  an optional output is the vector of coeficient a
+% describing the circle's equation
+%
+%   x^2+y^2+a(1)*x+a(2)*y+a(3)=0
+%
+%  By:  Izhak bucher 25/oct /1991,
+    x=x(:); y=y(:);
+    a=[x y ones(size(x))]\[-(x.^2+y.^2)];
+    xc = -.5*a(1);
+    yc = -.5*a(2);
+    R  =  sqrt((a(1)^2+a(2)^2)/4-a(3));

matlab/src/initializeDisturbances.m

@@ -10,15 +10,15 @@ arguments
     % Global parameter to enable or disable the disturbances
     args.enable logical {mustBeNumericOrLogical} = true
     % Ground Motion - X direction
-    args.Dwx logical {mustBeNumericOrLogical} = true
+    args.Dw_x logical {mustBeNumericOrLogical} = true
     % Ground Motion - Y direction
-    args.Dwy logical {mustBeNumericOrLogical} = true
+    args.Dw_y logical {mustBeNumericOrLogical} = true
     % Ground Motion - Z direction
-    args.Dwz logical {mustBeNumericOrLogical} = true
+    args.Dw_z logical {mustBeNumericOrLogical} = true
     % Translation Stage - X direction
-    args.Fty_x logical {mustBeNumericOrLogical} = true
+    args.Fdy_x logical {mustBeNumericOrLogical} = true
     % Translation Stage - Z direction
-    args.Fty_z logical {mustBeNumericOrLogical} = true
+    args.Fdy_z logical {mustBeNumericOrLogical} = true
     % Spindle - X direction
     args.Frz_x logical {mustBeNumericOrLogical} = true
     % Spindle - Y direction
@@ -31,175 +31,163 @@ end
 rng("shuffle");
 
 %% Ground Motion
-load('dist_psd.mat', 'dist_f');
+if args.enable
+    % Load the PSD of disturbance
+    load('ustation_disturbance_psd.mat', 'gm_dist')
 
-% Frequency Data
-Dw.f = dist_f.f(2:end);
-Dw.psd_x = dist_f.psd_gm(2:end);
-Dw.psd_y = dist_f.psd_gm(2:end);
-Dw.psd_z = dist_f.psd_gm(2:end);
+    % Frequency Data
+    Dw.f = gm_dist.f;
+    Dw.psd_x = gm_dist.pxx_x;
+    Dw.psd_y = gm_dist.pxx_y;
+    Dw.psd_z = gm_dist.pxx_z;
 
-% Time data
-Fs = 2*Dw.f(end);      % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz]
-N  = 2*length(Dw.f);   % Number of Samples match the one of the wanted PSD
-T0 = N/Fs;             % Signal Duration [s]
-Dw.t = linspace(0, T0, N+1)'; % Time Vector [s]
+    % Time data
+    Fs = 2*Dw.f(end);      % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz]
+    N  = 2*length(Dw.f);   % Number of Samples match the one of the wanted PSD
+    T0 = N/Fs;             % Signal Duration [s]
+    Dw.t = linspace(0, T0, N+1)'; % Time Vector [s]
 
-C = zeros(N/2,1);
-for i = 1:N/2
-    C(i) = sqrt(Dw.psd_x(i)/T0);
-end
-
-if args.Dwx && args.enable
-    theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
-    Cx = [0 ; C.*complex(cos(theta),sin(theta))];
-    Cx = [Cx; flipud(conj(Cx(2:end)))];;
-    Dw.x = N/sqrt(2)*ifft(Cx); % Ground Motion - x direction [m]
-else
-    Dw.x = zeros(length(Dw.t), 1);
-end
-
-if args.Dwy && args.enable
-    theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
-    Cx = [0 ; C.*complex(cos(theta),sin(theta))];
-    Cx = [Cx; flipud(conj(Cx(2:end)))];;
-    Dw.y = N/sqrt(2)*ifft(Cx); % Ground Motion - y direction [m]
-else
-    Dw.y = zeros(length(Dw.t), 1);
-end
-
-if args.Dwy && args.enable
-    theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
-    Cx = [0 ; C.*complex(cos(theta),sin(theta))];
-    Cx = [Cx; flipud(conj(Cx(2:end)))];;
-    Dw.z = N/sqrt(2)*ifft(Cx); % Ground Motion - z direction [m]
-else
-    Dw.z = zeros(length(Dw.t), 1);
-end
-
-load('dist_psd.mat', 'dist_f');
-dist_f.f = dist_f.f(2:end);
-dist_f.psd_gm = dist_f.psd_gm(2:end);
-dist_f.psd_ty = dist_f.psd_ty(2:end);
-dist_f.psd_rz = dist_f.psd_rz(2:end);
-
-%% Translation Stage
-load('dist_psd.mat', 'dist_f');
-
-% Frequency Data
-Ty.f = dist_f.f(2:end);
-Ty.psd_x = dist_f.psd_ty(2:end); % TODO - we take here the vertical direction which is wrong but approximate
-Ty.psd_z = dist_f.psd_ty(2:end);
-
-% Time data
-Fs = 2*Ty.f(end);      % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz]
-N  = 2*length(Ty.f);   % Number of Samples match the one of the wanted PSD
-T0 = N/Fs;             % Signal Duration [s]
-Ty.t = linspace(0, T0, N+1)'; % Time Vector [s]
-
-C = zeros(N/2,1);
-for i = 1:N/2
-    C(i) = sqrt(Ty.psd_x(i)/T0);
-end
-
-% Translation Stage - X
-if args.Fty_x && args.enable
-    phi = Ty.psd_x;
+    % ASD representation of the ground motion
     C = zeros(N/2,1);
     for i = 1:N/2
-        C(i) = sqrt(phi(i)/T0);
+        C(i) = sqrt(Dw.psd_x(i)/T0);
     end
-    theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
-    Cx = [0 ; C.*complex(cos(theta),sin(theta))];
-    Cx = [Cx; flipud(conj(Cx(2:end)))];;
-    u = N/sqrt(2)*ifft(Cx); % Disturbance Force Ty x [N]
-    Ty.x = u;
+
+    if args.Dw_x
+        theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
+        Cx = [0 ; C.*complex(cos(theta),sin(theta))];
+        Cx = [Cx; flipud(conj(Cx(2:end)))];;
+        Dw.x = N/sqrt(2)*ifft(Cx); % Ground Motion - x direction [m]
+    else
+        Dw.x = zeros(length(Dw.t), 1);
+    end
+
+    if args.Dw_y
+        theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
+        Cx = [0 ; C.*complex(cos(theta),sin(theta))];
+        Cx = [Cx; flipud(conj(Cx(2:end)))];;
+        Dw.y = N/sqrt(2)*ifft(Cx); % Ground Motion - y direction [m]
+    else
+        Dw.y = zeros(length(Dw.t), 1);
+    end
+
+    if args.Dw_y
+        theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
+        Cx = [0 ; C.*complex(cos(theta),sin(theta))];
+        Cx = [Cx; flipud(conj(Cx(2:end)))];;
+        Dw.z = N/sqrt(2)*ifft(Cx); % Ground Motion - z direction [m]
+    else
+        Dw.z = zeros(length(Dw.t), 1);
+    end
+
 else
-    Ty.x = zeros(length(Ty.t), 1);
+    Dw.t = [0,1]; % Time Vector [s]
+    Dw.x = [0,0]; % Ground Motion - X [m]
+    Dw.y = [0,0]; % Ground Motion - Y [m]
+    Dw.z = [0,0]; % Ground Motion - Z [m]
 end
 
-% Translation Stage - Z
-if args.Fty_z && args.enable
-    phi = Ty.psd_z;
+%% Translation stage
+if args.enable
+    % Load the PSD of disturbance
+    load('ustation_disturbance_psd.mat', 'dy_dist')
+
+    % Frequency Data
+    Dy.f = dy_dist.f;
+    Dy.psd_x = dy_dist.pxx_fx;
+    Dy.psd_z = dy_dist.pxx_fz;
+
+    % Time data
+    Fs = 2*Dy.f(end);      % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz]
+    N  = 2*length(Dy.f);   % Number of Samples match the one of the wanted PSD
+    T0 = N/Fs;             % Signal Duration [s]
+    Dy.t = linspace(0, T0, N+1)'; % Time Vector [s]
+
+    % ASD representation of the disturbance voice
     C = zeros(N/2,1);
     for i = 1:N/2
-        C(i) = sqrt(phi(i)/T0);
+        C(i) = sqrt(Dy.psd_x(i)/T0);
     end
-    theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
-    Cx = [0 ; C.*complex(cos(theta),sin(theta))];
-    Cx = [Cx; flipud(conj(Cx(2:end)))];;
-    u = N/sqrt(2)*ifft(Cx); % Disturbance Force Ty z [N]
-    Ty.z = u;
+
+    if args.Fdy_x
+        theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
+        Cx = [0 ; C.*complex(cos(theta),sin(theta))];
+        Cx = [Cx; flipud(conj(Cx(2:end)))];;
+        Dy.x = N/sqrt(2)*ifft(Cx); % Translation stage disturbances - X direction [N]
+    else
+        Dy.x = zeros(length(Dy.t), 1);
+    end
+
+    if args.Fdy_z
+        theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
+        Cx = [0 ; C.*complex(cos(theta),sin(theta))];
+        Cx = [Cx; flipud(conj(Cx(2:end)))];;
+        Dy.z = N/sqrt(2)*ifft(Cx); % Translation stage disturbances - Z direction [N]
+    else
+        Dy.z = zeros(length(Dy.t), 1);
+    end
+
 else
-    Ty.z = zeros(length(Ty.t), 1);
+    Dy.t = [0,1]; % Time Vector [s]
+    Dy.x = [0,0]; % Translation Stage disturbances - X [N]
+    Dy.z = [0,0]; % Translation Stage disturbances - Z [N]
 end
 
-%% Translation Stage
-load('dist_psd.mat', 'dist_f');
+%% Spindle
+if args.enable
+    % Load the PSD of disturbance
+    load('ustation_disturbance_psd.mat', 'rz_dist')
 
-% Frequency Data
-Rz.f = dist_f.f(2:end);
-Rz.psd_x = dist_f.psd_rz(2:end); % TODO - we take here the vertical direction which is wrong but approximate
-Rz.psd_y = dist_f.psd_rz(2:end); % TODO - we take here the vertical direction which is wrong but approximate
-Rz.psd_z = dist_f.psd_rz(2:end);
+    % Frequency Data
+    Rz.f = rz_dist.f;
+    Rz.psd_x = rz_dist.pxx_fx;
+    Rz.psd_y = rz_dist.pxx_fy;
+    Rz.psd_z = rz_dist.pxx_fz;
 
-% Time data
-Fs = 2*Rz.f(end);      % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz]
-N  = 2*length(Rz.f);   % Number of Samples match the one of the wanted PSD
-T0 = N/Fs;             % Signal Duration [s]
-Rz.t = linspace(0, T0, N+1)'; % Time Vector [s]
+    % Time data
+    Fs = 2*Rz.f(end);      % Sampling Frequency of data is twice the maximum frequency of the PSD vector [Hz]
+    N  = 2*length(Rz.f);   % Number of Samples match the one of the wanted PSD
+    T0 = N/Fs;             % Signal Duration [s]
+    Rz.t = linspace(0, T0, N+1)'; % Time Vector [s]
 
-C = zeros(N/2,1);
-for i = 1:N/2
-    C(i) = sqrt(Rz.psd_x(i)/T0);
-end
-
-% Translation Stage - X
-if args.Frz_x && args.enable
-    phi = Rz.psd_x;
+    % ASD representation of the disturbance voice
     C = zeros(N/2,1);
     for i = 1:N/2
-        C(i) = sqrt(phi(i)/T0);
+        C(i) = sqrt(Rz.psd_x(i)/T0);
     end
-    theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
-    Cx = [0 ; C.*complex(cos(theta),sin(theta))];
-    Cx = [Cx; flipud(conj(Cx(2:end)))];;
-    u = N/sqrt(2)*ifft(Cx); % Disturbance Force Rz x [N]
-    Rz.x = u;
-else
-    Rz.x = zeros(length(Rz.t), 1);
-end
 
-% Translation Stage - Y
-if args.Frz_y && args.enable
-    phi = Rz.psd_y;
-    C = zeros(N/2,1);
-    for i = 1:N/2
-        C(i) = sqrt(phi(i)/T0);
+    if args.Frz_x
+        theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
+        Cx = [0 ; C.*complex(cos(theta),sin(theta))];
+        Cx = [Cx; flipud(conj(Cx(2:end)))];;
+        Rz.x = N/sqrt(2)*ifft(Cx); % spindle disturbances - X direction [N]
+    else
+        Rz.x = zeros(length(Rz.t), 1);
     end
-    theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
-    Cx = [0 ; C.*complex(cos(theta),sin(theta))];
-    Cx = [Cx; flipud(conj(Cx(2:end)))];;
-    u = N/sqrt(2)*ifft(Cx); % Disturbance Force Rz y [N]
-    Rz.y = u;
-else
-    Rz.y = zeros(length(Rz.t), 1);
-end
 
-% Translation Stage - Z
-if args.Frz_z && args.enable
-    phi = Rz.psd_z;
-    C = zeros(N/2,1);
-    for i = 1:N/2
-        C(i) = sqrt(phi(i)/T0);
+    if args.Frz_y
+        theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
+        Cx = [0 ; C.*complex(cos(theta),sin(theta))];
+        Cx = [Cx; flipud(conj(Cx(2:end)))];;
+        Rz.y = N/sqrt(2)*ifft(Cx); % spindle disturbances - Y direction [N]
+    else
+        Rz.y = zeros(length(Rz.t), 1);
     end
-    theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
-    Cx = [0 ; C.*complex(cos(theta),sin(theta))];
-    Cx = [Cx; flipud(conj(Cx(2:end)))];;
-    u = N/sqrt(2)*ifft(Cx); % Disturbance Force Rz z [N]
-    Rz.z = u;
+
+    if args.Frz_z
+        theta = 2*pi*rand(N/2,1); % Generate random phase [rad]
+        Cx = [0 ; C.*complex(cos(theta),sin(theta))];
+        Cx = [Cx; flipud(conj(Cx(2:end)))];;
+        Rz.z = N/sqrt(2)*ifft(Cx); % spindle disturbances - Z direction [N]
+    else
+        Rz.z = zeros(length(Rz.t), 1);
+    end
+
 else
-    Rz.z = zeros(length(Rz.t), 1);
+    Rz.t = [0,1]; % Time Vector [s]
+    Rz.x = [0,0]; % Spindle disturbances - X [N]
+    Rz.y = [0,0]; % Spindle disturbances - X [N]
+    Rz.z = [0,0]; % Spindle disturbances - Z [N]
 end
 
 u = zeros(100, 6);
@@ -208,22 +196,16 @@ Fd = u;
 Dw.x   = Dw.x   - Dw.x(1);
 Dw.y   = Dw.y   - Dw.y(1);
 Dw.z   = Dw.z   - Dw.z(1);
-Ty.x   = Ty.x   - Ty.x(1);
-Ty.z   = Ty.z   - Ty.z(1);
+
+Dy.x   = Dy.x   - Dy.x(1);
+Dy.z   = Dy.z   - Dy.z(1);
+
 Rz.x   = Rz.x   - Rz.x(1);
 Rz.y   = Rz.y   - Rz.y(1);
 Rz.z   = Rz.z   - Rz.z(1);
 
 if exist('./mat', 'dir')
-    if exist('./mat/nass_disturbances.mat', 'file')
-        save('mat/nass_disturbances.mat', 'Dw', 'Ty', 'Rz', 'Fd', 'args', '-append');
-    else
-        save('mat/nass_disturbances.mat', 'Dw', 'Ty', 'Rz', 'Fd', 'args');
-    end
+    save('mat/nass_disturbances.mat', 'Dw', 'Dy', 'Rz', 'Fd', 'args');
 elseif exist('./matlab', 'dir')
-    if exist('./matlab/mat/nass_disturbances.mat', 'file')
-        save('matlab/mat/nass_disturbances.mat', 'Dw', 'Ty', 'Rz', 'Fd', 'args', '-append');
-    else
-        save('matlab/mat/nass_disturbances.mat', 'Dw', 'Ty', 'Rz', 'Fd', 'args');
-    end
+    save('matlab/mat/nass_disturbances.mat', 'Dw', 'Dy', 'Rz', 'Fd', 'args');
 end