基于MATLAB的2×2光纤定向耦合器设计

基于MATLAB的2×2光纤定向耦合器设计

1 设计原理

1.1 单模光纤的传导场

如图1，光纤的横截面有三层介质，分别是是芯层、包层和涂层，芯层折射率n1稍大于包层折射率n2，导波光由于全反射背包层约束在芯层中沿光纤延伸方向传播。假设光的传播方向为光纤中心轴方向。

图1 阶跃光纤横截面结构图

为简化讨论，只考虑基模的耦合。已知光纤中传导场表达式为

E(x,y,z,t)?a?z??e?x,y??ei?z?ei?t （1-1）

其中，a?z?为光纤中导波光场的场振幅，e?x,y?为光纤中导波光场的场分布，?为基模场的传播常数，?为角频率。

某时刻在光纤中的传导场的空间分布就与a?z?，e?x,y?和?为相关。

1.1.1 单模光纤的场分布

当给定波导（即光纤）的边界条件时，求解波导场方程可得本征解及相应的本征值?，即模式。模式是波导结构的固有电磁共振属性的表征。单模光纤中只能存在基模，其场分布是确定的，可由亥姆霍兹方程求得。在弱导光纤中的电磁波，其横向场分量Et、Ht远大于纵向场分量Ez、Hz ，而且横向场分量是线偏振的。于是我们把电场的横向分量取为y轴方向，即Et=Ey。

亥姆霍兹方程为

?Ey?k0n?r?Ey?0 （1-2）

222其中n?r????n1?n2r?ar?a，k0?2??为真空中的光波矢量。

利用分离变量法，将方程1-2在圆柱坐标系中求解，并结合电磁场的边界条件，可以解出电场的横向分量Ey：

Ey?Jm?Ura?Jm?U??Acosm???Km?Wra?Km?Wr?a?r?a （1-3）

其中，r是点到光纤中心轴的距离，m是整数，Jm和Km分别是m阶贝塞尔函数和m阶变态汉克尔函数；a是光纤芯层半径，一般单模光纤的a=2~5μm。?是极角。

横向磁场只包含Hx分量，由于横向电场与横向磁场的比等于波阻抗Z0=122?，故可由Ey计算出Hx：

Hx?Jm?Ura?Jm?U????AZ0??n?r??cosm????Km?Wra?Km?Wr?a?r?a （1-4）

再由麦克斯韦方程组可求出场的纵向分量Ez、Hz：

??Un1?[Jm?1?Ura?Jm?U??sin?m?1????Jm?1?Ura?Jm?U??sin?m?1??]?Ez??jA2k0a?????Wn1?[Km?1?Wra?Km?W??sin?m?1????Km?1?Wra?Km?W??sin?m?1??]??U[Jm?1?Ura?Jm?U??cos?m?1????Jm?1?Ura?Jm?U??cos?m?1??]????jA2k0aZ0????W[Km?1?Wra?Km?W??cos?m?1????Km?1?Wra?Km?W??cos?m?1??]?r?a （1-5）

r?aEzr?a（1-6）

r?a1.1.2 单模光纤的传播常数

由于在无耦合的情况下，光纤中导波光场的场振幅a?z??1。故当某时刻t的场分布只同传播常数有关系。

由式1-3可以得出光纤的本征方程本征方程，即r =0连续。

U?Jm?1?U?Jm?U??W?Km?1?W?Km?W?

（1-7）

同理，由横向磁场的连续性同样可以得出本征方程:

U?Jm?1?U?Jm?U??-W?Km?1?W?Km?W?

（1-8）

由于特征方程是超越方程，只能靠计算机数值求解U(或W)。 其中U与W定义为

U?an1k0??222 （1-9）

W?a?2?n2k0 （1-10）

22U与W是场的横向传播常数，U反映了导模在芯区中的驻波场的横向振荡频率，W值反映了导模在包层中的消逝场的衰减速度。W越大衰减越快。

式1-6、1-7与β相联系，因此本征方程实际是关于β的一个超越方程。对于不同的U值，可求得相应的β值。由于贝塞尔函数及其导数具有周期振荡性质，所以本征值方程可以有多个不同的解βmn(m=0,1,2,3..., n=1,2,3...),每一个βmn都对应于一个导模。单模光纤的本征值为β01 。

光纤的结构参数可由其归一化频率V表征：

