fdtd_1d_erand[1].m

fdtd一维到三维算法

%***********************************************************************
%     1-D FDTD code with E-field boundary condition
%***********************************************************************
%
%     Program author: Mats Gustafsson
%                     Department of Electroscience
%                     Lund Institute of Technology
%                     Sweden
%                     
%     Date of this version:  October 2002
%
%
%   Ein
%   --->                Eut                         PEC
%     ----------------------------------------------->
%    z=0               z=Lut                        z=Lz
%    r=1                                            r=Nz+1
%
%***********************************************************************
% Physical constants
mu0  = 4e-7 * pi;             % Permeability of vacuum
c0   = 299792458;             % Speed of light in vacuum
eps0 = 1/c0^2/mu0;            % Permittivity of vacuum
eta0 = sqrt(mu0/eps0);        % Wave impedance in vacuum  
GHz = 10^9;
mm = 10^(-3);

%***********************************************************************
% Parameter initiation
%***********************************************************************
Lz = 3;                     % Length in meters
Tmax = Lz/c0;                 % Time interval [0,Tmax]
Dz = 20*mm;                    % spatial grid size


R = 0.9;                  % stabillity number (grid cells per time step)
freq = 1*GHz;                  % center frequency of source excitation
fmax = 4*GHz;                  % Max frequency to plot

Lut = Lz/2;                        % Output point
n_pict = 20;                    % plot each n_pict step 
Nfft = 4*1024;                    % Number of points for the fft

Nz = round(Lz/Dz);                     % Number of cells in the z-direction
Nut = round(Lut/Dz);           % Coordinate for output
Dt = Dz*R/c0;                  % Time step
Nt = round(Tmax/Dt);           % Number of time steps
z = linspace(0,Lz,Nz+1);         % z-coordinates
t = Dt*[0:Nt-1]';               % time
lambda = c0/freq;               % center wavelength of source excitation
ppw = lambda/Dz                 % sample points per wave length at the center frequency
omega = 2.0*pi*freq;            % angular frequency


%***********************************************************************
% Allocate field vectors
%***********************************************************************
Ex = zeros(Nz+1  ,1);          % Electric field
Hy = zeros(Nz,1);          % Magnetic field

Eut = zeros(Nt,1);           % Output

res = zeros(Nz+1,2*Nt);
%***********************************************************************
%     Wave excitation, Incident puls from the left
%***********************************************************************
% Broad band puls
tau = 250.0e-12;
delay = 2.5*tau;
Ein = sin(omega*(t-delay)).*exp(-((t-delay).^2/tau^2));

% Incident time harmonic wave from the left
%Ein = sin(omega*(t));

Ehin = fft(Ein,Nfft)/Nt*2;
Df = 1/(Dt*Nfft);
freqN = round(fmax/Df);
f = linspace(0,Df*freqN,freqN);

% analytic solution
propdelay = (Nut-1)*Dz/c0;
Eut_ana = sin(omega*(t-delay-propdelay)).*exp(-((t-delay-propdelay).^2/tau^2)).*(sign(t-propdelay)+1)/2;
%Eut_ana = sin(omega*(t-propdelay)).*(sign(t-propdelay)+1)/2;
Ehut_ana = fft(Eut_ana,Nfft)/Nt*2;

%***********************************************************************
% Time stepping
%***********************************************************************
for n = 1:Nt;
 
  % Update H from n*Dt-Dt/2 to n*Dt+Dt/2
  Hy = Hy - (Dt/mu0/Dz) * diff(Ex,1);                     % main grid
  
  % Update E from n*Dt to (n+1)*Dt
  Ex(2:Nz) = Ex(2:Nz) - (Dt/eps0/Dz) * diff(Hy,1);    % main grid except boundary
  
  Ex(1) = Ein(n);          % E-field on the boundary  
  
  % Sample the electric field at chosen points
  if mod(n,n_pict) == 0
      figure(1);
      subplot(2,1,1);
      plot(z,Ex(1:Nz+1),'LineWidth',1);
      title(['time is ',num2str(n*Dt*10^9),'ns'])
      axis tight;
      ylabel('E-field');
      
      subplot(2,1,2);
      plot(z(1:Nz),Hy(1:Nz),'LineWidth',1);
      title(['time is ',num2str(n*Dt*10^9),'ns'])
      axis tight;
      ylabel('H-field');
      pause;%(0.05);
  end
  Eut(n) = Ex(Nut);
end

%***********************************************************************
%     Post processing
%***********************************************************************

Eerror = norm(Eut-Eut_ana)*Dt   % L2 error

figure(2);
clf;
subplot(3,1,1);
plot(t,Eut,t,Eut_ana,'r','LineWidth',1)
axis([0.9*propdelay 1.9*propdelay -1 1])
xlabel('time');
ylabel('E-field');

subplot(3,1,2);
warning off;
A=plotyy(f,abs(Ehin(1:freqN)),f,c0./f/Dz); 
warning on;
set(A(2),'yLim',[0 30],'yTick',[0:5:30],'xlim',[0 Df*freqN],'xGrid','on','yGrid','on');
set(A(1),'xlim',[0 Df*freqN]);
set(get(A(2),'Ylabel'),'String','points/wavelength');
set(get(A(1),'Ylabel'),'String','input');
set(get(A(1),'Xlabel'),'String','frequency');

subplot(3,1,3);
Ehut = fft(Eut,Nfft)/Nt*2;
%plot(f,20*log10(abs(abs(Ehut(1:freqN))-abs(Ehut_ana(1:freqN)))./abs(Ehut_ana(1:freqN))));
Eherror = abs(Ehut(1:freqN)-Ehut_ana(1:freqN))./abs(Ehut_ana(1:freqN));
Eherror1 = interp1(f,Eherror,freq)
plot(f,20*log10(Eherror),'LineWidth',1);
axis tight
ax=axis;
axis([ax(1:2) -120 20]);
xlabel('frequency');
ylabel('Error in dB');