fdtd_1d.m

FDTD的一维和三维i算法

%文件1

%======================================================================================
% function Ex_obs = FDTD_1D(Ex,dx,dt,srcIndex,srcJx,pf,nSteps,e_r,mu_r,sigma,sigma_m,k_obs)
%
% This function performs one-dimensional (Cartesian) finite-difference time-domain 
% simulations. It assumes the spatial grid extends from k = 1 to k = Nz and that the 
% simulation runs for 'nSteps' time steps. 
% The function should take the following input arguments, in this order:
%
%   Input Parameter         Size        Description
%   Ex                      1xNz        initial condition for Ex (V/m)
%   dx                      1x1         spatial step size (m)
%   dt                      1x1         time step size (s)
%   Src_indx                1x1         grid index for source location
%   Src_Jx                  1xnSteps    source current stimulus (A/m2)
%   pf                      1x1         plot flag - # steps between plots
%   nSteps                  1x1         # steps for simulation
%   e_r                     1xNz        relative permittivity
%   mu_r                    1x(Nz-1)    relative permeability
%   sigma                   1xNz        conductivity (S/m)
%   sigma_m                 1x(Nz-1)    magnetic conductivity (Wb/C穖)
%   k_obs                   1x1         index of observation point
% The one output parameter should be Ex_obs, a 1xnSteps vector
% of values at the observation point. If no observation point is desired,
% set k_obs = [ ] (an empty matrix) and return [ ] in Ex_obs.
% The function also generates a plot of Ex every 'pf' time steps
% so that we can observe what is going on
%
% Written by
% Aroh Barjatya
% For Utah State University class ECE 5830 Electromagnetics 2

function Ex_obs = FDTD_1D(Ex,dx,dt,srcIndex,srcJx,pf,nSteps,e_r,mu_r,sigma,sigma_m,k_obs)

c1 = ((2 .* e_r) - (sigma .* dt)) ./ ((2 .* e_r) + (sigma .* dt));
c2 = (2 .* dt) ./ (dx .* ((2 .* e_r) + (sigma .* dt)));
c3 = (2 .* dt) ./ ((2 .* e_r) + (sigma .* dt));
c4 = ((2 .* mu_r) - (sigma_m .* dt)) ./ ((2 .* mu_r) + (sigma_m .* dt));
c5 = (2 .* dt) ./ (dx .* ((2 .* mu_r) + (sigma_m .* dt)));

Ex(srcIndex) = srcJx(1);
Hy = zeros(1,length(Ex) - 1);
Ex_obs = zeros(1,nSteps);
if length(k_obs) > 0
    Ex_obs(1) = Ex(k_obs);
end

for i = 2:nSteps
    Ex(2:end-1) = c1(2:end-1) .* Ex(2:end-1) - c2(2:end-1) .* (Hy(2:end)-Hy(1:end-1));
    Ex(srcIndex) = Ex(srcIndex) - c3(srcIndex)*srcJx(i);
    Hy = (c4 .* Hy) - (c5 .* (Ex(2:end) - Ex(1:end-1)));
    if (mod(i, pf) == 0)
        plot(Ex);
        a = sprintf('step number = %d', i);
        title(a);
        set(gcf,'color','white');
        xlabel('Grid Points -->');
        ylabel('Magnitude -->');
%         axis([1 1001 -1 1]);
        pause(0.2);
    end
    if length(k_obs) > 0
        Ex_obs(i) = Ex(k_obs) ;
    end
end

%文件2

% function p=Gaussian_pulse(dt,f_c,len)
%
% The function creates a Gaussian pulse centered at frequency fc
% It assumes a sample frequency of fs = 1/dt
%
% the pulse vector will be of length len; specify too short of a value for len
% and an error will be generated; delay is needed before the start of the pulse
% to reduce the step function transient created at t=0
%
% the resulting pulse has a lowpass characteristic with a 3 dB bandwidth f_c 
% (i.e. extending from roughly DC to f_c)
%
% Originally Written by Michael Tompkins, Asst. Prof at USU Fall 2004
% Minor modifications by Aroh Barjatya