V?U2?W2?n1k0a2? （1-11）

其中，其中相对折射率差???n1?n2?n1，n1、n2（n1>n2）分别是芯层和包层的折射率，在阶跃折射率光纤中均为常数。实际的光纤的?很小，一般小于1%，因此也称之为弱导光纤。

归一化频率V描述了光纤的结构参数a，?，n1以及工作波长（包含在k0中），是一个重要的综合性参数。V值越大，允许存在的导模数就越多。除了基模之外，其他导模都可能在某一个V值以下截止，这时导模转化为辐射模。当某一模式截止时，它就不能沿光纤传输了。

当W2 > 0时，对于导波模，场在纤芯外是指数衰减的；当W2 < 0时，W为虚数，场在纤芯外不再衰减了，也就不被约束在芯层中传播了，这种波叫做辐射波。而W2 =0正好是临界状态，以此作为导波模被截止的条件，此时的W记做Wc，相应的U和V值记做Uc和Vc，并称Vc为归一化截止频率。

截止条件(W=0)下，特征方程可以简化为

Jm?1?Uc??Jm?1?Vc??0

（1-12）

对于m=0，n=1时，截止时由1-9式得Vc=Uc=0，说明这种模式永不截止，是光纤中的最低阶模，也称基模。第二个归一化截止频率较低的模式是当m=1，n=1时，称为二阶模，其归一化截止频率Vc=Uc=2.40483。当光纤中只有基模工作而其他模式均被截止时，就称为单模光纤。要保证单模传输，必须要求光纤的归一化频率小于二阶模LP11的截止 频率Vc，即让二阶模截止：

V?Vc?m?1,n?1??2.40483 （1-13）

此即单模传输条件。

1.1.3 LPmn模的功率密度以及功率

平均坡印廷矢量

???*?E?H?Sav的大小就代表功率密度，可由S?Re???2???1以及场解1-3~1-6

式，可以计算出光纤横截面上的功率密度(Sav)近似为：

Sz?A2?2Z0??n?r??cos22???????JUraJUr?a?mmm???2???????KWraKWr?am?m （1-14）

平均坡印廷矢量对横截面积分便得到光纤中的功率P：

P?n2?aA22?4Z0??[?V2U2?Km?1?W??Km?1?W?Km?W?]

2 （1-15）

令P=1 (归一化功率)，可算出以上各式的常数A：

A?[4Z0?n2?a2?]??U211V??Km?W?[Km?1?W?Km?1?W?]2 （1-16）

1.2 2×2光纤定向耦合器的场振幅

在两根单模光纤耦合的情况下，光纤中导波光场的场振幅a?z??1，而是随着轴向距离z变化的。故某时刻的场传导方程就还需求解场振幅。

2×2光纤定向耦合器的示意图如下图所示：

图(a) 2×2光纤定向耦合器的示意图 图(b) 2×2光纤耦合器横截面示意图 图(c) 耦合器纵向剖面示意图(即rz面)

如图所示，两平行直光纤的纤芯相互接近时，在其中传输的基模场分布就会互相渗透和交叠。这样1号光纤中传播的导模场将在2号光纤中产生极化作用，从而在2号光纤中激起传导模。设nc_01和 nc_02分别为1号光纤和2号光纤的芯层折射率，nc_clad为包层折射率。可将两光纤中的传导场写为：

（1-17） （1-18）

其中e1、e2是两根单模阶跃光纤中的横向电场，其表达式由（1-3）式给出。a1(z)、a2(z)为复数振幅，可由耦合波方程解得：

（1-19） （1-20）

其中，??=?1－?2为失配相位常数。k12、k21为耦合系数，取决于波导结构以及波长，且满足：

（1-21）

耦合系数可由下式计算：

（1-22） （1-23）

其中，n1(r)、n2(r)分别代表1号光纤和2号光纤的折射率分布。P0（=1，已经归一化）为输入光功率。（1-19）式两边对z求导，并利用(1-20)式，可得到：

???i???a1??k21a12a1?0

（1-24）

求解(1-24)式，解得

?a1(z)??c1e??j(???(??)?4k21222)?zj(???(??)?4k2122)?c2e2?z?? (20) （1-25） ??再利用(1-19)式解得

