一维等离子体FDTD的Matlab源代码(两种方法)

%***********************************************************************

%***********************************************************************

%

%

%

%

%

%***********************************************************************

function [Ex,Ey,Ez]=FDTD3D_Main(handles)

global SimRunStop

% if ~isdir('C:\\MATLAB7\\work\\cavity\\figures')

% end

%***********************************************************************

%***********************************************************************

p.ip = get(handles.XdirPar,'Value');

p.jp = get(handles.YdirPar,'Value');

p.PN = get(handles.PartNo,'Value');

%***********************************************************************

%***********************************************************************

kb=ke+1;

%***********************************************************************

%***********************************************************************

Ex=zeros(ie,jb,kb);

Ey=zeros(ib,je,kb);

Ez=zeros(ib,jb,ke);

Hx=zeros(ib,je,ke);

Hy=zeros(ie,jb,ke);

Hz=zeros(ie,je,kb);

%***********************************************************************

%***********************************************************************

param.epsz = 1.0/(param.cc*param.cc*param.muz); %permittivity of free space

%***********************************************************************

%***********************************************************************

%***********************************************************************

%***********************************************************************

param.rtau=get(handles.tau,'Value')*100e-12;

param.tau=param.rtau/param.dt;

param.ndelay=3*param.tau;

param.Amp = get(handles.sourceamp,'Value')*10e11;

param.srcconst=-param.dt*param.Amp;

%***********************************************************************

%***********************************************************************

param.eps= get(handles.epsilon,'Value');

%***********************************************************************

%***********************************************************************

param.ca=(1.0-(param.dt*param.sig)/(2.0*param.epsz*param.eps))/(1.0+(param.dt*param.sig)/(2.0*param.epsz*param.eps));

param.cb=(param.dt/param.epsz/param.eps/param.dx)/(1.0+(param.dt*param.sig)/(2.0*param.epsz*param.eps));

param.da=1.0;

param.db=param.dt/param.muz/param.dx;

%***********************************************************************

%***********************************************************************

ex=zeros(ib,jb,kb);

ey=zeros(ib,jb,kb);

ez=zeros(ib,jb,kb);

hx=zeros(ib,jb,kb);

hy=zeros(ib,jb,kb);

hz=zeros(ib,jb,kb);

[X,Y,Z] = meshgrid(1:ib,1:jb,1:kb); % Grid coordinates

Psim = zeros(param.nmax,1);

Panl = zeros(param.nmax,1);

if ((p.ip == 1)&(p.jp == 0))

    x = ceil(ie/p.PN)    

        for m1=1:1:p.PN                        

        end                

    end    

elseif ((p.ip == 0)&(p.jp == 1))

    x = ceil(ie/p.PN);    

    end    

                SimRunStop = get(handles.Stop_sim,'value');   

                    end                          

                pause(0.01);                

    end    

end

%文件2

% Cavity Field Viz. %

function [] = field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p)

%***********************************************************************

%***********************************************************************

tview = squeeze(ez(:,:,param.kobs));

sview = squeeze(ez(:,param.js,:));

ax1 = handles.axes1;

ax2 = handles.axes2;

ax3 = handles.axes3;

%***********************************************************************

%***********************************************************************

timestep=int2str(n);

ezview=squeeze(ez(:,:,param.kobs));

exview=squeeze(ex(:,:,param.kobs));

eyview=squeeze(ey(:,:,param.kobs));

xmin = min(X(:));

xmax = max(X(:));

ymin = min(Y(:));

ymax = max(Y(:));

zmin = min(Z(:));

daspect([2,2,1])

xrange = linspace(xmin,xmax,8);

yrange = linspace(ymin,ymax,8);

zrange = 3:4:15;

[cx cy cz] = meshgrid(xrange,yrange,zrange);

% sview=squeeze(ez(:,param.js,:));

rect = [-50 -35 360 210];

rectp = [-50 -40 350 260];

axes(ax1);

axis tight

set(gca,'nextplot','replacechildren');

E_total = sqrt(ex.^2 + ey.^2 + ez.^2);

