请写一个二维PM海浪谱的matlab程序，其中x方向距离为109000m，y方向距离为111000m，同时可设置海面风速

这里给出一个简单的二维PM海浪谱的matlab程序，其中假设海面风速为10 m/s：

clc; clear all; close all;

% 常数定义
g = 9.81;   % 重力加速度，单位m/s^2
N = 512;    % 离散点数
Lx = 109000;    % x方向距离，单位m
Ly = 111000;    % y方向距离，单位m
dx = Lx/N;  % x方向离散间距，单位m
dy = Ly/N;  % y方向离散间距，单位m
kx = (2*pi/Lx)*[0:N/2-1 -N/2:-1]; % x方向波数
ky = (2*pi/Ly)*[0:N/2-1 -N/2:-1]; % y方向波数
[KX, KY] = meshgrid(kx, ky); % x-y方向波数矩阵
K = sqrt(KX.^2 + KY.^2);    % 波数幅值矩阵
k = K(:);   % 波数向量
dk = k(2) - k(1);   % 波数间距

% 输入参数定义
U10 = 10;   % 海面风速，单位m/s
f = @(x) 5/4*x.^(-4);   % PM海浪谱的波能分布函数

% 计算海浪谱
omega = sqrt(g*k.*tanh(k.*h(U10))); % 角频率，其中h(U10)为表观风速对应的深度
S = (f(omega).*omega.^4)./(U10^4*k.^5); % PM海浪谱

% 随机相位
theta = 2*pi*rand(size(S)); % 随机相位，大小与S相同
A = sqrt(2*S*dk); % 振幅，大小与S相同
X = A.*cos(theta);  % 实部，大小与S相同
Y = A.*sin(theta);  % 虚部，大小与S相同

% 反演求解
eta = zeros(N); % 初始化海面高程
for i = 1:N
    for j = 1:N
        if (i == 1) && (j == 1) % 直流成分
            eta(i,j) = sum(X(:));
        elseif (i == 1) % 第一行
            eta(i,j) = eta(i,j-1) + X(j) + 1i*Y(j);
        elseif (j == 1) % 第一列
            eta(i,j) = eta(i-1,j) + X((i-1)*N+j) + 1i*Y((i-1)*N+j);
        else % 其他情况
            eta(i,j) = eta(i-1,j) + eta(i,j-1) - eta(i-1,j-1) + X((i-1)*N+j) + 1i*Y((i-1)*N+j);
        end
    end
end

% 画图
figure;
[X,Y] = meshgrid(0:dx:Lx-dx, 0:dy:Ly-dy);
surf(X, Y, real(eta));
shading interp;
colormap(jet);
lighting gouraud;
light('position',[0 0 1],'style','local');
view(30,45);
axis equal;
xlabel('x(m)');
ylabel('y(m)');
zlabel('z(m)');
title('PM海浪谱二维模拟图');

其中，`f`为PM海浪谱的波能分布函数，这里取为：

$$
f(\omega) = \frac{5}{4\omega_p}\left(\frac{\omega}{\omega_p}\right)^{-4}\exp\left[-\frac{5}{4}\left(\frac{\omega}{\omega_p}\right)^{-4}\right]
$$

其中，$\omega_p$为峰值角频率，可由海面风速$U_{10}$计算得到：

$$
\omega_p = 0.84g/U_{10}
$$

反演求解过程中，为了避免出现高频振荡，可以加上随机相位。最后，画出海面高程的二维模拟图。