jk12j(???(??)?4k21)2j(???2222j(????(??)?4k2122)a2(z)?[c1e2?z （1-26）

j(????(??)?4k2122 ?c2(??)?4k21)2)e2?z] (21)积分常数c1和c2由边界条件a1(0)?1; a2(0)?0 (20)，给出 则最终求解出

相位常数??=?1－?2=0。

图3 1号和2号光纤的耦合功率密度分布三维图

图4 1号和2号光纤的耦合功率密度分布随轴向变化比较图

由图可以看出，在没有损耗的情况下，两根光纤的功率密度随传播距离交替变化。当1号光纤在Z轴上的0.4404，1.3173，1.1981以及3.0751处时，1号光纤中的能量几乎全部耦合到2号光纤中。故而这些点到入射端的距离均为耦合比最大时的耦合距离。这个结论从下面的耦合效率—耦合距离曲线图中可以更清晰的获得。

图5 耦合效率—耦合距离曲线图

由上图可知，两根光纤的耦合效率在上述耦合长度处约为1，1号光纤中的能量几乎全部耦合到2号光纤中。通过相同仿真方法，在保持失配相位常数为零的情况下，耦合效率不随参数的变化而改变，始终为1。

图6 耦合距离虽纤芯半径变化曲线图

调试中将2号光纤的参数进行固定，对1号光纤的纤芯半径进行改变，获得如图6所示不同纤芯半径下耦合效率对耦合距离的图像。我们可以得出，耦合距离随着纤芯半径的减小也减小，即耦合达到最大值的距离随着光纤的纤芯半径的减小而缩短。

4.2.2 失配相位常数不为0时的耦合情况

通过MATLAB的调试，输入不同光纤参数并使失配相位常数不为零。在输入参数1号光纤的纤芯半径为4.8μm，2号光纤的纤芯半径4.5μm， 1号光纤和2号光纤的纤芯折射率1.5时，失配相位常数??=?1－?2=738.3652。仿真输出1号和2号光纤的耦合功率密度分布随轴向变化比较图与耦合效率—耦合距离曲线图分别如图7、图8所示。

图7 1号和2号光纤的耦合功率密度分布随轴向变化比较图

图8 耦合效率—耦合距离曲线图

由图7，图8可以得出，失配相位常数不为零时，耦合效率极低，只有10?4的数量级，即不到0.1%。光信号基本在直通臂中传输，在耦合臂中的耦合能量非常少。

4.2.3 失配相位常数同耦合效率的关系

通过仿真，输入不同的光纤参数，可以获得不同的??——耦合效率的对应数值，数据如下表所示。

表1 ??——耦合效率的对应数据表

?? -305.554 -105.682 -17.9916 -8.5928 -6.5234 -4.3285 0.035 0 1 0.1529 1.9763 0.8261 0.3382 2.1959 0.7938 0.4673 3.2941 0.6315 耦合效率 2.72E-04 0.0013 ?? -2.1835 -1.9385 0.8724 8.788 0.1955 耦合效率 0.7258 ?? 6.5899 10.9868 21.9916 110.6821 223.191 0.135 0.0381 0.0018 5.35E-04 耦合效率 0.3009 由上述数据获得的曲线图如下

1.20E+001.00E+008.00E-016.00E-014.00E-012.00E-010.00E+00-400-300-200-1000-2.00E-01100200300400

图9 ??——耦合效率曲线图

其中，横坐标为??，纵坐标为偶和效率。由图9可以得出，当??=0时，耦合效率获得最大。

5 心得体会

2×2光纤定向耦合器的课程设计是一次集合理论学习、程序设计、计算机仿真研究等项目的综合性设计，是具有较高难度和较全面锻炼个人综合能力的课程设计。通过这次课程设计，我不仅在对2×2光纤定向耦合器的相关理论知识有了全面而深入的学习，还对MATLAB编程仿真以及通过计算机研究物理理论等科研技能有了锻炼和提升，可以说这次课程设计对我的理论学习和程序设计都有很大作用。

5.1 深化理论学习

由于在平常的学习中只是对2×2光纤定向耦合器的性能和器件特点有了一定的学习，而对其理论原理的学习缺乏重视，故而在计算机设计仿真前需要做大量的理论学习。在这次设计中，我先后复习了电动力学、光纤光学、数学物理方法、数值方法等科目的知识，花费了较多时间在理论推导和器件原理的学习中，独立求解了耦合光纤中耦合波的振幅函数、耦合波场、能流密度、耦合效率（耦合比）的表达式。

