下载文档原格式
时域有限差分法对平面TE波的MATLAB仿真姓名：王云璐学号：2011302021（西北工业大学电子信息学院08041103，陕西西安，710072）摘要：本文分析FDTD算法的基本原理及两种典型边界条件的算法特点，通过MATLAB程序对TE波进行了仿真，模拟了高斯制下完全匹配层中磁场分量瞬态分布。
得到了相应的磁场幅值效果图。
最后得出用Matlab语言对FDTD算法编程的几点结论。
关键词：FDTD；MATLAB；PML；平面TE波1引言电磁场的有限差分一般是在时域上进行的，随着计算机技术的发展和应用，近年来，时域计算方法越来越受到重视。
时域有限差分法具有简单、结果直观、网格剖分简单等优点。
近些年FDTD发展的十分迅速，在许多方面都有重要应用，包括天线设计，微波电路设计，电磁兼容分析，电磁散射计算，光子学应用等等。
时域有限差分(FDTD)算法是K.S.Yee于1966年提出的直接对麦克斯韦方程作差分处理，用来解决电磁脉冲在电磁介质中传播和反射问题的算法。
基本思想是：FDTD计算域空间节点采用Yee元胞的方法，同时电场和磁场节点空间与时间上都采用交错抽样；把整个计算域划分成包括散射体的总场区以及只有反射波的散射场区，这两个区域是以连接边界相连接，最外边是采用特殊的吸收边界，同时在这两个边界之间有个输出边界，用于近、远场转换；在连接边界上采用连接边界条件加入入射波，从而使得入射波限制在总场区域；在吸收边界上采用吸收边界条件，尽量消除反射波在吸收边界上的非物理性反射波。
本文主要简述了FDTD算法的基本原理，解的稳定性以及边界条件特点，并用Matlab语言进行对平面TE波进行了编程计算。
2FDTD基本原理时域有限差分法的主要思想是把Maxwell方程在空间、时间上离散化，用差分方程代替一阶偏微分方程，求解差分方程组，从而得出各网格单元的场值。
FDTD空间网格单元上电场和磁场各分量的分布如图1所示。
图1.Yee 氏网格及其电磁场分量分布电场和磁场被交叉放置，电场分量位于网格单元每条棱的中心，磁场分量位于网格单元每个面的中心，每个磁场(电场)分量都有4个电场(磁场)分量环绕。
时域有限差分法
时域有限差分法的基本思想是用中心差商代替场量对时间和空间的一阶偏微商, 通过在时域的递推模拟波的传播过程, 从而得出场分布。
它最早由K.S.Yee 于1966 年提出，在此之后的20 年内，其研究进展缓慢，只是在电磁散射、电磁兼容领域有一些初步的应用。
自80 年代末，时域有限差分法成为电磁场数值计算的重要方法之一。
在声学数值计算中，时域有限差分法已应用于水声学、噪声控制及室内声学等方面的数值模拟。
时域有限差分法(Finite difference time domainmethod，FDTD)直接离散时域波动方程，不需要任何形式的导出方程，故不会因为数学模型而限制其应用范围。
它的差分格式中包含有介质的参量，只须赋予各网格相应的参量，就能模拟各种复杂的结构，这是时域有限差分法的一个突出优点。
另外，由于时域有限差分法采用步进法进行计算，故能很容易地实现各种复杂时域宽带信号的模拟，而且可以非常方便地获得空间某一点的时域信号波形。
时域有限差分法对平面TE波的MATLAB仿真摘要时域有限差分法是由有限差分法发展出来的数值计算方法。
自1966年Yee 在其论文中首次提出时域有限差分以来，时域有限差分法在电磁研究领域得到了广泛的应用。
主要有分析辐射条线、微波器件和导行波结构的研究、散射和雷达截面计算、分析周期结构、电子封装和电磁兼容的分析、核电磁脉冲的传播和散射以及在地面的反射及对电缆传输线的干扰、微光学元器件中光的传播和衍射特性等等。
由于电磁场是以场的形态存在的物质，具有独特的研究方法，采取重叠的研究方法是其重要的特点，即只有理论分析、测量、计算机模拟的结果相互佐证，才可以认为是获得了正确可信的结论。
时域有限差分法就是实现直接对电磁工程问题进行计算机模拟的基本方法。
在近年的研究电磁问题中，许多学者对时域脉冲源的传播和响应进行了大量的研究，主要是描述物体在瞬态电磁源作用下的理论。
另外，对于物体的电特性，理论上具有几乎所有的频率成分，但实际上，只有有限的频带内的频率成分在区主要作用。
文中主要谈到了关于高斯制下完全匹配层的差分公式的问题，通过MATLAB 程序对TE波进行了仿真，模拟了高斯制下完全匹配层中磁场分量瞬态分布。
得到了相应的磁场幅值效果图。
关键词：时域有限差分完全匹配层MATLAB 磁场幅值效果图目录摘要 (1)目录 (3)第一章绪论 (4)1.1 课题背景与意义 (4)1.2 时域有限差分法的发展与应用 (4)2.1 Maxwell方程和Yee氏算法 (7)2.2 FDTD的基本差分方程 (9)2.3 时域有限差分法相关技术 (11)2.3.1 数值稳定性问题 (11)2.3.2 数值色散 (12)2.3.3 离散网格的确定 (13)2.4 吸收边界条件 (13)2.4.1 一阶和二阶近似吸收边界条件 (14)2.4.2 二维棱边及角顶点的处理 (17)2.4.3 完全匹配层 (19)2.5 FDTD计算所需时间步的估计 (23)第三章MATLAB的仿真的程序及模拟 (25)3.1 MATLAB程序及相应说明 (25)3.2 出图及结果 (28)3.2.1程序部分 (28)3.2.2 所出的效果图 (29)第四章结论 (31)参考文献 (32)第一章绪论1.1 课题背景与意义20世纪60年代以来，随着计算机技术的发展，一些电磁场的数值计算方法逐步发展起来，并得到广泛应用，其中主要有：属于频域技术的有限元法（FEM）、矩量法(MM)和单矩法等；属于时域技术方面的时域有限差分法(FDTD)、传输线矩阵法(TLM)和时域积分方程法等。
时域有限差分法（FDTD 算法）时域有限差分法是1966年K.S.Yee 发表在AP 上的一篇论文建立起来的，后被称为Yee 网格空间离散方式。
这种方法通过将Maxwell 旋度方程转化为有限差分式而直接在时域求解, 通过建立时间离散的递进序列, 在相互交织的网格空间中交替计算电场和磁场。
FDTD 算法的基本思想是把带时间变量的Maxwell 旋度方程转化为差分形式，模拟出电子脉冲和理想导体作用的时域响应。
需要考虑的三点是差分格式、解的稳定性、吸收边界条件。
有限差分通常采用的步骤是：采用一定的网格划分方式离散化场域；对场内的偏微分方程及各种边界条件进行差分离散化处理，建立差分格式，得到差分方程组；结合选定的代数方程组的解法，编制程序，求边值问题的数值解。
1.FDTD 的基本原理FDTD 方法由Maxwell 旋度方程的微分形式出发，利用二阶精度的中心差分近似，直接将微分运算转换为差分运算，这样达到了在一定体积内和一段时间上对连续电磁场数据的抽样压缩。
Maxwell 方程的旋度方程组为：E E HtHH Emt（1）在直角坐标系中，（1）式可化为如下六个标量方程：zz x y y y z xx x yzE tE yH xH E t E x H z HE t E zHy H ，zmz x y y m y z x x mx y z H tH yE xE 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(),,,(xO xk j iF k j iF x t z y x F nnxi x 2),21,(),21,(),,,(yO yk j i F k ji F yt z y x F nnyj y 2)21,,()21,,(),,,(zO zk j i F k j i F zt z y x F nnzk z对时间离散：22121),,(),,(),,,(tO tk j i Fk j i Ftt z y x F n n tn t （4）Yee 把空间任一网格上的E 和H 的六个分量，如下图放置:oyxzEyHzExEzHxEyEyEzEx HyEzEx 图1 Yee 氏网格及其电磁场分量分布在FDTD 中，空间上连续分布的电磁场物理量离散的空间排布如图所示。
时域有限差分法的Matlab仿真
张通;孙晶
【期刊名称】《科技视界》
【年(卷),期】2017(000)003
【摘要】文章介绍了时域有限差分法的基本原理,利用matlab仿真,实现了用时域有限差分程序来计算二维问题空间中的电场分布.
【总页数】3页(P18-19,4)
【作者】张通;孙晶
【作者单位】吉首大学物理与机电工程学院,湖南吉首416000;吉首大学物理与机电工程学院,湖南吉首416000
【正文语种】中文
【相关文献】
1.基于高阶时域有限差分法与改进节点分析法混合求解复杂传输线网络瞬态响应[J], 王为;覃宇建;刘培国;周东明
2.时域有限差分法的Matlab仿真 [J], 郭春波
3.时域有限差分法和基于频域的有限元法模拟仿真在通r信车的天线布局设计的应用 [J], 李宏
4.时域有限差分法的Matlab仿真 [J], 黄明红;汪清泉
5.一种简化的分析光纤布拉格光栅的时域有限差分束传输法 [J], 戴劲草;王云明;张明德;孙小菡
因版权原因，仅展示原文概要，查看原文内容请购买。
一维有限时域差分方法(FDTD)计算中参数的选择
陈义万;李文兵;杜海霞;彭波勇
【期刊名称】《教育教学论坛》
【年(卷),期】2013(000)003
【摘要】在用一维时域有限差分方法计算波的传播时,为了防止出现发散的结果,要求时间步长△t＜Δx/v,但具体取多大的值,没有定论.通过编程试验,时间步长取△t＜0.75Δx/c,而时间步数Nt大于波从波源传到观察点所需要的时间,又小于反射波到达观察点的时间.而且计算的时间格点和位置格点全部都可以取整数.
【总页数】2页(P155-156)
【作者】陈义万;李文兵;杜海霞;彭波勇
【作者单位】湖北工业大学理学院,湖北武汉430068;湖北工业大学理学院,湖北武汉430068;第二炮兵工程大学,陕西西安710025;第二炮兵工程大学,陕西西安710025
【正文语种】中文
【中图分类】G642.4
【相关文献】
1.时域有限差分方法分析微带天线散射参数和方向图
2.快速时域有限差分方法计算矩形缺陷接地结构
3.区域分解时域有限差分方法(DD-FDTD)及其在散射问题中的应用
4.一种节约内存的局部一维减缩时域有限差分方法
5.利用时域有限差分方法计算雷电水平电场
因版权原因，仅展示原文概要，查看原文内容请购买。
3 时域有限差分法（FDTD ）1966年K.S. Yee 发表了时域有限差分法(Finite Difference -Time Domain ，简记FDTD)的奠基性论文[1]，之后在很长一段时间内，这一思想没有引起电磁理论界的足够重视。
直到七十年代末八十年代初，在A. Taflove [2]、K.S. Kunz [3]和R. Holland [4]等学者的推进下，这一方法才逐渐走向成熟并得到广泛的研究和应用。
时域有限差分法的原理非常简单，就是直接将时域Maxwell 方程组的两个旋度方程中关于空间变量和时间变量的偏导数用差商近似，从而转换为离散网格节点上的时域有限差分方程。
加入时域脉冲激励后，在时间上迭代就可直观地模拟出脉冲在求解区域上传播、反射和散射的过程，进而采用FFT 将时域响应变换到频域就可获得所希望的各种电参数，如无源电路的散射参数、天线的辐射方向图和输入阻抗、散射体的雷达散射截面(RCS)等。
随着FDTD 方法的迅猛发展，新的处理方法和技术不断涌现。
其中，子网格模型技术是用子网格或细网格划分薄片、裂缝和导线，其余部分用粗网格进行划分，以便在不显著增加计算时间的基础上提高计算精度；非正交和广义正交曲线网格技术适应于各种结构形状，可以模拟各种复杂的结构；非均匀正交网格技术在复杂结构区域或在场量快变化区域采用细网格，而在其它地方用粗网格，可以兼顾计算时间、存储量和计算精度；回路积分法从积分形式的Faraday 定律和推广的Ampere 定律出发导出回路积分表示的差分格式，使之适用于任意形状的网格结构；外推技术从前面有限时间步的瞬时响应外推以后瞬时响应以大量节省计算时间；网格压缩模型技术用于导波结构分析，通过解析处理，可将传播方向的网格压缩为零。
此外还有，超吸收边界条件技术、色散吸收边界条件技术、完全匹配层吸收边界条件技术、多分辨率技术、伪谱技术、及混合显-隐格式算法等。
新方法与技术的发展迅速扩大了时域有限差分法的应用范围，该方法不仅在目标电磁散射问题，而且在电磁兼容预测、微波电路分析、天线辐射特性计算和生物电磁学研究等方面中都获得了广泛的应用。
第十一章 时域有限差分方法自从1966年K. S. Yee 创建时域有限差分法 (Finite Difference Time Domain ，简称FDTD) 以来[1]，已经发展成为一种理论完整、应用广泛的数值方法，并且与矩量法和有限元法一起奠定了计算电磁学的基础。
本章将介绍时域有限差分的基本理论，数值模拟技术，若干相关的专题以及工程实例。
11-1 差分的基本概念时域有限差分法是对微分形式的Maxwell 方程进行差分求解的技术。
在详述其之前，首先简单回顾差分的基本概念。
已知分段连续函数()f x 在位置x 处的增量可表示为()()()f x f x x f x ∆=+∆- (11-1-1)其差商为()()()f x f x x f x xx∆+∆-=∆∆ (11-1-2)当x ∆→0时，()f x 的导数定义为差商的极限，即()()()()'limlimx x f x f x x f x f x xx∆→∆→∆+∆-==∆∆ (11-1-3)当x ∆足够小时，()f x 的导数可以近似为d d f fx x∆≈∆ (11-1-4) 根据导数取值位置的不同，差分格式分为前向差分、后向差分和中心差分。
前向差分定义为()()x f x x f x fx x+∆-∆=∆∆ (11-1-5) 后向差分定义为()()x f x f x x fx x--∆∆=∆∆ (11-1-6) 中心差分定义为()()22x f x x f x x fx x+∆--∆∆=∆∆ (11-1-7) 将()f x x +∆在点x 处展开为Taylor 级数，得()()()()()232323d d d 11d 2!d 3!d f x f x f x f x x f x x x x x x x+∆=+∆+∆+∆ (11-1-8)()()()()()232323d d d 11d 2!d 3!d f x f x f x f x x f x x x x x x x-∆=-∆+∆-∆ (11-1-9)将方程 (11-1-8) 和 (11-1-9) 代入 (11-1-5) ~ (11-1-7)后可以发现，前向和后向差分具有一阶精度，中心差分具有二阶精度。
专利名称：一种基于时域有限差分法的半确定性信道仿真方法专利类型：发明专利
发明人：周健义,董云扬,于志强,洪伟
申请号：CN202111431082.1
申请日：20211129
公开号：CN114117860A
公开日：
20220301
专利内容由知识产权出版社提供
摘要：本发明公开了一种基于时域有限差分法的半确定性信道仿真方法。
该方法采用时域有限差分法(FDTD)对全空间进行电磁仿真来获取收发之间的响应，并对散射物进行微调来引起信道的波动，从而获得信道的多径统计特性。
首先，在设置空间媒质分布的基础上，采用FDTD通过更新方程进行磁场和电场的不停迭代，直到达到预设的时间总步数后停止；然后，取出激励源信号s(t)和观测点的信号r(t)，将二者进行频率搬移并滤波取得信号的包络信号或称等效低通信号sl(t)和rl(t)，之后求观测信号rl(t)与激励信号sl(t)的相关，由相关峰的峰值和时延可得此信道每条路径的衰减和时延；最后，将计算空间内的散射物随机变动，多次之后求取各条路径上衰减与时延的统计特性。
本发明计算精确，适用性广泛，也可用于复杂介质下的信道计算。
又加入了散射物的随机性，得到了半确定性信道模型，体现出了信道的随机性特征。
申请人：东南大学
地址：210096 江苏省南京市玄武区四牌楼2号
国籍：CN
代理机构：南京瑞弘专利商标事务所(普通合伙)
代理人：沈廉
更多信息请下载全文后查看。
时域有限差分法仿真一维TE 波在分裂场完全匹配层吸收边界条件下的传输一、时域有限差分法 (FDTD, Finite-Difference Time-Domain)FDTD 是1966年K.S.Yee 发表在AP 上的一篇论文建立起来的，后被称为Yee 网格空间离散方式核心思想是把带时间变量的Maxwell 旋度方程转化为差分形式，模拟出电子脉冲和理想导体作用的时域响应号称目前计算电磁学界最受关注，最时髦的算法，但还在发展完善之中国外已有多种基于FDTD 算法的电磁场计算的软件:XFDTD 等。
二、差分算法x h ∆= ()()f f x h f x ∆=+-0()()()()lim ()()()()=2x df f x f x f x h f x dx x x hf x f x h hf x h f x h h∆→∆∆+-=≈=∆∆--=+--前向差分 后向差分中心差分 222()1()()()2!df x d f x f x h f x h h dx dx +=+++222()1()()()2!df x d f x f x h f x h h dx dx -=-++333()2()()()23!df x d f x f x h f x h h h dx dx +--=++三、时域有限差分算法步骤(1)采用一定的网格划分方式离散化场域；(2)对场内的偏微分方程及各种边界条件进行差分离散化处理，建立差分格式，得到差分方程组；(3)结合选定的代数方程组的解法，编制程序，求边值问题的数值解。
四、吸收边界条件1、问题的提出在电磁场的辐射和散射问题中，边界总是开放的，电磁场占据无限大空间，而计算机内存是有限的，所以只能模拟有限空间。
即：时域有限差分网格将在某处被截断。
这要求在网格截断处不能引起波的明显反射，因而对向外传播的波而言，就像在无限大的空间传播一样，一种行之有效的方法是在截断处设置一种吸收边界条件。
使传播到截断出的波被边界吸收而不产生反射。
2、良匹配层基本原理1994年Beernger 提出一种新颖的由吸收媒质构成的良匹配层(PML)的概念,这种人工设计的良匹配层由有耗,导电,导磁媒质组成,可吸收任意入射角,任意频率,任意偏振态的入射电磁波,且透射波以原来波速传播,且有相同的特征阻抗,在垂直入射方向衰减。
根据推导，只要这种媒质的结构参数满足一定条件,这种媒质就可以起到完全吸收入射波的作用。
五、计算所用到的电磁场公式根据麦克斯韦方程组和时域差分法，得到以下公式：()()()()1/21/21/21/2nn n n y y x x H k H k E k E k t x +-+---=∆∆()()()()11/21/201/21/211n n n n y y x x H k H k E k E k txμ++++--+-=-⋅∆∆12=1/21/21/2()()[(1/2)(1/2)]n n n nx x y y rE k E k H k H k ε++=++--1/21/21(1/2)=(1/2)-[(1)()]2n nn n y y x x H k H k E k E k +++++-六、仿真结果timeE -f i e l di n p u tfrequencyp o i n t s /w a v e l e n g t h-4time [ns]R e f l e c t e d f i e l dfrequency [GHz]R e f l e c t i o n i n d Btime [ns]T r a n s m i t t e d f i e l d-3frequency [GHz]T r a n s m i s s i o n i n d B七、源代码%*********************************************************************** % 时域有限差分法仿真一维TE波在分裂场完全匹配层吸收边界条件下的传输%*********************************************************************** %%% Eref Ein Erans% PEC PML 介质PML PEC% +-----+---+------+---------------------+-----+------+-->% z=0 z=Lz% r=1 Nref Nin Ntrans r=Nz+1%%*********************************************************************** % 物理常量mu0 = 4e-7 * pi; % 磁导率c0 = 299792458; % 光速eps0 = 1/c0^2/mu0; % 电导率eta0 = sqrt(mu0/eps0); % 波导率GHz = 10^9;mm = 10^(-3);%*********************************************************************** % 参数初始化%*********************************************************************** Lz = 1; % 长度单位：米Tmax = 4*Lz/c0; % 时间最大值Dz = 1*mm; % 空间网格尺寸R = 1/sqrt(3); %每一单位时间计算的网格单元数freq = 5*GHz; % 激励源的中心频率fmax = 9*GHz; % 画图的最大频率fmin = 1*GHz;Nz = round(Lz/Dz); % z轴上的单元格数量（round 四舍五入）n_pict = round(Nz/20); % 每n_pict步画一次图linew = 2; % 画图的线宽Dt = Dz*R/c0; % 时间步长Nt = round(Tmax/Dt); % 时间步长数量Nfft = max(2^(round(log2(Nt)+1)),2*1024); %做傅里叶变换点的数量z = linspace(0,Lz,Nz+1); % Z轴坐标t = Dt*[0:Nt-1]'; % 时间lambda = c0/freq; % 中心波长的激励源ppw = lambda/Dz % 中心频率处每一段波长的采样点omega = 2.0*pi*freq; % 角频率Nabc = 1*10; % PML边界左边和右边的层数Nin = 4; % 表面散射场Nref = 2; % 由近到远场表面Ntrans = Nz-2;Labc = Nabc*Dz; % PML边界长度%*********************************************************************** % 材料参数%*********************************************************************** Nmedia=3; % 材料数量media_par = [1.0 0.0; % 真空：电容率和电导率3.4 1*10^(-12); %3.0 0];objects = 0; % 介质的数量obj_z = Dz/2+[0.500 0.504; % Z1和Z2之间的介质10.624 0.628;0.24 0.26];obj_m = [2 2 1]; % 介质的材料obj_c = 'rrb'; % 介质的颜色epsr = ones(Nz-1,1);sig = zeros(Nz-1,1);obj_ind = [];obj_ind(:,1) = ceil(obj_z(:,1)/Dz);obj_ind(:,2) = floor(obj_z(:,2)/Dz);obj_frac(:,1) = obj_ind(:,1)-obj_z(:,1)/Dz;obj_frac(:,2) = -obj_ind(:,2)+obj_z(:,2)/Dz;for n=1:objectsepsr(obj_ind(n,1):obj_ind(n,2)) = media_par(obj_m(n),1);sig(obj_ind(n,1):obj_ind(n,2)) = media_par(obj_m(n),2);end%*********************************************************************** % 分配场矢量%*********************************************************************** Ex = zeros(Nz+1 ,1);Hy = zeros(Nz,1);Eref = zeros(Nt,1);Etrans = zeros(Nt,1);% 吸收边界条件E1abc = zeros(Nabc+1,1); % 左端E2abc = zeros(Nabc+1,1); % 右端H1abc = zeros(Nabc,1);H2abc = zeros(Nabc,1);zpmlH = linspace(0,Labc,Nabc)'/Labc;zpmlE = linspace(Dz/2,Labc-Dz/2,Nabc-1)'/Labc;abc_pow = 3;sigma_max = -(abc_pow+1)*log(10^(-4))/(2*eta0*Labc);sigmaH2 = sigma_max*zpmlH.^abc_pow*mu0/eps0;sigmaE2 = sigma_max*zpmlE.^abc_pow;sigmaE1 = sigmaE2((Nabc-1):-1:1);sigmaH1 = sigmaH2(Nabc:-1:1);figure(4)plot(sigmaH2)%***********************************************************************% 波激励%***********************************************************************Hdelay = Dt/2-Dz/2/c0;% 从左侧加入射脉冲tau = 80.0e-12;delay = 2.1*tau;Ein = sin(omega*(t-delay)).*exp(-((t-delay).^2/tau^2));Hin = 1/eta0*sin(omega*(t-delay+Hdelay)).*exp(-((t-delay+Hdelay).^2/tau^2));% 谐波从左侧入射的时间Einmax = max(abs(Ein));%***********************************************************************% 时间步进%***********************************************************************for n = 1:Nt;% 从-Dt/2到Dt/2更新磁场Hy = Hy - (Dt/mu0/Dz) * diff(Ex,1); % 主网格Hy(Nin) = Hy(Nin) - (Dt/mu0/Dz) * Ein(n);H1abc = (H1abc.*(mu0/Dt-sigmaH1/2) - diff(E1abc)/Dz) ./ (mu0/Dt+sigmaH1/2); % 从ABC向左H2abc = (H2abc.*(mu0/Dt-sigmaH2/2) - diff(E2abc)/Dz) ./ (mu0/Dt+sigmaH2/2); % 从ABC向右% 从0到Dt更新EEx(2:Nz) = (Ex(2:Nz).*(eps0*epsr/Dt-sig/2) - diff(Hy,1)/Dz)./(eps0*epsr/Dt+sig/2); % 除去边界的主网格Ex(Nin) = Ex(Nin) - Hin(n)/(Dz*eps0/Dt);Ex(1) = Ex(1) - (Dt /eps0) * (Hy(1)-H1abc(Nabc))/Dz; %左边ABC和主网格之间的边界E1abc(Nabc+1) = Ex(1); % 补充一点，以简化方案E1abc(2:Nabc) = (E1abc(2:Nabc).*(eps0/Dt-sigmaE1/2) - diff(H1abc,1)/Dz) ./ (eps0/Dt+sigmaE1/2);Ex(Nz+1) = Ex(Nz+1) - (Dt /eps0) * (H2abc(1)-Hy(Nz))/Dz;E2abc(1) = Ex(Nz+1);E2abc(2:Nabc) = (E2abc(2:Nabc).*(eps0/Dt-sigmaE2/2) - diff(H2abc,1)/Dz) ./ (eps0/Dt+sigmaE2/2);% 在选定点采样电场if mod(n,n_pict) == 0figure(1);subplot(2,1,1);plot(z,Ex,'LineWidth',1);title(['time is ',num2str(n*Dt*10^9),'ns'])axis tight;ax=axis;axis([ax(1:2) -1.2*Einmax 1.2*Einmax]);ax=axis;ylabel('E-field');for p=1:objectsline([obj_z(p,1) obj_z(p,1)], ax(3:4),'color',obj_c(p));line([obj_z(p,2) obj_z(p,2)], ax(3:4),'color',obj_c(p));line([obj_z(p,1) obj_z(p,2)], [ax(3) ax(3)],'color',obj_c(p),'LineWidth',4);line([obj_z(p,1) obj_z(p,2)], [ax(4) ax(4)],'color',obj_c(p),'LineWidth',4);endsubplot(2,1,2);plot([E1abc(:)' Ex' E2abc(:)'],'LineWidth',1)title(['time is ',num2str(n*Dt*10^9),'ns'])ax = axis;line([Nabc Nabc], ax(3:4),'color','r');line(Nin+[Nabc Nabc], ax(3:4),'color','b');line(Ntrans+2+[Nabc Nabc], ax(3:4),'color','b');line([Nz+2+Nabc Nz+2+Nabc], ax(3:4),'color','r');axis tight;ylabel('E-field');%plot(Eabc(2:Nabc,1:2));pause(0.01);endEref(n) = Ex(Nref);Etrans(n) = Ex(Ntrans);end%***********************************************************************% 后处理%***********************************************************************Ehin = fft(Ein,Nfft)/Nt*2;Df = 1/(Dt*Nfft);freqN = round(fmax/Df);freqN1 = round(fmin/Df);f = linspace(Df*freqN1,Df*freqN,freqN-freqN1+1)/10^9;t = t*10^9; % 纳米秒（ns）figure(2);% 入射波clf;subplot(2,3,1);plot(t,Ein,'r','LineWidth',linew)axis tight;xlabel('time');ylabel('E-field');subplot(2,3,4);warning off;A=plotyy(f,abs(Ehin(freqN1:freqN)),f,c0./f/Dz/10^9);warning on;set(A(2),'yLim',[0 40],'yTick',[0:5:40],'xlim',[fmin fmax]/10^9,'xGrid','on','yGrid','on','LineWidth',linew); set(A(1),'xlim',[fmin fmax]/10^9,'LineWidth',linew);set(get(A(2),'Ylabel'),'String','points/wavelength');set(get(A(1),'Ylabel'),'String','input');set(get(A(1),'Xlabel'),'String','frequency');% 反射subplot(2,3,2);plot(t,Eref,'LineWidth',linew); axis tight;xlabel('time [ns]');ylabel('Reflected field');subplot(2,3,5);Ehref = fft(Eref,Nfft)/Nt*2;plot(f,20*log10(abs(Ehref(freqN1:freqN)./Ehin(freqN1:freqN))),'LineWidth',linew); axis tight;ax=axis;axis([fmin/10^9 fmax/10^9 ax(3:4)]);axis tight;xlabel('frequency [GHz]');ylabel('Reflection in dB');% 传播subplot(2,3,3);plot(t,Etrans,'LineWidth',linew); axis tight;xlabel('time [ns]');ylabel('Transmitted field');subplot(2,3,6);Ehtrans = fft(Etrans,Nfft)/Nt*2;plot(f,20*log10(abs(Ehtrans(freqN1:freqN)./Ehin(freqN1:freqN))),'LineWidth',linew); axis tight;ax=axis;axis([fmin/10^9 fmax/10^9 ax(3:4)]);xlabel('frequency [GHz]');ylabel('Transmission in dB');tr=interp1(f,abs(Ehtrans(freqN1:freqN)./Ehin(freqN1:freqN)),5) re=interp1(f,abs(Ehref(freqN1:freqN)./Ehin(freqN1:freqN)),5) tr^2+re^2。
