瑞利波模拟程序，用于面波入射模拟

%% ---- Function Rayleigh_Waves --- %% %% ---- 本程序是显示瑞雷面波的程序 --- %% clc clear close filename= 'forward.gif'; % 生成波场快照的文件名 %%%%%% set_delay = 1/24; % 每一幅图和每一幅图的时间延迟，注意人眼识别每分钟24幅图画 f = 0.5; T = 1/f; omega = 2*pi*f; vel = 3000; ncol = 7; % Set up array of depths x = [0:300:15000]; z = [0:-100:-2000]; Ax = 200; Az = 200; nx = length(x); nz = length(z); k = 2*pi*f/vel; lamda = 750; %Start time step loop n=40; M=moviein(n); for it=1:n t=T*(it)/n; for ix=1:nx dx(ix) = 0.001*(x(ix)+Ax*exp(z(1)/lamda)*sin(k*x(ix)-omega*t)); dz(ix) = 0.001*(z(1)+Az*exp(z(1)/lamda)*cos(k*x(ix)-omega*t)); end figure plot(dx,dz,'o-') hold on plot(dx(20),dz(20),'or') axis([-0.2,15.2,-2.05,1]) title('Rayleigh wave f=0.5 Hz') ylabel ('depth(km)') xlabel ('(km)') for iz=2:ncol-1 for ix=1:nx dx(ix) = 0.001*(x(ix)+Ax*exp(z(iz)/lamda)*sin(k*x(ix)-omega*t)); dz(ix) = 0.001*(z(iz)+Az*exp(z(iz)/lamda)*cos(k*x(ix)-omega*t)); end plot(dx,dz,'o') end for ix=1:nx dx(ix) = 0.001*(x(ix)+Ax*exp(z(ncol)/lamda)*sin(k*x(ix)-omega*t)); dz(ix) = 0.001*(z(ncol)+Az*exp(z(ncol)/lamda)*cos(k*x(ix)-omega*t)); end plot(dx,dz,'ro') for iz=ncol+1:nz for ix=1:nx dx(ix) = 0.001*(x(ix)+Ax*exp(z(iz)/lamda)*sin(k*x(ix)-omega*t)); dz(ix) = 0.001*(z(iz)+Az*exp(z(iz)/lamda)*cos(k*x(ix)-omega*t)); end plot(dx,dz,'o') end set(gcf,'color',[1 1 1]) saveas(gcf,strcat('第',num2str(it)),'bmp') frame=getframe; % 把每一幅图当做一帧 % im=frame2im(frame); % 读入 % [I,map]=rgb2ind(im,256); % 颜色变为数值 % if it==1; imwrite(I,map,filename,'gif','Loopcount',inf,... 'DelayTime',set_delay); % loopcount第一个波场快照格式? else imwrite(I,map,filename,'gif','WriteMode','append',... 'DelayTime',set_delay); % l后面波场快照格式 end hold off M(:,it)=getframe; end movie(M,10) close all