首先，在这些理论学习过程中，我更加认识到做理论研究的艰难和繁琐，不论是基本的概念学习还是理论过程的推导，都需要细心、认真、耐心的一步一步去做，任何一个符号的错误就可能使前后结论天壤之别。其次，理论过程的推导是需要扎实的数学功底的，这就要求我在平时学习数学知识的同时，注重计算技巧和成规方法的练习和记忆，这样才能有助于理论知识的学习和巩固。

5.2 加强编程练习

今后的科学研究是要求学生具有一定的程序编程基础的，而这一点在2×2光纤定向耦合器的课程设计中体现的至关重要。这次课程设计要求在MATLAB中对理论结果进行计算机仿真实现，其中需要用到许多诸如变量输入，特殊函数求解，三维曲面绘制等初级MATLAB语言，这些语言如果没有一定的基础，短时间学习根本不能达到灵活应用，对于较好实现理论结果，高效完成设计有很大影响。

MATLAB语言是当今国际上自动控制领域的首选计算机语言，也是众多理工科专业最适合的计算机数学语言。所以，系统、全面掌握MATLAB在数学运算问题中的应用，对我们今后的发展显得更加重要。通过这次课程设计，虽然只是对MATLAB的一些基础有了学习，但这对今后的深入学习仍然显得很重要。

5.3 通过计算机仿真研究

在MATLAB上通过程序实现理论结果后，我们看到了直观的曲线图像，这是以前通过手工计算几乎无法实现的。这一点在求解由贝塞尔函数和变态汉克尔函数构成的特征方程时体现的尤为突出。特征方程本身就是一个超越方程，通过数值求解这样大运算量、计算复杂才能获得结果的方法只有计算机才得以胜任。在设计中通过改变输入值，经过几秒便可以得出结果，这若是通过手工计算是难以想象的困难的。而且通过改变输入可以很快获得??——耦合效率的关系式，这样获得的结果不但可信度很高，而且高效方便。

6 参考文献

[1]刘德明.光纤光学[M].北京:科学出版社,2008.35-38;132-139. [2]王元明. 数学物理方法[M].北京:高等教育出版社，2004.18-28

[3]谢处方，饶克谨. 电磁场与电磁波[M]. 北京:高等教育出版社,2006.205-208 [4]薛定宇，陈阳泉. 高等应用数学问题的MATLAB求解[M].北京:清华大学出版社,2004.9-41.

附录

%定义光纤耦合器的基本参数 m=0;n=1;%定义光纤模式 n_clad=1.4955;%包层折射率 lamda_0=1.55*1e-06;%波导光波长 k0=2*pi/lamda_0;%真空中波矢量 z0=122*pi;%真空波阻抗 a1_0=1;a2_0=0;%输入边界条件 eption_0=1/36/pi*1e-09;%真空介电常数 omega=2*pi/lamda_0*3*1e8;%光波角频率

%定义1号光纤的参数

r01_core=input('输入1号光纤的纤芯半径(单位为μm)=')*1e-06; %定义1号光纤纤芯半径 r01_clad=3*r01_core; %定义1号光纤的包层半径

nc01=input('输入1号光纤的纤芯折射率='); %定义1号光纤的纤芯折射率 delta1=(nc01-n_clad)/nc01;%相对折射率差 V1=nc01*k0*r01_core*sqrt(2*delta1);%归一化频率

while V1>2.40483,disp('1号光纤不是单模光纤,请重新输入');%判断1号光纤是否满足单模条件

r01_core=input('输入1号光纤的纤芯半径(单位为μm)=')*1e-06; %定义1号光纤纤芯半径 r01_clad=3*r01_core; %定义1号光纤的包层半径

nc01=input('输入1号光纤的纤芯折射率='); %定义1号光纤的纤芯折射率 delta1=(nc01-n_clad)/nc01; %相对折射率差 V1=nc01*k0*r01_core*sqrt(2*delta1); %归一化频率 end;

disp('1号光纤是单模光纤,输入正确') %1号光纤定义结束

%由特征方程求1号光纤的归一化径向相位常数和归一化径向衰减常数 fun1=@(U1)

U1.*besselj(m+1,U1)/besselj(m,U1)-sqrt(V1^2-U1^2).*besselk(m+1,sqrt(V1^2-U1^2))./besselk(

m,sqrt(V1^2-U1^2)); %特征方程

