matlab 2d pml,二维TM波，PML边界，matlab程序。但是仿真结果有问题，求高手。

该楼层疑似违规已被系统折叠 隐藏此楼查看此楼

clc;

clear;

ft=1.5;

fs=20;

NX=100;%X方向网格数

NY=150;%Y方向网格数

NA=20;%吸收层厚度

bxrhx=1:100;

bxirhx=2:150;

byrhy=2:100;

byirhy=1:150;

bxrez=1:100;

byrez=1:150;

hx=zeros(100,151);

hy=zeros(101,150);

ezx=zeros(100,150);

ezy=zeros(100,150);

sigmamy=hx;

sigmamx=hy;

sigmax=ezx;

sigmay=ezy;

sigmamax=0.2;

sigma=sigmamax*(NA:-1:1)/NA;

sigmam=sigmamax*(2*NA-1:-2:1)/2/NA;

sigmamx([1:NA,end:-1:end-NA+1],:)=repmat([sigmam,sigmam]',1,150);

sigmax([1:NA,end:-1:end-NA+1],:)=repmat([sigma,sigma]',1,150);

sigmamy(:,[1:NA,end:-1:end-NA+1])=repmat([sigmam,sigmam],100,1);

sigmay(:,[1:NA,end:-1:end-NA+1])=repmat([sigma,sigma],100,1);

f1hx=exp(-sigmamy(bxrhx,bxirhx));

f1hy=exp(-sigmamx(byrhy,byirhy));

f2hx=zeros(size(hx));

f2hy=zeros(size(hy));

f2hx(sigmamy == 0)=1/ft;

f2hy(sigmamx==0)=1/ft;

f2hx(sigmamy~=0)=(1-exp(-sigmamy(sigmamy~=0)))./sigmamy(sigmamy~=0)/ft;

f2hy(sigmamx~=0)=(1-exp(-sigmamx(sigmamx~=0)))./sigmamx(sigmamx~=0)/ft;

f2hx=f2hx(bxrhx,bxirhx);

f2hy=f2hy(byrhy,byirhy);

f1ezx=exp(-sigmax);

f1ezy=exp(-sigmay);

f2ezx=zeros(size(ezx));

f2ezy=zeros(size(ezy));

f2ezx(sigmax==0)=1/ft;

f2ezy(sigmay==0)=1/ft;

f2ezx(sigmax~=0)=(1-exp(-sigmax(sigmax~=0)))./sigmax(sigmax~=0)/ft;

f2ezy(sigmay~=0)=(1-exp(-sigmay(sigmay~=0)))./sigmay(sigmay~=0)/ft;

centerx=floor((1+NX)/2);

centery=floor((1+NY)/2);

for nn=1:150

ezx(centerx,centery)=ezx(centerx,centery)+sin(2*pi*nn/ft/fs);

ezy(centerx,centery)=ezy(centerx,centery)+sin(2*pi*nn/ft/fs);

hx(bxrhx,bxirhx)=f1hx.*hx(bxrhx,bxirhx)-f2hx.*(ezx(bxrhx,bxirhx)-ezx(bxrhx,bxirhx-1)+ezy(bxrhx,bxirhx)-ezy(bxrhx,bxirhx-1));

hy(byrhy,byirhy)=f1hy.*hx(byrhy,byirhy)+f2hy.*(ezx(byrhy,byirhy)-ezx(byrhy-1,byirhy)+ezy(byrhy,byirhy)-ezy(byrhy-1,byirhy));

ezx=f1ezx.*ezx+f2ezx.*(hy(bxrez+1,:)-hy(bxrez,:));

ezy=f1ezy.*ezy-f2ezy.*(hx(:,byrez+1)-hx(:,byrez));

end