etview=squeeze(E_total(:,:,param.kobs));

sview = squeeze(E_total(:,param.js,:));

surf(etview');

shading interp;

caxis([-1.0 1.0]);

colorbar;

axis image;

title(['Total E-Field, time step = ',timestep],'fontname','Times New Roman','fontsize',12,'FontWeight','BOLD');

xlabel('i coordinate');

ylabel('j coordinate');

set(gca,'fontname','Times New Roman','fontsize',10);

% F1 = getframe(gca,rect);

% M1 = frame2im(F1);

% imwrite(M1,filename,'jpeg')

axes(ax2);

axis tight

set(gca,'nextplot','replacechildren');

surf(sview');

shading interp;

caxis([-1.0 1.0]);

colorbar;

axis image;

title(['Ez(i,j=13,k), time step = ',timestep],'fontname','Times New Roman','fontsize',12,'FontWeight','BOLD');

xlabel('i coordinate');

ylabel('k coordinate');

set(gca,'fontname','Times New Roman','fontsize',10);

% F2 = getframe(gca,rect);

% M2 = frame2im(F2);

% imwrite(M2,filename,'jpeg')

%***********************************************************************

% Cavity Power - Analitic expression 

%***********************************************************************

axes(ax3);

% axis tight

% set(gca,'nextplot','replacechildren');

t = [1:1:param.nmax];

Psim = 1e3*Psim./max(Psim);

Panl = 1e3*Panl./max(Panl);

semilogy(t,Psim,'b.-',t,Panl,'rx-');

str(1) = {'Normalized Analytic Vs. Simulated Power'};

str(2) = {'as function of time-steps'};

title(str,'fontname','Times New Roman','fontsize',12,'FontWeight','BOLD' );

xlabel('Time Step');

ylabel('Log(Power)');

legend('Simulation','Analytic','location','SouthEast');

set(gca,'fontname','Times New Roman','fontsize',10);

% F3 = getframe(gca,rectp);

% M3 = frame2im(F3);

% imwrite(M3,filename,'jpeg')

%文件3

% Source waveform function

function [source]=waveform(param,handles,n)

ax1 = handles.axes1;

type = get(handles.sourcetype,'value');

amp = get(handles.sourceamp,'value')*1e4;

f = get(handles.Frequency,'value')*1e9;

switch type

end

%文件4

function [hx,hy,hz] = Hfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,x,y,p)

a = param.a(x);

b = param.b(x);

c = param.c(y);

d = param.d(y);

hx(a+1:b,c:d,1:ke)=...

hy(a:b,c+1:d,1:ke)=...

hz(a:b,c:d,2:ke)=...

%文件5

function [ex,ey,ez]=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,x,y,p)

a = param.a(x);

b = param.b(x);

c = param.c(y);

d = param.d(y);

if (param.is>=a)&(param.is<=b)|(param.js>=c)&(param.js<=d)

else

end

ex(a:b,c+1:d,2:ke)=...

ey(a+1:b,c:d,2:ke)=...

ez(a+1:b,c+1:d,1:ke)=...

ez(param.is,param.js,1:ke)=ez(param.is,param.js,1:ke)+source;

function FDTDonedimensionpipei(L,d,T)

%version1.0  终端匹配

%FDTDonedimensionpipei(6,0.18,0.5e-9)

t0=3*T;

c=3e8;

u=4*pi*1e-7;

e=8.8541878e-12;

dz=T*c/10;

Nz=L/dz;

Nz=fix(Nz);

dt=dz/2/c;

Ex=zeros(1,Nz+1);

B=zeros(1,Nz+1);

Hy=zeros(1,Nz);

Nt=2*Nz;

