?? untitled.m
字號:
for i = 1:ihT
for j = 1:jhT
CA(i,j) = (1.0 - sig(i,j)*dt / (2.0*epsi(i,j))) / (1.0 + sig(i,j) * dt / (2.0*epsi(i,j)));
CB(i,j) = (dt/(epsi(i,j))) / (1.0 + sig(i,j)*dt / (2.0*epsi(i,j)));
end
end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% FILL IN DENOMINATORS FOR FIELD UPDATES
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% CPML, Psi_w,v
% (7.109)
ii = iBC-1;
for i = 1:ieT
if i <= iBC-1
F_hx(i) = 1.0/(K_HxBC(i)*dxyz);
elseif i >= ihT+1-iBC
F_hx(i) = 1.0/(K_HxBC(ii)*dxyz);
ii = ii-1;
else
F_hx(i) = 1.0/dxyz;
end
end
jj = jBC-1;
for j = 1:jeT;
if j <= jBC-1
F_hy(j) = 1.0/(K_HyBC(j)*dxyz);
elseif j >= jhT+1-jBC
F_hy(j) = 1.0/(K_HyBC(jj)*dxyz);
jj = jj-1;
else
F_hy(j) = 1.0/dxyz;
end
end
ii = iBC;
for i = 1:ieT
if i <= iBC
F_ex(i) = 1.0/(K_ExBC(i)*dxyz);
elseif i >= ihT+1-iBC
F_ex(i) = 1.0/(K_ExBC(ii)*dxyz);
ii = ii-1;
else
F_ex(i) = 1.0/dxyz;
end
end
jj = jBC;
for j = 1:jeT
if j <= jBC
F_ey(j) = 1.0/(K_EyBC(j)*dxyz);
elseif j >= jhT+1-jBC
F_ey(j) = 1.0/(K_EyBC(jj)*dxyz);
jj = jj-1;
else
F_ey(j) = 1.0/dxyz;
end
end
%.:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:. .:.
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% BEGIN TIME STEP
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
for nTimeStep = 1:MaxTime
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% UPDATE Ex
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
for i = 1:ieT
for j = 2:jeT
Ex(i,j) = CA(i,j) * Ex(i,j) + CB(i,j) * ( (Hz(i,j) - Hz(i,j-1))*F_ey(j) ); % The first four terms in (7.106)
end
end
for i = 1:ieT
%.....................................................................
% CPML for bottom Ex, j-direction
%.....................................................................
for j = 2:jBC
psi_Exy1(i,j) = bEy(j)*psi_Exy1(i,j) + cEy(j) *(Hz(i,j) - Hz(i,j-1))/dxyz; %(7.105)
Ex(i,j) = Ex(i,j) + CB(i,j)*psi_Exy1(i,j); %The last two terms in (7.106)
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% UPDATE Ey
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
for i = 2:ieT
for j = 1:jeT
Ey(i,j) = CA(i,j) * Ey(i,j) + CB(i,j) * ( (Hz(i-1,j) - Hz(i,j))*F_ex(i) );
end
end
for j = 1:jeT
%.....................................................................
% CPML for bottom Ey, i-direction
%.....................................................................
for i = 2:iBC
psi_Eyx1(i,j) = bEx(i)*psi_Eyx1(i,j) + cEx(i)*(Hz(i-1,j) - Hz(i,j))/dxyz;
Ey(i,j) = Ey(i,j) + CB(i,j)*psi_Eyx1(i,j);
end
%.....................................................................
% CPML for top Ey, i-direction
%.....................................................................
ii = iBC;
for i = ihT+1-iBC:ieT
psi_Eyx2(ii,j) = bEx(ii)*psi_Eyx2(ii,j) + cEx(ii)*(Hz(i-1,j) - Hz(i,j))/dxyz;
Ey(i,j) = Ey(i,j) + CB(i,j)*psi_Eyx2(ii,j);
ii = ii-1;
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% UPDATE Ez
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
for i = 2:ieT
for j = 2:jeT
Ez(i,j) = CA(i,j) * Ez(i,j) + CB(i,j) * ((Hy(i,j) - Hy(i-1,j))*F_ex(i) + (Hx(i,j-1) - Hx(i,j))*F_ey(j));
end
end
for j = 2:jeT
%.....................................................................
% CPML for bottom Ez, x-direction
%.....................................................................
for i = 2:iBC
psi_Ezx1(i,j) = bEx(i)*psi_Ezx1(i,j) + cEx(i) *(Hy(i,j) - Hy(i-1,j))/dxyz;
Ez(i,j) = Ez(i,j) + CB(i,j)*psi_Ezx1(i,j);
end
%.....................................................................
% CPML for top Ez, x-direction
%.....................................................................
ii = iBC;
for i = ihT+1-iBC:ieT
psi_Ezx2(ii,j) = bEx(ii)*psi_Ezx2(ii,j) + cEx(ii) *(Hy(i,j) - Hy(i-1,j))/dxyz;
Ez(i,j) = Ez(i,j) + CB(i,j)*psi_Ezx2(ii,j);
ii = ii-1;
end
end
for i = 2:ieT
%.....................................................................
% CPML for bottom Ez, y-direction
%.....................................................................
for j = 2:jBC
psi_Ezy1(i,j) = bEy(j)*psi_Ezy1(i,j) + cEy(j)*(Hx(i,j-1) - Hx(i,j))/dxyz;
Ez(i,j) = Ez(i,j) + CB(i,j)*psi_Ezy1(i,j);
end
%.....................................................................
% CPML for top Ez, y-direction
%.....................................................................
jj = jBC;
for j = jhT+1-jBC:jeT
psi_Ezy2(i,jj) = bEy(jj)*psi_Ezy2(i,jj) + cEy(jj)*(Hx(i,j-1) - Hx(i,j))/dxyz;
Ez(i,j) = Ez(i,j) + CB(i,j)*psi_Ezy2(i,jj);
jj = jj-1;
end
end
%-----------------------------------------------------------------------
% SOURCE
%-----------------------------------------------------------------------
% i = isource;
% j = jsource;
% k = ksource;
% source = -2.0*((nTimeStep*dt-tO)/tw) *
% exp(-((nTimeStep*dt-tO)/tw)^2.0);
%source =exp(-((nTimeStep*dt-t0)/tw)^2);
%if nTimeStep<=floor(t_periodic)
% source = sin(omega*dt*nTimeStep);
%else
% source = 0;
%end
source =sin(omega*nTimeStep*dt);
Ez(isrc,30) = Ez(isrc,30)+source;
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% UPDATE Hx
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
for i = 2:ieT
for j = 1:jeT
Hx(i,j) = DA * Hx(i,j) + DB * ( (Ez(i,j) - Ez(i,j+1))*F_hy(j) ) ;
end
end
for i = 2:ieT
%.....................................................................
% CPML for bottom Hx, j-direction
%.....................................................................
for j = 1:jBC-1
psi_Hxy1(i,j) = bHy(j)*psi_Hxy1(i,j) + cHy(j) *(Ez(i,j) - Ez(i,j+1))/dxyz;
Hx(i,j) = Hx(i,j) + DB*psi_Hxy1(i,j);
end
%.....................................................................
% CPML for top Hx, j-direction
%.....................................................................
jj = jBC-1;
for j = jhT+1-jBC:jeT
psi_Hxy2(i,jj) = bHy(jj)*psi_Hxy2(i,jj) + cHy(jj) *(Ez(i,j) - Ez(i,j+1))/dxyz ;
Hx(i,j) = Hx(i,j) + DB*psi_Hxy2(i,jj) ;
jj = jj-1 ;
end
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% UPDATE Hy
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
for i = 1:ieT
for j = 2:jeT
Hy(i,j) = DA * Hy(i,j) + DB * ( (Ez(i+1,j) - Ez(i,j))*F_hx(i) );
end
end
for j = 2:jeT
%.....................................................................
% CPML for bottom Hy, i-direction
%.....................................................................
for i = 1:iBC-1
psi_Hyx1(i,j) = bHx(i)*psi_Hyx1(i,j) + cHx(i)*(Ez(i+1,j) - Ez(i,j))/dxyz;
Hy(i,j) = Hy(i,j) + DB*psi_Hyx1(i,j);
end
%.....................................................................
% CPML for top Hy, i-direction
%.....................................................................
ii = iBC-1;
for i = ihT+1-iBC:ieT
psi_Hyx2(ii,j) = bHx(ii)*psi_Hyx2(ii,j) + cHx(ii)*(Ez(i+1,j) - Ez(i,j))/dxyz;
Hy(i,j) = Hy(i,j) + DB*psi_Hyx2(ii,j);
ii = ii-1;
end
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% UPDATE Hz
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
for i = 1:ieT
for j = 1:jeT
Hz(i,j) = DA * Hz(i,j) + DB * ((Ey(i,j) - Ey(i+1,j))*F_hx(i) + (Ex(i,j+1) - Ex(i,j))*F_hy(j));
end
end
for j = 1:jeT
%.....................................................................
% CPML for bottom Hz, x-direction
%.....................................................................
for i = 1:iBC-1
psi_Hzx1(i,j) = bHx(i)*psi_Hzx1(i,j) + cHx(i) *(Ey(i,j) - Ey(i+1,j))/dxyz;
Hz(i,j) = Hz(i,j) + DB*psi_Hzx1(i,j);
end
%.....................................................................
% CPML for top Hz, x-direction
%.....................................................................
ii = iBC-1;
for i = ihT+1-iBC:ieT
psi_Hzx2(ii,j) = bHx(ii)*psi_Hzx2(ii,j) + cHx(ii) *(Ey(i,j) - Ey(i+1,j))/dxyz;
Hz(i,j) = Hz(i,j) + DB*psi_Hzx2(ii,j);
ii = ii-1;
end
end
for i = 1:ieT
%.....................................................................
% CPML for bottom Hz, y-direction
%.....................................................................
for j = 1:jBC-1
psi_Hzy1(i,j) = bHy(j)*psi_Hzy1(i,j) + cHy(j)*(Ex(i,j+1) - Ex(i,j))/dxyz;
Hz(i,j) = Hz(i,j) + DB*psi_Hzy1(i,j);
end
%.....................................................................
% CPML for top Hz, y-direction
%.....................................................................
jj = jBC-1;
for j = jhT+1-jBC:jeT
psi_Hzy2(i,jj) = bHy(jj)*psi_Hzy2(i,jj) + cHy(jj)*(Ex(i,j+1) - Ex(i,j))/dxyz;
Hz(i,j) = Hz(i,j) + DB*psi_Hzy2(i,jj);
jj = jj-1;
end
end
%-----------------------------------------------
%-----------------------------------------------
%-----------------------------------------------
n=nTimeStep;
EzPMLpower(n)=0;
%{
for ja=jBC+1:jeT-jBC
for ia=iBC+1:2*iBC
EzPMLpower(n)= EzPMLpower(n) + abs(Ez(ia,ja)) + abs(Ez(ie-iBC+ia,ja)) ; %right,left power
end
end
%}
% キА?? 快捷鍵說明
復(fù)制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -