Last report for Lab of Physics

同轴光子晶体实验报告

learning_rate 未央软-11 2022.12.12

摘要:

本实验使用信号发生器和示波器测试电路传输效率和相位，推导同轴电缆的特征阻抗、传播常数、衰减系数参数，理解相速度和群速度的概念。使用同一种特征阻抗的同轴电缆交替连接组成同轴光子晶体，测量并计算传输效率、等效折射率、群速度等物理量。在同轴晶体中引入缺陷后，观察透射峰并同时测试特性。使用 MATLAB 编程计算上述物理量的理论值。实验结果与理论吻合度较好。

关键词:

同轴光子晶体; 能带理论; 色散关系; 衰减关系; "超光速"现象

目录

实验仪器

信号发生器、示波器、测试电缆（约 \(1\text{m}\) 长）、同轴电缆（约 \(5\text{m}\) 长）、三通接头、\(50\Omega\) 电阻。

实验原理

同轴电缆电学参数

对于两种特征阻抗的同轴电览交替连接构成的同轴光子晶体，当两种电缆的特征阻抗之比 \(Z_{0 H} / Z_{0 L}\) 越大、连接的电缆的数量越多，电磁带隙特性以及快慢光效应越明显。

常用同轴电缆的特征阻抗为 \(75 \Omega\) 和 \(50 \Omega\) 。使用这两种电缆的话阻抗比偏小，且两种电缆内外导体间介质的介电常数 \(\varepsilon_r\) 通不同，对应的半波长以及电缆长度也不同。

本实验只使用一种特征阻抗的同轴电缆。通过并联两根电缆的方法获得另一特征阻抗以及更高的阻抗比。

图1 单根同轴电缆的等效电路模型