for n=0:Nt

    t=n*dt;

    F=exp(-(t-t0).^2./T^2);

    Ex(1)=F;

    for k=1:Nz

        Hy(k)=Hy(k)+dt./u.*(Ex(k)-Ex(k+1))./dz;

    end

    for k=1:Nz-1

        Ex(k+1)=Ex(k+1)+dt./e.*(Hy(k)-Hy(k+1))./dz;

    end

    Ex(1)=B(2)+(c*dt-dz)./(c*dt+dz).*(Ex(2)-B(1));

    Ex(Nz+1)=B(Nz)+(c*dt-dz)./(c*dt+dz).*(Ex(Nz)-B(Nz+1));

    Vref1=d.*Ex(Nz-300);

    Vref2=d.*Ex(Nz-100);

    plot(t,Vref1,'s');

    hold on;

    plot(t,Vref2,'rx');

    hold on;

    B=Ex;

end

function FDTD_debug

%constants

c_0 = 3E8;                  % Speed of light in free space

Nz = 100;                    % Number of cells in z-direction

Nt = 200;                    % Number of time steps

N = Nz;                      % Global number of cells

E = zeros(Nz,1);            % E component

E2 = zeros(Nz,1);

E1 = zeros(Nz,1);

Es = zeros(Nz,Nz);          % Es is the output for surf function

%H = zeros(Nz,1);            % H component

pulse = 0;        % Pulse

%to be set

mu_0 = 4.0*pi*1.0e-7;      % Permeability of free space

eps_0 = 1.0/(c_0*c_0*mu_0);  % Permittivty of free space

c_ref = c_0;                % Reference velocity

eps_ref = eps_0;

mu_ref = mu_0;

f_0 = 10e9;                  % Excitation frequency

f_ref = f_0;                % Reference frequency 

omega_0 = 2.0*pi*f_0;        % Excitation circular frequency

omega_ref = omega_0;

lambda_ref = c_ref/f_ref;    % Reference wavelength

Dx_ref = lambda_ref/20;      % Reference cells width

Dz = Dx_ref;

nDz = Dz/Dx_ref;

Z = Nz*Dz;

r = 1;                      % Normalized time step

Dt_ref = r*Dx_ref/c_ref;    % Reference time step

Dt = Dt_ref;

% Source position

Nz_Source = 0.5*Nz;

N_Source = Nz_Source;

for n = 1:Nt-1

    t = Dt_ref*r*(n-1);    % Actual time

    %Source function series

    Source_type = 1;

    switch Source_type 

        case 1  % modified source function

            ncycles = 2;

            if t < ncycles*2.0*pi/(omega_0)  

                pulse = -0.5*( 1.0 - cos(omega_0*t/ncycles) ) * sin(omega_0*t);

            else 

                pulse = 0;

            end

        case 2  % sigle cos source function

            if t < 2.0*pi/(omega_0)

                pulse = 8*c_0^2*Dt_ref^2*mu_ref*omega_0*cos(omega_0*t);

            else 

                pulse = 0;

            end

        case 3  % Gaussian pulse

            if t < Dt_ref*r*50

                pulse = -40*c_0*(t-t*25/(n-1))*exp(-(t-t*25/(n-1))^2/2/(50/2.3548)^2/(t/(n-1))^2);

            else

                pulse = 0;

            end      

    end

    E2(N_Source) = E2(N_Source) - r*pulse;

    E1 = E2;

    E2 = E;

    m = 2:Nz-1;

    E(m) = r^2*(E2(m+1) - 2*E2(m) + E2(m-1)) + 2*E2(m) - E1(m);

    %Boundary

    E(1) = 0; E(Nz) = 0;                                % Dirichlet

    %E(1) = E(2);E(Nz) = E(Nz-1);                      % Neumann

    %E(Nz) = E2(Nz-1) + (r-1)/(r+1)*(E(Nz-1)-E2(Nz));    % Mur ABC @ right side z = Nz

    %E(1) = E2(2) + (r-1)/(r+1)*(E(2)-E2(1));            % Mur ABC @ left side z = 0

    %display

    %choice***********************************************

    display = 3; % choice: '1' = line plot; '2' = imagesc; '3' = surf

    switch display

        case 1

            i = 1:Nz;

            plot(i, E(i));

            axis([0 Nz -4 4]);

        case 2

            A = rot90(E);

            imagesc(A);

        case 3

            i = 1:Nz;

            for j = 1:Nz

                Es(i,j) = E(i);

            end

            Es = rot90(Es);

            j = 1:Nz;

            surf(i,j,Es);

            axis([0 Nz 0 Nz -4 4])    

        otherwise

            A = rot90(E);

            imagesc(A);

    end 

    pause(0.005);  

