这个程序在matlab中有错误啊
循环参数有问题，自己修正一下就行了 发信人: mygodest (AAA), 信区: Electromagnetics cs]
标题: 2-D FDTD TE code
发信站: 飘渺水云间 (Sun Oct 10 10:01:31 2004), 转信
请不要侵犯版权！谢谢！
%***********************************************************************
% 2-D FDTD TE code with PML absorbing boundary conditions
%***********************************************************************
%
% Program author: Susan C. Hagness
% Department of Electrical and Computer Engineering
% University of Wisconsin-Madison
% 1415 Engineering Drive
% Madison, WI 53706-1691
% 608-265-5739
% hagness@engr.wisc.edu
%
% Date of this version:February 2000
%
% 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 metal cylindrical
% scatterer in free space is modeled. The source excitation is
% a Gaussian pulse with a carrier frequency of 5 GHz.
%
% The grid resolution (dx = 3 mm) was chosen to provide 20 samples
% per wavelength at the center frequency of the pulse (which in turn
% provides approximately 10 samples per wavelength at the high end
% of the excitation spectrum, around 10 GHz).
%
% The computational domain is truncated using the perfectly matched
% layer (PML) absorbing boundary conditions.The formulation used
% in this code is based on the original split-field Berenger PML. The
% PML regions are labeled as shown in the following diagram:
%
% ----------------------------------------------
% || BACK PML ||
% ----------------------------------------------
% |L | /| R|
% |E | (ib,jb) | I|
% |F | | G|
% |T | | H|
% || MAIN GRID | T|
% |P | ||
% |M | | P|
% |L | (1,1) | M|
% ||/ | L|
% ----------------------------------------------
% || FRONT PML ||
% ----------------------------------------------
%
% To execute this M-file, type "fdtd2D" at the MATLAB prompt.
% This M-file displays the FDTD-computed Ex, Ey, and Hz fields at
% every 4th time step, and records those frames in a movie matrix,
% M, which is played at the end of the simulation using the "movie"
% command.
%
%***********************************************************************
clear
%***********************************************************************
% Fundamental constants
%***********************************************************************
cc=2.99792458e8; %speed of light in free space
muz=4.0*pi*1.0e-7; %permeability of free space
epsz=1.0/(cc*cc*muz); %permittivity of free space
freq=5.0e+9; %center frequency of source excitation
lambda=cc/freq; %center wavelength of source excitation
omega=2.0*pi*freq;
%***********************************************************************
% Grid parameters
%***********************************************************************
ie=100; %number of grid cells in x-direction
je=50; %number of grid cells in y-direction
ib=ie+1;
jb=je+1;
is=15; %location of z-directed hard source
js=je/2; %location of z-directed hard source
dx=3.0e-3; %space increment of square lattice
dt=dx/(2.0*cc); %time step
nmax=300; %total number of time steps
iebc=8; %thickness of left and right PML region
jebc=8; %thickness of front and back PML region
rmax=0.00001;
orderbc=2;
ibbc=iebc+1;
jbbc=jebc+1;
iefbc=ie+2*iebc;
jefbc=je+2*jebc;
ibfbc=iefbc+1;
jbfbc=jefbc+1;
%***********************************************************************
% Material parameters
%***********************************************************************
media=2;
eps=;
sig=;
mur=;
sim=;
%***********************************************************************
% Wave excitation
%***********************************************************************
rtau=160.0e-12;
tau=rtau/dt;
delay=3*tau;
source=zeros(1,nmax);
for n=1:7.0*tau
source(n)=sin(omega*(n-delay)*dt)*exp(-((n-delay)^2/tau^2));
end
%***********************************************************************
% Field arrays
%***********************************************************************
ex=zeros(ie,jb); %fields in main grid
ey=zeros(ib,je);
hz=zeros(ie,je);
exbcf=zeros(iefbc,jebc); %fields in front PML region
eybcf=zeros(ibfbc,jebc);
hzxbcf=zeros(iefbc,jebc);
hzybcf=zeros(iefbc,jebc);
exbcb=zeros(iefbc,jbbc); %fields in back PML region
eybcb=zeros(ibfbc,jebc);
hzxbcb=zeros(iefbc,jebc);
hzybcb=zeros(iefbc,jebc);
exbcl=zeros(iebc,jb); %fields in left PML region
eybcl=zeros(iebc,je);
hzxbcl=zeros(iebc,je);
hzybcl=zeros(iebc,je);
exbcr=zeros(iebc,jb); %fields in right PML region
eybcr=zeros(ibbc,je);
hzxbcr=zeros(iebc,je);
hzybcr=zeros(iebc,je);
%***********************************************************************
% 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
%***********************************************************************
% Geometry specification (main grid)
%***********************************************************************
% Initialize entire main grid to free space
caex(1:ie,1:jb)=ca(1);
cbex(1:ie,1:jb)=cb(1);
caey(1:ib,1:je)=ca(1);
cbey(1:ib,1:je)=cb(1);
dahz(1:ie,1:je)=da(1);
dbhz(1:ie,1:je)=db(1);
% Add metal cylinder
diam=20; % diameter of cylinder: 6 cm
rad=diam/2.0; % radius of cylinder: 3 cm
icenter=4*ie/5; % i-coordinate of cylinder's center
jcenter=je/2; % j-coordinate of cylinder's center
for i=1:ie
for j=1:je
dist2=(i+0.5-icenter)^2 + (j-jcenter)^2;
if dist2 <= rad^2
caex(i,j)=ca(2);
cbex(i,j)=cb(2);
end
dist2=(i-icenter)^2 + (j+0.5-jcenter)^2;
if dist2 <= rad^2
caey(i,j)=ca(2);
cbey(i,j)=cb(2);
end
end
end
%***********************************************************************
% Fill the PML regions
%***********************************************************************
delbc=iebc*dx;
sigmam=-log(rmax/100.0)*epsz*cc*(orderbc+1)/(2*delbc);
bcfactor=eps(1)*sigmam/(dx*(delbc^orderbc)*(orderbc+1));
% FRONT region
caexbcf(1:iefbc,1)=1.0;
cbexbcf(1:iefbc,1)=0.0;
for j=2:jebc
y1=(jebc-j+1.5)*dx;
y2=(jebc-j+0.5)*dx;
sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
ca1=exp(-sigmay*dt/(epsz*eps(1)));
cb1=(1.0-ca1)/(sigmay*dx);
caexbcf(1:iefbc,j)=ca1;
cbexbcf(1:iefbc,j)=cb1;
end
sigmay = bcfactor*(0.5*dx)^(orderbc+1);
ca1=exp(-sigmay*dt/(epsz*eps(1)));
cb1=(1-ca1)/(sigmay*dx);
caex(1:ie,1)=ca1;
cbex(1:ie,1)=cb1;
caexbcl(1:iebc,1)=ca1;
cbexbcl(1:iebc,1)=cb1;
caexbcr(1:iebc,1)=ca1;
cbexbcr(1:iebc,1)=cb1;
for j=1:jebc
y1=(jebc-j+1)*dx;
y2=(jebc-j)*dx;
sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
sigmays=sigmay*(muz/(epsz*eps(1)));
da1=exp(-sigmays*dt/muz);
db1=(1-da1)/(sigmays*dx);
dahzybcf(1:iefbc,j)=da1;
dbhzybcf(1:iefbc,j)=db1;
caeybcf(1:ibfbc,j)=ca(1);
cbeybcf(1:ibfbc,j)=cb(1);
dahzxbcf(1:iefbc,j)=da(1);
dbhzxbcf(1:iefbc,j)=db(1);
end
% BACK region
caexbcb(1:iefbc,jbbc)=1.0;
cbexbcb(1:iefbc,jbbc)=0.0;
for j=2:jebc
y1=(j-0.5)*dx;
y2=(j-1.5)*dx;
sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
ca1=exp(-sigmay*dt/(epsz*eps(1)));
cb1=(1-ca1)/(sigmay*dx);
caexbcb(1:iefbc,j)=ca1;
cbexbcb(1:iefbc,j)=cb1;
end
sigmay = bcfactor*(0.5*dx)^(orderbc+1);
ca1=exp(-sigmay*dt/(epsz*eps(1)));
cb1=(1-ca1)/(sigmay*dx);
caex(1:ie,jb)=ca1;
cbex(1:ie,jb)=cb1;
caexbcl(1:iebc,jb)=ca1;
cbexbcl(1:iebc,jb)=cb1;
caexbcr(1:iebc,jb)=ca1;
cbexbcr(1:iebc,jb)=cb1;
for j=1:jebc
y1=j*dx;
y2=(j-1)*dx;
sigmay=bcfactor*(y1^(orderbc+1)-y2^(orderbc+1));
sigmays=sigmay*(muz/(epsz*eps(1)));
da1=exp(-sigmays*dt/muz);
db1=(1-da1)/(sigmays*dx);
dahzybcb(1:iefbc,j)=da1;
dbhzybcb(1:iefbc,j)=db1;
caeybcb(1:ibfbc,j)=ca(1);
cbeybcb(1:ibfbc,j)=cb(1);
dahzxbcb(1:iefbc,j)=da(1);
dbhzxbcb(1:iefbc,j)=db(1);
end
% LEFT region
caeybcl(1,1:je)=1.0;
cbeybcl(1,1:je)=0.0;
for i=2:iebc
x1=(iebc-i+1.5)*dx;
x2=(iebc-i+0.5)*dx;
sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
ca1=exp(-sigmax*dt/(epsz*eps(1)));
cb1=(1-ca1)/(sigmax*dx);
caeybcl(i,1:je)=ca1;
cbeybcl(i,1:je)=cb1;
caeybcf(i,1:jebc)=ca1;
cbeybcf(i,1:jebc)=cb1;
caeybcb(i,1:jebc)=ca1;
cbeybcb(i,1:jebc)=cb1;
end
sigmax=bcfactor*(0.5*dx)^(orderbc+1);
ca1=exp(-sigmax*dt/(epsz*eps(1)));
cb1=(1-ca1)/(sigmax*dx);
caey(1,1:je)=ca1;
cbey(1,1:je)=cb1;
caeybcf(iebc+1,1:jebc)=ca1;
cbeybcf(iebc+1,1:jebc)=cb1;
caeybcb(iebc+1,1:jebc)=ca1;
cbeybcb(iebc+1,1:jebc)=cb1;
for i=1:iebc
x1=(iebc-i+1)*dx;
x2=(iebc-i)*dx;
sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
sigmaxs=sigmax*(muz/(epsz*eps(1)));
da1=exp(-sigmaxs*dt/muz);
db1=(1-da1)/(sigmaxs*dx);
dahzxbcl(i,1:je)=da1;
dbhzxbcl(i,1:je)=db1;
dahzxbcf(i,1:jebc)=da1;
dbhzxbcf(i,1:jebc)=db1;
dahzxbcb(i,1:jebc)=da1;
dbhzxbcb(i,1:jebc)=db1;
caexbcl(i,2:je)=ca(1);
cbexbcl(i,2:je)=cb(1);
dahzybcl(i,1:je)=da(1);
dbhzybcl(i,1:je)=db(1);
end
% RIGHT region
caeybcr(ibbc,1:je)=1.0;
cbeybcr(ibbc,1:je)=0.0;
for i=2:iebc
x1=(i-0.5)*dx;
x2=(i-1.5)*dx;
sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
ca1=exp(-sigmax*dt/(epsz*eps(1)));
cb1=(1-ca1)/(sigmax*dx);
caeybcr(i,1:je)=ca1;
cbeybcr(i,1:je)=cb1;
caeybcf(i+iebc+ie,1:jebc)=ca1;
cbeybcf(i+iebc+ie,1:jebc)=cb1;
caeybcb(i+iebc+ie,1:jebc)=ca1;
cbeybcb(i+iebc+ie,1:jebc)=cb1;
end
sigmax=bcfactor*(0.5*dx)^(orderbc+1);
ca1=exp(-sigmax*dt/(epsz*eps(1)));
cb1=(1-ca1)/(sigmax*dx);
caey(ib,1:je)=ca1;
cbey(ib,1:je)=cb1;
caeybcf(iebc+ib,1:jebc)=ca1;
cbeybcf(iebc+ib,1:jebc)=cb1;
caeybcb(iebc+ib,1:jebc)=ca1;
cbeybcb(iebc+ib,1:jebc)=cb1;
for i=1:iebc
x1=i*dx;
x2=(i-1)*dx;
sigmax=bcfactor*(x1^(orderbc+1)-x2^(orderbc+1));
sigmaxs=sigmax*(muz/(epsz*eps(1)));
da1=exp(-sigmaxs*dt/muz);
db1=(1-da1)/(sigmaxs*dx);
dahzxbcr(i,1:je) = da1;
dbhzxbcr(i,1:je) = db1;
dahzxbcf(i+ie+iebc,1:jebc)=da1;
dbhzxbcf(i+ie+iebc,1:jebc)=db1;
dahzxbcb(i+ie+iebc,1:jebc)=da1;
dbhzxbcb(i+ie+iebc,1:jebc)=db1;
caexbcr(i,2:je)=ca(1);
cbexbcr(i,2:je)=cb(1);
dahzybcr(i,1:je)=da(1);
dbhzybcr(i,1:je)=db(1);
end
%***********************************************************************
% Movie initialization
%***********************************************************************
subplot(3,1,1),pcolor(ex');
shading flat;
caxis([-80.0 80.0]);
axis();
colorbar;
axis image;
axis off;
title(['Ex at time step = 0']);
subplot(3,1,2),pcolor(ey');
shading flat;
caxis([-80.0 80.0]);
axis();
colorbar;
axis image;
axis off;
title(['Ey at time step = 0']);
subplot(3,1,3),pcolor(hz');
shading flat;
caxis([-0.2 0.2]);
axis();
colorbar;
axis image;
axis off;
title(['Hz at time step = 0']);
rect=get(gcf,'Position');
rect(1:2)=;
M=moviein(nmax/4,gcf,rect);
%***********************************************************************
% BEGIN TIME-STEPPING LOOP
%***********************************************************************
for n=1:nmax
%***********************************************************************
% Update electric fields (EX and EY) in main grid
%***********************************************************************
ex(:,2:je)=caex(:,2:je).*ex(:,2:je)+...
cbex(:,2:je).*(hz(:,2:je)-hz(:,1:je-1));
ey(2:ie,:)=caey(2:ie,:).*ey(2:ie,:)+...
cbey(2:ie,:).*(hz(1:ie-1,:)-hz(2:ie,:));
%***********************************************************************
% Update EX in PML regions
%***********************************************************************
% FRONT
exbcf(:,2:jebc)=caexbcf(:,2:jebc).*exbcf(:,2:jebc)-...
cbexbcf(:,2:jebc).*(hzxbcf(:,1:jebc-1)+hzybcf(:,1:jebc-1)-...
hzxbcf(:,2:jebc)-hzybcf(:,2:jebc));
ex(1:ie,1)=caex(1:ie,1).*ex(1:ie,1)-...
cbex(1:ie,1).*(hzxbcf(ibbc:iebc+ie,jebc)+...
hzybcf(ibbc:iebc+ie,jebc)-hz(1:ie,1));
% BACK
exbcb(:,2:jebc-1)=caexbcb(:,2:jebc-1).*exbcb(:,2:jebc-1)-...
cbexbcb(:,2:jebc-1).*(hzxbcb(:,1:jebc-2)+hzybcb(:,1:jebc-2)-...
hzxbcb(:,2:jebc-1)-hzybcb(:,2:jebc-1));
ex(1:ie,jb)=caex(1:ie,jb).*ex(1:ie,jb)-...
cbex(1:ie,jb).*(hz(1:ie,jb-1)-hzxbcb(ibbc:iebc+ie,1)-...
hzybcb(ibbc:iebc+ie,1));
% LEFT
exbcl(:,2:je)=caexbcl(:,2:je).*exbcl(:,2:je)-...
cbexbcl(:,2:je).*(hzxbcl(:,1:je-1)+hzybcl(:,1:je-1)-...
hzxbcl(:,2:je)-hzybcl(:,2:je));
exbcl(:,1)=caexbcl(:,1).*exbcl(:,1)-...
cbexbcl(:,1).*(hzxbcf(1:iebc,jebc)+hzybcf(1:iebc,jebc)-...
hzxbcl(:,1)-hzybcl(:,1));
exbcl(:,jb)=caexbcl(:,jb).*exbcl(:,jb)-...
cbexbcl(:,jb).*(hzxbcl(:,je)+hzybcl(:,je)-...
hzxbcb(1:iebc,1)-hzybcb(1:iebc,1));
% RIGHT
exbcr(:,2:je)=caexbcr(:,2:je).*exbcr(:,2:je)-...
cbexbcr(:,2:je).*(hzxbcr(:,1:je-1)+hzybcr(:,1:je-1)-...
hzxbcr(:,2:je)-hzybcr(:,2:je));
exbcr(:,1)=caexbcr(:,1).*exbcr(:,1)-...
cbexbcr(:,1).*(hzxbcf(1+iebc+ie:iefbc,jebc)+...
hzybcf(1+iebc+ie:iefbc,jebc)-...
hzxbcr(:,1)-hzybcr(:,1));
exbcr(:,jb)=caexbcr(:,jb).*exbcr(:,jb)-...
cbexbcr(:,jb).*(hzxbcr(:,je)+hzybcr(:,je)-...
hzxbcb(1+iebc+ie:iefbc,1)-...
hzybcb(1+iebc+ie:iefbc,1));
%***********************************************************************
% Update EY in PML regions
%***********************************************************************
% FRONT
eybcf(2:iefbc,:)=caeybcf(2:iefbc,:).*eybcf(2:iefbc,:)-...
cbeybcf(2:iefbc,:).*(hzxbcf(2:iefbc,:)+hzybcf(2:iefbc,:)-...
hzxbcf(1:iefbc-1,:)-hzybcf(1:iefbc-1,:));
% BACK
eybcb(2:iefbc,:)=caeybcb(2:iefbc,:).*eybcb(2:iefbc,:)-...
cbeybcb(2:iefbc,:).*(hzxbcb(2:iefbc,:)+hzybcb(2:iefbc,:)-...
hzxbcb(1:iefbc-1,:)-hzybcb(1:iefbc-1,:));
% LEFT
eybcl(2:iebc,:)=caeybcl(2:iebc,:).*eybcl(2:iebc,:)-...
cbeybcl(2:iebc,:).*(hzxbcl(2:iebc,:)+hzybcl(2:iebc,:)-...
hzxbcl(1:iebc-1,:)-hzybcl(1:iebc-1,:));
ey(1,:)=caey(1,:).*ey(1,:)-...
cbey(1,:).*(hz(1,:)-hzxbcl(iebc,:)-hzybcl(iebc,:));
% RIGHT
eybcr(2:iebc,:)=caeybcr(2:iebc,:).*eybcr(2:iebc,:)-...
cbeybcr(2:iebc,:).*(hzxbcr(2:iebc,:)+hzybcr(2:iebc,:)-...
hzxbcr(1:iebc-1,:)-hzybcr(1:iebc-1,:));
ey(ib,:)=caey(ib,:).*ey(ib,:)-...
cbey(ib,:).*(hzxbcr(1,:)+hzybcr(1,:)- hz(ie,:));
%***********************************************************************
% Update magnetic fields (HZ) in main grid
%***********************************************************************
hz(1:ie,1:je)=dahz(1:ie,1:je).*hz(1:ie,1:je)+...
dbhz(1:ie,1:je).*(ex(1:ie,2:jb)-ex(1:ie,1:je)+...
ey(1:ie,1:je)-ey(2:ib,1:je));
hz(is,js)=source(n);
%***********************************************************************
% Update HZX in PML regions
%***********************************************************************
% FRONT
hzxbcf(1:iefbc,:)=dahzxbcf(1:iefbc,:).*hzxbcf(1:iefbc,:)-...
dbhzxbcf(1:iefbc,:).*(eybcf(2:ibfbc,:)-eybcf(1:iefbc,:));
% BACK
hzxbcb(1:iefbc,:)=dahzxbcb(1:iefbc,:).*hzxbcb(1:iefbc,:)-...
dbhzxbcb(1:iefbc,:).*(eybcb(2:ibfbc,:)-eybcb(1:iefbc,:));
% LEFT
hzxbcl(1:iebc-1,:)=dahzxbcl(1:iebc-1,:).*hzxbcl(1:iebc-1,:)-...
dbhzxbcl(1:iebc-1,:).*(eybcl(2:iebc,:)-eybcl(1:iebc-1,:));
hzxbcl(iebc,:)=dahzxbcl(iebc,:).*hzxbcl(iebc,:)-...
dbhzxbcl(iebc,:).*(ey(1,:)-eybcl(iebc,:));
% RIGHT
hzxbcr(2:iebc,:)=dahzxbcr(2:iebc,:).*hzxbcr(2:iebc,:)-...
dbhzxbcr(2:iebc,:).*(eybcr(3:ibbc,:)-eybcr(2:iebc,:));
hzxbcr(1,:)=dahzxbcr(1,:).*hzxbcr(1,:)-...
dbhzxbcr(1,:).*(eybcr(2,:)-ey(ib,:));
%***********************************************************************
% Update HZY in PML regions
%***********************************************************************
% FRONT
hzybcf(:,1:jebc-1)=dahzybcf(:,1:jebc-1).*hzybcf(:,1:jebc-1)-...
dbhzybcf(:,1:jebc-1).*(exbcf(:,1:jebc-1)-exbcf(:,2:jebc));
hzybcf(1:iebc,jebc)=dahzybcf(1:iebc,jebc).*hzybcf(1:iebc,jebc)-...
dbhzybcf(1:iebc,jebc).*(exbcf(1:iebc,jebc)-exbcl(1:iebc,1));
hzybcf(iebc+1:iebc+ie,jebc)=...
dahzybcf(iebc+1:iebc+ie,jebc).*hzybcf(iebc+1:iebc+ie,jebc)-...
dbhzybcf(iebc+1:iebc+ie,jebc).*(exbcf(iebc+1:iebc+ie,jebc)-...
ex(1:ie,1));
hzybcf(iebc+ie+1:iefbc,jebc)=...
dahzybcf(iebc+ie+1:iefbc,jebc).*hzybcf(iebc+ie+1:iefbc,jebc)-...
dbhzybcf(iebc+ie+1:iefbc,jebc).*(exbcf(iebc+ie+1:iefbc,jebc)-...
exbcr(1:iebc,1));
% BACK
hzybcb(1:iefbc,2:jebc)=dahzybcb(1:iefbc,2:jebc).*hzybcb(1:iefbc,2:jebc)-...
dbhzybcb(1:iefbc,2:jebc).*(exbcb(1:iefbc,2:jebc)-exbcb(1:iefbc,3:jbbc));
hzybcb(1:iebc,1)=dahzybcb(1:iebc,1).*hzybcb(1:iebc,1)-...
dbhzybcb(1:iebc,1).*(exbcl(1:iebc,jb)-exbcb(1:iebc,2));
hzybcb(iebc+1:iebc+ie,1)=...
dahzybcb(iebc+1:iebc+ie,1).*hzybcb(iebc+1:iebc+ie,1)-...
dbhzybcb(iebc+1:iebc+ie,1).*(ex(1:ie,jb)-exbcb(iebc+1:iebc+ie,2));
hzybcb(iebc+ie+1:iefbc,1)=...
dahzybcb(iebc+ie+1:iefbc,1).*hzybcb(iebc+ie+1:iefbc,1)-...
dbhzybcb(iebc+ie+1:iefbc,1).*(exbcr(1:iebc,jb)-...
exbcb(iebc+ie+1:iefbc,2));
% LEFT
hzybcl(:,1:je)=dahzybcl(:,1:je).*hzybcl(:,1:je)-...
dbhzybcl(:,1:je).*(exbcl(:,1:je)-exbcl(:,2:jb));
% RIGHT
hzybcr(:,1:je)=dahzybcr(:,1:je).*hzybcr(:,1:je)-...
dbhzybcr(:,1:je).*(exbcr(:,1:je)-exbcr(:,2:jb));
%***********************************************************************
% Visualize fields
%***********************************************************************
if mod(n,4)==0;
timestep=int2str(n);
subplot(3,1,1),pcolor(ex');
shading flat;
caxis([-80.0 80.0]);
axis();
colorbar;
axis image;
axis off;
title(['Ex at time step = ',timestep]);
subplot(3,1,2),pcolor(ey');
shading flat;
caxis([-80.0 80.0]);
axis();
colorbar;
axis image;
axis off;
title(['Ey at time step = ',timestep]);
subplot(3,1,3),pcolor(hz');
shading flat;
caxis([-0.2 0.2]);
axis();
colorbar;
axis image;
axis off;
title(['Hz at time step = ',timestep]);
nn=n/4;
M(:,nn)=getframe(gcf,rect);
end;
%***********************************************************************
% END TIME-STEPPING LOOP
nn=n/4;
M(:,nn)=getframe(gcf,rect);
end;
%***********************************************************************
% END TIME-STEPPING LOOP
%***********************************************************************
end
movie(gcf,M,0,10,rect);
--
·╔══╗╭═══╮╔═╗╔╗╔═══╗ ╭╮ ·
·║ ║║╭╮║║╰╯║║ ═╣ ╱╱ ( ︶ ) ·
·║ ╚╗║╰╯║║ ║║ ═╣ ∠╱╔╮╮ ╲╱ ·
·╚═══╝╰═══╯╰═══╯╚═══╝〆 ╠╣║╭╦═╮╔╮╮ ·
· 用心笔，写下心语： ╰╩═╣╠╣║╠╣║ ·
· 爱你，不长，就一生。。。 ╰══╯╰╩═╯╰╩═╰ ·
※ 来源:·飘渺水云间 freecity.cn·
一维的我倒是写了，二维和三维却没有 不如传上去吧，是一个计算一维光子晶体透射系数的程序，ＭＡＴＬＡＢ的
function =photostal_1D_simpson(N,D1,epsilon1,D2,epsilon2,f_width,forDB)
% PHOTOSTAL_1D calculates the transmission spectrum
% of 1D photonic crystal,which can not be only single layer.
% In calculating the Fourier Transform,we use Simpson's Rule
%
% INPUTS: photostal_1D(N,D1,epsilon1,D2,epsilon2,f_width,forDB)
%
% N:Number of the layers,which can be either odd or even.
% D1:Thickness of first material
% epsilon1:Dielectric constant of first material
% D2:Thickness of second material
% epsilon2:Dielectric constant of second material
% f_width:Specifies the frequencies we are interest in,
% it must be a 1x2 or 2x1 vector ,the first element
% specifies the low frequency,and the second element specifies
% the up frequency.they are given through the form below
% a1<=(frequency*D/c)<=a2,where D=D1+D2 is the lattice constant
% and c is the light speed in vacuum
% forDB:true if we will use the dB plot
%
% OUTPUTS:
%
% freq:frequency,which is a1*c/D<=f<=a2*c/D
% transmission:transmission spectrum of the PC over frequency f
% D:Lattice constant
%
c=3e8;
D=D1+D2;
epsilon=max(epsilon1,epsilon2);
% Argument check
if nargin<5
error('Input parameters must be at least 5!!!');
end
if nargin<6
f_width=;
end
if nargin<7
forDB=0;
end
=size(f_width);
if ~((row==1&col==2)|(row==2&col==1))
error('The sixth argument must be 1x2 or 2x1 vector');
end
f_low=min(f_width(1),f_width(2));
f_up=max(f_width(1),f_width(2));
f_b=f_up-f_low;
if f_low==0
f_low=1e-4*f_up;
end
Nf=ceil(70*N*f_b);
if Nf<150
Nf=150;
elseif Nf>800
Nf=800;
end
freq=linspace(f_low*c/D,f_up*c/D,Nf);
freq_band=max(freq)-min(freq);
% Define the FDTD cells dimensions
lambda_min=c/sqrt(epsilon)/max(freq);
dx=lambda_min/15;
k0=50;
k_1=ceil(D1/dx);% Spatial steps of 1st material
dx1=D1/k_1;
k_2=ceil(D2/dx);% Spatial steps of 2nd material
dx2=D2/k_2;
dx=max(dx1,dx2);
dt=dx/2/c;
k_end=k_1+k_2;
% Specify the incident pulse
freq_in=(max(freq)+min(freq))/2; % Center frequency of the incident pulse
T=0.5/freq_band;
T0=3*T/dt;% Center of Gaussian pulse at T0 time step
% The Epsilon Matrix
Epsilon=ones(1,k_end);
Epsilon(1:k_1)=epsilon1;
Epsilon(k_1+1:k_end)=epsilon2;
% The coefficient of c*dt/dx
ca=0.5*ones(1,k_end);
ca(1:k_1)=c*dt/dx1;
ca(k_1+1:k_end)=c*dt/dx2;
if N==1
N=2;
end
N0=floor(N/2);
Epsilon=repmat(Epsilon,1,N0);
ca=repmat(ca,1,N0);
if mod(N,2)==1 % N is an odd number
Epsilon=;
ca=;
else
Epsilon=;
ca=;
end
ca=ca./Epsilon;
r=0.5;
factor=(1-r)/(1+r);
KE=length(Epsilon);
k_source=10;
k_sample_i=k0;% Point where we samples the incident pulse
k_sample_t=KE-k0;% Point where we samples the transmit pulse
% Initialize the E and H fields
Ex=zeros(1,KE);
Hy=zeros(1,KE-1);
Ex0=zeros(1,KE);
Hy0=zeros(1,KE-1);
% Parameters for the ABCs
Ex_Low_1=Ex(1,1);
Ex_Low_2=Ex(1,2);
Ex_High_1=Ex(1,end);
Ex_High_2=Ex(1,end-1);
Ex0_Low_1=Ex0(1,1);
Ex0_Low_2=Ex0(1,2);
Ex0_High_1=Ex0(1,end);
Ex0_High_2=Ex0(1,end-1);
% Total time steps estimation
T_Des=2*k0+(N0+1)*(sqrt(epsilon1)*k_1+sqrt(epsilon2)*k_2);
T_Des=6*round(sqrt(epsilon)*T_Des);
% Vectors to record the incident and transmit pulses
L_Ei=floor(8*T0);
real_in=zeros(1,Nf);
imag_in=zeros(1,Nf);
real_tr=zeros(1,Nf);
imag_tr=zeros(1,Nf);
n=1;
disp('Computing the Transmission spectrum...')
disp(['Estimate time steps is ',num2str(T_Des)])
disp('This process may take several seconds or several minutes')
disp('or even tens of minutes !!!')
disp('Please wait patiently')
disp('......')
% Before FDTD iteration,clear the variables
% no longer needed
clear N D1 D2 epsilon1 epsilon2 f_width epsilon
clear row col f_b freq_band
clear dx Nf N0 KE k0 k_1 k_2 k_end T T_Des
clear Epsilon lambda_min
sc=2; % Coefficient for Simpson sum
while 1
% Update the Ex fields
Ex(2:end-1)=Ex(2:end-1)+ca(2:end-1).*(Hy(1:end-1)-Hy(2:end));
if n<=L_Ei
% Update the Ex0 fields
Ex0(2:end-1)=Ex0(2:end-1)+r*(Hy0(1:end-1)-Hy0(2:end));
pulse=exp(-(n-T0)^2*9/T0^2)*sin(2*pi*freq_in*dt*n);
Ex(1,k_source)=Ex(1,k_source)+pulse;
Ex0(1,k_source)=Ex0(1,k_source)+pulse;
Ex0(1)=Ex0_Low_2+factor*(Ex0_Low_1-Ex0(2));
Ex0_Low_1=Ex0(1);
Ex0_Low_2=Ex0(2);
Ex0(end)=Ex0_High_2+factor*(Ex0_High_1-Ex0(end-1));
Ex0_High_1=Ex0(end);
Ex0_High_2=Ex0(end-1);
% Records the pulses and calculate the FT
real_in=real_in+sc*Ex0(1,k_sample_i)*cos(2*pi*freq*dt*n);
imag_in=imag_in+sc*Ex0(1,k_sample_i)*sin(2*pi*freq*dt*n);
% Update the Hy0 fields
Hy0=Hy0+r*(Ex0(1:end-1)-Ex0(2:end));
end
% Mur Absorbing Boundary Conditions
Ex(1)=Ex_Low_2+factor*(Ex_Low_1-Ex(2));
Ex_Low_1=Ex(1);
Ex_Low_2=Ex(2);
Ex(end)=Ex_High_2+factor*(Ex_High_1-Ex(end-1));
Ex_High_1=Ex(end);
Ex_High_2=Ex(end-1);
% Records the pulses and calculate the FT
real_tr=real_tr+sc*Ex(1,k_sample_t)*cos(2*pi*freq*dt*n);
imag_tr=imag_tr+sc*Ex(1,k_sample_t)*sin(2*pi*freq*dt*n);
% Change the coefficient sc
if sc==2
sc=4;
else
sc=2;
end
% Update the Hy fields
Hy=Hy+r*(Ex(1:end-1)-Ex(2:end));
n=n+1;
if (n>3*T0)&max(abs(Ex))<1e-4
break
end
end
% At this time ,there must be too many variables in
% the memory,so before calculate the Fourier transform,
% clear the ones no longer needed
clear Ex Ex0 Hy Hy0 ca
clear Ex_Low_1 Ex_Low_2 Ex0_Low_1 Ex0_Low_2
clear Ex_High_1 Ex_High_2 Ex0_High_1 Ex0_High_2
clear r pulse n T0 freq_in factor L_Ei
clear k_sample_i k_sample_t k_source
disp(blanks(1))
disp('FDTD iteration has complete')
disp('Now plots the transmission spectrum figure')
disp(blanks(1))
pause(0.01)
amp_in=abs(real_in+i*imag_in);
amp_tr=abs(real_tr+i*imag_tr);
transmission=amp_tr./amp_in;
if forDB
transmission=20*log10(transmission);
end
% Another time clear no longer used variables
clear real_in imag_in real_tr imag_tr dt forDB
% Now handle the output
t_low=min(transmission);
t_low=t_low/1.1;
t_up=max(transmission);
t_up=t_up*1.1;
W=300;H=200;
fig=figure('position',);
set(fig,'Name','Photostal_1D--->HuaiMin.Peng','NumberTitle','off');
plot(freq*D/c,transmission)
xlabel('\omegaa/2\pic','fontweight','bold')
ylabel('Author: Pheigenbaum','fontweight','bold')
title('Transmission Spectrum by Simpson','fontweight','bold')
axis()
grid on 其实到现在，上面的程序已经经过了很多次修改，那个版本的程序已经相当的落后了 原帖由 pheigenbau 于 2007-7-20 19:59 发表 http://www.chinavib.com/forum/images/common/back.gif
其实到现在，上面的程序已经经过了很多次修改，那个版本的程序已经相当的落后了
这个很正常 请教楼主一些问题，虽然已经都是两年前的帖了，还是占用楼主几分钟的时间，谢谢了，呵呵！
第一就是：PML边界条件中为什么选择8层，这个是根据什么定的，另外可以给出PML边界条件的迭代式么，呵呵
第二就是：十七楼的程序中，里面有好些参变量不知道表示什么意思，比如：ib,jb，iebc,jebc,orderbc,ibbc,iefbc,rtau,想问下这些命名是不是有规律，我看到好多程序里面貌似有相同之处，但是我不太清楚它们所代表的意思，因为是新手，初学编程，望给与指导，谢谢了，呵呵
页:
1
[2]