发帖回复

862阅读
0回复

[求助]用matlab模拟Sullivan那本书里的三维介质球

离线kangyel

UID ：101723

注册：2012-11-16
登录：2013-10-25
发帖：27
等级：仿真新人

0楼发表于: 2012-12-17 16:58:29

有大侠用matlab模拟过Sullivan那本书里的三维介质球吗？小弟的程序不能随时间推进，真心求助！我的程序如下。

clc;clear;

%**********************************************************************

% fdtd4-2c program

%**********************************************************************

%Three dimensional FDTD program of a planewave impinging on a

%dielectric sphere

%The source is in position at (ic,jc,kc) cellof the problem space with

%PML absorbing boundary conditions.

%**********************************************************************

%***********************************************************************

% Fundamental constants

%***********************************************************************

pi =3.14159;

ddx =0.01 ; % cell size

epsz = 8.85e-12; %permittivity of free space

muz =4.0*pi*1.0e-7; %permeability offree space

%***********************************************************************

% Grid parameters

%***********************************************************************

ie =40; % number of gridcells in X direction.

je =40; % number of gridcells in Y-direction.

ke =40; % number of gridcells in Z-direction.

nmax=65;

ddx =0.01;

dt =ddx/(2*3e8); %The center ofthe incident pulse

%Total/Scattered field boundaries.

%********************************

ia=7;

ja=7;

ka=7;

ib=ie-ia-1;

jb=je-ja-1;

kb=ke-ka-1;

npml = 9; %The number of PML cells

%***********************************************************************

% Initialize Field arrays

%***********************************************************************

ex(1:ie,1:je,1:ke) = 0.0;

ey(1:ie,1:je,1:ke) = 0.0;

ez(1:ie,1:je,1:ke) = 0.0;

dx(1:ie,1:je,1:ke) = 0.0;

dy(1:ie,1:je,1:ke) = 0.0;

dz(1:ie,1:je,1:ke) = 0.0;

hx(1:ie,1:je,1:ke) = 0.0;

hy(1:ie,1:je,1:ke) = 0.0;

hz(1:ie,1:je,1:ke) = 0.0;

gax(1:ie,1:je,1:ke) = 1.0;

gay(1:ie,1:je,1:ke) = 1.0;

gaz(1:ie,1:je,1:ke) = 1.0;

gbx(1:ie,1:je,1:ke) = 0.0;

gby(1:ie,1:je,1:ke) = 0.0;

gbz(1:ie,1:je,1:ke) = 0.0;

ix(1:ie,1:je,1:ke) = 0.0;

iy(1:ie,1:je,1:ke) = 0.0;

iz(1:ie,1:je,1:ke) = 0.0;

ez_inc(1:je+1) =0.0;

hx_inc(1:je+1) =0.0;

ez_low_m2 =0.0;

ez_high_m2 =0.0;

ez_low_m1 =0.0;

ez_high_m1 =0.0;

%************************************************************************

%Fourier analysis parameters

%************************************************************************

NFREQS=3;

real_in(1:NFREQS)=0.0;

imag_in(1:NFREQS)=0.0;

real_pt(1:NFREQS,1:ie,1:je)=0.0;

imag_pt(1:NFREQS,1:ie,1:je)=0.0;

amp(1:ie,1:je)=0.0;

phase(1:ie,1:je)=0.0;

amp_in(1:NFREQS)=0.0;

phase_in(1:NFREQS)=0.0;

freq(1)=10e6;

freq(2)=100e6;

freq(3)=433e6;

for n=1:NFREQS

arg(n)=2*pi*freq(n)*dt;

end

idxl(1:ia+1,1:je,1:ke)=0.0;idxh(1:ia+1,1:je,1:ke)=0.0; %review

ihxl(1:ia+1,1:je,1:ke)=0.0; ihxh(1:ia+1,1:je,1:ke)=0.0;

idyl(1:ie,1:ja+1,1:ke)=0.0;idyh(1:ie,1:ja+1,1:ke)=0.0;

ihyl(1:ie,1:ja+1,1:ke)=0.0;ihyh(1:ie,1:ja+1,1:ke)=0.0;

idzl(1:ie,1:je,1:ka+1)=0.0;idzh(1:ie,1:je,1:ka+1)=0.0;

ihzl(1:ie,1:je,1:ka+1)=0.0;ihzh(1:ie,1:je,1:ka+1)=0.0;