end

%***********************************************************************

% 3-D FDTD code with PEC boundaries

%***********************************************************************

%

% Program author: Susan C. Hagness

% Department of Electrical and Computer Engineering

% University of Wisconsin-Madison

% 1415 Engineering Drive

% Madison, WI 53706-1691

% 608-265-5739

% *****************.edu

%

% Date of this version: February 2000

%

% This MATLAB M-file implements the finite-difference time-domain

% solution of Maxwell's curl equations over a three-dimensional

% Cartesian space lattice comprised of uniform cubic grid cells.

%

% To illustrate the algorithm, an air-filled rectangular cavity

% resonator is modeled. The length, width, and height of the

% cavity are 10.0 cm (x-direction), 4.8 cm (y-direction), and

% 2.0 cm (z-direction), respectively.

%

% The computational domain is truncated using PEC boundary

% conditions:

% ex(i,j,k)=0 on the j=1, j=jb, k=1, and k=kb planes

% ey(i,j,k)=0 on the i=1, i=ib, k=1, and k=kb planes

% ez(i,j,k)=0 on the i=1, i=ib, j=1, and j=jb planes

% These PEC boundaries form the outer lossless walls of the cavity.

%

% The cavity is excited by an additive current source oriented

% along the z-direction. The source waveform is a differentiated

% Gaussian pulse given by

% J(t)=-J0*(t-t0)*exp(-(t-t0)^2/tau^2),

% where tau=50 ps. The FWHM spectral bandwidth of this zero-dc-

% content pulse is approximately 7 GHz. The grid resolution

% (dx = 2 mm) was chosen to provide at least 10 samples per

% wavelength up through 15 GHz.

%

% To execute this M-file, type "fdtd3D" at the MATLAB prompt.

% This M-file displays the FDTD-computed Ez fields at every other

% time step, and records those frames in a movie matrix, M, which

% is played at the end of the simulation using the "movie" command.

%

%***********************************************************************

clear

%***********************************************************************

% Fundamental constants

%***********************************************************************

cc=2.99792458e8; %speed of light in free space

muz=4.0*pi*1.0e-7; %permeability of free space

epsz=1.0/(cc*cc*muz); %permittivity of free space

%***********************************************************************

% Grid parameters

%***********************************************************************

ie=50; %number of grid cells in x-direction

je=24; %number of grid cells in y-direction

ke=10; %number of grid cells in z-direction

ib=ie+1;

jb=je+1;

kb=ke+1;

is=26; %location of z-directed current source

js=13; %location of z-directed current source

kobs=5;

dx=0.002; %space increment of cubic lattice

dt=dx/(2.0*cc); %time step

nmax=500; %total number of time steps

%***********************************************************************

% Differentiated Gaussian pulse excitation

%***********************************************************************

rtau=50.0e-12;

tau=rtau/dt;

ndelay=3*tau;

srcconst=-dt*3.0e+11;

%***********************************************************************

% Material parameters

%***********************************************************************

eps=1.0;

sig=0.0;

%***********************************************************************

% Updating coefficients

%***********************************************************************

ca=(1.0-(dt*sig)/(2.0*epsz*eps))/(1.0+(dt*sig)/(2.0*epsz*eps));

cb=(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps));

da=1.0;

db=dt/muz/dx;

%***********************************************************************

% Field arrays

%***********************************************************************

ex=zeros(ie,jb,kb);

ey=zeros(ib,je,kb);

ez=zeros(ib,jb,ke);

hx=zeros(ib,je,ke);

hy=zeros(ie,jb,ke);

hz=zeros(ie,je,kb);

%***********************************************************************

% Movie initialization