U1=fzero(fun1,1);%求解归一化径向相位常数U W1=sqrt(V1^2-U1^2); %求解归一化径向衰减常数W

beta01=sqrt(nc01^2*k0^2-(U1/r01_core)^2); %求解1号光纤的传播常数beta

A1=U1/V1*besselk(m,W1)/sqrt(besselk(m+1,W1)*besselk(m+1,W1))*sqrt(4*z0/n_clad/pi/r01_core^2); %求解1号光纤的A

%求1号单模光纤中的横向场分量Ey1

Ey01=@(r,theta) A1*cos(m*theta).*besselj(m,U1*r/r01_core)./besselj(m,U1); Ey02=@(r,theta) A1*cos(m*theta).*besselk(m,W1*r/r01_core)./besselk(m,W1);

%定义2号光纤的参数

r02_core=input('输入2号光纤的纤芯半径(单位为μm)=')*1e-06; %定义2号光纤纤芯半径 r02_clad=3*r02_core; %定义2号光纤的包层半径

nc02=input('输入2号光纤的纤芯折射率='); %定义2号光纤的纤芯折射率 delta2=(nc02-n_clad)/nc02; %相对折射率差 V2=nc02*k0*r02_core*sqrt(2*delta2); %归一化频率

while V2>2.40483,disp('2号光纤不是单模光纤,请重新输入');%判断2号光纤是否满足单模条件

r02_core=input('输入2号光纤的纤芯半径(单位为μm)=')*1e-06; %定义2号光纤纤芯半径 r02_clad=3*r02_core; %定义2号光纤的包层半径

nc02=input('输入2号光纤的纤芯折射率='); %定义2号光纤的纤芯折射率 delta2=(nc02-n_clad)/nc02; %相对折射率差 V2=nc02*k0*r02_core*sqrt(2*delta2); %归一化频率 end;

disp('2号光纤是单模光纤,输入正确') disp('开始计算，输出仿真图')

%由特征方程求2号光纤的归一化径向相位常数和归一化径向衰减常数 fun2=@(U2)

U2.*besselj(m+1,U2)/besselj(m,U2)-sqrt(V2^2-U2^2).*besselk(m+1,sqrt(V2^2-U2^2))./besselk(

m,sqrt(V2^2-U2^2));

U2=fzero(fun2,1); %求解归一化径向相位常数U W2=sqrt(V2^2-U2^2); %求解归一化径向衰减常数W

beta02=sqrt(nc02^2*k0^2-(U2/r02_core)^2); %求解2号光纤的传播常数beta %求解2号光纤的A

A2=U2/V2*besselk(m,W2)/sqrt(besselk(m+1,W2)*besselk(m+1,W2))*sqrt(4*z0/n_clad/pi/r02_core^2);

%求2号单模光纤中的横向场分量Ey2

Ey11=@(r,theta) A2*cos(m*theta).*besselj(m,U2*r/r02_core)./besselj(m,U2); Ey12=@(r,theta) A2*cos(m*theta).*besselk(m,W2*r/r02_core)./besselk(m,W2);

d=r01_clad+r02_clad;

delta_beta=beta01-beta02;%定义失配相位常数 F=@(r,theta)

conj(Ey01(r,theta)).*Ey12(sqrt(d^2+r.^2-2*d*r.*cos(theta)),theta)*(nc01^2-n_clad^2).*r*omega*eption_0/4;

K12=dblquad(F,0,r01_core,0,2*pi); 求解耦合系数K12 K21=conj(K12); %求解耦合系数K21

%求耦合光纤中的功率密度

c2=((delta_beta+sqrt(delta_beta^2+4*abs(K21)^2))*a1_0+2*K12*a2_0)/(2*sqrt(delta_beta^2+4*abs(K21)^2)); c1=a1_0-c2;

%求解1号光纤中的耦合计算

r_core1=linspace(-r01_core,r01_core,400); %rz平面上纤芯内某点的纵坐标 z1=linspace(0,1,1000); %耦合距离 [Z1,R_core1]=meshgrid(z1,r_core1);

theta0=0; %沿零度方位角取光纤的轴向剖面（rz平面）

r_clad1=linspace(r01_core,(d+r02_clad),500); %rz包层内某点的纵坐标

[Z11,R_clad1]=meshgrid(z1,r_clad1);