%***********************************************************************

% Specify the dipole Source

%***********************************************************************

ic =ie/2; jc=je/2; kc=ke/2; %Theposition of the source

to = 40;

spread = 10; %The width of the incident pulse

T = 0 ;

%*************************************************************************

% Specify the dielectric sphere

%*************************************************************************

radius=10;

epsilon=[1.0 30];

sigma =[0.0 0.3];

%Calculate the gax and gbx

%*************************

for i=ia:ib-1

forj=ja:jb-1

for k=ka:kb-1

eps=epsilon(1);

cond=sigma(1);

xdist=(ic-i-0.5);

ydist=(jc-j);

zdist=(kc-k);

dist=sqrt(power(xdist,2)+power(ydist,2)+power(zdist,2));

if dist<=radius

eps=epsilon(2);

cond=sigma(2);

end

gax(i,j,k)=1.0/(eps+(cond*dt/epsz));

gbx(i,j,k)=cond*dt/epsz;

end

%Calculate the gay and gby

%*************************

for i=ia:ib-1

forj=ja:jb-1

for k=ka:kb-1

eps=epsilon(1);

cond=sigma(1);

xdist=(ic-i);

ydist=(jc-j-0.5);

zdist=(kc-k);

dist=sqrt(power(xdist,2)+power(ydist,2)+power(zdist,2));

if dist<=radius

eps=epsilon(2);

cond=sigma(2);

end

gay(i,j,k)=1.0/(eps+(cond*dt/epsz));

gby(i,j,k)=cond*dt/epsz;

end

%Calculate the gaz and gbz

%*************************

for i=ia:ib-1

forj=ja:jb-1

for k=ka:kb-1

eps=epsilon(1);

cond=sigma(1);

xdist=(ic-i);

ydist=(jc-j);

zdist=(kc-k-0.5);

dist=sqrt(power(xdist,2)+power(ydist,2)+power(zdist,2));

if dist<=radius

eps=epsilon(2);

cond=sigma(2);

end

gaz(i,j,k)=1.0/(eps+(cond*dt/epsz));

gbz(i,j,k)=cond*dt/epsz;

end

%***********************************************************************

% Calculate the PML parameters

%***********************************************************************

for i=1:ie

gi1(i)=0.0;

fi1(i)=0.0;

gi2(i)=1.0;

fi2(i)=1.0;

gi3(i)=1.0;

fi3(i)=1.0;

end

for j=1:je

gj1(j)=0.0;

fj1(j)=0.0;

gj2(j)=1.0;

fj2(j)=1.0;

gj3(j)=1.0;

fj3(j)=1.0;

end

for k=1:ke

gk1(k)=0.0;

fk1(k)=0.0;

gk2(k)=1.0;

fk2(k)=1.0;

gk3(k)=1.0;

fk3(k)=1.0;

end

for i=1:npml

xxn=(npml-i)/npml;

xn=0.33*power(xxn,3);

fi1(i) =xn;

fi1(ie-1-i)=xn;

gi2(i) =1.0/(1.0+xn);

gi2(ie-1-i)=1.0/(1.0+xn);

gi3(i) =(1.0-xn)/(1.0+xn);

gi3(ie-1-i)=(1.0-xn)/(1.0+xn);

xxn=(npml-i-0.5)/npml;

xn=0.25*power(xxn,3);

gi1(i) =xn;

gi1(ie-2-i)=xn;

fi2(i) =1.0/(1.0+xn);

fi2(ie-2-i)=1.0/(1.0+xn);

fi3(i) =(1.0-xn)/(1.0+xn);

fi3(ie-2-i)=(1.0-xn)/(1.0+xn);

end

for j=1:npml

xnum=npml-j;

xd=npml;

xxn=xnum/xd;

xn=0.33*power(xxn,3);

fj1(j) =xn;

fj1(je-1-j)=xn;

gj2(j) =1.0/(1.0+xn);

gj2(je-1-j)=1.0/(1.0+xn);

gj3(j) =(1.0-xn)/(1.0+xn);

gj3(je-1-j)=(1.0-xn)/(1.0+xn);

xxn=(xnum-0.5)/xd;

xn=0.33*power(xxn,3);

gj1(j) =xn;

gj1(je-2-j)=xn;

fj2(j) =1.0/(1.0+xn);