%***********************************************************************

tview(:,:)=ez(:,:,kobs);

sview(:,:)=ez(:,js,:);

subplot('position',[0.15 0.45 0.7 0.45]),pcolor(tview');

shading flat;

caxis([-1.0 1.0]);

colorbar;

axis image;

title(['Ez(i,j,k=5), time step = 0']);

xlabel('i coordinate');

ylabel('j coordinate');

subplot('position',[0.15 0.10 0.7 0.25]),pcolor(sview');

shading flat;

caxis([-1.0 1.0]);

colorbar;

axis image;

title(['Ez(i,j=13,k), time step = 0']);

xlabel('i coordinate');

ylabel('k coordinate');

rect=get(gcf,'Position');

rect(1:2)=[0 0];

M=moviein(nmax/2,gcf,rect);

%***********************************************************************

% BEGIN TIME-STEPPING LOOP

%***********************************************************************

for n=1:nmax

%***********************************************************************

% Update electric fields

%***********************************************************************

ex(1:ie,2:je,2:ke)=ca*ex(1:ie,2:je,2:ke)+...

cb*(hz(1:ie,2:je,2:ke)-hz(1:ie,1:je-1,2:ke)+...

hy(1:ie,2:je,1:ke-1)-hy(1:ie,2:je,2:ke));

ey(2:ie,1:je,2:ke)=ca*ey(2:ie,1:je,2:ke)+...

cb*(hx(2:ie,1:je,2:ke)-hx(2:ie,1:je,1:ke-1)+...

hz(1:ie-1,1:je,2:ke)-hz(2:ie,1:je,2:ke));

ez(2:ie,2:je,1:ke)=ca*ez(2:ie,2:je,1:ke)+...

cb*(hx(2:ie,1:je-1,1:ke)-hx(2:ie,2:je,1:ke)+...

hy(2:ie,2:je,1:ke)-hy(1:ie-1,2:je,1:ke));

ez(is,js,1:ke)=ez(is,js,1:ke)+...

srcconst*(n-ndelay)*exp(-((n-ndelay)^2/tau^2));

%***********************************************************************

% Update magnetic fields

%***********************************************************************

hx(2:ie,1:je,1:ke)=hx(2:ie,1:je,1:ke)+...

db*(ey(2:ie,1:je,2:kb)-ey(2:ie,1:je,1:ke)+...

ez(2:ie,1:je,1:ke)-ez(2:ie,2:jb,1:ke));

hy(1:ie,2:je,1:ke)=hy(1:ie,2:je,1:ke)+...

db*(ex(1:ie,2:je,1:ke)-ex(1:ie,2:je,2:kb)+...

ez(2:ib,2:je,1:ke)-ez(1:ie,2:je,1:ke));

hz(1:ie,1:je,2:ke)=hz(1:ie,1:je,2:ke)+...

db*(ex(1:ie,2:jb,2:ke)-ex(1:ie,1:je,2:ke)+...

ey(1:ie,1:je,2:ke)-ey(2:ib,1:je,2:ke));

%***********************************************************************

% Visualize fields

%***********************************************************************

if mod(n,2)==0;

timestep=int2str(n);

tview(:,:)=ez(:,:,kobs);

sview(:,:)=ez(:,js,:);

subplot('position',[0.15 0.45 0.7 0.45]),pcolor(tview');

shading flat;

caxis([-1.0 1.0]);

colorbar;

axis image;

title(['Ez(i,j,k=5), time step = ',timestep]);

xlabel('i coordinate');

ylabel('j coordinate');

subplot('position',[0.15 0.10 0.7 0.25]),pcolor(sview');

shading flat;

caxis([-1.0 1.0]);

colorbar;

axis image;

title(['Ez(i,j=13,k), time step = ',timestep]);

xlabel('i coordinate');

ylabel('k coordinate');

nn=n/2;

M(:,nn)=getframe(gcf,rect);

end;

%***********************************************************************

% END TIME-STEPPING LOOP

%***********************************************************************

end

movie(gcf,M,0,10,rect);