一维).txt

FDTD一维程序代码

        case 1
            surf(i,j,Ez0);
            axis([0 Nx 0 Ny -0.03 0.03]);
            set(gca, 'XTick',[1 Nx/4 Nx/2 3*Nx/4 Nx],'FontSize',8);
            set(gca, 'XTickLabel',[0 X/4 X/2 3*X/4 X],'FontSize',8);
            xlabel('x in m');  
            set(gca, 'YTick',[1 Ny/4 Ny/2 3*Ny/4 Ny],'FontSize',8);
            set(gca, 'YTickLabel',[0 Y/4 Y/2 3*Y/4 Y],'FontSize',8);
            ylabel('y in m');
            zlabel('Amplitude of Ez');
        case 2
            surf(i,j,Hx0);
            axis([0 Nx 0 Ny -1e-4 1e-4]);
            set(gca, 'XTick',[1 Nx/4 Nx/2 3*Nx/4 Nx],'FontSize',8);
            set(gca, 'XTickLabel',[0 X/4 X/2 3*X/4 X],'FontSize',8);
            xlabel('x in m');  
            set(gca, 'YTick',[1 Ny/4 Ny/2 3*Ny/4 Ny],'FontSize',8);
            set(gca, 'YTickLabel',[0 Y/4 Y/2 3*Y/4 Y],'FontSize',8);
            ylabel('y in m');
            zlabel('Amplitude of Hx');
        case 3
            surf(i,j,Hy0);
            axis([0 Nx 0 Ny -1e-4 1e-4]);
            set(gca, 'XTick',[1 Nx/4 Nx/2 3*Nx/4 Nx],'FontSize',8);
            set(gca, 'XTickLabel',[0 X/4 X/2 3*X/4 X],'FontSize',8);
            xlabel('x in m');  
            set(gca, 'YTick',[1 Ny/4 Ny/2 3*Ny/4 Ny],'FontSize',8);
            set(gca, 'YTickLabel',[0 Y/4 Y/2 3*Y/4 Y],'FontSize',8);
            ylabel('y in m');
            zlabel('Amplitude of Hy');
        otherwise 
            subplot(2,2,1);
            surf(i,j,Ez0);
            title('Ex');
            axis([0 Nx 0 Ny -0.03 0.03]);
            set(gca, 'XTick',[1 Nx/4 Nx/2 3*Nx/4 Nx],'FontSize',8);
            set(gca, 'XTickLabel',[0 X/4 X/2 3*X/4 X],'FontSize',8);
            xlabel('x in m');  
            set(gca, 'YTick',[1 Ny/4 Ny/2 3*Ny/4 Ny],'FontSize',8);
            set(gca, 'YTickLabel',[0 Y/4 Y/2 3*Y/4 Y],'FontSize',8);
            ylabel('y in m');
            zlabel('Amplitude of Ez');
            
            subplot(2,2,2);
            surf(i,j,Hx0);
            title('Hx');
            axis([0 Nx 0 Ny -1e-4 1e-4]);
            set(gca, 'XTick',[1 Nx/4 Nx/2 3*Nx/4 Nx],'FontSize',8);
            set(gca, 'XTickLabel',[0 X/4 X/2 3*X/4 X],'FontSize',8);
            xlabel('x in m');  
            set(gca, 'YTick',[1 Ny/4 Ny/2 3*Ny/4 Ny],'FontSize',8);
            set(gca, 'YTickLabel',[0 Y/4 Y/2 3*Y/4 Y],'FontSize',8);
            ylabel('y in m');
            zlabel('Amplitude of Hx');
            
            subplot(2,2,3);
            surf(i,j,Hy0);
            title('Hy');
            axis([0 Nx 0 Ny -1e-4 1e-4]);
            set(gca, 'XTick',[1 Nx/4 Nx/2 3*Nx/4 Nx],'FontSize',8);
            set(gca, 'XTickLabel',[0 X/4 X/2 3*X/4 X],'FontSize',8);
            xlabel('x in m');  
            set(gca, 'YTick',[1 Ny/4 Ny/2 3*Ny/4 Ny],'FontSize',8);
            set(gca, 'YTickLabel',[0 Y/4 Y/2 3*Y/4 Y],'FontSize',8);
            ylabel('y in m');
            zlabel('Amplitude of Hy');
    end

    %surf(i,j,Ez0);
    %axis([0 Nx 0 Ny -0.03 0.03])
    %A = rot90(Ez0);
    %imagesc(A);
    pause(0.005);
end


3D FDTD for hexahedral cavity with conducting walls

% Physical constants
eps0 = 8.8541878e-12;        % Permittivity of vacuum
mu0  = 4e-7 * pi;            % Permeability of vacuum
c0  = 299792458;            % Speed of light in vacuum

