fdtd1[1].m

fdtd一维到三维算法

function [y1,y1max,y2,y2max,x,dt]=fdtd1(L,xnum,tnum,tfactor)
%
% Function file for the example problem of a 1-dimensional
%  plane wave propagating in the x direction.  This file calculates
%  the finite difference solution for the electric and magnetic fields
%  Ez and Hy using the FDTD.  

%
%  The input to the file includes "L", "xnum", "tnum, "tfactor"
%
%     It is assumed that our solution space is x=0 to x=L
%
%    "xnum" is the number of points at which the solution is to be
%     calculated.  Note that dz=L/(xnum-1)
%
%    "tnum" is the number of time iterations desired.
%
%    "tfactor" is used to pick the propagation speed of simulated wave.
%        dt=tfactor*(dz/c)
%
%
%  There are six outputs:  
%     (1) matrix of values of the electric field Ez, which
%         is stored in y1 (which is tnum x znum in size)
%     (2) the maximum value of all the Ez values - stored as y1max.
%     (3) matrix of values of the magnetic field Hy, which
%         is stored in y2 (which is tnum x znum in size)
%     (4) the maximum value of all the Hy values - stored as y2max.
%     (5) the vector of spatial points, x (which is 1 x znum).
%     (6) the time increment dt (1x1)
%
%     This version assumes:
%
%         (1) a delta function source of magnitude eta0 which is located
%             in the middle of the solution space.  It is turned "on" for
%             one time step, and then turned off.
%
%         (2) perfectly conducting walls (i.e., Etan=0) at x=0 and x=L.
%        
%         (3) When tangential electric fields are specified on the boundary,
%             the value of Ez also needs to be calculated at one point outside
%             of the solution space.  In particular, at x=L+dx.
%
c=3e8;
dx=L/(xnum-1);
dt=tfactor*(dx/c);
x=linspace(0,L,xnum);
Ez=zeros(tnum,xnum+1);
Hy=zeros(tnum,xnum);


%  Note that in order to calculate the field components at future
%  time steps, the electric field also need to be specifid at one point "outside'
%  the solution space. This "extra point" will be taken out later for
%  the matrix passed back.

mu0=4e-7*pi;
eps0=8.854e-12;
eta0=sqrt(mu0/eps0);
isrce=round(L/2);

for n=1:tnum

	if (n==1) 
 
 %   n=1 and n-2 establish the initial source excitation,
 %    which is assumed to be a point source of magnitude eta0
 %    in the middle of the solution space.  It is on for one time step.
 
		Ez(1,isrce)=eta0;
       
        for i=1:xnum
            Hy(1,i)=0.0+1./eta0.*(Ez(1,i+1)-Ez(1,i));
 
 %             assuming that Hy=0 for t<0.
 
        end
                 
    elseif (n==2)
        	    Ez(2,isrce)=0;   % turns off the source.
         
                for i=2:(isrce-1)
                    
  %              start at i=2, since Ez(1,i)=0 for conductor
                    
                   Ez(2,i)=Ez(1,i)+eta0.*(Hy(1,i)-Hy(1,i-1));
               end
               
            for i=(isrce+1):(xnum-1)
                Ez(2,i)=Ez(1,i)+eta0.*(Hy(1,i)-Hy(1,i-1));
            end
            
                Ez(2,xnum)=0.;    % Etan=0 on conductor;
                Ez(2,xnum+1)=0.;  % Again, Etan=0 in conductor.  Used to
   %                              calculate Hy.             
            for i=1:xnum
               Hy(2,i)=Hy(1,i)+1./eta0.*(Ez(2,i+1)-Ez(2,i));
           end
              
    else  
              
		for i=2:xnum-1
                Ez(n,i)=Ez(n-1,i)+eta0.*(Hy(n-1,i)-Hy(n-1,i-1));
        end
             Ez(n,xnum)=0.;  % These are already a 0 in the matrix, but
             Ez(n,xnum+1)=0.;  % put here as a reminder.
             
        for i=1:xnum
               Hy(n,i)=Hy(n-1,i)+1./eta0.*(Ez(n,i+1)-Ez(n,i));
        end

	end

	
end

% Since we added an point outside the solution space for Ez, we
%   need to take that point out before passing back the answer.

for i=1:xnum
    y1(:,i)=Ez(:,i);
end

y1max=max(max(Ez));
y2=Hy;
y2max=max(max(Hy));