Finite Difference Time Domain Program... Adding 2 Dimensions to my 1 Dimension Program

2 次查看（过去 30 天）

cestcandice 2017-11-16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/367306-finite-difference-time-domain-program-adding-2-dimensions-to-my-1-dimension-program

I am writing a program about Finite Difference Time Domain. This program shows E&M waves in 1 dimension. I'm trying to figure out how I can add 2 dimensions to this particular program.

close all;
clear all;
%define 
NN = 200; %# of iterations to perform
N = 500; %# spatial samples in domain
L = 5; %domain length in meters
fs = 300e6; %source frequency in Hz
ds = L/N; %spatial step in meters
x = linspace(0,L,N); %x coordinate of spatial samples
dt = ds/300e6; %"magic time step"
eps0 = 8.854e-12; %permittivity of free space
mu0 = pi*4e-7; %permeability of free space
%scale factors for E and H
ae = ones(N,1)*dt/(ds*eps0);
am = ones(N,1)*dt/(ds*mu0);
as = ones(N,1);
epsr = ones(N,1);
mur= ones(N,1);
sigma = zeros(N,1);
p = 1;
for i=1:N
   epsr(i) = 1;
   mur(i) = 1;
   w1 = 0.5;
   w2 = 1.5;
   if (p==1) %dielectric window
      if (abs(x(i)-L/2)<0.5) epsr(i)=4; end
   end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ae = ae./epsr;
am = am./mur;
ae = ae./(1+dt*(sigma./epsr)/(2*eps0));
as = (1-dt*(sigma./epsr)/(2*eps0))./(1+dt*(sigma./epsr)/(2*eps0));
%initialize fields 
Hz = zeros(N,1);
Ey = zeros(N,1);
for iter=1:NN
   %Yee Updataed Algorithm 
   Ey(3) = Ey(3)+2*(1-exp(-((iter-1)/50)^2))*sin(2*pi*fs*dt*iter); 
     %%%%%%%%%%%%%%%%%
     Hz(1) = Hz(2); %absorbing boundary conditions for H
     for i=2:N-1 %H field Update
        Hz(i) = Hz(i)-am(i)*(Ey(i+1)-Ey(i));
     end
     %%%%%%%%%%%%%%%
     Ey(N) = Ey(N-1); %absorbing boundary conditions for E
     for i=2:N-1 % E field Update
        Ey(i) = as(i)*Ey(i)-ae(i)*(Hz(i)-Hz(i-1));
     end
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     figure(1)
     hold off
     plot(x,Ey,'b');
     axis([3*ds L -2 2]);
     grid on;
     title('FDTD Gaussian Pulse');
     hold on
     plot(x,380*Hz,'r');%multiplied by 380 so that the amplitude will reach same as E.
     legend('E Field','H-Field')
     xlabel('Time (Mintues)');
     ylabel('Amplitude') 
     pause(0);
  end