fj2(je-2-j)=1.0/(1.0+xn);

fj3(j) =(1.0-xn)/(1.0+xn);

fj3(je-2-j)=(1.0-xn)/(1.0+xn);

end

for k=1:npml

xnum=npml-k;

xd=npml;

xxn=xnum/xd;

xn=0.33*power(xxn,3);

fk1(k) =xn;

fk1(ke-1-k)=xn;

gk2(k) =1.0/(1.0+xn);

gk2(ke-1-k)=1.0/(1.0+xn);

gk3(k) =(1.0-xn)/(1.0+xn);

gk3(ke-1-k)=(1.0-xn)/(1.0+xn);

xxn=(xnum-0.5)/xd;

xn=0.33*power(xxn,3);

gk1(k) =xn;

gk1(ke-2-k)=xn;

fk2(k) =1.0/(1.0+xn);

fk2(ke-2-k)=1.0/(1.0+xn);

fk3(k) =(1.0-xn)/(1.0+xn);

fk3(ke-2-k)=(1.0-xn)/(1.0+xn);

end

%***********************************************************************

% BEGIN TIME-STEPPING LOOP

%***********************************************************************

for n=1:nmax

T=T+1

%Calculate the incident buffer

%*****************************

forj=2:je

ez_inc(j)=ez_inc(j)+0.5*(hx_inc(j-1)-hx_inc(j));

end

%Fourier transform of the incident field

%****************************************

form=1:NFREQS

real_in(m)=real_in(m)+cos(arg(m)*T)*ez_inc(ja-1);

imag_in(m)=imag_in(m)-sin(arg(m)*T)*ez_inc(ja-1);

end

%Putsthe sinusoidal source in the middle

% pulse=sin(2*pi*400*1e6*dt*T);

pulse=exp(-0.5*(power((to-T)/spread,2)));

ez_inc(3)=pulse;

%ABCfor the incident buffer

%&&&&&&&&&&&&&&&&&&&&&&&&&&&&

ez_inc(1)=ez_low_m2;

ez_low_m2=ez_low_m1;

ez_low_m1=ez_inc(2);

ez_inc(je)=ez_high_m2;

ez_high_m2=ez_high_m1;

ez_high_m1=ez_inc(je-1);

%***********************************************************************

% Update Dx field

%***********************************************************************

for i=2:ia-1

forj=2:je

for k=2:ke

curl_h=hz(i,j,k)-hz(i,j-1,k)-hy(i,j,k)+hy(i,j,k-1);

idxl(i,j,k)=idxl(i,j,k)+curl_h;

dx(i,j,k)=gj3(j)*gk3(k)*dx(i,j,k)+gj2(j)*gk2(k)*0.5*(curl_h+gi1(i)*idxl(i,j,k));

end

for i=ia:ib

forj=2:je

for k=2:ke

curl_h=hz(i,j,k)-hz(i,j-1,k)-hy(i,j,k)+hy(i,j,k-1);

dx(i,j,k)=gj3(j)*gk3(k)*dx(i,j,k)+gj2(j)*gk2(k)*0.5*curl_h; %review

end

for i=ib+1:ie

ixh=i-ib; %review

forj=2:je

for k=2:ke

curl_h=hz(i,j,k)-hz(i,j-1,k)-hy(i,j,k)+hy(i,j,k-1);

idxh(ixh,j,k)=idxh(ixh,j,k)+curl_h;

dx(i,j,k)=gj3(j)*gk3(k)*dx(i,j,k)+gj2(j)*gk2(k)*0.5*(curl_h+gi1(i)*idxh(ixh,j,k));

end

%***********************************************************************

% Update Dy field

%***********************************************************************

for i=2:ie

forj=2:ja-1

for k=2:ke

curl_h=hx(i,j,k)-hx(i,j,k-1)-hz(i,j,k)+hz(i-1,j,k);

idyl(i,j,k)=idyl(i,j,k)+curl_h;

dy(i,j,k)=gi3(i)*gk3(k)*dy(i,j,k)+gi2(i)*gk2(k)*0.5*(curl_h+gj1(j)*idyl(i,j,k));

end

for i=2:ie

forj=ja:jb

for k=2:ke

curl_h=hx(i,j,k)-hx(i,j,k-1)-hz(i,j,k)+hz(i-1,j,k);

