1D sea surface generation using Pierson-Moskowitz spectrum

Hi everyone,

I am working on a 1D sea surface generation. I used Pierson-Moskowitz spectrum to generate sea surface. I get stuck in somewhere in the code. Sea surface heights as output should be in meter (at least tens of centimeter corresponding wind speed) range, but they are in mm ranges. I think my fault is in spectrum vector. Code is developed from attached file (between page 22-26). Code is also shown in question. Thank you for your helps

clear all; close all;
%Generating sea surface with Pierson-Moskowitz variance spectrum
% S(k) is PM spectrum
% alpha is constant
% beta is constant
% V is wind speed (m/s) at 19.5 m above sea
% k is spatial frequency (rad/m)

%definition
L=100;         %surface length
N=128;         %fft points
delx= L/N;     %shortest wavelength that can be resolved
kf= 2*pi/L;    %fundamental spatial freq.
kny= (N/2)*kf; %the highest freq. sampled

k= linspace(-kny/2,kny/2,N); 

alpha = 0.0081;   %constant
beta = 0.74;      %constant
V = 10;       % wind speed m/s
g = 9.81;     % m/s^2

S = [];   %preallocation for PM spectrum
z0= [];   %preallocation for wave profile
z = []; 

%loop for PM spectrum and wave profile
for i=1:length(k)
    S(i)= (alpha/(2*k(i)^3))*exp(-beta*(g/k(i))*(g/k(i))*1/(V^4));  %PM spectrum
    rho = randn;
    sig = randn;
    z0(i) = (1/sqrt(2))*(rho + 1i*sig).*sqrt(0.5*S(i)*kf);    %wave profile
    z(i) = (1/sqrt(2))*(z0(i) + conj(z0(i)));   
end

Z=ifft(z);     %inverse dft for extracting surface elevations

%plot random surface wave 
figure,
plot(linspace(0,L,N),real(Z))
axis([0,L,-0.015,0.015]);
xlabel('distance (m)');
ylabel('surface elevation (m)');
out= sprintf('1D random sea surface V= %g m/sec at 19.5m and L= %g m',V,L);  
title(out);

%plot PM spectrum
figure,
plot(k,S);
%axis([0,2,0,1.5*max(S)]);
grid;
xlabel('spatial frequency (rad/m)');
ylabel('PM spectral density (m^2/rad/m)');
out1=sprintf('Pierson-Moskowitz Spectrum  V = %g m/sec at 19.5 m',V);
title(out1);

trapz(S)

[Jian Sun]

please read the equation 2.6 before ifft by MATLAB.

[Yi Ling Tan]

There is some scaling required after using fft/ifft. See https://www.mathworks.com/matlabcentral/answers/15770-scaling-the-fft-and-the-ifft

I have attached my matlab code.