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

*z__$!LR  
clear; ]JlM/  
clc; ddEV@2F  
e0=1e-9/(36*pi);%epsilon0 Rw|P$dbu  
m0=1e-7*(4*pi);%mu0 ~T9wx  
c=3e8;%light speed Xj$'i/=-+c  
f=10e9; U,V+qnS  
lambda=c/f;%wavelength S'=}eeG  
omega=2*pi*f;%omega i}v3MO\X  
nd=20;%space grid ratio to lambda w\ddC DZ  
d=lambda/nd;%grid size !JbWxGN`jn  
s=0.5;%courantsteadycoefficient #,;Q|)AD:e  
dt=s*d/c;%time step LUEZqIf  
nx=100;%grid x number %18%T{|$e  
ny=100;%grid y number /|8/C40aY  
na=10;%absorber layers thickness G=&nwSL  
ez=zeros(nx+1,ny+1); :_pn|  
dz=zeros(nx-1,ny-1); }r|$\ms  
dz0=zeros(nx-1,ny-1); zE?@_p1gei  
hx=zeros(nx-1,ny); } ^WmCX2a  
hy=zeros(nx,ny-1); QW2SFpE  
bx=zeros(nx-1,ny); y>]Yq-  
by=zeros(nx,ny-1); g1&q6wCg|  
bx0=zeros(nx-1,ny); '6GW.;  
by0=zeros(nx,ny-1); \I7,
1I  
sigmax=zeros(nx,ny); 3zzl|+# 6  
sigmay=zeros(nx,ny); 0?=a$0_C  
m=2; "e
d A  
sigmamax=(m+1)/(150*pi*d);%lost tangent max =gHUY&sPu8  
sigmaa=sigmamax*((na-1:-1:0)/(na-1)).^m;%linearly distribution sigma la>H&  
sigmax(1:na,:)=bsxfun(@plus,sigmax(1:na,:),sigmaa'); $t.M`:G  
sigmax(end-na+1:end,:)=bsxfun(@plus,sigmax(end-na+1:end,:),fliplr(sigmaa)'); !rff/0/x"  
sigmay(:,1:na)=bsxfun(@plus,sigmay(:,1:na),sigmaa); :~'R|l  
sigmay(:,end-na+1:end)=bsxfun(@plus,sigmay(:,end-na+1:end),fliplr(sigmaa)); #pk  
kappa=1; EU.!/'<  
sigmaxEz=(sigmax(1:end-1,1:end-1)+sigmax(1:end-1,2:end)+sigmax(2:end,1:end-1)+sigmax(2:end,2:end))/4; "f>`ZFp^  
sigmayEz=(sigmay(1:end-1,1:end-1)+sigmay(1:end-1,2:end)+sigmay(2:end,1:end-1)+sigmay(2:end,2:end))/4; "X\6tl7a|  
sigmaxHx=(sigmax(1:end-1,:)+sigmax(2:end,:))/2; )F4BVPI  
sigmayHx=(sigmay(1:end-1,:)+sigmay(2:end,:))/2; ^AC2
 zC  
sigmaxHy=(sigmax(:,1:end-1)+sigmax(:,2:end))/2; y0,>_MS  
sigmayHy=(sigmay(:,1:end-1)+sigmay(:,2:end))/2; 1B~[L 5p9  
caDz=(kappa/dt-sigmaxEz/e0/2)./(kappa/dt+sigmaxEz/e0/2); GxA[N

 
cbDz=1/d./(kappa/dt+sigmaxEz/e0/2); MGH2z:  
caEz=(kappa/dt-sigmayEz/e0/2)./(kappa/dt+sigmayEz/e0/2); B:(a?X-7  
cbEz=1/dt./(kappa/dt+sigmayEz/e0/2)/e0; @j=rSS  
cpBx=(kappa/dt-sigmayHx/e0/2)./(kappa/dt+sigmayHx/e0/2); P9qZjBS  
cqBx=1/d./(kappa/dt+sigmayHx/e0/2); f}7/UGd  
cpHx=(kappa/dt+sigmaxHx/e0/2)./(m0/dt); R+tQvxp#  
cqHx=(kappa/dt-sigmaxHx/e0/2)./(m0/dt); TnJNs  
cpBy=(kappa/dt-sigmaxHy/e0/2)./(kappa/dt+sigmaxHy/e0/2); *V{Y.`\  
cqBy=1/d./(kappa/dt+sigmaxHy/e0/2); gq050Bl)  
cpHy=(kappa/dt+sigmayHy/e0/2)./(m0/dt); URj2 evYW  
cqHy=(kappa/dt-sigmayHy/e0/2)./(m0/dt); u|75r%p>  
h=imagesc(ez,0.1*[-1,1]); ,&s%^I+CC  
axis equal tight; [m(n-MuF  
for n=0:400 #dkSAS  
    ez(nx/2+1,ny/2+1)=ez(nx/2+1,ny/2+1)+sin(omega*n*dt); l]S%k&  
    dz=caDz.*dz+cbDz.*(diff(hy,1,1)-diff(hx,1,2)); 0D&-BAzi  
    ez(2:end-1,2:end-1)=caEz.*ez(2:end-1,2:end-1)+cbEz.*(dz-dz0); hSG1f`  
    dz0=dz; 5?8jj  
    bx=cpBx.*bx-cqBx.*diff(ez(2:end-1,:),1,2); X;CRy,  
    hx=hx+cpHx.*bx-cqHx.*bx0; 1'b}Y8YO  
    bx0=bx; U|+c&TY
 
NB3ar&.$S  

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