1314阅读
0回复

[RFEDA原创]FDTD迷你代码-2D BPML

UID ：90389

注册：2012-03-16
登录：2014-09-11
发帖：17
等级：仿真新人

0楼发表于: 2014-06-11 08:58:30

clear;
clc;
e0=1e-9/(36*pi);%epsilon0
m0=1e-7*(4*pi);%mu0
c=3e8;%light speed
f=1e9;
lambda=c/f;%wavelength
omega=2*pi*f;%omega
nd=20;%space grid ratio to lambda
d=lambda/nd;%grid size
s=0.5;%courantsteadycoefficient
dt=s*d/c;%time step
nx=100;%grid x number
ny=100;%grid y number
na=10;%absorber layers thickness
ezx=zeros(nx+1,ny+1);%Hz
ezy=zeros(nx+1,ny+1);%Hz
hx=zeros(nx-1,ny);%Ex,derive to y
hy=zeros(nx,ny-1);%Hy,derive to x
m=2;
sigmaMax=(m+1)/(150*pi*d);%lost tangent max
sigma=sigmaMax*((na-1:-1:0)/(na-1)).^m;%linearly distribution sigma
sigmam=sigma*m0/e0;% sigma m
sigma=(sigma(1:end-1)+sigma(2:end))/2;
sigmax=zeros(nx-1,ny-1);
sigmay=zeros(nx-1,ny-1);
sigmamx=zeros(nx,ny-1);
sigmamy=zeros(nx-1,ny);
sigmax(1:na-1,:)=repmat(sigma',1,ny-1);
sigmax(end-na+2:end,:)=repmat(fliplr(sigma)',1,ny-1);
sigmay(:,1:na-1)=repmat(sigma,nx-1,1);
sigmay(:,end-na+2:end)=repmat(fliplr(sigma),nx-1,1);
sigmamx(1:na,:)=repmat(sigmam',1,ny-1);
sigmamx(end-na+1:end,:)=repmat(fliplr(sigmam)',1,ny-1);
sigmamy(:,1:na)=repmat(sigmam,nx-1,1);
sigmamy(:,end-na+1:end)=repmat(fliplr(sigmam),nx-1,1);
caEzx=(e0/dt-sigmax/2)./(e0/dt+sigmax/2);
cbEzx=1/d./(e0/dt+sigmax/2);
caEzy=(e0/dt-sigmay/2)./(e0/dt+sigmay/2);
cbEzy=1/d./(e0/dt+sigmay/2);
cpHx=(m0/dt-sigmamy/2)./(m0/dt+sigmamy/2);
cqHx=1/d./(m0/dt+sigmamy/2);
cpHy=(m0/dt-sigmamx/2)./(m0/dt+sigmamx/2);
cqHy=1/d./(m0/dt+sigmamx/2);
h=imagesc(ezx+ezy,0.1*[-1,1]);
axis equal tight;
for n=0:300
    ezx(nx/2+1,ny/2+1)=ezx(nx/2+1,ny/2+1)+sin(omega*n*dt);
    ezy(nx/2+1,ny/2+1)=ezy(nx/2+1,ny/2+1)+sin(omega*n*dt);
    ezx(2:end-1,2:end-1)=caEzx.*ezx(2:end-1,2:end-1)+cbEzx.*diff(hy,1,1);
    ezy(2:end-1,2:end-1)=caEzy.*ezy(2:end-1,2:end-1)-cbEzy.*diff(hx,1,2);
    hx=cpHx.*hx-cqHx.*diff(ezx(2:end-1,:)+ezy(2:en ..

未注册仅能浏览部分内容，查看全部内容及附件请先登录或注册

描述:请输入描述

图片:untitled.jpg

发帖回复