fdtd_2d_tm_3.m

fdtd一维到三维算法

% fdtd_2d_TM_3.m
% Excercise ....
% 2D FDTD free space propagation striking dielectric slab on top right corner
% Coded by Pradip Mukhopadhyay
% 7 March 2005
%
% For EEE355F AJW
%

clear all;
cc=2.99792458e8;			%speed of light in free space
mu_0=4.0*pi*1.0e-7;			%permeability of free space
eps_0=1.0/(cc*cc*mu_0);			%permittivity of free space


ie = 150;				% Grid pixels in Y direction
je = 150;				% Grid pixels in X direction

is = floor(ie/2) ;
js = floor(je/2) ;

% Electric and Magnetic field

Ez = zeros(ie,je);
Hx = zeros(ie,je);
Hy = zeros(ie,je);

eps_r = ones(ie,je);

id = is-15;				% dielectric interface
jd = js+15;
eps_r(1:id,jd:je) = 40.0;		% dielectric constant


nmax = 120;				% number of time step
ddx = 1.0e-3;				% X grid size
ddy = ddx;				% Y grid size

dt = 0.98/(cc*sqrt( (1/ddx)^2 + (1/ddy)^2 ));

%***********************************************************************
tw = 26.53e-12;
t0 = 4.0*tw;
T = (0:1:nmax-1)'.*dt;

source = -2.0*((T-t0)./tw).*exp(-1.0*((T-t0)./tw).^2.0);

[Emax] = max(source);
[Emin] = min(source);

C1 = dt./(eps_r.*eps_0);
C2 = dt/mu_0;


for n = 1:nmax

    %-----------------------------------update Ez-----------------

    for jj = 2:je
        for ii = 2:ie

            Ez(ii,jj) = Ez(ii,jj) + C1(ii,jj).*((Hy(ii,jj) - Hy(ii-1,jj))./ddx -(Hx(ii,jj) - Hx(ii,jj-1))./ddy);

        end
    end

    %-----------------------------------source---------------------

    Ez(is,js) =  source(n);

    %-----------------------------------update Hx------------------

    for jj = 1:je-1
        for ii = 1:ie

            Hx(ii,jj) = Hx(ii,jj) - C2*((Ez(ii,jj+1) - Ez(ii,jj))./ddy);

        end
    end

    %-----------------------------------update Hy------------------------

    for jj = 1:je
        for ii = 1:ie-1

            Hy(ii,jj) = Hy(ii,jj) + C2*((Ez(ii+1,jj) - Ez(ii,jj))./ddx);

        end
    end



    imagesc(Ez, 1.0*[Emin Emax]);
    patch([jd  jd],[1 id],'w');
    patch([jd  je],[id id],'w');
    title('2D EM propagation --striking a dielectric slab on the right corner');
    colormap(hot);

    pause(0.0002);
end