fdtd_1d.m

1D的有限差分时域法计算电磁波散射c语言模拟程序

%======================================================================================
% 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