function p=Gaussian_pulse(dt,f_c,len)
if 0
   Gaussian_pulse(1e-7,2e5,15000);
   Gaussian_pulse(1.6e-10,30e6,15000);
end

spread = 1/2/pi/f_c/dt; % this is the standard deviation of our Gaussian function

% must wait a while to launch the pulse so that we don't get too big of a 
% "jump" (step function) at t=0; this step has higher frequency harmonics
% in it that will not meet our requirement that dx << lambda

delay = 5*spread; 
if delay>len
   error('Need longer simulation to accommodate desired pulse');
end

n = 0:len-1;
p = exp(-((n-delay)/spread).^2); % Gaussian pulse

%文件3

%=====================================================================================
% This is a test script for FDTD_1D
% We test it by giving it a Gaussian-pulse
% plane-wave source in the middle of a 1001-point grid.
% Assumed free space everywhere.
%
% Written by
% Aroh Barjatya
% For Utah State University class ECE 5830 Electromagnetics 2

clear all
close all

% Physical Constants
e_0 = 8.85e-12;
mu_0 = 4e-7 * pi;
c = 2.99e8;
eta_0 = sqrt(mu_0/e_0);

% Simulation Constants
Nz = 1001;
nSteps = 3500;
Ex = zeros(1,Nz);
dx = 0.1;
dt = dx / (2*c);
e_r = ones(1,Nz) .* e_0;
mu_r = ones(1,Nz - 1) .* mu_0;
pf = 20;

% Define the source and its location
pulse = Gaussian_pulse(dt,30e6,nSteps);
srcJx = (-2/(eta_0*dx))*pulse;
srcIndex = (Nz - 1)/2;

sigma = zeros(1,Nz);
sigma_m = zeros(1,Nz-1);

FDTD_1D(Ex, dx, dt, srcIndex, srcJx, pf, nSteps, e_r, mu_r, sigma, sigma_m, []);


%文件4

%====================================================================================
% This is the second test script for FDTD_1D
% In addition to what was done in test cript 1, we add a 50 unit think perfectly 
% matched layer (PML) on both ends of the simulation space.
%
% Written by
% Aroh Barjatya
% For Utah State University class ECE 5830 Electromagnetics 2


clear all
close all

% Physical Constants
e_0 = 8.85e-12;
mu_0 = 4e-7 * pi;
c = 2.99e8;
eta_0 = sqrt(mu_0/e_0);

% Simulation Constants
Nz = 1001;
nSteps = 3500;
Ex = zeros(1,Nz);
dx = 0.1;
dt = dx / (2*c);
e_r = ones(1,Nz) .* e_0;
mu_r = ones(1,Nz - 1) .* mu_0;
pf = 20;

% Define the source and its location
pulse = Gaussian_pulse(dt,30e6,nSteps);
srcJx = (-2/(eta_0*dx))*pulse;
srcIndex = (Nz - 1)/2;

sigma = zeros(1,Nz);
sigma_m = zeros(1,Nz-1);

% for PML the following lines define sigma and sigma_m
% if no PML is required, then comment the following lines

% The Left Edge
k = 50;
k0 = 50;
for i = k:-1:1
    sigma(i) = ((2 * e_r(i)) / dt) * 0.33 * ((k0 - i + 1) / k)^3;
    sigma_m(i) = sigma(i) * mu_r(i) / e_r(i);
end

% The Right Edge
k0 = Nz-k;
for i = Nz-k:Nz
    sigma(i) = ((2 * e_r(i)) / dt) * 0.33 * ((i - k0 + 1) / k)^3;
    sigma_m(i-1) = sigma(i) * mu_r(i-1) / e_r(i);
end

FDTD_1D(Ex, dx, dt, srcIndex, srcJx, pf, nSteps, e_r, mu_r, sigma, sigma_m, []);