a1=c1*exp(i*(delta_beta+sqrt(delta_beta^2+4*abs(K21)^2))/2*z1)+c2*exp(i*(delta_beta-sqrt(delta_beta^2+4*abs(K21)^2))/2*z1); %求解1号光纤中场振幅

Sz_core1=A1^2/2/z0*(cos(m*theta0))^2*nc01.*(besselj(m,U1*abs(r_core1)'/r01_core)./besselj(m,U1)).^2*abs(a1).^2; %求解纤芯功率密度

Sz_clad1=A1^2/2/z0*(cos(m*theta0))^2*n_clad.*(besselk(m,W1*r_clad1'/r01_core)/besselk(m,W1)).^2*abs(a1).^2; %求解包层功率密度

%求解2号光纤中的耦合计算

r_core2=linspace(d-r02_core,d+r02_core,500); %rz平面上纤芯内某点的纵坐标 z2=linspace(0,1,1000); %耦合距离 [Z2,R_core2]=meshgrid(z2,r_core2); r_clad2=linspace(0,d-r02_core,500); [Z22,R_clad2]=meshgrid(z2,r_clad2);

a2=-(c1*(delta_beta+sqrt(delta_beta^2+4*abs(K21)^2))/2*exp(i*(-delta_beta+sqrt(delta_beta^2+4*abs(K21)^2))/2*z2)+c2*(delta_beta-sqrt(delta_beta^2+4*abs(K21)^2))/2*exp(-i*(delta_beta+sqrt(delta_beta^2+4*abs(K21)^2))/2*z2))/K12; %求解2号光纤中场振幅

Sz_core2=A2^2/2/z0*(cos(m*theta0))^2*nc02.*(besselj(m,U2*abs(r_core2-d)'/r02_core)./besselj(m,U2)).^2*abs(a2).^2; %求解2号光纤纤芯功率密度

Sz_clad2=A2^2/2/z0*(cos(m*theta0))^2*n_clad.*(besselk(m,W2*(d-r_clad2)'/r02_core)/besselk(m,W2)).^2*abs(a2).^2; %求解2号光纤包层功率密度

%耦合功率密度绘图

cmap=[linspace(1,0,256);linspace(1,0,256);zeros(1,256)]'; colormap(cmap); subplot(1,1,1);

surf(Z1,R_core1,Sz_core1),view(-10,60); hold on;

surf(Z11,R_clad1,Sz_clad1),view(-10,60);

subplot(1,1,1);

surf(Z2,R_core2,Sz_core2),view(-10,60); hold on;

surf(Z22,R_clad2,Sz_clad2),view(-10,60); shading flat;colorbar;axis tight;

xlabel('耦合距离z（m）');ylabel('r','Fontsize',13,'FontName','Times'); title('1号和2号光纤的耦合功率密度分布');

%绘制耦合光纤功率密度对比图 figure

cmap=[linspace(1,0,256);linspace(1,0,256);zeros(1,256)]'; colormap(cmap); subplot(2,1,2);

mesh(Z1,R_core1,Sz_core1),view(0,0); hold on;

mesh(Z11,R_clad1,Sz_clad1); shading flat;colorbar;axis tight;

xlabel('耦合距离z(m)');ylabel('r','Fontsize',13,'Fontname','Times'); title('1号光纤的耦合功率密度分布'); subplot(2,1,1);

mesh(Z2,R_core2,Sz_core2); hold on;

mesh(Z22,R_clad2,Sz_clad2),view(0,0); shading flat;colorbar;axis tight;

xlabel('耦合距离z（m）');ylabel('r','Fontsize',13,'FontName','Times'); title('2号光纤的耦合功率密度分布');

%绘制耦合效率曲线图 figure

Pout_1=abs(a1).^2; Pout_2=abs(a2).^2; enta=Pout_2; subplot(3,1,1);

plot(z1,Pout_1);xlabel('耦合距离z(m)');ylabel('直通臂功率Pout_1'); subplot(3,1,2);

plot(z1,Pout_2);xlabel('耦合距离z(m)');ylabel('耦合臂功率Pout_2'); subplot(3,1,3);

plot(z1,enta);xlabel('耦合距离z(m)');ylabel('耦合效率'); title('耦合效率——耦合距离曲线图');

本文来源：https://www.bwwdw.com/article/hv26.html