fdtd_1d_abc[1].m

fdtd一维到三维算法

%***********************************************************************
%     1-D FDTD code with PML
%***********************************************************************
%
%     Program author: Mats Gustafsson
%                     Department of Electroscience
%                     Lund Institute of Technology
%                     Sweden
%                     
%     Date of this version:  October 2002
%
%
%             Eref   Ein                  Erans
%   PEC  PML                objects                 PML  PEC
%    +-----+---+------+---------------------+-----+------+-->
%         z=0                                    z=Lz
%         r=1 Nref   Nin                 Ntrans  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 = 1;                     % Length in meters
Tmax = 4*Lz/c0;               % Time interval [0,Tmax]
Dz = 1*mm;                    % spatial grid size

R = 1/sqrt(3);                  % stabillity number, grid cells per time step
freq = 5*GHz;                  % center frequency of source excitation
fmax = 9*GHz;                  % Max frequency to plot
fmin = 1*GHz;

Nz = round(Lz/Dz);                     % Number of cells in the z-direction
n_pict = round(Nz/20);                    % plot each n_pict step 
linew = 2;                      % linewidth in plot

Dt = Dz*R/c0;                  % Time step
Nt = round(Tmax/Dt);           % Number of time steps
Nfft = max(2^(round(log2(Nt)+1)),2*1024);   % Number of points for the fft
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


Nabc = 1*10;                    % Number of PML layers to the left and to the right
Nin = 4;                      % scattered field surface
Nref = 2;                      % near to far field surface
Ntrans = Nz-2;
Labc = Nabc*Dz;               % PML length

%***********************************************************************
% Material parameters
%***********************************************************************
Nmedia=3;                   % Number of material
media_par = [1.0 0.0;       % vacuum: permitivity and conductivity 
             3.4 1*10^(-12);  %       
             3.0 0];   

objects = 0;               % Number of objects
obj_z = Dz/2+[0.500 0.504;        % object 1 between z1 and z2
         0.624 0.628;
         0.24 0.26];
obj_m = [2 2 1];           % Material of objects
obj_c = 'rrb';             % Color of object

epsr = ones(Nz-1,1);
sig = zeros(Nz-1,1);
obj_ind = [];
obj_ind(:,1) = ceil(obj_z(:,1)/Dz);
obj_ind(:,2) = floor(obj_z(:,2)/Dz);
obj_frac(:,1) = obj_ind(:,1)-obj_z(:,1)/Dz;
obj_frac(:,2) = -obj_ind(:,2)+obj_z(:,2)/Dz;
%break

%epsr(64:73) = 3;
for n=1:objects
    epsr(obj_ind(n,1):obj_ind(n,2)) = media_par(obj_m(n),1);
    %epsr(obj_ind(n,1)) = obj_frac(n,1) * (media_par(obj_m(n),1)-1)+1;
    %epsr(obj_ind(n,2)+1) = obj_frac(n,2) * (media_par(obj_m(n),1)-1)+1;
    sig(obj_ind(n,1):obj_ind(n,2)) = media_par(obj_m(n),2);
    %sig(obj_ind(n,1)) = obj_frac(n,1) * media_par(obj_m(n),2);
    %sig(obj_ind(n,2)+1) = obj_frac(n,2) * media_par(obj_m(n),2);
end


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

Eref = zeros(Nt,1);
Etrans = zeros(Nt,1);

% Absorbing boundary conditions
E1abc = zeros(Nabc+1,1);        % left
E2abc = zeros(Nabc+1,1);        % right
H1abc = zeros(Nabc,1);
H2abc = zeros(Nabc,1);

zpmlH = linspace(0,Labc,Nabc)'/Labc;
zpmlE = linspace(Dz/2,Labc-Dz/2,Nabc-1)'/Labc;
%sigmaH2 = linspace(1,Nabc,Nabc)'/Nabc;
abc_pow = 3;
sigma_max = -(abc_pow+1)*log(10^(-4))/(2*eta0*Labc);
sigmaH2 = sigma_max*zpmlH.^abc_pow*mu0/eps0;
sigmaE2 = sigma_max*zpmlE.^abc_pow;
sigmaE1 = sigmaE2((Nabc-1):-1:1);
sigmaH1 = sigmaH2(Nabc:-1:1);

figure(4)
plot(sigmaH2)

%break
%***********************************************************************
%     Wave excitation
%***********************************************************************
Hdelay = Dt/2-Dz/2/c0;

% Incident puls from the left
tau = 80.0e-12;
delay = 2.1*tau;
Ein = sin(omega*(t-delay)).*exp(-((t-delay).^2/tau^2));
Hin = 1/eta0*sin(omega*(t-delay+Hdelay)).*exp(-((t-delay+Hdelay).^2/tau^2));

% Incident time harmonic wave from the left
%Ein = sin(omega*(t));
%Hin = 1/eta0*sin(omega*(t+Hdelay));

Einmax = max(abs(Ein));

