线性波传播至垂直壁面反射形成驻波的动画MATLAB程序

最近在写博士论文非常繁忙，很久没有过啦更新，这次带来个之前写的好玩的程序，可以帮助大家理解线性波传递以及驻波的形成。

%% linear waves add
%% code by JN_CUI
%%   Variables
%   eta     - water surface profile
%   u       - horizontal velocity of water particle
%   w       - vortical velocity of water particle
%   zeta    - horizontal displacement of water particle
%   epsilon - vertical displacement of water particle
%   t       - time
%   x       - horizontal coordinate
%   z       - vertical coordinate
%   k       - wavenumber
%   H       - wave height
%   T       - wave period
%   sigma   - phase angle
clear;clc;
%% Function
%   Vedio settings
Path        = 'C:\Users\JunnanCui\Desktop\Sandingwave\';
Folder      = 'gif\';
V           = VideoWriter([Path,Folder,'linear_wave_add_position3.avi']); % 输出AVI文件名字
V.FrameRate = 50; % 帧率，50fps
V.Quality   = 90;   % 质量，百分比，90%
GIF         = [Path,Folder,'linear_wave_add_position3.gif'];
open(V);
%   Wave settings
H     = 0.5;
g     = 9.8;
T     = 1.96134;
d     = 10;
K     = g*T^2/(2*pi);
fun   = @(L1) L1/tanh(2*pi/L1*d)-K;
L     = Division(fun,0.0001,0,K);
k     = (2*pi/L);
sigma = L/T;
c     = sigma/k;
% k=1;
x     = 0:0.01:12;
x1    = 0:0.5:12;
z1    = -1:-0.5:-10;
[X,Z] = meshgrid(x1,z1);
tt= 0:0.1:20;
i=1;

for t = 0:0.1:20
    
    set(0,'DefaultFigureVisible', 'off');
    eta1     = H/2*(x<=c*t).*cos(k*x-sigma*t+pi/2);
    eta2     = H/2*((4*L-x)<=c*t).*cos(k*(-x)-sigma*t+pi/2);
    Eta      = eta1+eta2;
    u1       = pi*H*(X<=c*t)/T.*cosh(k*(d+Z))/sinh(k*d).*cos(k*X-sigma*t+pi/2);
    w1       = pi*H*(X<=c*t)/T.*sinh(k*(d+Z))/sinh(k*d).*sin(k*X-sigma*t+pi/2);
    zeta1    = -H*(X<=c*t)/2.*cosh(k*(d+Z))/sinh(k*d).*sin(k*X-sigma*t+pi/2);
    epsilon1 = H*(X<=c*t)/2.*sinh(k*(d+Z))/sinh(k*d).*cos(k*X-sigma*t+pi/2);
    
    u2       = pi*0.5*H*((4*L-X)<=c*t)/T.*cosh(k*(d+Z))/sinh(k*d).*cos(-k*X-sigma*t+pi/2);
    w2       = pi*0.5*H*((4*L-X)<=c*t)/T.*sinh(k*(d+Z))/sinh(k*d).*sin(-k*X-sigma*t+pi/2);
    zeta2    = -0.5*H*((4*L-X)<=c*t)/2.*cosh(k*(d+Z))/sinh(k*d).*sin(-k*X-sigma*t+pi/2);
    epsilon2 = 0.5*H*((4*L-X)<=c*t)/2.*sinh(k*(d+Z))/sinh(k*d).*cos(-k*X-sigma*t+pi/2);
%     if i==1
    XX{i}   = X+zeta1-zeta2;
    ZZ{i}   = Z+epsilon1+epsilon2;
    U       = u1-u2;
    W       = w1+w2;
    plot(XX{i}(1:2:end,1:end),ZZ{i}(1:2:end,1:end),'r.','markersize',15);hold on
    for j = 2:2:20
        if i>=j
            plot(XX{i-j+1}(1:2:end,1:end),ZZ{i-j+1}(1:2:end,1:end),'r.');hold on
        end
    end
    plot(x,Eta,'b','linewidth',2);hold on
    plot(x1,x1*0,'r--');hold on
    fill([2*L,2*L,12.2,12.2],[-10,2,2,-10],[0.8,0.8,0.8]);hold on
    quiver(X,Z,U,W,'k','linewidth',2,'autoscale','off');
    title(['t = ',num2str(t,'%.1f'),' s']);
    xlim([0,13]);
    ylim([-10,2]);
    xlabel('X (m)');
    ylabel('Depth (m)');
    set(gcf,'position',[100,100,1500,700])
    set(gca,'fontsize',25);
    img = gcf;
    % 输出GIF文件
    CurrFrame = getframe(img);   % 获取像素，否则无法显示动画
    imshow(CurrFrame.cdata(1:5:end,1:5:end));
    writeVideo(V, CurrFrame.cdata); %%% 输出视频文件
    im = frame2im(CurrFrame);
    [A,map] = rgb2ind(im,256);  % 将RGB图像转换为索引图像
    if t == 0
        imwrite(A,map,GIF,'gif','LoopCount',Inf,'DelayTime',0.005);  % DelayTime表示写入的时间间隔
    else
        imwrite(A,map,GIF,'gif','WriteMode','append','DelayTime',0.005);
    end
    close
    i=i+1;
end
close(V); %%% 关闭视频输出

Division子函数跟之前随机波里面的那个函数一样就不贴了

最后输出的视频和GIF是这样婶儿的：