periodic 3d tongzhou resonanter.m

FDTD方法计算三维结构的程序

%COXIAL
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 
ii=sqrt(-1)

%*********************************************************************** 
%     Grid parameters 
%*********************************************************************** 

ie=34;      %number of grid cells in x-direction  r=17mm  rout=34mm
je=20;
ke=15;       %number of grid cells in z-direction  l=7.5mm
  
ib=ie+1;
jb=je+1;
kb=ke+1; 

dtheta=2*pi/je;
dx=0.0005;          %space increment of cubic lattice 0.5mm
dt=dx/(2.0*cc);     %time step 
nmax=2000;          %total number of time steps 
cb=dt/epsz;
db=dt/muz/dx;

er=zeros(ie,jb,kb);
et=zeros(ib,jb,kb);
ez=zeros(ib,jb,ke);
hr=zeros(ib,jb,ke);
ht=zeros(ie,jb,ke);
hz=zeros(ie,jb,kb);
r=zeros(ib,jb,kb);
ez0=zeros(ib,jb);
Nf=1000
fmax=100e+9
df=fmax/Nf

f0=30e+9


ndecay=1/(2*fmax*dt);
n0=1/(2*f0*dt);
t=3.123/dx/ie/2;
%compute the r(i)
for i=1:ib
    r(i,:,:)=dx*ie+dx*(i-1);
end
    for i=1:ib
        ez0(i,:)=besselj(0,t*(ie+i-1)*dx)*bessely(0,t*dx*2*ie)-bessely(0,t*(ie+i-1)*dx)*besselj(0,t*2*dx*ie); 
    end
   
    %*********************************************************************** 
%     BEGIN TIME-STEPPING LOOP 
%*********************************************************************** 
  
for n=1:nmax 
  
%*********************************************************************** 
%     Update electric fields 
%*********************************************************************** 
er(1:ie,2:jb,2:ke)=er(1:ie,2:jb,2:ke)+...
                    dt/dtheta/epsz/(r(1:ie,2:jb,2:ke)+dx/2).*(hz(1:ie,2:jb,2:ke)-hz(1:ie,1:je,2:ke))-...
                    dt/epsz/dx*(ht(1:ie,2:jb,2:ke)-ht(1:ie,2:jb,1:ke-1)); 

ez(2:ie,2:jb,1:ke)=ez(2:ie,2:jb,1:ke)+... 
                   dt/epsz/dx*(ht(2:ie,2:jb,1:ke)-ht(1:ie-1,2:jb,1:ke))+dt/epsz/2*(ht(2:ie,2:jb,1:ke)+ht(1:ie-1,2:jb,1:ke))./r(2:ie,2:jb,1:ke)-... 
                   dt/dtheta/epsz*(hr(2:ie,2:jb,1:ke)-hr(2:ie,1:je,1:ke))./r(2:ie,2:jb,1:ke);

et(2:ie,1:jb,2:ke)=et(2:ie,1:jb,2:ke)+dt/epsz/dx*(hr(2:ie,1:jb,2:ke)-hr(2:ie,1:jb,1:ke-1)-hz(2:ie,1:jb,2:ke)+hz(1:ie-1,1:jb,2:ke));
                   

   %add excite**************************depend upon the mode which compute.
         ez(5,5,3)=exp(-(((n-n0)/ndecay)*(n-n0)/ndecay))+ez(5,5,3);
        %ez(:,:,3)=ez0(:,:)*exp(-(((n-n0)/ndecay)*(n-n0)/ndecay))+ez(:,:,3);
   %************* periodic  boundary condition   for theta direction *****************************

 er(1:ie,1,2:ke)=er(1:ie,1,2:ke)+...
                    dt/dtheta/epsz/(r(1:ie,1,2:ke)+dx/2).*(hz(1:ie,1,2:ke)-hz(1:ie,je,2:ke))-...
                   dt/epsz/dx*(ht(1:ie,1,2:ke)-ht(1:ie,1,1:ke-1)); 

 ez(2:ie,1,1:ke)=ez(2:ie,1,1:ke)+... 
                   dt/epsz/dx*(ht(2:ie,1,1:ke)-ht(1:ie-1,1,1:ke))+dt/epsz/2*(ht(2:ie,1,1:ke)+ht(1:ie-1,1,1:ke))./r(2:ie,1,1:ke)-... 
                   dt/dtheta/epsz*(hr(2:ie,1,1:ke)-hr(2:ie,je,1:ke))./r(2:ie,1,1:ke);
 

%*********************************************************************** 
%     Update magnetic fields 
%*********************************************************************** 
hr(2:ie,1:je,1:ke)=hr(2:ie,1:je,1:ke)-...
                    dt/muz/dtheta*(ez(2:ie,2:jb,1:ke)-ez(2:ie,1:je,1:ke))./r(2:ie,1:je,1:ke)+dt/muz/dx*(et(2:ie,1:je,2:kb)-et(2:ie,1:je,1:ke));
  
ht(1:ie,1:jb,1:ke)=ht(1:ie,1:jb,1:ke)+... 
                   dt/muz/dx*(er(1:ie,1:jb,1:ke)-er(1:ie,1:jb,2:kb)+ez(2:ib,1:jb,1:ke)-ez(1:ie,1:jb,1:ke));

hz(1:ie,1:je,2:ke)=hz(1:ie,1:je,2:ke)-dt/muz/2/(r(1:ie,1:je,2:ke)+dx/2).*(et(1:ie,1:je,2:ke)+et(2:ib,1:je,2:ke))-...
                     dt/muz/dx*(et(2:ib,1:je,2:ke)-et(1:ie,1:je,2:ke))+dt/muz/dtheta/(r(1:ie,1:je,2:ke)+dx/2).*(er(1:ie,2:jb,2:ke)-er(1:ie,1:je,2:ke));  
                     
%************* periodic  boundary condition   for theta direction *****************************
hr(2:ie,jb,1:ke)=hr(2:ie,jb,1:ke)-...
                    dt/muz/dtheta*(ez(2:ie,2,1:ke)-ez(2:ie,jb,1:ke))./r(2:ie,jb,1:ke)+dt/muz/dx*(et(2:ie,jb,2:kb)-et(2:ie,jb,1:ke));
                
hz(1:ie,jb,2:ke)=hz(1:ie,jb,2:ke)-dt/muz/2/(r(1:ie,jb,2:ke)+dx/2).*(et(1:ie,jb,2:ke)+et(2:ib,jb,2:ke))-...
                     dt/muz/dx*(et(2:ib,jb,2:ke)-et(1:ie,jb,2:ke))+dt/muz/dtheta/(r(1:ie,jb,2:ke)+dx/2).*(er(1:ie,2,2:ke)-er(1:ie,jb,2:ke)); 
               
%*********************************************************************** 
%% restroe the imformation at the picked point
%*********************************************************************** 
gVer1(n)=ez(10,9,10);

 end
 
ff=abs(gVer1);
%plot(ff)