?? untitled.m
字號:
%***********************************************************************
% 2-D FDTD code with CPML absorbing boundary conditions
%***********************************************************************
%***********************************************************************
% fdtd2D_CPML
% ..................................
clear;
clear all;
clc;
%function [TP,eps0,dt,dxyz,freq,lambda,EzPMLpower]=fdtd2D_cpml_e01_newxx(ax)
% Input Fundamental Constants (MKS units)
c0 = 2.99792458E8;
mu0 = 4.0 * pi * 1.0E-7;
eps0 = 1.0/(c0*c0*mu0);
% ..................................
% Specify Material Relative Permittivity and Conductivity
epsR = 1.0;
sigM1 = 0.0 ; % free space
% ..................................
% Specify Number of Time Steps and Grid Size Parameters
MaxTime = 500; % total number of time steps
% ..................................
ie=300; % Size of main grid
je=300;
ke=300;
ih=ie+1;
jh=je+1;
kh=ke+1;
% ..................................
% Specify the CPML Thickness in Each Direction (Value of Zero
% Corresponds to No PML, and the Grid is Terminated with a PEC)
% PML thickness in each direction
iBC = 21;
jBC = iBC;
kBC = iBC;
%
ieT=ie+2*iBC; % Size of total computational domain
jeT=je+2*jBC;
keT=ke+2*kBC;
ihT=ieT+1;
jhT=jeT+1;
khT=keT+1;
% number of Ez field components %
Ex =zeros(ieT, jhT) ;
Ey =zeros(ihT, jeT);
Ez =zeros(ihT, jhT);
EzPMLpower=zeros(1,ie);
EzPMLpowert=zeros(1,ie);
%
Hx =zeros(ihT, jeT) ;
Hy =zeros(ieT, jhT);
Hz =zeros(ieT, jeT);
CA =zeros(ihT, jhT);
CB =zeros(ihT, jhT);
sig =zeros(ihT, jhT);
epsi=zeros(ihT, jhT);
% ..................................
% Specify Grid Cell Size in Each Direction and Calculate the
% Resulting Courant-Stable Time Step
lambda=5.5e-7; %wavelength bandwidth of source excitation
freq=c0/lambda; %frequency bandwidth of source excitation
omega=2.0*pi*freq;
dxyz = lambda/20; % cell size in each direction
dt = 0.99 / (c0*(1.0/dxyz^2+1.0/dxyz^2 )^0.5);
t_periodic=1/(freq*dt);
ha=1/freq;
ta=zeros(MaxTime,1); % time step increment
EzPMLpower=zeros(MaxTime,1);
for n=1:MaxTime
ta(n)=(n-1)*dt;
end
% ..................................
% Specify the Impulsive Source (See Equation 7.134)
tw = 53.0E-12; t0 = 4.0*tw ;
rtau=50.0e-12;
tau=rtau/dt;
delay=tau;
srcconst=-dt*3.0e+11;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ..................................
% Specify the Time Step at which the Grid is Recorded for Visualization
record_grid = 300 ;
NoFig=99;
Nplt_jmp=10;
% ..................................
% Source point.
isrc=round(ihT/2); % Location of source
jsrc=round(jhT/2);
%
iview= floor(ieT/2);
jview= floor(jeT/2);
% ..................................
% Specify the CPML Order and Other Parameters
m = 3; ma = 1 ;
sigXmax = (0.8*(m+1)/(dxyz*(mu0/eps0*epsR)^0.5));
sigYmax = (0.8*(m+1)/(dxyz*(mu0/eps0*epsR)^0.5));
sigZmax = (0.8*(m+1)/(dxyz*(mu0/eps0*epsR)^0.5));
aXmax = 0.0003;
aYmax = aXmax;
aZmax = aXmax;
kXmax = 15.0;
kYmax = kXmax;
kZmax = kXmax;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CPML (7.100)
psi_Exy1=zeros(ieT,jBC) ;
psi_Exy2=zeros(ieT,jBC);
psi_Exz1=zeros(ieT,jhT);
psi_Exz2=zeros(ieT,jhT) ;
psi_Eyx1=zeros(iBC,jeT);
psi_Eyx2=zeros(iBC,jeT);
psi_Eyz1=zeros(ihT,jeT);
psi_Eyz2=zeros(ihT,jeT) ;
psi_Ezx1=zeros(iBC,jhT);
psi_Ezx2=zeros(iBC,jhT) ;
psi_Ezy1=zeros(ihT,jBC);
psi_Ezy2=zeros(ihT,jBC);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
psi_Hxy1=zeros(ihT,jBC-1) ;
psi_Hxy2=zeros(ihT,jBC-1) ;
psi_Hxz1=zeros(ihT,jeT);
psi_Hxz2=zeros(ihT,jeT) ;
psi_Hyz1=zeros(ieT,jhT) ;
psi_Hyz2=zeros(ieT,jhT);
psi_Hyx1=zeros(iBC-1,jhT);
psi_Hyx2=zeros(iBC-1,jhT) ;
psi_Hzx1=zeros(iBC-1,jeT);
psi_Hzx2=zeros(iBC-1,jeT);
psi_Hzy1=zeros(ieT,jBC-1) ;
psi_Hzy2=zeros(ieT,jBC-1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
bEx=zeros(iBC,1);
cEx=zeros(iBC,1);
A_ExBC=zeros(iBC,1);
sigExBC=zeros(iBC,1);
K_ExBC=zeros(iBC,1);
bEy=zeros(jBC,1);
cEy=zeros(jBC,1);
A_EyBC=zeros(jBC,1);
sigEyBC=zeros(jBC,1);
K_EyBC=zeros(jBC,1);
bEz=zeros(kBC,1) ;
cEz=zeros(kBC,1) ;
A_EzBC=zeros(kBC,1) ;
sigEzBC=zeros(kBC,1) ;
K_EzBC=zeros(kBC,1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (7.105)
bHx=zeros(iBC-1,1);
cHx=zeros(iBC-1,1);
A_HxBC=zeros(iBC-1,1);
sigHx_BC=zeros(iBC-1,1);
K_HxBC=zeros(iBC-1,1);
bHy=zeros(jBC-1,1) ;
cHy=zeros(jBC-1,1) ;
A_HyBC=zeros(jBC-1,1) ;
sigHy_BC=zeros(jBC-1,1) ;
K_HyBC=zeros(jBC-1,1) ;
bHz=zeros(kBC-1,1) ;
cHz=zeros(kBC-1,1) ;
A_HzBC=zeros(kBC-1,1) ;
sigHz_BC=zeros(kBC-1) ;
K_HzBC=zeros(kBC-1,1) ;
% denominators for the update equations
% (7.109)
F_ex=zeros(ieT,1);
F_hx=zeros(ieT,1);
F_ey=zeros(jeT,1);
F_hy=zeros(jeT,1);
F_ez=zeros(keT,1);
F_hz=zeros(keT,1);
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% INITIALIZE VARIABLES
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
sig(:,:) = sigM1;
sig(1:isrc-18,jBC+90:jBC+120)=2e+12;
sig(isrc+18:ihT,jBC+90:jBC+120)=2e+12;
sig(1:isrc-18,jBC+160:jBC+190)=2e+12;
sig(isrc+18:ihT,jBC+160:jBC+190)=2e+12;
epsi(:,:)=epsR*eps0;
%{
neff=1.5;
b_periodic=floor(lambda/(2*neff*dxyz));
delta_n=neff*(0.0001);
for ia=1:ihT
for ja=1:jhT
if ia>=40 && ia<=80 && ja>=30 && ja<=90
epsi(ia,ja)=eps0*(neff+delta_n*(1+cos(2*pi*ia/b_periodic)))^2;
elseif ia>=26 && ia<40 && ja>=30 && ja<=90
epsi(ia,ja)=(neff*(1-0.003))^2*eps0;
elseif ia>80 && ia<=94 && ja>=30 && ja<=90
epsi(ia,ja)=(neff*(1-0.003))^2*eps0;
else
epsi(ia,ja)=epsR*eps0;
end
end
end
%}
% write(*,*)"ihT: ", ihT
% write(*,*)"jhT: ", jhT
% write(*,*)"dt: ", dt
% write(*,*)"MaxTime : ", MaxTime
% write(*,*)"max time: ", MaxTime *dt
% write(*,*)"record grid after ", record_grid, "dt"
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% SET CPML PARAMETERS IN EACH DIRECTION
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% x-dir
%
for i = 1:iBC
sigExBC(i) = sigXmax * ( (iBC - i) / (iBC - 1.0) )^m ;
K_ExBC(i) = 1.0+(kXmax-1.0)*((iBC - i) / (iBC - 1.0))^m ;
A_ExBC(i) = aXmax*((i-1.0)/(iBC-1.0))^ma ;%(7.79)
bEx(i) = exp(-(sigExBC(i) / K_ExBC(i) + A_ExBC(i))*dt/eps0) ;
if ( (sigExBC(i) == 0.0) && (A_ExBC(i) == 0.0) && (i == iBC) )
cEx(i) = 0.0;
else
cEx(i) = sigExBC(i)*(bEx(i)-1.0)/(sigExBC(i)+K_ExBC(i)*A_ExBC(i)) / K_ExBC(i);
end
end
for i = 1:iBC-1
sigHx_BC(i) = sigXmax * ( (iBC - i - 0.5)/(iBC-1.0))^m;
K_HxBC(i) = 1.0+(kXmax-1.0)*((iBC - i - 0.5) / (iBC - 1.0))^m ;
A_HxBC(i) = aXmax*((i-0.5)/(iBC-1.0))^ma;
bHx(i) = exp(-(sigHx_BC(i) / K_HxBC(i) + A_HxBC(i))*dt/eps0) ;
cHx(i) = sigHx_BC(i)*(bHx(i)-1.0)/ (sigHx_BC(i)+K_HxBC(i)*A_HxBC(i)) / K_HxBC(i);
end
% y-dir
%
for j = 1:jBC
sigEyBC(j) = sigYmax * ( (jBC - j ) / (jBC - 1.0) )^m;
K_EyBC(j) = 1.0+(kYmax-1.0)* ((jBC - j) / (jBC - 1.0))^m;
A_EyBC(j) = aYmax*((j-1)/(jBC-1.0))^ma;
bEy(j) = exp(-(sigEyBC(j) / K_EyBC(j) + A_EyBC(j))*dt/eps0);
if ( (sigEyBC(j) == 0.0) && (A_EyBC(j) == 0.0) && (j == jBC))
cEy(j) = 0.0;
else
cEy(j) = sigEyBC(j)*(bEy(j)-1.0)/(sigEyBC(j)+K_EyBC(j)*A_EyBC(j)) / K_EyBC(j);
end
end
for j = 1:jBC-1
sigHy_BC(j) = sigYmax * ( (jBC - j - 0.5)/(jBC-1.0))^m;
K_HyBC(j) = 1.0+(kYmax-1.0)*((jBC - j - 0.5) / (jBC - 1.0))^m;
A_HyBC(j) = aYmax*((j-0.5)/(jBC-1.0))^ma;
bHy(j) = exp(-(sigHy_BC(j) / K_HyBC(j) + A_HyBC(j))*dt/eps0);
cHy(j) = sigHy_BC(j)*(bHy(j)-1.0)/ (sigHy_BC(j)+K_HyBC(j)*A_HyBC(j)) / K_HyBC(j);
end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% FILL IN UPDATING COEFFICIENTS
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
DA = 1.0; %(7-108)
DB = dt/mu0;
% (7.107)
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -