Finite Difference Time Domain Program... Adding 2 Dimensions to my 1 Dimension Program
    3 vues (au cours des 30 derniers jours)
  
       Afficher commentaires plus anciens
    
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
0 commentaires
Réponses (0)
Voir également
Catégories
				En savoir plus sur Analysis dans Help Center et File Exchange
			
	Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!