24小时热门版块排行榜

返回列表

wjdwin

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 41.5
帖子: 7
在线: 9.6小时
虫号: 1138733
注册: 2010-11-04
专业: 建筑物理

[求助] 二维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

程序可以运行。但是自己画出的等位线图不对，contour（ezx+ezy）,contour(hx),contour(hy),都不对，求解决。

回复此楼