% Parameter initiation
Lx = .05; Ly = .04; Lz = .03; % Cavity dimensions in meters
Nx =  25; Ny =  20; Nz =  15; % Number of cells in each direction
Cx = Nx / Lx;                % Inverse cell dimensions
Cy = Ny / Ly;
Cz = Nz / Lz;
Nt = 1024;                    % Number of time steps
Dt = 1/(c0*norm([Cx Cy Cz])); % Time step

% Allocate field matrices
Ex = zeros(Nx  , Ny+1, Nz+1);
Ey = zeros(Nx+1, Ny  , Nz+1);
Ez = zeros(Nx+1, Ny+1, Nz  );
Hx = zeros(Nx+1, Ny  , Nz  );
Hy = zeros(Nx  , Ny+1, Nz  );
Hz = zeros(Nx  , Ny  , Nz+1);

% Allocate time signals
Et = zeros(Nt,3);

% Initiate fields (near but not on the boundary)
Ex(1,2,2) = 1;
Ey(2,1,2) = 2;
Ez(2,2,1) = 3;

% Time stepping
for n = 1:Nt;

  % Update H everywhere
  Hx = Hx + (Dt/mu0)*(diff(Ey,1,3)*Cz - diff(Ez,1,2)*Cy);
  Hy = Hy + (Dt/mu0)*(diff(Ez,1,1)*Cx - diff(Ex,1,3)*Cz);
  Hz = Hz + (Dt/mu0)*(diff(Ex,1,2)*Cy - diff(Ey,1,1)*Cx);

  % Update E everywhere except on boundary
  Ex(:,2:Ny,2:Nz) = Ex(:,2:Ny,2:Nz) + (Dt /eps0) * ...
      (diff(Hz(:,:,2:Nz),1,2)*Cy - diff(Hy(:,2:Ny,:),1,3)*Cz);
  Ey(2:Nx,:,2:Nz) = Ey(2:Nx,:,2:Nz) + (Dt /eps0) * ...
      (diff(Hx(2:Nx,:,:),1,3)*Cz - diff(Hz(:,:,2:Nz),1,1)*Cx);
  Ez(2:Nx,2:Ny,:) = Ez(2:Nx,2:Ny,:) + (Dt /eps0) * ...
      (diff(Hy(:,2:Ny,:),1,1)*Cx - diff(Hx(2:Nx,:,:),1,2)*Cy);

  % Sample the electric field at chosen points
  Eto(n,:) = [Ex(4,4,4) Ey(4,4,4) Ez(4,4,4)];
end



Version with for-loops for update of Ex and Hx.

% Physical constants
eps0 = 8.8541878e-12;        % Permittivity of vacuum
mu0  = 4e-7 * pi;            % Permeability of vacuum
c0  = 299792458;            % Speed of light in vacuum

% Parameter initiation
Lx = .05; Ly = .04; Lz = .03; % Cavity dimensions in meters
Nx =  25; Ny =  20; Nz =  15; % Number of cells in each direction
Cx = Nx / Lx;                % Inverse cell dimensions
Cy = Ny / Ly;
Cz = Nz / Lz;
Nt = 32;                    % Number of time steps
Dt = 1/(c0*norm([Cx Cy Cz])); % Time step

% Allocate field matrices
Ex = zeros(Nx  , Ny+1, Nz+1);
Ey = zeros(Nx+1, Ny  , Nz+1);
Ez = zeros(Nx+1, Ny+1, Nz  );
Hx = zeros(Nx+1, Ny  , Nz  );
Hy = zeros(Nx  , Ny+1, Nz  );
Hz = zeros(Nx  , Ny  , Nz+1);

% Allocate time signals
Et = zeros(Nt,3);

% Initiate fields (near but not on the boundary)
Ex(1,2,2) = 1;
Ey(2,1,2) = 2;
Ez(2,2,1) = 3;

% Time stepping
for n = 1:Nt;
  % Update H everywhere

  for i = 1:Nx+1
    for j = 1:Ny
        for k = 1:Nz
          Hx(i,j,k) = Hx(i,j,k) + (Dt/mu0)* ...
          ((Ey(i,j,k+1)-Ey(i,j,k))*Cz - (Ez(i,j+1,k)-Ez(i,j,k))*Cy);
        end
    end
  end
  
  Hy = Hy + (Dt/mu0)*((Ez(2:Nx+1,:,:)-Ez(1:Nx,:,:))*Cx ...
                    - (Ex(:,:,2:Nz+1)-Ex(:,:,1:Nz))*Cz);
  Hz = Hz + (Dt/mu0)*((Ex(:,2:Ny+1,:)-Ex(:,1:Ny,:))*Cy ...
                    - (Ey(2:Nx+1,:,:)-Ey(1:Nx,:,:))*Cx);

