%***********************************************************************
% 2-D FDTD TM code with UPML absorbing
% boundary conditions.
%***********************************************************************
%
%
% J.Y.Zhu,2013/09/05
%
% This MATLAB M-file implements the finite-difference time-domain
% solution of Maxwell's curl equations over a two-dimensional
% Cartesian space lattice comprised of uniform square grid cells.
%
% To illustrate the algorithm, a 6-cm-diameter pec cylindrical
% scatterer in free space is modeled. The source excitation is
% a modulated Gaussian pulse.
%
% The computational domain is truncated using Mur's second order and
% UPML absorbing boundary conditions.
%
% To execute this M-file, type "fdtd_TM" at the MATLAB prompt.
% This M-file displays the FDTD-computed Hx, Hy, and Ez fields at
% every 5th time step, and records those frames in a movie matrix,
% M, which is played at the end of the simulation using the "movie"
% command.
%
%***********************************************************************
close all;
clear all;
%***********************************************************************
%% Calculation parameters
%***********************************************************************
% the angle between incident direction and x-axis, for example,
% if the incident direction is \hat{x}, Phi_i = 0
Phi_i = 45 /180*pi;
ki = [cos(Phi_i), sin(Phi_i)];
Phi_s = Phi_i+pi;
ks = [cos(Phi_s), sin(Phi_s)];
%***********************************************************************
%% Fundamental constants
%***********************************************************************
muz=4.0*pi*1.0e-7; %permeability of free space
epsz=1.0/(35950207149.4727056*pi); %permittivity of free space
cc=1/sqrt(muz*epsz); %speed of light in free space
z0 = sqrt(muz/epsz);
freq = 20e9; %max frequency of plane wave excitation
freq0 = 12e9;
lambda=cc/freq; %minimum wavelength of plane wave excitation
omega=2.0*pi*freq;
%***********************************************************************
%% Grid parameters
%***********************************************************************
% Size of main grid
Nx= 350; %number of grid cells in x-direction of fdtd region
Ny= 350; %number of grid cells in y-direction of fdtd region
TF_SF = 50;
TF_OB = 40;
%define the boundary of the total-field and scatter field
Nx_SF_min = TF_SF;
Nx_SF_max = Nx - TF_SF;
Ny_SF_min = TF_SF;
Ny_SF_max = Ny - TF_SF;
i_sample_sf = linspace(Nx_SF_min,Nx_SF_max,Nx_SF_max-Nx_SF_min+1);
j_sample_sf = linspace(Ny_SF_min,Ny_SF_max,Ny_SF_max-Ny_SF_min+1);
%define the output boundary
Nx_OB_min = TF_OB;
Nx_OB_max = Nx - TF_OB;
Ny_OB_min = TF_OB;
Ny_OB_max = Ny - TF_OB;
OBxsize = Nx_OB_max-Nx_OB_min+1;
OBysize = Ny_OB_max-Ny_OB_min+1;
i_sample_ob = linspace(Nx_OB_min,Nx_OB_max,OBxsize);
j_sample_ob = linspace(Ny_OB_min,Ny_OB_max,OBysize);
Nxp1=Nx+1;
Nyp1=Ny+1;
upml = 10; % Thickness of UPML boundaries
iu_abc = upml+1;
ju_abc = upml+1;
% Size of total computational domain
Nxabc=Nx+2*upml;
Nyabc=Ny+2*upml;
Nxp1abc=Nxabc+1;
Nyp1abc=Nyabc+1;
%***********************************************************************
%% FDTD computation parameters
%***********************************************************************
dx=lambda/20; %space increment of square lattice
dt=dx/(2.0*cc); %time step
TratioN = 3;
nmax=Nx*TratioN; %total number of time steps
%***********************************************************************
%% Material parameters
%***********************************************************************
media=2;
eps=[1.0 1.0];
sig=[0.0 1.0e+24];
mur=[1.0 1.0];
sim=[0.0 0.0];
%***********************************************************************
%% Wave excitation
%***********************************************************************
tau = 1.7 / (freq-freq0);
delay=0.8*tau;
%choose the striking corner by the plane wave
if(Phi_i >= 0 && Phi_i < pi/2)
istrike = Nx_SF_min;
jstrike = Ny_SF_min;
obxmin = Nx_OB_min;
obymin = Ny_OB_min;
elseif (Phi_i >= pi/2 && Phi_i < pi)
istrike = Nx_SF_max;
jstrike = Ny_SF_min;
obxmin = Nx_OB_max;
obymin = Ny_OB_min;
elseif (Phi_i >= pi && Phi_i < 1.5*pi)
istrike = Nx_SF_max;
jstrike = Ny_SF_max;
obxmin = Nx_OB_max;
obymin = Ny_OB_max;
else
istrike = Nx_SF_min;
jstrike = Ny_SF_max;
obxmin = Nx_OB_min;
obymin = Ny_OB_max;
end
Wavesize = ceil(sqrt((Nx-TF_SF*2)^2+(Ny-TF_SF*2)^2)*TratioN)+4+1;
ez_inc = zeros(1,Wavesize);
hi_inc = zeros(1,Wavesize-1);
hx_inc = zeros(1,Wavesize-1);
hy_inc = zeros(1,Wavesize-1);
ez_inc_coff = dt/epsz/dx;
hi_inc_coff = dt/muz/dx;
%add plane wave at left boundary
y_inc_lm = (Nx_SF_min-0.5-istrike)*cos(Phi_i)+(j_sample_sf-jstrike)*sin(Phi_i);
index_inc_lm = floor(y_inc_lm)+3;
w_lm = y_inc_lm - index_inc_lm+3;
%add plane wave at right boundary
y_inc_rm = (Nx_SF_max+0.5-istrike)*cos(Phi_i)+(j_sample_sf-jstrike)*sin(Phi_i);
index_inc_rm = floor(y_inc_rm)+3;
w_rm = y_inc_rm - index_inc_rm+3;
%add plane wave at buttom boundary
y_inc_bm = (i_sample_sf-istrike)*cos(Phi_i)+(Ny_SF_min-0.5-jstrike)*sin(Phi_i);
index_inc_bm = floor(y_inc_bm)+3;
w_bm = y_inc_bm - index_inc_bm+3;
%add plane wave at top boundary
y_inc_tm = (i_sample_sf-istrike)*cos(Phi_i)+(Ny_SF_max+0.5-jstrike)*sin(Phi_i);
index_inc_tm = floor(y_inc_tm)+3;
w_tm = y_inc_tm - index_inc_tm+3;
%a plane wave at left boundary electric field
y_inc_le = (Nx_SF_min-istrike)*cos(Phi_i)+(j_sample_sf-jstrike)*sin(Phi_i);
index_inc_le = floor(y_inc_le)+3;
w_le = y_inc_le - index_inc_le+3;
%a plane wave at left boundary electric field
y_inc_re = (Nx_SF_max-istrike)*cos(Phi_i)+(j_sample_sf-jstrike)*sin(Phi_i);
index_inc_re = floor(y_inc_re)+3;
w_re = y_inc_re - index_inc_re+3;
% plane wave at bottum boundary electric field
y_inc_be = (i_sample_sf-istrike)*cos(Phi_i)+(Ny_SF_min-jstrike)*sin(Phi_i);
index_inc_be = floor(y_inc_be)+3;
w_be = y_inc_be - index_inc_be+3;
% plane wave at bottum boundary electric field
y_inc_te = (i_sample_sf-istrike)*cos(Phi_i)+(Ny_SF_max-jstrike)*sin(Phi_i);
index_inc_te = floor(y_inc_te)+3;
w_te = y_inc_te - index_inc_te+3;
%***********************************************************************
%% Field arrays
%***********************************************************************
hx=zeros(Nxp1abc,Nyabc); %fields in main grid
bx=zeros(Nxp1abc,Nyabc);
hy=zeros(Nxabc,Nyp1abc);
by=zeros(Nxabc,Nyp1abc);
ez=zeros(Nxp1abc,Nyp1abc);
dz=zeros(Nxp1abc,Nyp1abc);
ezabcb=zeros(Nxabc-1,1); %buffer for 2 order Mur ABC
ezabct=zeros(Nxabc-1,1);
ezabcl=zeros(1,Nyabc-1);
ezabcr=zeros(1,Nyabc-1);
%***********************************************************************
%% Updating coefficients
%***********************************************************************
for i=1:media
eaf =dt*sig(i)/(2.0*epsz*eps(i));
ca(i)=(1.0-eaf)/(1.0+eaf);
cb(i)=dt/epsz/eps(i)/dx/(1.0+eaf);
haf =dt*sim(i)/(2.0*muz*mur(i));
da(i)=(1.0-haf)/(1.0+haf);
db(i)=dt/muz/mur(i)/dx/(1.0+haf);
end
% Initialize update coefficient arrays of 2 order Mur
ratio1 = cc*dt/dx;%一阶mur吸收边界系数
coeabc1 = (ratio1-1)/(ratio1+1);
coeabc2 = cc*cc*muz*dt/(2*(cc*dt+dx));
coeabc_c = (cc*dt-sqrt(2)*dx)/(cc*dt+sqrt(2)*dx);
% Initialize update coefficient arrays of UPML
C1ez=zeros(size(ez));
C2ez=zeros(size(ez));
C3ez=zeros(size(ez));
C4ez=zeros(size(ez));
D1hx=zeros(size(hx));
D2hx=zeros(size(hx));
D5hx=zeros(size(hx));
D6hx=zeros(size(hx));
D1hy=zeros(size(hy));
D2hy=zer
评论1