close all
clear all

fs = 10;  %Output timeseries frequency (Hz)
T = 300;  %Final time
t = 0:1/fs:T;

%Compute Bretshnieder PSD for positive frequencies.
Hs = 1.9;
Tp = 8.8;
fp = 1/Tp;
df = .5/500;
f = 0:df:.5;

S(1) = 0;
for i = 2:length(f)
  S(i) = (5/16)*fp^4*Hs^2*exp(-1.25*(fp/f(i))^4)/(f(i)^5);
end
std_dev = sqrt(trapz(f,S));


%Compute Timeseries
for n = 1:length(f)
  A(n) = sqrt(2*df*S(n));
  phi(n) = rand*2*pi;
end

for i = 1:length(t)
    z(i) = 0;
  for n = 1:length(f)
      z(i) = z(i) + A(n)*cos(2*pi*f(n)*t(i)+phi(n));
  end
end

fh1 = figure
plot(t,z)
xlabel('time [s]');
ylabel('z [m]');
grid on;
hold all
legend(['4*std(z) = ' num2str(4*std(z))]);


fh2 = figure
plot(f,S);
xlabel('f [Hz]');
ylabel('S [m^2/Hz]');
title(['Bretschieder: Hs = ' num2str(Hs) 'm  Tp = ' num2str(Tp) 's']);
legend(['4*sqrt(Area under S) = ' num2str(4*std_dev)]);
grid on;

% Make spectrum double-sided (Sxx shoud be your input)
  Sxx_positiveFreq = S/2;
  Sxx_negativeFreq = flip(conj(Sxx_positiveFreq(2:end-1)));
  Sxx = [Sxx_positiveFreq Sxx_negativeFreq]'; % [m^2/Hz]


plot_on = 1;
[timeseries, time] = TimeseriesFromPSD(Sxx, fs, T, plot_on);


figure(fh1)

plot(time, timeseries);
