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

+zFEx%3^  
clear; l~`JFWur]  
clc; \ ]h$8JwV  
e0=1e-9/(36*pi);%epsilon0 U%2{PbL  
m0=1e-7*(4*pi);%mu0 kdm@1x  
c=3e8;%light speed ,+g0#8?p^x  
f=10e9; n{F&GE="  
lambda=c/f;%wavelength 4,6?sTuX  
omega=2*pi*f;%omega 3V/|"R2s  
nd=20;%space grid ratio to lambda aOW~! f/M  
d=lambda/nd;%grid size w4&-9[@Y  
s=0.5;%courantsteadycoefficient KU0;}GSNX}  
dt=s*d/c;%time step jAFJ?L(  
nx=100;%grid x number o, qBMo^.  
ny=100;%grid y number QRY7ck:N  
na=10;%absorber layers thickness SQ`ec95',  
ez=zeros(nx+1,ny+1); 5sMyH[5zY  
dz=zeros(nx-1,ny-1); .V^h<d{  
dz0=zeros(nx-1,ny-1); '^t(=02J  
hx=zeros(nx-1,ny); H!g9~a  
hy=zeros(nx,ny-1); >3ASrM+>w  
bx=zeros(nx-1,ny); )%?SWuS?N  
by=zeros(nx,ny-1); Q mz3GH@wg  
bx0=zeros(nx-1,ny); hniT
MO  
by0=zeros(nx,ny-1);
6FI`0j=~  
sigmax=zeros(nx,ny); Su`] ku'  
sigmay=zeros(nx,ny); w>#.id[k  
m=2; FwSV \N+#'  
sigmamax=(m+1)/(150*pi*d);%lost tangent max y{qKb:~wv  
sigmaa=sigmamax*((na-1:-1:0)/(na-1)).^m;%linearly distribution sigma
cT^x^%  
sigmax(1:na,:)=bsxfun(@plus,sigmax(1:na,:),sigmaa'); 1b"3]?  
sigmax(end-na+1:end,:)=bsxfun(@plus,sigmax(end-na+1:end,:),fliplr(sigmaa)'); +3;[1dpgf  
sigmay(:,1:na)=bsxfun(@plus,sigmay(:,1:na),sigmaa); h6gtO$A|p=  
sigmay(:,end-na+1:end)=bsxfun(@plus,sigmay(:,end-na+1:end),fliplr(sigmaa)); D|5Fo'O^AV  
kappa=1; $-]PD`wmY  
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; g>Kh?(  
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; 771r(X?Fa  
sigmaxHx=(sigmax(1:end-1,:)+sigmax(2:end,:))/2; [%7oq;^J  
sigmayHx=(sigmay(1:end-1,:)+sigmay(2:end,:))/2; w3oe.hWP3N  
sigmaxHy=(sigmax(:,1:end-1)+sigmax(:,2:end))/2; x@"`KiEUs  
sigmayHy=(sigmay(:,1:end-1)+sigmay(:,2:end))/2; !%yd'"6Dl  
caDz=(kappa/dt-sigmaxEz/e0/2)./(kappa/dt+sigmaxEz/e0/2); U[l{cRT   
cbDz=1/d./(kappa/dt+sigmaxEz/e0/2); *&yt;|y  
caEz=(kappa/dt-sigmayEz/e0/2)./(kappa/dt+sigmayEz/e0/2); 1A
9
Gf  
cbEz=1/dt./(kappa/dt+sigmayEz/e0/2)/e0; E@ !~q  
cpBx=(kappa/dt-sigmayHx/e0/2)./(kappa/dt+sigmayHx/e0/2); o7 X5{  
cqBx=1/d./(kappa/dt+sigmayHx/e0/2); "C%* 'k  
cpHx=(kappa/dt+sigmaxHx/e0/2)./(m0/dt); ;hU~nj+{  
cqHx=(kappa/dt-sigmaxHx/e0/2)./(m0/dt); cX-)]D  
cpBy=(kappa/dt-sigmaxHy/e0/2)./(kappa/dt+sigmaxHy/e0/2); :j!N7c{  
cqBy=1/d./(kappa/dt+sigmaxHy/e0/2);
AQz&u  
cpHy=(kappa/dt+sigmayHy/e0/2)./(m0/dt); 
ehYGw2  
cqHy=(kappa/dt-sigmayHy/e0/2)./(m0/dt); vn=0
=(  
h=imagesc(ez,0.1*[-1,1]); so\8.(7n  
axis equal tight; RjQdlr6*  
for n=0:400 gNo}\ lm4V  
    ez(nx/2+1,ny/2+1)=ez(nx/2+1,ny/2+1)+sin(omega*n*dt); qgLj^{  
    dz=caDz.*dz+cbDz.*(diff(hy,1,1)-diff(hx,1,2)); {nmBIk2v  
    ez(2:end-1,2:end-1)=caEz.*ez(2:end-1,2:end-1)+cbEz.*(dz-dz0); QW"BGg~6c  
    dz0=dz; 0Z~G:$O/i  
    bx=cpBx.*bx-cqBx.*diff(ez(2:end-1,:),1,2); e;)&Hc:Z  
    hx=hx+cpHx.*bx-cqHx.*bx0; V\><6v  
    bx0=bx; ?t];GNU`l  
5-X(K 'Q  

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