同轴电缆的等效电路如图 \(5\) 所示，其中 \(R,L,G\) 分别为单位长度同轴电缆的电阻、电感、电容，且结合电磁学知识可知 \[ \begin{equation} \begin{cases} C & =\dfrac{2 \pi \varepsilon_0 \varepsilon_r}{\ln (b / a)} \\ L & =\dfrac{\mu_0 \mu_r}{2 \pi} \ln \dfrac{b}{a} \end{cases} \end{equation} \] 其中 \(a,b\) 为同轴电缆的内外径，\(\epsilon_r,\mu_r\) 分别为同轴电缆内外的相对介电常数以及相对磁导率，沿同轴电缆方向共有两个方向叠加，对应通解为 \[ \begin{equation} \begin{aligned} & v=v^{+}(x)+v^{-}(x)=V_0^{+} e^{-\gamma x}+V_0^{-} e^{\gamma x} \\ & i=i^{+}(x)+i^{-}(x)=\frac{1}{Z_0}\left(V_0^{+} e^{-\gamma x}-V_0^{-} e^{\gamma x}\right) \end{aligned} \end{equation} \] 对应电压波，其中 \(\gamma=\alpha+j\beta\) 为电缆的传播常数， 其表达式如下 \[ \begin{equation} \alpha+j \beta=\sqrt{(R+j \omega L)\left(G+j \omega C\right)} \end{equation} \] 且定义特征阻抗 \[ \begin{equation} Z_0=\sqrt{\dfrac{R+j \omega L}{G+j \omega C}} \end{equation} \] 本实验中，电缆的损耗很小，有 \(R\ll\omega L,G\ll\omega C\)，做近似有 \[ \begin{equation} \left\{\begin{array}{l} Z_0=\sqrt{\dfrac{L}{C}}=\dfrac{\ln (b / a)}{2 \pi} \sqrt{\dfrac{\mu_0 \mu_r}{\epsilon_0 \epsilon_r}} \\ \beta=\omega \sqrt{L C}=\omega \sqrt{\epsilon_0 \mu_0 \epsilon_r \mu_r} \end{array}\right. \end{equation} \]

透射与反射

对于半无限长的同轴电缆，信号从同轴电缆始端输入后，只有沿 \(+\) 方向传输的行波。实际同轴电缆都是有限长的，此时需考虑同轴电缆末端对信号的影响。

两根特征阻抗分别为 \(Z_0\) 和 \(Z_1\) 的同轴电缆连接在一起，两根同轴电缆之间的连接线很短，远小于信号波长和传输线长度。输入信号通过阻抗为 \(Z_0\) 的同轴电缆到达连接处的电压为 \(V_i\)，透射信号在连接处的电压为 \(V_t\) ，反射信号在连接处的电压为 \(V_r\) ，满足如下方程 \[ \begin{equation} \left\{\begin{array}{l} V_i+V_r=V_t \\ V_i / Z_0=V_r / Z_0+V_t / Z_1 \end{array}\right. \end{equation} \]

测量同轴光子晶体特性

由群速度公式 \[ \begin{equation} v_g=\dfrac{c} {n+\omega \dfrac{d n}{d \omega}} \end{equation} \]

可测得某一频率下示波器 \(\mathrm{Ch} 2\) 相对于 \(\mathrm{Ch} 1\) 的相位差 \(\varphi\)，投输入待测电路的信号为 \(\cos \omega t\)，初始相位为 0。输出信号为 \(\cos (\omega t+\varphi)\)，相位变化为 \(\varphi\) \[ \begin{equation} \cos (\omega t+\varphi)=\cos \omega\left(t-t_p\right) \Rightarrow \varphi=-\omega t_p \end{equation} \] 将同轴光子晶体效为一根折射率为 \(n\) 的均匀同轴电线，则延时 \(t_p\) 与等效折射率 \(n\) 的关系为 \[ \begin{equation} t_p=\dfrac{Ln}{c} \end{equation} \] 公式 (10)(11) 联立计算得到 \[ \begin{equation} n=-\frac{\varphi c}{\omega L}=-\frac{\varphi c}{2 \pi f L} \end{equation} \]

电学参数理论计算\((\,\mathbf{A}\,\)部分\()\)

衰减系数和幅度

由电压幅度值每米减少 \(\mathbf{0.571\%}\ll1\)，从而由衰减幅度公式 \[ \begin{equation} V(x)=V_0 e^{-\alpha x} \end{equation} \] 取近似 \(V_0e^{-\alpha x}\approx V_0(1-\alpha x)\)，故 \(\alpha \approx 0.00571\)，或者严格计算有 \[ \alpha =-\ln (1-0.571\%)=\textcolor{blue}{0.0057264}\color{blue}\text{m}^{-1} \] 代入公式 (1) 计算 \(10\text{m},20\text{m},30\text{m}\) 处的幅度列表如下（保留三位小数）

距离/m	0	10	20	30
幅度/V	1.000	0.944	0.892	0.842

表1 信号幅度随距离的变化

电缆并联

由图 1 可知，记每一小段 \(\Delta x\) 对应的复阻抗为 \(\Delta Z\)，总复阻抗为 \(N\Delta Z\)，则始末并联和每一小段处并联是等效的 \[ \sum_{i}\dfrac{\Delta Z\cdot \Delta Z}{\Delta Z+\Delta Z}=\dfrac{1}{2}\sum_{i}\Delta Z=\dfrac{\displaystyle \sum_{i}\Delta Z\cdot \displaystyle \sum_{i}\Delta Z}{\displaystyle \sum_{i}\Delta Z+\displaystyle \sum_{i}\Delta Z} \] 由并联性质 \[ R^{\prime}=\textcolor{blue}{\dfrac{R}{2}},L^{\prime}=\textcolor{blue}{\dfrac{L}{2}},G^{\prime}=\textcolor{blue}{2G},C^{\prime}=\textcolor{blue}{2C} \] 代入公式 (5)(6) 得 \[ \begin{equation} Z_0^{\prime}=\textcolor{blue}{\frac{1}{2} Z_0}, \alpha^{\prime}=\textcolor{blue}{\alpha}, \beta^{\prime}=\textcolor{blue}{\beta} \end{equation} \] 一种更简便的想法是电压相同时，电流变成原来的两倍，从而复阻抗折半，而每个电缆上的电流衰减幅度以及相位变化与原来相同，总比例也不变，从而传播系数不变。

高低阻抗交替连接

由上问可知，高特征阻抗使用单根电缆，低特征阻抗使用并联电缆，且有 \[ \begin{equation} \frac{Z_{0 H}}{Z_{0 L}}=\textcolor{blue}{2} \end{equation} \] 可以通过多次并联的方法得到更高的阻抗比。

频率计算

由电磁波公式 \[ \begin{equation} f=\frac{v}{\lambda}=\frac{c}{\sqrt{\epsilon_r \mu_r} \lambda} \end{equation} \] 代入计算得 \[ f_0=\frac{c}{\sqrt{\epsilon_r \mu_r} \lambda} \approx \frac{299792458 \mathrm{~m} / s}{\sqrt{2.354} \cdot 10 \mathrm{~m}}=\textcolor{blue}{19.540 \mathrm{~}\mathrm{MHz}},\dfrac{1}{2}f_0=\textcolor{blue}{9.770 \mathrm{~}\mathrm{MHz}} \]

间接测量色散关系\((\,\mathbf{B}\,\)部分\()\)

正常光子晶体

测量数据

以 \(1 \mathrm{MHz}\) 为间隔，在 \(1 \mathrm{MHz}\) 至 \(20 \mathrm{MHz}\) 范围内改变频率，观察示波器 \(\mathrm{Ch} 1\) 和 \(\mathrm{Ch} 2\) 的峰峰值及其变化情况，测试电路图如下

图2 同轴光子晶体测试电路

测量数据整理如下

f (MHz)	V1pp (V)	V2pp (V)	φ (°)
1	10.250	9.220	-72.10
2	10.260	9.380	-139.80
3	10.270	8.750	147.70
4	10.280	8.730	80.80
5	10.270	8.250	2.70
6	10.320	7.590	-66.40
7	10.480	5.890	-166.90
8	10.510	2.320	138.10
9	10.570	1.450	112.60
10	10.640	1.360	92.10
11	10.750	1.850	68.20
12	10.860	4.160	27.00
13	10.980	8.000	-63.80
14	11.110	8.160	-145.60
15	11.280	8.770	140.40
16	11.450	8.700	71.10
17	11.630	9.190	2.40
18	11.810	8.670	-67.30
19	12.060	9.340	-131.90
20	12.310	9.030	155.50

表2 正常光子晶体Ch1 和 Ch2 峰峰值和相位差随频率变化

传输效率和等效折射率

同轴晶体未接入、阻抗匹配状态下 \(50\Omega\) 电阻上消耗的功率为 \(𝑃_1\)，接入同轴光子晶体、阻抗匹配状态下 \(50\Omega\) 电阻上消耗的有功功率为 \(P_2\)，故计算公式如下 \[ \begin{equation} P_i=\frac{V_i^2}{R_0} \quad i=1,2 \end{equation} \] 则传输效率 \[ \begin{equation} \eta=\frac{P_2}{P_1}=\left(\frac{V_2}{V_1}\right)^2 \end{equation} \] 代入公式 (18)(12) 计算效率以及折射率，其中相位差相对小频率来说需要减去 \(360^{\circ}\) 的整数倍，保证相位差是单调递减的，计算以及列表如下

f (MHz)	φ' (°)	η	n
1	-72.10	80.91%	1.71548
2	-139.80	83.58%	1.66313
3	-212.30	72.59%	1.68375
4	-279.20	72.12%	1.66076
5	-357.30	64.53%	1.70025
6	-426.40	54.09%	1.69089
7	-526.90	31.59%	1.79094
8	-581.90	4.87%	1.73065
9	-607.40	1.88%	1.60577
10	-627.90	1.63%	1.49397
11	-651.80	2.96%	1.40985
12	-693.00	14.67%	1.37405
13	-783.80	53.09%	1.43454
14	-865.60	53.95%	1.47109
15	-939.60	60.45%	1.49040
16	-1008.90	57.73%	1.50030
17	-1077.60	62.44%	1.50820
18	-1147.30	53.89%	1.51654
19	-1211.90	59.98%	1.51762
20	-1284.50	53.81%	1.52811

表3 正常光子晶体传输效率与等效折射率计算表格

其效率和频率关系如下

图3 正常光子晶体传输效率与频率关系曲线

等效折射率与频率的变化如下

图4 正常光子晶体等效折射率与频率关系曲线

由图表可得 \(f=\textcolor{blue}{1,2,3,4,5,13,14,15,16,17,18,19,20 \mathrm{MHz}}\) 为通带, \(f=\textcolor{blue}{8,9,10,11 \mathrm{MHz}}\) 为带隙，在 \(f=\textcolor{blue}{7,8,9,10,11 M H Z}\) 出现反常色散。

群速度

代入公式 (9) ，其中离散化处理 \(\dfrac{dn}{d\omega}\approx\dfrac{n_{k+1}-n_k}{\omega_{k+1}-\omega_{k}},k=1,\cdots,19\)，列表如下（取 \(\dfrac{v_g}{c}\) 为纵坐标）

f (MHz)	ωdn/dω	vg/c
1	-0.05234	0.60127
2	0.04124	0.58673
3	-0.06900	0.61929
4	0.15799	0.54983
5	-0.04679	0.60479
6	0.60026	0.43646
7	-0.42203	0.73051
8	-0.99904	1.36686
9	-1.00621	1.66789
10	-0.84119	1.53193
11	-0.39378	0.98418
12	0.72587	0.47621
13	0.47518	0.52364
14	0.27029	0.57426
15	0.14856	0.61015
16	0.12638	0.61475
17	0.14183	0.60605
18	0.01941	0.65106
19	0.19927	0.58245
20	N/A	N/A

表4 正常光子晶体群速度计算表格

将群速度和频率作图如下

图5 正常光子晶体群速度与频率关系曲线

可以在处于反常色散频率时群速度可能超光速（对应 \(f=8 ,9 ,10\text{MHz}\)）。

缺陷光子晶体

测量数据

在同轴光子晶体中增加一根同轴电缆，形成有缺陷的同轴电缆，按照 4.1 节内容测量数据结果如下

f (MHz)	V1pp (V)	V2pp (V)	φ (°)
1	10.250	9.280	-81.60
2	10.260	9.210	-163.30
3	10.270	7.860	119.80
4	10.280	8.570	44.90
5	10.260	8.540	-45.00
6	10.310	6.110	-148.40
7	10.380	2.660	150.00
8	10.430	1.980	118.80
9	10.490	2.950	84.30
10	10.560	4.990	-34.40
11	10.670	2.340	-90.10
12	10.790	2.340	-122.70
13	10.910	4.550	-172.50
14	11.030	7.790	96.50
15	11.190	8.470	4.30
16	11.370	7.760	-76.20
17	11.530	8.570	-151.80
18	11.720	8.490	126.90
19	11.960	8.930	51.80
20	12.220	8.670	-27.30

表5 缺陷光子晶体Ch1 和 Ch2 峰峰值和相位差随频率变化

传输效率和等效折射率

同理代入公式 (16)(12) 列表如下

f (MHz)	φ' (°)	η	n
1	-81.60	81.97%	1.94151
2	-163.30	80.58%	1.94270
3	-240.20	58.57%	1.90503
4	-315.10	69.50%	1.87430
5	-405.00	69.28%	1.92724
6	-508.40	35.12%	2.01606
7	-570.00	6.57%	1.93743
8	-601.20	3.60%	1.78805
9	-635.70	7.91%	1.68058
10	-754.40	22.33%	1.79495
11	-810.10	4.81%	1.75225
12	-842.70	4.70%	1.67087
13	-892.50	17.39%	1.63348
14	-983.50	49.88%	1.67146
15	-1075.70	57.29%	1.70628
16	-1156.20	46.58%	1.71935
17	-1231.80	55.25%	1.72402
18	-1313.10	52.48%	1.73570
19	-1388.20	55.75%	1.73840
20	-1467.30	50.34%	1.74558

表6 缺陷光子晶体传输效率与等效折射率计算表格

其效率和频率关系如下，其中透射峰使用 \(\textcolor{red}{\swarrow}\) 标出，两侧由于多项式拟合的边界效应其峰可以排除

图6 缺陷光子晶体传输效率与频率关系曲线

图7 缺陷光子晶体等效折射率与频率关系曲线

计算表 6和表 3计算的传输效率之比，列表如下

f (MHz)	η1	η2	η2/η1
1	80.91%	81.97%	1.013
2	83.58%	80.58%	0.964
3	72.59%	58.57%	0.807
4	72.12%	69.50%	0.964
5	64.53%	69.28%	1.074
6	54.09%	35.12%	0.649
7	31.59%	6.57%	0.208
8	4.87%	3.60%	0.740
9	1.88%	7.91%	4.202
10	1.63%	22.33%	13.667
11	2.96%	4.81%	1.624
12	14.67%	4.70%	0.321
13	53.09%	17.39%	0.328
14	53.95%	49.88%	0.925
15	60.45%	57.29%	0.948
16	57.73%	46.58%	0.807
17	62.44%	55.25%	0.885
18	53.89%	52.48%	0.974
19	59.98%	55.75%	0.929
20	53.81%	50.34%	0.935

表7 缺陷光子晶体传输效率与正常光子晶体传输效率之比

将两者的传输效率之比绘图如下

图8 两种光子晶体传输效率之比随频率的变化曲线

可以看出在 \(f=9,10,11\text{MHz}\) 处的效率提升非常显著，其余基本在 \(1\) 附近波动。

群速度

同理离散化代入公式 (9) 计算群速度，列表如下

f (MHz)	ωdn/dω	vg/c
1	0.00119	0.51475
2	-0.07534	0.53552
3	-0.09220	0.55162
4	0.21176	0.47937
5	0.44414	0.42170
6	-0.47178	0.64755
7	-1.04570	1.12142
8	-0.85972	1.07721
9	1.02929	0.36902
10	-0.42698	0.73101
11	-0.89521	1.16681
12	-0.44859	0.81815
13	0.49371	0.47010
14	0.48744	0.46320
15	0.19600	0.52569
16	0.07474	0.55739
17	0.19867	0.52011
18	0.04846	0.56049
19	0.13645	0.53338
20	N/A	N/A

表8 缺陷光子晶体群速度计算表格

其群速度和频率的曲线绘制如下

图9 正常光子晶体和缺陷光子晶体群速度对比曲线

由图像可知，在 \(f=9\text{MHz}\) 处出现强正常色散，该处的 \(v_g\) 显著小于 \(c\)，为慢光速现象。

色散关系直接测量\((\,\mathbf{C}\,\)部分\()\)

特征频率

在 50 kHz 至 150 kHz 范围内改变信号发生器频率，观察脉冲信号的波形、幅度和延时变化情况。从上到下分别为 \(114、130、160 \text{ kHz}\) 时的波形，前两者为幅度最小的两个频率，后者为频率之间信号幅度最大且延时最长的频率。

图10 频率为 114kHz 脉冲信号波形

图11 频率为 130kHz 脉冲信号波形

图12 频率为 160kHz 脉冲信号波形

群速度计算

对图片截图计算中心波腹的延迟像素与示波器背景单位格点像素比例，例如对第一张图两者像素为 333 和 543，故其延迟时间为 \[ t_1=\frac{333 \mathrm{~pix}}{543 \mathrm{~pix}} \cdot 500 \mathrm{~ns}=306.6\text{~ns} \] 同理对三幅图片处理列表如下，其中由于设置 \(100T\) 脉冲波形，高频正弦波信号的频率为原始频率的 \(100\) 倍，使用 \(\text{MHz}\) 单位方便统一

f (MHz)	基准像素	波腹像素	Δt/ns
11.4	543	333	306.6
13.0	493	99	100.4
16.0	498	164	164.7

表9 像素法延时长度计算

由群速度代表波包的运动，\(v_{g}=\dfrac{l}{\Delta t}\)，计算 \(\dfrac{v_g}{c}=\dfrac{l}{c\Delta t}\) 列表如下

f (MHz)	vg/c
11.4	0.3807
13.0	1.1628
16.0	0.7090

表10 群速度直接测量结果

将该结果标注在图5中如下，其中 \(\textcolor{red}{\bullet}\) 代表测量点

图13 利用脉冲直接测量和间接测量群速度比对

可以发现，除中间点外，测量结果基本符合预期。

传输特性理论计算\((\,\mathbf{D.1-D.4}\,\)部分\()\)

叠加电压求解

由公式 (6) 的两个方程反解得 \[ \begin{equation} \begin{cases} V_0^{+} =\dfrac{1}{2}\left(v_2+Z_0 i_2\right) e^{\gamma l} \\ V_0^{-} =\dfrac{1}{2}\left(v_2-Z_0 i_2\right) e^{-\gamma l} \end{cases} \end{equation} \]

传递矩阵计算

取 \(x=0\) 的点对应 \(v_1,i_1\) \[ \begin{equation} \begin{cases} v_1=v^{+}(0)+v^{-}(0)=V_0^{+}+V_0^{-} \\ i_1=i^{+}(0)+i^{-}(0)=\dfrac{1}{Z_0}\left(V_0^{+}-V_0^{-}\right) \end{cases} \end{equation} \] 代入公式 (17) 得 \[ \begin{equation} \begin{aligned} v_1 & =\frac{1}{2}\left(v_2+Z_0 i_2\right) e^{\gamma l}+\frac{1}{2}\left(v_2-Z_0 i_2\right) e^{-\gamma l} \\ i_1 & =\frac{1}{Z_0}\left(\frac{1}{2}\left(v_2+Z_0 i_2\right) e^{\gamma l}-\frac{1}{2}\left(v_2-Z_0 i_2\right) e^{-\gamma l}\right) \end{aligned} \end{equation} \] 写成矩阵的形式如下 \[ \begin{equation} \left[\begin{array}{l} v_1 \\ i_1 \end{array}\right]=\left[\begin{array}{cl} \cosh (\gamma l) & Z_0 \sinh (\gamma l) \\ \dfrac{1}{Z_0} \sinh (\gamma l) & \cosh (\gamma l) \end{array}\right] \cdot\left[\begin{array}{l} v_2 \\ i_2 \end{array}\right] \end{equation} \] 则其中 \(\left[\begin{array}{cl} \cosh (\gamma l) & Z_0 \sinh (\gamma l) \\ \dfrac{1}{Z_0} \sinh (\gamma l) & \cosh (\gamma l) \end{array}\right]\) 为传递矩阵。

传递矩阵性质

将两个同轴电缆串联连接如下

图14 同轴光子晶体的其中两段级联同轴电缆

由上节传递矩阵的性质 \[ \begin{equation} \left[\begin{array}{l} v_1 \\ i_1 \end{array}\right]=A_1\left[\begin{array}{l} v_2 \\ i_2 \end{array}\right],\left[\begin{array}{l} v_2 \\ i_2 \end{array}\right]=A_2\left[\begin{array}{l} v_3 \\ i_3 \end{array}\right] \end{equation} \] 后者代入前者可得 \[ \begin{equation} \left[\begin{array}{l} v_1 \\ i_1 \end{array}\right]=A_1 A_2\left[\begin{array}{l} v_3 \\ i_3 \end{array}\right] \end{equation} \] 故总的传递矩阵为 \(A=A_1 A_2\)，从而通过传递矩阵相乘的方法可以得到若干个同轴光子晶体连接形成的等效传递矩阵。

负载电压表达式

令总传递矩阵 \(A_T=\left[\begin{array}{cc}A & B \\ C & D\end{array}\right]\)，信号源电动势 \(e\) 与内阻 \(Z_0\)，联立方程 \[ \begin{equation} \left\{\begin{array}{l} {\left[\begin{array}{l} v_1 \\ i_1 \end{array}\right]=\left[\begin{array}{cc}A & B \\ C & D\end{array}\right]\cdot\left[\begin{array}{l} v_L \\ i_L \end{array}\right]} \\ v_L=i_L \cdot Z_0 \\ e=v_1+i_1 \cdot Z_0 \end{array}\right. \end{equation} \] 解得负载电压的表达式 \[ \begin{equation} V_L=\frac{e}{A+\dfrac{1}{Z_0} B+Z_0 C+D} \end{equation} \]

编程计算\((\,\mathbf{D.5}\,\)部分\()\)

使用 matlab 对给定数据进行模拟，已知单根同轴电缆的特征阻抗 \(Z_0=50 \Omega\)，介质参数 \(\varepsilon_r=2.354 、 \mu_r \approx 1\)，衰减系数 \(\alpha=\) \(1.810 \times 10^{-6} \sqrt{f} m^{-1}\)，其中 \(f\) 的单位为 \(\mathrm{Hz}\)。

由连接图可知，令 \(A_1\) 为原始同轴电缆，\(A_2\) 为并联之后的同轴电缆，其总传递矩阵如下 \[ \begin{equation} A_T=A_2 A_1 A_2 A_1 A_2 A_1 A_2 \end{equation} \]

程序设计

将转移矩阵对应的元素设计在函数中，例如第一个元素获取函数 A1() 如下

function [a] = A1(f)
    A1=[cosh(ga(f)*5) 50*sinh(ga(f)*5);1/50*sinh(ga(f)*5) cosh(ga(f)*5)];
    A2=[cosh(ga(f)*5) 25*sinh(ga(f)*5);1/25*sinh(ga(f)*5) cosh(ga(f)*5)];
    % A=A2*A1*A2*A1*A2*A1*A2; % 缺陷光子晶体
    A=A2*A2*A2*A1*A2*A1*A2;   % 正常光子晶体
    a=A(1,1) ;
end

主函数针对传输效率、等效折射率、群速度计算，纵坐标标注不同结果

f = 0:0.1:80;
x = f;
angs = f;
% 计算效率
for i=1:801
    x(i) = eta(f(i) * 1e6);
end
% 计算相位差以及使用 unwrap() 函数
for i=1:801
   angs(i) = ang(f(i)*1e6);
end
angs = unwrap(angs);
% 计算等效等效折射率
for i=1:801
   x(i) = n(f(i)*1e6, angs(i));
end
% 计算群速度
vg = x(1:800);
for i=1:800
   vg(i) = 1 / (x(i) + f(i)/0.1*(x(i+1)-x(i)));
end
% 对不同情况绘制图像，选定不同的 y 轴标注
% plot(f, x);
plot(f(1:800), vg);
xlabel('$f/MHz$', 'interpreter', 'latex', 'FontSize', 15);
% ylabel('$\eta$', 'interpreter', 'latex', 'FontSize', 20);
% ylabel('$n$', 'interpreter', 'latex', 'FontSize', 20);
ylabel('$\frac{v_g}{c}$', 'interpreter', 'latex', 'FontSize', 20);

其中使用 unwrap() 函数对相位差有跃变的情况加以处理，自动增加 \(n\pi\)。

正常光子晶体

取定 \(\textcolor{blue}{f\in [0,8]\text{~MHz}}\) 步长为 \(\textcolor{blue}{0.1\text{~MHZ}}\) 进行程序模拟，绘制图像如下

传输效率

图15 正常光子晶体效率与频率关系曲线（程序模拟）

从左往右依次为通带、带隙、通带交替产生，并且随着频率地增加，多次反射和透射产生的损耗增加，从而引起效率极值点下降。

等效折射率

图16 正常光子晶体等效折射率与频率关系曲线（程序模拟）

可以看出在频率逐渐增大的过程中存在多个反常色散区域，这与反常色散曲线吻合。

群速度

图17 正常光子晶体群速度与频率关系曲线（程序模拟）

可以看出在反常色散区域群速度超过光速 \(c\)，与预期吻合。

缺陷光子晶体

假定增加的电缆位于左侧，其总传递矩阵更改为 \[ \begin{equation} A_T^{\prime}=A_2 A_2 A_2 A_1 A_2 A_1 A_2 \end{equation} \]

传输效率

仅仅更改矩阵，其余代码不变，传输效率关系图如下

图18 缺陷光子晶体效率与频率关系曲线（程序模拟）

与图 6对比，在每个通带和带隙之间增加了若干个透射峰（确切来说有 \(4\) 个），与实测结论完全相符。

等效折射率

图19 缺陷光子晶体等效折射率与频率关系曲线（程序模拟）

群速度

图20 缺陷光子晶体群速度与频率关系曲线（程序模拟）

有缺陷的情形均与实际情况吻合，同时符合能带理论。

原始数据

图21 原始数据截图1

图22 原始数据截图2