%***********************************************************************
% Time stepping
%***********************************************************************
for n = 1:Nt;
 
  % Update H from -Dt/2 to Dt/2
  Hy = Hy - (Dt/mu0/Dz) * diff(Ex,1);                                       % main grid
  
  Hy(Nin) = Hy(Nin) - (Dt/mu0/Dz) * Ein(n);   
  
  H1abc = (H1abc.*(mu0/Dt-sigmaH1/2) - diff(E1abc)/Dz) ./ (mu0/Dt+sigmaH1/2); % ABC to the left
  H2abc = (H2abc.*(mu0/Dt-sigmaH2/2) - diff(E2abc)/Dz) ./ (mu0/Dt+sigmaH2/2); % ABC to the right
  
  % Update E from 0 to Dt
  Ex(2:Nz) = (Ex(2:Nz).*(eps0*epsr/Dt-sig/2) - diff(Hy,1)/Dz)./(eps0*epsr/Dt+sig/2);    % main grid except boundary
  
  Ex(Nin) = Ex(Nin) - Hin(n)/(Dz*eps0/Dt);
  
  Ex(1) = Ex(1) - (Dt /eps0) * (Hy(1)-H1abc(Nabc))/Dz;          % boundary between left ABC and main grid 
  E1abc(Nabc+1) = Ex(1);                                        % add one point to simplify the scheme
  E1abc(2:Nabc) = (E1abc(2:Nabc).*(eps0/Dt-sigmaE1/2) - diff(H1abc,1)/Dz) ./ (eps0/Dt+sigmaE1/2);
  
  Ex(Nz+1) = Ex(Nz+1) - (Dt /eps0) * (H2abc(1)-Hy(Nz))/Dz;
  E2abc(1) = Ex(Nz+1);
  E2abc(2:Nabc) = (E2abc(2:Nabc).*(eps0/Dt-sigmaE2/2) - diff(H2abc,1)/Dz) ./ (eps0/Dt+sigmaE2/2);
  
  % Sample the electric field at chosen points
  if mod(n,n_pict) == 0
      figure(1);
      subplot(2,1,1);
      plot(z,Ex,'LineWidth',1);
      title(['time is ',num2str(n*Dt*10^9),'ns'])
      axis tight;
      ax=axis;
      axis([ax(1:2) -1.2*Einmax 1.2*Einmax]);
      ax=axis;
      ylabel('E-field');
      for p=1:objects
          line([obj_z(p,1) obj_z(p,1)], ax(3:4),'color',obj_c(p));
          line([obj_z(p,2) obj_z(p,2)], ax(3:4),'color',obj_c(p));
          line([obj_z(p,1) obj_z(p,2)], [ax(3) ax(3)],'color',obj_c(p),'LineWidth',4);
          line([obj_z(p,1) obj_z(p,2)], [ax(4) ax(4)],'color',obj_c(p),'LineWidth',4);
      end
      subplot(2,1,2);
      plot([E1abc(:)' Ex' E2abc(:)'],'LineWidth',1)
      title(['time is ',num2str(n*Dt*10^9),'ns'])
      ax = axis;
      line([Nabc Nabc], ax(3:4),'color','r');
      line(Nin+[Nabc Nabc], ax(3:4),'color','b');
      line(Ntrans+2+[Nabc Nabc], ax(3:4),'color','b');
      line([Nz+2+Nabc Nz+2+Nabc], ax(3:4),'color','r');
      axis tight;
      ylabel('E-field');
      %plot(Eabc(2:Nabc,1:2));
      pause(0.01);
  end
  Eref(n) = Ex(Nref);
  Etrans(n) = Ex(Ntrans);
end

%***********************************************************************
%     Post processing
%***********************************************************************
Ehin = fft(Ein,Nfft)/Nt*2;
Df = 1/(Dt*Nfft);
freqN = round(fmax/Df);
freqN1 = round(fmin/Df);
f = linspace(Df*freqN1,Df*freqN,freqN-freqN1+1)/10^9;
t = t*10^9;   % nano seconds [ns]

figure(2);
% incident wave 
clf;
subplot(2,3,1);
plot(t,Ein,'r','LineWidth',linew)
axis tight;
xlabel('time');
ylabel('E-field');

subplot(2,3,4);
warning off;
A=plotyy(f,abs(Ehin(freqN1:freqN)),f,c0./f/Dz/10^9); 
warning on;
set(A(2),'yLim',[0 40],'yTick',[0:5:40],'xlim',[fmin fmax]/10^9,'xGrid','on','yGrid','on','LineWidth',linew);
set(A(1),'xlim',[fmin fmax]/10^9,'LineWidth',linew);
set(get(A(2),'Ylabel'),'String','points/wavelength');
set(get(A(1),'Ylabel'),'String','input');
set(get(A(1),'Xlabel'),'String','frequency');

% reflection
subplot(2,3,2);
plot(t,Eref,'LineWidth',linew); axis tight;
xlabel('time [ns]');
ylabel('Reflected field');
subplot(2,3,5);
Ehref = fft(Eref,Nfft)/Nt*2;
plot(f,20*log10(abs(Ehref(freqN1:freqN)./Ehin(freqN1:freqN))),'LineWidth',linew); axis tight;
ax=axis;
axis([fmin/10^9 fmax/10^9 ax(3:4)]);
axis tight;
xlabel('frequency [GHz]');
ylabel('Reflection in dB');

% transmission
subplot(2,3,3);
plot(t,Etrans,'LineWidth',linew); axis tight;
xlabel('time [ns]');
ylabel('Transmitted field');
subplot(2,3,6);
Ehtrans = fft(Etrans,Nfft)/Nt*2;
plot(f,20*log10(abs(Ehtrans(freqN1:freqN)./Ehin(freqN1:freqN))),'LineWidth',linew); axis tight;
ax=axis;
axis([fmin/10^9 fmax/10^9 ax(3:4)]);
xlabel('frequency [GHz]');
ylabel('Transmission in dB');

tr=interp1(f,abs(Ehtrans(freqN1:freqN)./Ehin(freqN1:freqN)),5)
re=interp1(f,abs(Ehref(freqN1:freqN)./Ehin(freqN1:freqN)),5)
tr^2+re^2