% Update E everywhere except on boundary

  for i = 1:Nx
    for j = 2:Ny
        for k = 2:Nz
          Ex(i,j,k) = Ex(i,j,k) + (Dt /eps0) * ...
              ((Hz(i,j,k)-Hz(i,j-1,k))*Cy-(Hy(i,j,k)-Hy(i,j,k-1))*Cz);
        end
    end
  end

  Ey(2:Nx,:,2:Nz) = Ey(2:Nx,:,2:Nz) + (Dt /eps0) * ...
      ((Hx(2:Nx,:,2:Nz)-Hx(2:Nx,:,1:Nz-1))*Cz ...
    - (Hz(2:Nx,:,2:Nz)-Hz(1:Nx-1,:,2:Nz))*Cx);
  Ez(2:Nx,2:Ny,:) = Ez(2:Nx,2:Ny,:) + (Dt /eps0) * ...
      ((Hy(2:Nx,2:Ny,:)-Hy(1:Nx-1,2:Ny,:))*Cx ...
    - (Hx(2:Nx,2:Ny,:)-Hx(2:Nx,1:Ny-1,:))*Cy);

  % Sample the electric field at chosen points
  Etss(n,:) = [Ex(4,4,4) Ey(4,4,4) Ez(4,4,4)];
end



Version with diff written out as plain Matlab. 

% Physical constants
eps0 = 8.8541878e-12;        % Permittivity of vacuum
mu0  = 4e-7 * pi;            % Permeability of vacuum
c0  = 299792458;            % Speed of light in vacuum

% Parameter initiation
Lx = .05; Ly = .04; Lz = .03; % Cavity dimensions in meters
Nx =  25; Ny =  20; Nz =  15; % Number of cells in each direction
Cx = Nx / Lx;                % Inverse cell dimensions
Cy = Ny / Ly;
Cz = Nz / Lz;
Nt = 1024;                    % Number of time steps
Dt = 1/(c0*norm([Cx Cy Cz])); % Time step

% Allocate field matrices
Ex = zeros(Nx  , Ny+1, Nz+1);
Ey = zeros(Nx+1, Ny  , Nz+1);
Ez = zeros(Nx+1, Ny+1, Nz  );
Hx = zeros(Nx+1, Ny  , Nz  );
Hy = zeros(Nx  , Ny+1, Nz  );
Hz = zeros(Nx  , Ny  , Nz+1);

% Allocate time signals
Et = zeros(Nt,3);

% Initiate fields (near but not on the boundary)
Ex(1,2,2) = 1;
Ey(2,1,2) = 2;
Ez(2,2,1) = 3;

% Time stepping
for n = 1:Nt;
  % Update H everywhere
  Hx = Hx + (Dt/mu0)*((Ey(:,:,2:Nz+1)-Ey(:,:,1:Nz))*Cz ...
                    - (Ez(:,2:Ny+1,:)-Ez(:,1:Ny,:))*Cy);
  Hy = Hy + (Dt/mu0)*((Ez(2:Nx+1,:,:)-Ez(1:Nx,:,:))*Cx ...
                    - (Ex(:,:,2:Nz+1)-Ex(:,:,1:Nz))*Cz);
  Hz = Hz + (Dt/mu0)*((Ex(:,2:Ny+1,:)-Ex(:,1:Ny,:))*Cy ...
                    - (Ey(2:Nx+1,:,:)-Ey(1:Nx,:,:))*Cx);

  % Update E everywhere except on boundary
  Ex(:,2:Ny,2:Nz) = Ex(:,2:Ny,2:Nz) + (Dt /eps0) * ...
      ((Hz(:,2:Ny,2:Nz)-Hz(:,1:Ny-1,2:Nz))*Cy ...
    - (Hy(:,2:Ny,2:Nz)-Hy(:,2:Ny,1:Nz-1))*Cz);
  Ey(2:Nx,:,2:Nz) = Ey(2:Nx,:,2:Nz) + (Dt /eps0) * ...
      ((Hx(2:Nx,:,2:Nz)-Hx(2:Nx,:,1:Nz-1))*Cz ...
    - (Hz(2:Nx,:,2:Nz)-Hz(1:Nx-1,:,2:Nz))*Cx);
  Ez(2:Nx,2:Ny,:) = Ez(2:Nx,2:Ny,:) + (Dt /eps0) * ...
      ((Hy(2:Nx,2:Ny,:)-Hy(1:Nx-1,2:Ny,:))*Cx ...
    - (Hx(2:Nx,2:Ny,:)-Hx(2:Nx,1:Ny-1,:))*Cy);

  % Sample the electric field at chosen points
  Ets(n,:) = [Ex(4,4,4) Ey(4,4,4) Ez(4,4,4)];
end