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

-EU~ %/=m+  
clear; ?eR^\-e  
clc; 'p'nAB''!  
e0=1e-9/(36*pi);%epsilon0 3],[6%w  
m0=1e-7*(4*pi);%mu0 C5#$NV99p  
c=3e8;%light speed 55x.Q  
lambda=1;%dx, dy Y{m1\s/o  
f=c/lambda;%frequency NVo=5  
omega=2*pi*f; gO!:WD  
nd=20;%spatial grid ratio \OF"hPq  
d=lambda/nd;%mesh size E`3[62C  
s=0.5; w
Q9fPOm  
dt=s*d/c;
7:S)J~s*O  
mur1=(s-1)/(s+1); %coefficient for Mur absorbing 5CM]-qbf@  
murc=(s-sqrt(2))/(s+sqrt(2));%coefficient for corner l8!n!sC[
,  
mur2=0.5*c*m0*s/(s+1);%coeffient for Mur2 absorbing x$V[xX  
cb=dt/(d*e0);%coefficient for E from H ,tOc+3Qz$  
cq=dt/(d*m0);%coefficient for H from E MqswYK-s  
nx=100;%网格数 p\,PY  
ny=100; sX
=_|<[  
ez=zeros(nx+1,ny+1);%field cz*Z/5XH  
hx=zeros(nx+1,ny); ,W;2A0A?X  
hy=zeros(nx,ny+1); [Q20c<,  
fh=@(t)(heaviside(t).*sin(t).*((heaviside(t)-heaviside(t-6*pi)).*(1-cos(pi*t/(6*pi)))/2+heaviside(t-6*pi))); ljj}XJQ  
fh=@(t)(fh(omega.*t)); ("@ih]zYf  
h=imagesc(ez,0.05*[-1,1]);%handle for image j}DG +M  
axis equal tight; N6S}u@{J~N  
nt=1000; q%}54E80  
for n=1:nt%time iteration J.npv1F  
    ez0x=ez([2,end-1],:); -B#>Jn#F  
    ez0y=ez(:,[2,end-1]); ]4oF!S%F  
    ez0c=ez([2,end-1],[2,end-1]); rIF6^?  
    ez(2:end-1,2:end-1)=ez(2:end-1,2:end-1)+cb*(diff(hy(:,2:end-1),1,1)-diff(hx(2:end-1,:),1,2));%更新E XiUae{j`  
    ez([1,end],2:end-1)=ez0x(:,2:end-1)+mur1*(ez([2,end-1],2:end-1)-ez([1,end],2:end-1))-mur2*(diff(hx([1,end],:),1,2)+diff(hx([2,end-1],:),1,2));%边界条件，加内 I!,FxOM|$  
    ez(2:end-1,[1,end])=ez0y(2:end-1,:)+mur1*(ez(2:end-1,[2,end-1])-ez(2:end-1,[1,end]))+mur2*(diff(hy(:,[1,end]),1,1)+diff(hy(:,[2,end-1]),1,1));%后面是mur2条件 v!?>90a  
    ez([1,end],[1,end])=ez0c+murc*(ez([2,end-1],[2,end-1])-ez([1,end],[1,end]));%边界条件，角点 .AHww7  
    ez(nx/2+1,ny/2+1)=ez(nx/2+1,ny/2+1)+fh(n*dt); p< jM%fbZk  
    hx=hx-cq*diff .. nSV OS6  
Tk0Senq,  

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