手搓超声聚焦[MATLAB原创]

clc;clear;close
a=0.00005;   %波速
dx=0.002; dy=0.002;   %空间步长c
dt=0.01;  %时间步长
Lx=0.7; Ly=1; %求解域
Ti=150;     %传播时间
U0=0.5;       %声强
x=0:dx:Lx;y=0:dy:Ly;t=0:dt:Ti;
[X,Y]=meshgrid(x,y);
nx=length(x); ny=length(y);
U=zeros(length(x),length(y),2); %存放结果
i=1;j=1;k=2;   %初始化变量
U(:,1,1)=U0;U(:,1,2)=U0;
beta=0.005;K=1/(1+beta*dt);%衰减系数设置

%找到左边的弧面==============================
Ucx=Lx/2; Ucy=Ly/2;
M=(((X-Ucx).^2+(Y-Ucy).^2)<(Ly/2-dx)^2);
M=M+M([end,1:end-1],:)
[Xbw,Ybw]=find((Y<Ly/2).*(M==1));

 %主循环==============================
for it=dt:dt:Ti
    i=1;k=k+1;      %下标变量
    for ix=dx:dx:Lx-dx
        i=i+1; j=1;      %下标变量
        for mk=1:length(Ybw); U(mk,1:Xbw(mk),k)=0;end %左边界
        U(i,ny,k)=U(i,ny-1,k); %右边界
        for iy=dy:dy:Ly-dy
            j=j+1;
            U(1,j,k)=0*U(3,j,k); U(nx,j,k)=0*U(nx-2,j,k); %上下边界
            U(i,j,k)=K*(a*(U(i,j-1,k-1)-2*U(i,j,k-1)+U(i,j+1,k-1))/dx^2*dt^2+…
                     K*(a*(U(i-1,j,k-1)-2*U(i,j,k-1)+U(i+1,j,k-1))/dy^2*dt^2)+…
                     +(2+beta*dt)*U(i,j,k-1)-U(i,j,k-2));
        end 
        
        for mk=1:length(Ybw) %波动源
         U(Ybw(mk),Xbw(mk),k)=sin(it);
        end
 end
% 后处理==============================
if mod(k,50)==0;Ut=U(:,:,k);imagesc(Ut,[-5,5]);
colormap jet; axis equal; axis off;%后处理
drawnow;k
end
end