dy(i,j,k)=gi3(i)*gk3(k)*dy(i,j,k)+gi2(i)*gk2(k)*0.5*curl_h;

end

for i=2:ie

forj=jb+1:je

jyh=j-jb;

for k=2:ke

curl_h=hx(i,j,k)-hx(i,j,k-1)-hz(i,j,k)+hz(i-1,j,k);

idyh(i,jyh,k)=idyh(i,jyh,k)+curl_h;

dy(i,j,k)=gi3(i)*gk3(k)*dy(i,j,k)+gi2(i)*gk2(k)*0.5*(curl_h+gj1(j)*idyh(i,jyh,k));

end

%Incident Dy

%***********

for i=ia:ib

forj=ja:jb-1

dy(i,j,ka)=dy(i,j,ka)-0.5*hx_inc(j);

dy(i,j,kb+1)=dy(i,j,kb+1)+0.5*hx_inc(j);

end

%***********************************************************************

% Update DZ field

%***********************************************************************

for i=2:ie

forj=2:je

for k=1:ka-1

curl_h=hy(i,j,k)-hy(i-1,j,k)-hx(i,j,k)+hx(i,j-1,k);

idzl(i,j,k)=idzl(i,j,k)+curl_h;

dz(i,j,k)=gi3(i)*gj3(j)*dz(i,j,k)+gi2(i)*gj2(j)*0.5*(curl_h+gk1(k)*idzl(i,j,k));

end

for i=2:ie

forj=2:je

for k=ka:kb

curl_h=hy(i,j,k)-hy(i-1,j,k)-hx(i,j,k)+hx(i,j-1,k);

dz(i,j,k)=gi3(i)*gj3(j)*dz(i,j,k)+gi2(i)*gj2(j)*0.5*curl_h;

end

for i=2:ie

forj=2:je

for k=kb+1:ke

kzh=k-kb;

curl_h=hy(i,j,k)-hy(i-1,j,k)-hx(i,j,k)+hx(i,j-1,k);

idzh(i,j,kzh)=idzh(i,j,kzh)+curl_h;

dz(i,j,k)=gi3(i)*gj3(j)*dz(i,j,k)+gi2(i)*gj2(j)*0.5*(curl_h+gk1(k)*idzh(i,j,kzh));

end

%Incident DZ

%***********

for i=ia:ib

fork=ka:kb

dz(i,ja,k)=dz(i,ja,k)+0.5*hx_inc(ja-1);

dz(i,jb,k)=dz(i,jb,k)-0.5*hx_inc(jb);

end

% Thesource

%pulse=sin(2*pi*400*1e6*dt*T);

%fork=kc-6:kc+6

% dz(ic,jc,k)=0;

%end

%pulse=exp(-0.5*(power((to-T)/spread,2)));

%dz(ic,jc,kc)=pulse;

%***********************************************************************

% Caluculate E fields from D fields

%***********************************************************************

for i=2:ie-1

forj=2:je-1

for k=2:ke-1

ex(i,j,k)=gax(i,j,k)*(dx(i,j,k)-ix(i,j,k));

ix(i,j,k)=ix(i,j,k)+gbx(i,j,k)*ex(i,j,k);

ey(i,j,k)=gay(i,j,k)*(dy(i,j,k)-iy(i,j,k));

iy(i,j,k)=iy(i,j,k)+gby(i,j,k)*ey(i,j,k);

ez(i,j,k)=gaz(i,j,k)*(dz(i,j,k)-iz(i,j,k));

iz(i,j,k)=iz(i,j,k)+gbz(i,j,k)*ez(i,j,k);

end

%Calculate the fourier transform of Ex

%***************************************

for j=1:je

fori=1:ie

for m=1:NFREQS

real_pt(m,i,j)=real_pt(m,i,j)+cos(arg(m)*T)*ez(i,j,kc);

imag_pt(m,i,j)=imag_pt(m,i,j)+sin(arg(m)*T)*ez(i,j,kc);

end

% Calculate the incident field

%******************************

for j=1:je-1

hx_inc(j)=hx_inc(j)+0.5*(ez_inc(j)-ez_inc(j+1));

end

%***********************************************************************

% Update magnetic fields

%***********************************************************************

%1- Culculate Hx

%*****************

for i=1:ia-1

forj=1:je-1

for k=1:ke-1

