function [] = initializeDisturbances(args)
% initializeDisturbances - Initialize the disturbances
%
% Syntax: [] = initializeDisturbances(args)
%
% Inputs:
%    - args -

arguments
    % Global parameter to enable or disable the disturbances
    args.enable logical {mustBeNumericOrLogical} = true
    % Ground Motion - X direction
    args.Dwx logical {mustBeNumericOrLogical} = true
    % Ground Motion - Y direction
    args.Dwy logical {mustBeNumericOrLogical} = true
    % Ground Motion - Z direction
    args.Dwz logical {mustBeNumericOrLogical} = true
    % Translation Stage - X direction
    args.Fty_x logical {mustBeNumericOrLogical} = true
    % Translation Stage - Z direction
    args.Fty_z logical {mustBeNumericOrLogical} = true
    % Spindle - X direction
    args.Frz_x logical {mustBeNumericOrLogical} = true
    % Spindle - Y direction
    args.Frz_y logical {mustBeNumericOrLogical} = true
    % Spindle - Z direction
    args.Frz_z logical {mustBeNumericOrLogical} = true
end

% Initialization of random numbers
rng("shuffle");

%% Ground Motion
load('dist_psd.mat', 'dist_f');

% 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);

% 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;
    C = zeros(N/2,1);
    for i = 1:N/2
        C(i) = sqrt(phi(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;
else
    Ty.x = zeros(length(Ty.t), 1);
end

% Translation Stage - Z
if args.Fty_z && args.enable
    phi = Ty.psd_z;
    C = zeros(N/2,1);
    for i = 1:N/2
        C(i) = sqrt(phi(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;
else
    Ty.z = zeros(length(Ty.t), 1);
end

%% Translation Stage
load('dist_psd.mat', 'dist_f');

% 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);

% 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;
    C = zeros(N/2,1);
    for i = 1:N/2
        C(i) = sqrt(phi(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);
    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);
    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;
else
    Rz.z = zeros(length(Rz.t), 1);
end

u = zeros(100, 6);
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);
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
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
end