FDTD(时域有限差分法)算法的Matlab源程序
时域有限差分法(FDTD算法)的基本原理及仿真
时域有限差分法(FDTD算法)的基本原理及仿真时域有限差分法（FDTD 算法）时域有限差分法是1966年K.S.Yee 发表在AP 上的一篇论文建立起来的，后被称为Yee 网格空间离散方式。
这种方法通过将Maxwell 旋度方程转化为有限差分式而直接在时域求解, 通过建立时间离散的递进序列, 在相互交织的网格空间中交替计算电场和磁场。
FDTD 算法的基本思想是把带时间变量的Maxwell 旋度方程转化为差分形式，模拟出电子脉冲和理想导体作用的时域响应。
需要考虑的三点是差分格式、解的稳定性、吸收边界条件。
有限差分通常采用的步骤是：采用一定的网格划分方式离散化场域；对场内的偏微分方程及各种边界条件进行差分离散化处理，建立差分格式，得到差分方程组；结合选定的代数方程组的解法，编制程序，求边值问题的数值解。
1.FDTD 的基本原理FDTD 方法由Maxwell 旋度方程的微分形式出发，利用二阶精度的中心差分近似，直接将微分运算转换为差分运算，这样达到了在一定体积内和一段时间上对连续电磁场数据的抽样压缩。
Maxwell 方程的旋度方程组为：E E H σε+∂∂=⨯∇t H HE m tσμ-∂∂-=⨯∇ （1） 在直角坐标系中，（1）式可化为如下六个标量方程：⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫+∂∂=∂∂-∂∂+∂∂=∂∂-∂∂+∂∂=∂∂-∂∂z z x y y y z x x x yz E t E y H x H E t E x H z H E t E z H y H σεσεσε，⎪⎪⎪⎪⎭⎪⎪⎪⎪⎬⎫-∂∂-=∂∂-∂∂-∂∂-=∂∂-∂∂-∂∂-=∂∂-∂∂z m zx y y m y z x x m x y z H t H y E x E H t H x E z E H t H z E y E σμσμσμ （2）上面的六个偏微分方程是FDTD 算法的基础。
Yee 首先在空间上建立矩形差分网格，在时刻t n ∆时刻，F(x,y,z)可以写成),,(),,,(),,,(k j i F t n z k y j x i F t z y x F n =∆∆∆∆= （3）用中心差分取二阶精度： 对空间离散：()[]2),,21(),,21(),,,(x O xk j i F k j i F x t z y x F n n xi x ∆+∆--+≈∂∂∆= ()[]2),21,(),21,(),,,(y O yk j i F k j i F y t z y x F n n yj y ∆+∆--+≈∂∂∆= ()[]2)21,,()21,,(),,,(z O zk j i F k j i F z t z y x F n n zk z ∆+∆--+≈∂∂∆=对时间离散：()[]22121),,(),,(),,,(t O tk j i F k j i F t t z y x F n n tn t ∆+∆-≈∂∂-+∆= （4） Yee 把空间任一网格上的E 和H 的六个分量，如下图放置:oyxzEyHzExEzHxEyEyEzEx HyEzEx图1 Yee 氏网格及其电磁场分量分布在FDTD 中，空间上连续分布的电磁场物理量离散的空间排布如图所示。
时域有限差分法的Matlab仿真
时域有限差分法的Matlab仿真关键词: Matlab 矩形波导时域有限差分法摘要：介绍了时域有限差分法的基本原理，并利用Matlab仿真，对矩形波导谐振腔中的电磁场作了模拟和分析。
关键词：时域有限差分法；Matlab；矩形波导；谐振腔目前，电磁场的时域计算方法越来越引人注目。
时域有限差分（Finite Difference Time Domain，FDTD）法［1］作为一种主要的电磁场时域计算方法，最早是在1966年由K. S. Yee提出的。
这种方法通过将Maxwell旋度方程转化为有限差分式而直接在时域求解，通过建立时间离散的递进序列，在相互交织的网格空间中交替计算电场和磁场。
经过三十多年的发展，这种方法已经广泛应用到各种电磁问题的分析之中。
Matlab作为一种工程仿真工具得到了广泛应用［2］。
用于时域有限差分法，可以简化编程，使研究者的研究重心放在FDTD法本身上，而不必在编程上花费过多的时间。
下面将采用FDTD法，利用Matlab仿真来分析矩形波导谐振腔的电磁场，说明了将二者结合起来的优越性。
1FDTD法基本原理时域有限差分法的主要思想是把Maxwell方程在空间、时间上离散化，用差分方程代替一阶偏微分方程，求解差分方程组，从而得出各网格单元的场值。
FDTD 空间网格单元上电场和磁场各分量的分布如图1所示。
电场和磁场被交叉放置，电场分量位于网格单元每条棱的中心，磁场分量位于网格单元每个面的中心，每个磁场（电场）分量都有4个电场（磁场）分量环绕。
这样不仅保证了介质分界面上切向场分量的连续性条件得到自然满足，而且还允许旋度方程在空间上进行中心差分运算，同时也满足了法拉第电磁感应定律和安培环路积分定律，也可以很恰当地模拟电磁波的实际传播过程。
1.1Maxwell方程的差分形式旋度方程为：将其标量化，并将问题空间沿3个轴向分成若干网格单元，用Δx，Δy和Δz 分别表示每个网格单元沿3个轴向的长度，用Δt表示时间步长。
基于MATLAB的FDTD算法编程
FD TD ar ithm etic programm ing ba sed on M ATLAB language
ZHAO J ia (Com pu ter Engineering D epartm en t, Guangx i U n iversity
of T echno logy, L iuzhou 545006, Ch ina) Abstract: T he basic p rincip le of FD TD w as in troduced and the 22D TM m ode FD TD arithm etic exp ression w as derived. T he p rogramm ing m ethod based on M A TLAB w as illu strated w ith exam p le. Key words: M A TLAB; FD TD; Yee A rithm etic; P rogramm ing
自由空间的光速。∃x、∃y 和 ∃ t 的选取必须使它们的关系满足 Cournan t 稳定性条件。尽管时间和空间步长越
小, FD TD 运算的精度会越高, 但太细的网格会使计算机运算时间太长, 同时, 在网格分辨率 (即波长与空间
步长之比) 小到一定程度时, 数值误差基本就不变了。对于采用M A TLAB 编程, 需要特别注意时间步长和空
211 算例
考察区域为正方形的二维自由空间, 区域的边
界为理想电导体, 为了激励外向辐射波, 设在区域的
中心有电场分量 E z 随时间按高斯变化, 求解该区
域内电场随时间的变化, 场域如图 2 所示。
212 M ATLAB 编程方法及结果
基于M A TLAB 对本算例 FD TD 算法进行编程 求 解可归纳为以下几个步骤: (1) 正确建立二维
matlab有限差分法
matlab有限差分法一、前言Matlab是一种广泛应用于科学计算和工程领域的计算机软件，它具有简单易学、功能强大、易于编程等优点。
有限差分法（Finite Difference Method）是一种常用的数值解法，它将微分方程转化为差分方程，通过对差分方程进行离散化求解，得到微分方程的数值解。
本文将介绍如何使用Matlab实现有限差分法。
二、有限差分法基础1. 有限差分法原理有限差分法是一种通过将微分方程转化为离散形式来求解微分方程的数值方法。
其基本思想是将求解区域进行网格划分，然后在每个网格点上进行逼近。
假设要求解一个二阶常微分方程：$$y''(x)=f(x,y(x),y'(x))$$则可以将其转化为离散形式：$$\frac{y_{i+1}-2y_i+y_{i-1}}{h^2}=f(x_i,y_i,y'_i)$$其中$h$为网格步长，$y_i$表示在$x_i$处的函数值。
2. 一维情况下的有限差分法对于一维情况下的常微分方程：$$\frac{d^2 y}{dx^2}=f(x,y,y')$$可以使用中心差分法进行离散化：$$\frac{y_{i+1}-2y_i+y_{i-1}}{h^2}=f(x_i,y_i,y'_i)$$这个方程可以写成矩阵形式：$$A\vec{y}=\vec{b}$$其中$A$为系数矩阵，$\vec{y}$为函数值向量，$\vec{b}$为右端项向量。
三、Matlab实现有限差分法1. 一维情况下的有限差分法假设要求解的方程为：$$\frac{d^2 y}{dx^2}=-\sin(x)$$首先需要确定求解区域和网格步长。
在本例中，我们将求解区域设为$[0,2\pi]$，网格步长$h=0.01$。
则可以通过以下代码生成网格：```matlabx = 0:0.01:2*pi;```接下来需要构造系数矩阵和右端项向量。
根据上面的公式，系数矩阵应该是一个三对角矩阵，可以通过以下代码生成：```matlabn = length(x)-2;A = spdiags([-ones(n,1), 2*ones(n,1), -ones(n,1)], [-1 0 1], n, n); ```其中`spdiags`函数用于生成一个稀疏矩阵。
基于MATLAB的时域有限差分法_FDTD_的研究
,
j)
−
0.5 ×
Ez
_ inc
(
j)
H
y
(ib
+
1 2
,
j)
=
H
y
(ib
+
1 2
,
j)
+
0.5 ×
E z _ inc
(
j)
当问题空域中不存在散射体时，总场中不存在散射
场，因此在总场区中的总场就是平面波，可用此时平面
波的质量来衡量连接条件的精度。同理，若不存在散射
体就没有散射波的存在，平面波不进入散射场区，这时
n+ 1
I Hy 2
(i
+
1 2
,
j)
=
n−1
I Hy 2
(i
+
1 2
,
j)
+
curl
_
e
H
n+1 y
(i
+
1, 2
j)
=
fi3(i +
1 2
)
⋅
H
n y
(i
+
1, 2
j) −
fi2(i
+
1 ) ⋅ 0.5 ⋅ curl 2
_e+
n+ 1
fj1( j) ⋅ I Hy 2 (i +
1, 2
j)
n+ 1
有三个区域的点的场值需要重新设定：
1．在 j = ja, j = jb 处设置 Dz
Dz
(i,
ja)
=
Dz
(i,
ja)
+
MATLAB编程在FDTD算法中的应用
MATLAB编程在FDTD算法中的应用摘要：文章介绍了时域有限差分（FDTD）法的基本原理，推导了二维TM 波的FDTD表达式。
给出了MATLAB语言编程的步骤和应注意的问题，并结合算例阐述了基于MATLAB仿真的基本方法,最后得出用MATLAB语言对FDTD 算法仿真的几点结论。
关键词：MATLAB；FDTD；编程时域有限差分（FDTD）法是在21世纪60年代由K.S.Yee首先提出并用于求解电磁场散射问题，其主要思路是在空间轴和时间轴上对场量进行离散，并用中心差分代替偏微分，这就将麦克斯韦方程组转化为了差分方程，通过在时间轴和空间轴上采取蛙跳法（leapfrog）逐步推进地求解，最终求得在一定边值与初值条件下的空间场解。
随着计算机技术的发展，FDTD的应用越来越多，对于FDTD算法的编程求解，最常用的程序语言有VC和FORTRUN，而MATLAB 作为一种可视化效果好的软件，在FDTD计算中可视化程度较高，并具有能显示动态场效果的特点。
文章采用FDTD法对高压开关柜内超高频电磁波的辐射和传播特性进行仿真，仿真中将对二维TM电磁波进行FDTD表达式的推导，并结合FDTD算法边界条件的特点，用MATLAB语言进行编程。
1二维TM电磁波FDTD算法1.1算式推导在自由空间中，对于二维TM电磁波，，MAXWELL的两个旋度方程可以分解为：＝-0 （1）＝0（2）０=-（3）构造二维TM波FDTD cell如图1所示：按照FDTD元胞对以上三式中的偏导用中心差商代替，分别可得：=-0 （4）=0 （5）０=-（6）由于方程中出现了半网格和半时间步，为了便于编程，可以将上面差分式改为如下FDTD算式：Hx（i,j,k+1）=Hx（i,j,k）-[EZ（i,j+1,n）-EZ（i,j,n）]（7）Hy（i,j,n+1）=Hy（i,j,n）+[EZ（i+1,j,n）-EZ（i,j,n）]（8）Ez（i,j,n+1）=Ez（i,j,n）+[-] （9）1.2数值稳定性的条件FDTD方法是以一组有限差分方程来替代MAXWELL旋度方程来进行数值计算的方法，在执行形如FDTD算法时，随着时间步长的增长，如何保证该算法的稳定性是一个重要问题。
matlab实现一维fdtd仿真
麦克斯韦方程组为
(2—1—1)
（2—1—2）
各项同性线性介质中的本构关系为
（2—1—3）
在直角坐标系中，（2—1—1）、（2—1—2）式写为
（2—1—4）
以及
（2—1—5）下面我们考虑（2—1—4）、（2—1—5）式的FDTD差分离散。令 代表 或 直角坐标系中某一分量，在时间和空间域中的离散取以下符号表示：
FDTD在电磁研究的多个领域获得广泛应用，其中有：辐射天线的分析，微波器件和导行波结构的研究，散射和雷达截面计算，周期机构分析，电子封装和电磁兼容分析，核电磁脉冲的传播和散射，微光学元器件中光的传播和颜射特性。
二、FDTD模拟电磁波传播的原理
麦克斯韦方程组是支配宏观电磁现象的一组基本方程。这组方程既可以写成微分形式，又可以写成积分形式。FDTD方法是由微分形式的麦克斯韦旋度方程出发进行差分离散从而得到一组时域推进公式。首先给出麦克斯韦旋度方程及其在直角坐标系中的FDTD离散形式，由三维情况到一维情况，即一维电磁波传播的理论基础。
c = 1/sqrt(mu_0*eps_0); %光速
%定义问题空间和参数
domain_size = 1; %空间问题的长度
dx = 1e-3; %单元网格大小
dt = 3e-12; %步长时间的持续
number_of_time_steps = 2000; %迭代次数
nx = round(domain_size/dx); %在问题空间的单元数目
3.2参数设置
程序设计了电流屏在z轴上传播，产生的电场和磁场，此电流屏放置在问题空间的中心，与两个相互平行的完善导体平面平行，整个区间充满空气。激励源在两板之间传播和反射。
无论是简单目标还是复杂目标，在进行FDTD离散时网格尺寸的确定，除了受计算资源的限制不可能取得很小外，还需要考虑以下几个因素：
FDTD(时域有限差分法)算法的Matlab源程序
rect=get(gcf,'Position'); rect(1:2)=[0 0];
M=moviein(nmax/2,gcf,rect);
%*********************************************************************** % BEGIN TIME-STEPPING LOOP %***********************************************************************
%
%***********************************************************************
clear
%*********************************************************************** % Fundamental constants
% Gaussian pulse given by
matlab模拟的电磁学时域有限差分法 pdf
matlab模拟的电磁学时域有限差分法 pdf电磁学时域有限差分法（FDTD）是一种基于数值模拟的电磁场计算方法，它使用有限差分来近似微分方程。
该方法广泛用于电磁学、电波传播、微波技术、光学等领域，以求解电磁场分布和场的辐射、散射等问题。
而在这个领域中，MATLAB是非常流行的工具之一。
本文将围绕“MATLAB模拟的电磁学时域有限差分法”这一主题，从以下几个方面进行阐述：1.时域有限差分法的基础概念在FDTD方法中，将时域中的Maxwell方程组转化为差分形式，使得可以在计算机上进行数值解法。
通过在空间和时间上的离散，可以得到电磁场在时域内的各种分布，进而求得特定情况下的电磁场变化。
2.MATLAB中的FDTD仿真在MATLAB中，我们可以使用PDE工具箱中的电磁学模块来实现FDTD仿真。
通过选择适当的几何形状和边界条件，可以利用该工具箱演示电磁场的传输、反射、折射、透射等现象。
同时，MATLAB中还提供了不同的场分量计算和可视化工具，以便用户可以更好地理解电磁场分布。
3.MATLAB代码实现以下是一些MATLAB代码示例，展示了FDTD模拟的基础实现方法。
代码中的示例模拟了平面波在一个矩形和圆形障碍物上的传播情况。
% 1. Square obstaclegridSize = 200; % Grid sizemaxTime = 600; % Maximum time (in steps)imp0 = 377.0; % Impedance of free spacecourantNumber = 0.5; % Courant numbercdtds = ones(gridSize,gridSize); % Courant number in space% (not variable in this example)Ez = zeros(gridSize, gridSize); % Define EzHy = zeros(gridSize, gridSize); % Define Hy% Simulation loopfor n = 1:maxTime% Update magnetic fieldHy(:,1:end-1) = Hy(:,1:end-1) + ...(Ez(:,2:end) - Ez(:,1:end-1)) .*cdtds(:,1:end-1) / imp0;% Update electric fieldEz(2:end-1,2:end-1) = Ez(2:end-1,2:end-1) + ...(Hy(2:end-1,2:end-1) - Hy(1:end-2,2:end-1)) .* cdtds(2:end-1,2:end-1) .* imp0;end% 2. Circular obstacleradius = 50;xAxis = [-100:99];[X,Y] = meshgrid(xAxis);obstacle = sqrt((X-50).^2 + (Y).^2) < radius;gridSize = length(xAxis); % Grid sizemaxTime = 500; % Maximum time (in steps)imp0 = 377.0; % Impedance of free space courantNumber = 0.5; % Courant numbercdtds = ones(gridSize,gridSize); % Courant number in space% (not variable in this example)Ez = zeros(gridSize, gridSize); % Define EzHy = zeros(gridSize, gridSize); % Define Hy% Simulation loopfor n = 1:maxTime% Update magnetic fieldHy(:,1:end-1) = Hy(:,1:end-1) + ...(Ez(:,2:end) - Ez(:,1:end-1)) .*cdtds(:,1:end-1) / imp0;% Update electric field, with obstacleEz(2:end-1,2:end-1) = Ez(2:end-1,2:end-1) + ...(Hy(2:end-1,2:end-1) - Hy(1:end-2,2:end-1)) .* cdtds(2:end-1,2:end-1) .* imp0;Ez(obstacle) = 0;end以上代码仅供参考，不同条件下的模拟需要适当修改，以便获得特定的模拟结果。
计算电磁学之FDTD算法的MATLAB语言实现
South China Normal University课程设计实验报告课程名称：计算电磁学指导老师：专业班级： 2014级电路与系统姓名：学号：FDTD算法的MATLAB语言实现摘要：时域有限差分(FDTD)算法是K.S.Yee于1966年提出的直接对麦克斯韦方程作差分处理,用来解决电磁脉冲在电磁介质中传播和反射问题的算法。
其基本思想是:FDTD计算域空间节点采用Yee元胞的方法,同时电场和磁场节点空间与时间上都采用交错抽样；把整个计算域划分成包括散射体的总场区以及只有反射波的散射场区,这两个区域是以连接边界相连接,最外边是采用特殊的吸收边界,同时在这两个边界之间有个输出边界,用于近、远场转换;在连接边界上采用连接边界条件加入入射波,从而使得入射波限制在总场区域；在吸收边界上采用吸收边界条件,尽量消除反射波在吸收边界上的非物理性反射波。
本文主要结合FDTD算法边界条件特点,在特定的参数设置下，用MATLAB语言进行编程，在二维自由空间TEz网格中，实现脉冲平面波。
关键词：FDTD；MATLAB；算法1 绪论1.1 课程设计背景与意义20世纪60年代以来，随着计算机技术的发展，一些电磁场的数值计算方法逐步发展起来，并得到广泛应用，其中主要有：属于频域技术的有限元法（FEM）、矩量法(MM)和单矩法等；属于时域技术方面的时域有限差分法(FDTD)、传输线矩阵法(TLM)和时域积分方程法等。
其中FDTD是一种已经获得广泛应用并且有很大发展前景的时域数值计算方法。
时域有限差分（FDTD）方法于1966年由K.S.Yee提出并迅速发展，且获得广泛应用。
K.S.Yee用后来被称作Yee氏网格的空间离散方式，把含时间变量的Maxwell旋度方程转化为差分方程，并成功地模拟了电磁脉冲与理想导体作用的时域响应。
但是由于当时理论的不成熟和计算机软硬件条件的限制，该方法并未得到相应的发展。
20世纪80年代中期以后，随着上述两个条件限制的逐步解除，FDTD便凭借其特有的优势得以迅速发展。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
% 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%%% 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 spacemuz=4.0*pi*1.0e-7; %permeability of free spaceepsz=1.0/(cc*cc*muz); %permittivity of free space%*********************************************************************** % Grid parameters%***********************************************************************ie=50; %number of grid cells in x-directionje=24; %number of grid cells in y-directionke=10; %number of grid cells in z-directionib=ie+1;jb=je+1;kb=ke+1;is=26; %location of z-directed current sourcejs=13; %location of z-directed current sourcekobs=5;dx=0.002; %space increment of cubic latticedt=dx/(2.0*cc); %time stepnmax=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');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)+...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;%***********************************************************************%*********************************************************************** endmovie(gcf,M,0,10,rect);。