curl_e=ey(i,j,k+1)-ey(i,j,k)-ez(i,j+1,k)+ez(i,j,k);

ihxl(i,j,k)=ihxl(i,j,k)+curl_e;

hx(i,j,k)=fj3(j)*fk3(k)*hx(i,j,k)+fj2(j)*fk2(k)*0.5*(curl_e+fi1(i)*ihxl(i,j,k));

end

fori=ia:ib

forj=1:je-1

for k=1:ke-1

curl_e=ey(i,j,k+1)-ey(i,j,k)-ez(i,j+1,k)+ez(i,j,k);

hx(i,j,k)=fj3(j)*fk3(k)*hx(i,j,k)+fj2(j)*fk2(k)*0.5*curl_e;

end

for i=ib+1:ie

ixh=i-ib;

forj=1:je-1

for k=1:ke-1

curl_e=ey(i,j,k+1)-ey(i,j,k)-ez(i,j+1,k)+ez(i,j,k);

ihxh(ixh,j,k)=ihxh(ixh,j,k)+curl_e;

hx(i,j,k)=fj3(j)*fk3(k)*hx(i,j,k)+fj2(j)*fk2(k)*0.5*(curl_e+fi1(i)*ihxh(ixh,j,k));

end

%Incident Hx value

%&&&&&&&&&&&&&&&&&&

for i=ia:ib

for k=ka:kb

hx(i,ja-1,k)= hx(i,ja-1,k)+0.5*ez_inc(ja);

hx(i,jb,k) = hx(i,jb,k)-0.5*ez_inc(jb);

end

%***********************************************************************

% Update Hy field

%***********************************************************************

for i=1:ie-1

forj=1:ja-1

for k=1:ke-1

curl_e=ez(i+1,j,k)-ez(i,j,k)-ex(i,j,k+1)+ex(i,j,k);

ihyl(i,j,k)=ihyl(i,j,k)+curl_e;

hy(i,j,k)=fi3(i)*fk3(k)*hy(i,j,k)+fi2(i)*fk3(k)*0.5*(curl_e+fj1(j)*ihyl(i,j,k));

end

for i=1:ie-1

forj=ja:jb

for k=1:ke-1

curl_e=ez(i+1,j,k)-ez(i,j,k)-ex(i,j,k+1)+ex(i,j,k);

hy(i,j,k)=fi3(i)*fk3(k)*hy(i,j,k)+fi2(i)*fk3(k)*0.5*curl_e;

end

for i=1:ie-1

forj=jb+1:je

jyh=j-jb;

for k=1:ke-1

curl_e=ez(i+1,j,k)-ez(i,j,k)-ex(i,j,k+1)+ex(i,j,k);

ihyh(i,jyh,k)=ihyh(i,jyh,k)+curl_e;

hy(i,j,k)=fi3(i)*fk3(k)*hy(i,j,k)+fi2(i)*fk3(k)*0.5*(curl_e+fj1(j)*ihyh(i,jyh,k));

end

%Incident Hy

%***********

for j=ja:jb

for k=ka:kb

hy(ia-1,j,k)=hy(ia-1,j,k)-0.5*ez_inc(j);

hy(ib,j,k)=hy(ib,j,k)+0.5*ez_inc(j);

end

%***********************************************************************

% Update HZ field

%***********************************************************************

for i=1:ie-1

forj=1:je-1

for k=1:ka-1

curl_e=ex(i,j+1,k)-ex(i,j,k)-ey(i+1,j,k)+ey(i,j,k);

ihzl(i,j,k)=ihzl(i,j,k)+curl_e;

hz(i,j,k)=fi3(i)*fj3(j)*hz(i,j,k)+fi2(i)*fj2(j)*0.5*(curl_e+fk1(k)*ihzl(i,j,k));

end

for i=1:ie-1

forj=1:je-1

for k=ka:kb

curl_e=ex(i,j+1,k)-ex(i,j,k)-ey(i+1,j,k)+ey(i,j,k);

hz(i,j,k)=fi3(i)*fj3(j)*hz(i,j,k)+fi2(i)*fj2(j)*0.5*curl_e;

end

for i=1:ie-1

forj=1:je-1

for k=kb+1:ke

kzh=k-kb;

curl_e=ex(i,j+1,k)-ex(i,j,k)-ey(i+1,j,k)+ey(i,j,k);

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

发帖回复