时域有限差分法仿真一维TE波在分裂场完全匹配层【含源码】
时域有限差分法
引言
时域有限差分法的软件
• • FDTDA，三维时域有限差分法的软件，源程序用FORTRAN语言 编写（1993年） XFDTD，具有多种功能，包含有瞬态近—远场外推，亚网格技 术，介质可以是有耗介质、磁化铁氧体，可用以分析生物体对电 磁波的吸收特性（SAR），螺旋及微带天线，天线阻抗的频率特 性，移动电话场强分布，细导线及复杂物体电磁散射和RCS （1996年） EMA3D，分析核电磁脉冲（NEMP）及雷电耦合，高功率微波， 宽带RCS，天线，屏蔽特性，印刷电路板的电磁兼容。软件具有 多种边界条件，亚网格剖分，适用于有耗介质、平面波源及电压 电流源（1997年）
其中E为电场强度，单位为伏特/米 D为电通量密度，单位为库仑/米2 H为磁场强度，单位为安培/米 B为磁通量密度，单位为韦伯/米2 J为电流密度，单位为安培/米2 Jm为磁流密度，单位为伏特/米2
麦克斯韦方程
各向同性线性介质中的本构关系为
B = μH
D = εE
其中 ε 为介质介电系数，单位为法拉/米 μ 为磁导系数，单位为亨利/米 σ 为电导率，单位为西门子/米 σ m 为导磁率，单位为欧姆/米 σ 和 σ m 分别为介质的电损耗和磁损耗 在真空中， σ = 0 ， σ = 0 ， ε = ε = 8.85 ×10−12 法拉/米
引言
时域有限差分法的产生与发展
• 1989年，Britt首次给出时域远场的结果，但未给出外 推的具体方法 • 1989年，Larson、Perlik和Taflove等人提出研究适用于 时域有限差分法的专用计算机，以便用于计算电磁波 与电大尺寸物体的相互作用 • 1990年，Maloney等人用柱坐标系下的时域有限差分法 分析了柱状和锥状天线位于理想导体平面上的辐射， 得到宽带天线的输入阻抗及瞬态辐射场的直观可视化 显示
时域有限差分法对平面TE波的MATLAB仿真
时域有限差分法对平面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个电场(磁场)分量环绕。
时域有限差分法
j 2 c 2 jk 2
(1-8)
（1-9） 可见，相速与频率无关，称为非色散。非色散意味着对于具有任意 调制的包络或脉冲形状的波传播任意距离后波形保持不变。进一步 由(1-8)可以得到群速关系 d vg c (1-10) dk 这种情况下，群速也是与频率无关。
1.2 数值色散关系(2)
(200空间格)，数值模拟传播了199.378格，相位误差为11.1960，也 减少了4倍。误差减少了4倍反映了差分算法是二阶精度的。
1.3 数值相速（2）
• 情况1：非常细网格 t 0, x 0
2 cos x 1 x 2 , 当 x 0 ，数值色散关系(1-12)变为 根据 ~ 2 2 2 t k x ct 1 1 1 1 2 2 x ~ ~ v 。 2 2 ~2 即， c k ，最后得 k k ，于是有 v p p c
定义数值相速为 （1-14） • 情况1 非常细网格 利用正弦函数的一阶Taylor展开，可得 ~ 2 2~ （1-15） c t k x c k ~ v g x0 x t t 0 所以，群速与相速一样，在细网格条件下趋近精确解。这证明了 当空间步长和时间步长趋于零时，数值解变得精确。 • 情况2 魔时间步 将魔时间步条件和波数代入(1-14)，得 2 （1-16） c t sin c ct ~ vg c v g ct sint 再次验证了魔时间步下数值解等于精确解。
时域有限差分法
第1讲 一维标量波动方程
引言（1）
• 1966年,K.S. Yee（美籍香港人）首先提出了FiniteDifference Time-Domain Method,并用于柱形金属柱 电磁散射分析。由于当时计算机技术还比较落后，这 一方法并未引起重视。
时域有限差分法
时域有限差分法
时域有限差分法的基本思想是用中心差商代替场量对时间和空间的一阶偏微商, 通过在时域的递推模拟波的传播过程, 从而得出场分布。
它最早由K.S.Yee 于1966 年提出，在此之后的20 年内，其研究进展缓慢，只是在电磁散射、电磁兼容领域有一些初步的应用。
自80 年代末，时域有限差分法成为电磁场数值计算的重要方法之一。
在声学数值计算中，时域有限差分法已应用于水声学、噪声控制及室内声学等方面的数值模拟。
时域有限差分法(Finite difference time domainmethod，FDTD)直接离散时域波动方程，不需要任何形式的导出方程，故不会因为数学模型而限制其应用范围。
它的差分格式中包含有介质的参量，只须赋予各网格相应的参量，就能模拟各种复杂的结构，这是时域有限差分法的一个突出优点。
另外，由于时域有限差分法采用步进法进行计算，故能很容易地实现各种复杂时域宽带信号的模拟，而且可以非常方便地获得空间某一点的时域信号波形。
时域有限差分
z
z
Ezn (i, j 1/ 2, k 1/ 2) y
Ezn (i, j 1, k 1/ 2) Ezn (i, j, k 1/ 2) y
偏微分方程离散
作业报告9，推导其他5个偏微分方程 的差分格式。
整理后得
Hx
n
1 2
n 1 1 1 1 i, j , k H x 2 i, j , k 2 2 2 2
H x t
n
i , j 1/2, k 1/2
1 E y Ez z y i , j 1/2,k 1/2
Ezn
n
Eyn Hxn
z y
(i,j+1/2,k+1/2)
方程离散点
Ezn
Eyn
偏微分方程离散
H x t
n i , j 1/2, k 1/2
t
E n 1 x, y, z E n x, y, z
Eyn Hxn
n
1 E y Ez z y i , j 1/2,k 1/2
Ezn
方程离散点
z y
(i,j+1/2,k+1/2)
Ezn
写成取样离散格式
Eyn
n H xn (i, j 1/ 2, k 1/ 2) 1 E y (i, j 1/ 2, k 1/ 2) Ezn (i, j 1/ 2, k 1/ 2) t z y
Ampère’s Equation
D x, y, z, t H x, y , z , t t
Maxwell方程组
B x, y, z, t E x, y, z, t t
第二章 时域有限差分法_II-一维FDTD
2017/5/2
4
2017/5/2
5
进一步，得到迭代公式
E
n 1/2 x
k E
n 1/2 x
t n n H k 1/ 2 H k y y k 1/ 2 0 x t n 1/2 n 1/2 E k 1 E k x x 0 x
n
n n 1 D 1
0
t n E 0
关于频域依赖媒质的迭代方程及代码 D(k) = D(k)+ eta *( H(k-1)-H(k) ); E(m)=ga(m)*(D(m)-I(m)); I(m)=I(m)+gb(m)*E(m); H(j)=H(j)+eta*(E(j)-E(j+1))； ga(m) = 1/(epsilon+(sigma*dt/epszero));
ct n n 2 r n 1/2 0.5 n 1/2 H y k 1/ 2 H y k 1/ 2 Ex k E k ct x c t 1 r 1 2 r 2 r
t
x 2c0
This value of η motivates Sullivan's choice of boundary conditions at the left boundary given by
n Ex 1 Exn2 2
Similarly, for the right boundary conditions we use
上面两方程的迭代方程
c t 1 取 x 2
ct n n 2 r n 1/2 0.5 n 1/2 H y k 1/ 2 H y k 1/ 2 Ex k Ex k ct c t 1 r 1 2 r 2 r
时域有限差分法对平面TE波的MATLAB仿真
时域有限差分法对平面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算法）的基本原理及仿真
时域有限差分法（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 中，空间上连续分布的电磁场物理量离散的空间排布如图所示。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
时域有限差分法仿真一维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。