实验一数据的可视化：光纤横截面上的场分布模拟

1、其他的二维绘图函数如ezplot，可用于绘制隐函数图、参数绘图，其调用格式:实验一数据的可视化：光纤横截面上的场分布模拟一、实验目的MATLAB 不仅在数值计算方面独占鳌头，而且在数据可视化方面也是功能强大。MATLAB 可以给出数据的二维、三维甚至四维的图形表现。通过对图形的线型、立面、色 彩、渲染、光线、视角等的控制，可以把数据的特征表现的淋漓尽致。本次实验拟通过对光 纤横截面上的模场分布进行模拟，使大家熟悉MATLAB常用的二维、三维绘图函数以及和绘图有关的命令，学会如何用不同的色彩来表示数值的大小，同时对单模、多模光纤的模场分布规律建立感性认识，为更好的学习后续专业课程打基础。二、MA

2、TLAB的常用绘图函数1. 二维绘图函数MATLAB中最常用的二维绘图函数是plot函数，其调用格式如下： plot ( X, Y, s )若X、Y为同维向量，则分别以X、Y中的元素为横、纵坐标绘制曲线；若X为向量，Y是有一维和X等维的矩阵，则绘制多个不同色彩的曲线，曲线数 等于Y的另一维；若X为矩阵，Y为向量，情况和上相同；若X，Y为同维矩阵，则以 X，Y对应的列元素为横、纵坐标绘制曲线族，曲线条 数等于矩阵的列数；s是字串，是用来指定线型、色彩的选项。各种可选项如下表所示曲裁线型曲拔濒色|标记符号选项选项童义选项意义选项童义iS项童义实线蓝色y蓝録色星号pen七曙rasp五角星3 -_ ?

3、虚霾k黑色点号W圆圈).J点线虹紫色红色叉号匚J点划线白色黄色叶Vflianond'0无线用一个打3向量任意指定 幻时虹绿蓝三原色1 * 9AA六角星< plot ( X1, Y1, si: X2, Y2, s2,)同时绘制多个曲线，每个绘线三元组( X, Y, s)的结构和作用与 plot(X, Y, s)相同，不 同的是三元组之间可以互不相关。例 1: t=0:0.01:2*pi; y仁t.*sin(t92) ; y2=exp(t)+cos(t42) ; plot (t, y1 ; y2, -r );% 或者 plot(t, y1, -r: t, y2, -.m )函数说明e

4、zplot( f',a, b)绘制隐函数f(x,y)=0的图形，横坐标范围a,bezplot( ',xmin,xmax,ymin,ymax)绘制隐函数 f(x,y)=0的图形，横坐标范围xmin, xmax, 纵坐标范围 ymin, ymaxezplot( fx ' fy ',tmin,tmax)参数绘图，绘出fx(t), fy(t), t的范围tmi n,tmax2 2 2例 2：画圆：x +y =R ( R=5 pm )a=5*1e-06; h1= ezplot('xA2 + yA2 - (5*1e-6)A2',-a,a); %hl 返回图形

5、的句柄(标识)2. 三维绘图函数 画三维曲线图一plot3函数调用格式：plot3(X1, Y1, Z1, 's1',X2, Y2, Z2, 's2'.J，除包含第三维之外，用法与plot函数相同。 画三维网格图一 mesh函数MATLAB的三维网格图的形成原理：在x-y平面上指定一个矩形区域，采用与坐标轴平行的直线将其分格，计算网格点上的函数值，即z轴的值，得到三维空间中的数据点；将这些数据点分别用处于x-z面或其平行平面内的曲线和处于y-z面或其平行平面内的曲线连接，即形成三维网格图。调用格式：mesh(X, Y, Z, C)X, Y,乙C为矩阵，X，丫存放

6、格点的横、纵坐标，可由meshgrid函数生成：X, Y=meshgrid(x,y), x,y为向量，指定画图的矩形区域。C是指定颜色的数据矩阵。mesh(X, Y, Z)(常用)(即 z轴的的值)成正比。)不着色，如例3图所示；而surf同上，只是颜色矩阵C=Z，故网格图的颜色和网格的高度例3：x=li nspace(-2,2,30);y=li nspace(-2,2,30); xx,yy=meshgrid(x,y);zz=xx.*exp(-xx.A2-yy.A2);mesh(xx,yy,zz)meshgrid函数将向量x,y指定的矩形区转化为二维 矩阵xx, yy(二维网格点坐标)，以供三

7、维绘图所需 画三维曲面图一surf函数调用格式与mesh函数一样，只是 mesh函数 所绘的网格图其格线着色，补片(格线间的空白区，网格面0所绘曲面图，其格线都为黑色，补片着不同的颜色(具体有颜色矩阵 C指定)，如下图所示。例 4： (1-4 条命令同例 3) ; surf(xx,yy,zz); colorbar;其中，colorbar命令的作用是在图形的右边加一竖直的颜色条，从中可见曲面图补片的颜色和z轴值成正比。对于surf图，可用shading命令进行着色的细 化处理：shading flat曲面图的某整个补片都着一种颜 色，颜色取自该补片四顶点数据中下标最小的那点 的颜色；shadin

8、g faceted 在flat用色基础上，再在补片 的四周勾画黑色网格线，是MATLAB的默认设置。 如上图所示。shadi ng in terp曲面图补片上各点的颜色由该补片四顶点的颜色经二次线性插值得到，这种方法用色细腻，但计算量大。如下图所示。由此可此，shadi ng flat或in terp都将去掉曲面图中 的黑色网格线。3. 矢量场图MATLAB 中绘制矢量图的函数有矢量场图(又称速度图)在工程中的应用很广泛。compass、feather、quiver、quiver3 这里仅介绍 quiver 禾口 quiver3。quiver函数绘制二维矢量场图，使用格式：quiver(X,

9、Y , U, V, S, 'SC '其中，X，Y，U, V为大小相同的向量或矩阵，X,Y指定每个矢量的起始点坐标，U，V指定每个矢量的x和y分量。S是指定所画箭头大小的参数，默认S=1，若S=2，则箭头大小加倍，若S=0.5,则缩小一倍，若 S=0，则箭头大小不可调整。第6个参数LSC是个字串，用于指定合法的线型和颜色。例5:绘制标量函数sin(x)*cos(y)的梯度(矢量)图xx,yy=meshgrid(-2:0.2:2,-2:0.2:2);zz=si n( xx).*cos(yy);u,v=gradie nt(zz);%计算梯度，输出u、v是二维梯度矢量的x、y分量qui

10、ver(xx,yy,u,v)结果如右图所示。quiver3函数quiver函数相同:用于三维矢量图绘制，除了包含第三维的信息之外，用法和quiver3（X, Y,乙 U, V , W, S, 'SC 'X, Y, Z指定矢量的起始点坐标，U, V, W指定矢量的x, y, z分量。4. 图形对象的颜色映像在MATLAB中，图形的颜色是作为第四维引入的，即在x,y,z轴的基础上引入一个颜色轴（caix），因此MATLAB图形的颜色其实是伪色彩 （pseudocolor）。每一种颜色可以使用R G B 基色三元行数组来表示。数值元素R、G、B在0, 1区间取值，分别表示红、绿、蓝三

11、基色的相对亮度。通过对 R、G、B大小的设置，可以调配出不同的颜色，例如：红色1 0 0;绿色0 1 0;蓝色0 0 1;灰色0.5 0.5 0.5;黑色0 0 0;白色1 1 1等。由m个不同的R G B三元行数组构成的 mX3的矩阵就定义了一种 色图（colormap，也 称调色板或 颜色表）。例如矩阵1 0 0;0 1 0;0 0 1定义了一个颜色表，它由红、绿、蓝三种颜色构成。MATLAB提供了一些现成的色图矩阵，如jet '代表的是一种蓝头红尾饱和色值的颜色表，hot '代表一种黑红白黄浓淡色的颜色表，这些色图矩阵可以在帮助系统中查到（help colormap命令）

12、。需要注意的是，代表颜色表的色图矩阵的维数必须是mX3,且元素值在01之间，m是矩阵的行数，矩阵的每一行定义了一种颜色，因此m代表了颜色表中颜色的数目，第 k（1 < kw m）行则对应着颜色表的第k个颜色。MATLAB的每个图形窗口（使用绘图命令之后自动弹出）只能有一个颜色表（色图）用于绘图。使用colormap命令可以指定图形窗口所使用的颜色表，其使用格式：colormap（cmap）cmap是一个mX3的矩阵。例如>> colormap（jet）%使用MATLAB提供的色图矩阵jet作为当前图形窗口的颜色表再如：>> colormap（ra nd（128,3

13、）%使用rand函数生成色图矩阵，定义了 128色的颜色表那么，在MATLAB中是如何使用一个颜色表中的颜色来进行绘图、着色？这就是所谓颜色映射示意图的颜色映射（Mapping）问题了。在许多绘图函数 中都提供了一个用来指定颜色的输入参数C,如前面的surf （X, Y,乙C）, C是一个和 X , Y , Z 同维的数值矩阵（称CData）,其作用是通过映射的方式来指定图形各个部分使用的颜色。注意C不是上面的色图矩阵 cmap，因为C矩阵的大 小不受mX3的限制而且元素值也不受01区间的限制。颜色映射的原理可以用右图说明：其中，cmin,cmax定义了颜色轴的范围 （如同xmin,xmax定

14、义了 x轴的范围），等于（或小于）cmin的数值映射到颜色表中的第一个颜色，等于（或大于）cmax的数值映射到颜色表的最后一个颜色，处于cmin和cmax之间的数值 C按照如下规律线性映射到颜色表中的其他颜色：k=fix（C-cmin）/（cmax-cmin）*m）+1 （不包括 C=cmax） （*）则映射到颜色表的第 k个颜色。函数fix（x） 取x的整数部分。除非特别指定，MATLAB 默认CData（指定颜色的数据矩阵，即参数C）的最小值为cmin，最大值为 cmax。也可通过 caxis（cmin,cmax）命令改变颜色轴的范围 cmin,cmax。因 此CData中元素值等于（或小

15、于）cmin者被映射为颜色表的第一个颜色，而等于（或大于）cmax者则被映射为颜色表的最后一种颜色，其他的则按上述规律映射。因此，图形的颜色实际上是利用绘图命令中输入的CData数据，按照其数值的大小以映射的方式来选择使用颜色表中的颜色。例6:三维球面绘制。X,Y,Z = sphere; %C = Z; colormap(jet) %绘制单位球的函数，输出为球面点的坐标指定当前图形窗口的颜色表surf（X,Y,Z,C）; axis equal %使 x,y,z图形如左下图所示。若在surf 命令后添加：轴刻度相同，否则球体看起来像椭球caxis(-1 0) %将颜色轴的范围由默认的-1,1 修

16、改为-1,0则上半球着单一的暗红色（右下图），因为z>0的值均被映射到颜色表jet的最后一个颜色0.5 0 0 了。plot、plot3 以及con tour和co ntour3（分别画二维和三维等高线函数）不使用颜色映像，即不接受Cdata参数，也不使用当前的颜色表，它们只使用列在plot颜色和线型表中的颜色（参见1页的表格）例7：>> t=0:0.01:2*pi;>> plot(t,si n( t),'l in ewidth',3,'color',1 0 0)%使用线宽为3像素的红色曲线描图5.二维伪彩色绘图一pcolor函数（

17、或非矩形）pcolor是专门绘制二维伪色彩图形的函数。它首先在二维平面内绘制出矩形单元阵列，然后使用输入的一般的使用格式：CData数据以映射的方式对各个矩形单兀指定颜色。 pcolor ( C )C为矩阵。首先利用 C矩阵中每相邻的4个元素绘制矩形单元，再将 C矩阵映射到颜 色表以指定每个单元的颜色。映射时，C的最后一行和最后一列不使用，元素C(i,j)的值指定第i行j列的矩形单元的颜色。例 8：>> colormap(1,0,0;0,1,0;0,0,1)%定义了只有 RGB三基色的颜色表>> C=1.5,1,-2.5,-1;0,-1.5,3,2;0.5,2.5,2,

18、-1;-2,-1,0,1>> pcolor(C)>> axis ij运行结果为图中右上角部分。映射规则是将 C中的最小值-2.5映射到颜色表的 第一个色红色，最大值 +3映射为第3个色蓝色。将C阵的最后1行和列划掉后，第 1行的3个元 素值分别指定第1行的3个矩形单元的颜色，按 映射规则，相应的单元色是B、G、R,以此类推。(例如，C(1,1)=1.5,按(*)式k=3，故映射为蓝色) pcolor (X, Y, C )X, Y, C是同维的矩阵，X, Y指定单元的顶点坐标。 例9:在极坐标系中绘制半径为5 pm的圆形伪色彩图>> a=5*1e-06;r=l

19、i nspace(0,a,8);>> theta=li nspace(0,2*pi,12);>> X=cos(theta ' )*r;Y= sin (theta%转化为直角坐标，X,>> C=randn(size(X); % C>> pcolor(X,Y,C);axis equal说明：可将极坐标下的%CDataG A I 0叩B吐丄'*r;Y扇形顶点坐标要和X, Y同维r维(径向)和theta维(角向)看成x维和y维，则C(i,j)方向指theta增加方向；指定第i行j列的扇形颜色，此时的行方向指r增加的方向，列MATLAB 中

20、，默认的着色方式是shad ing faceted，即每个单元着单一的颜色，而且勾画黑色格线，因此若在程序中使用shading flated命令，着色方式不变，但能去掉黑色格线。6.图形的控制与修饰(1)坐标轴的控制函数axis，调用格式如下axis(xmi n, xmax, ymi n, ymax, zmi n, zmax)用此命令可以控制坐标轴的范围。与axis相关的几条常用命令还有:axis equal严格控制各坐标的刻度使其相等axis square使绘图区为正方形axis on恢复对坐标轴的一切设置axis off取消对坐标轴的一切设置（轴不可见）axis tight使图形的绘图区紧

21、贴图形窗口（2）grid on在图形中绘制坐标网格grid off取消坐标网格.（3）图形的标注xlabel, ylabel, zlabel分别为x轴,y轴,z轴添加标注。title为图形添加标题。以上函数的调用格式大同小异，我们以 xlabel为例进行介绍：xlabel（'标注文本','属性1'，属性值1','属性2','属性值2',）这里的属性是标注文本的属性，包括字体大小（FontSize）、字体名（FontName）、字体粗细（Fo ntWeight）等。如：xlabel（'时间 t（单位秒）',&

22、#39;FontSize ' ,12, ' FontName', .'隶书','FontWeight' , ' Bold '）;（4）hold on在当前坐标轴设置的基础上，往当前图形窗口中添加多幅图形；MatLAB在绘图时，后画的图形在显示时会自动擦除原来的图形以显示新的图形，采用该命令后，将保留原有图形。常用于多幅图形的比较。hold off关闭上述功能。subplot图形窗口的分割函数格式：subplot（m,n,p），将图形窗口分割成 mx n的子窗口，并在第p个窗口绘制图形， 利用该命令可以在一个图形窗口中同时显

23、示多幅小图。三、实验原理如图所示，光纤的横截面有两层介质，里面的是芯层，外面的是包层，芯层折射率稍大于包层折射率，由于全反射，光主要约束在芯层中传播。为叙述方便，约定光的传播方向（i.e.光纤的轴向）为z轴方向，并先定义几个光纤参数1. 相对折射率差:g ni其中n1、n2分别是芯层和包层的折射率，这里只考虑n1、n2均为常数的阶跃折射率光纤。实际的光纤的厶很小，一般小于1%，因此也称为弱导光纤。2. 归一化径向相位常数 U和归一化径向衰减常数 WU uaJnkO2 - P2* （2）W = a2 -n；k0其中，a是光纤芯层半径，一般单模光纤的a=2 5 m。是光纤中电磁波的传播常数。k0=

24、2/ 0是真空中的波数。3.归一化频率VV = U2 W2 = n1 k0a(3)它概括了光纤的结构参数（a, , ni）以及工作波长（包含在 ko中），是一个重要的综合性参 数，光纤的很多特性都与归一化频率V有关。4.光纤横截面上的场的标量解在弱导光纤中的电磁波类似于一个横电磁波（TEM波），它的横向场分量 Et、Ht要远大于纵向场分量Ez、Hz，而且横向场分量是线偏振的，因此总是可以把电场的横向分量取为 y轴方向，即Et=Ey，它满足标量赫姆霍兹方程：2 2 2I Ey k°n Ey =0（4）其中，n是介质折射率（芯层中取ni，包层中取 门2）。在圆柱坐标系中利用分离变量法求解

25、（2）式，并使用电磁场的边界条件，可以解出电场的横向分量Ey：Ey 二 AcosimJm(Ur/a)/Jm(U)Km(Wr/a)/Km(W)其中，r是到光纤中心（z轴）的距离，m是整数,Jm和Km分别是m阶贝塞尔函数和 m阶第类修正的贝塞尔函数。7是极角。横向磁场只包含 Hx分量，由于横向电场与横向磁场的比等于波阻抗，故可由Ey计算出Hx：Hx十A/Z°)cosmnJm(Ur/a)/Jm(U)gKm(Wr/a)/Km(W)其中Z°=120二为真空波阻抗。再由麦氏方程组可求出场的纵向分量Ez、Hz：(U/nJm卅(Ur/a)/Jm(U ) sin(m+1)日匸 z-a/o.、

26、+JmUr/a)/Jm(U) sin(m-1)町 r"Ez 二(jA/ 2k0a)(W/nJKm勺(Wr/a)/Km(W) sin(m+1)。.-Km(Wr/a)/Km(W) sin(m1)日MaUJm+(Ur/a)/Jm(U) 9os(m+1)8出=(-jA/2k°aZ°) “-Jm(Ur/a)/Jm(U) cos(mT)日r" WKm出(Wr/a)/Km(W) cos(m + 1)日(.+Kmj(Wr/a)/Km(W) cos(m-1)日ra5.特征方程与LPmn模场解（5）（8）式中还有两个参量U和W未知，它们可由光纤的本征方程求出。利用电磁场的边

27、界条件Ez1=Ez2可以得到本征值方程:U Jmi(U)/Jm(U)=W Kmi(W)/Km(W)(9)or U Jm(U)/Jm(U) = -W Km(W)/Km(W)(10)（9）、（10）两式等价。利用（9）或（10）式可以求出U（或W），再由 式求出W（或U），然后代入（5）（8）式计算截面上的场分布。特征方程是超越方程，只能靠计算机数值求解。显然特征方程的解 U（或W），从而横向场解，与贝塞尔函数的阶数m有关。此外，特征方程的解也和贝塞尔函数的根有关。比如在大V值的情况下，特征方程可以简化为Jm（U ）=0，故此时U即是贝塞尔函数的根，但贝塞尔函数的根通常有若干个，如表2所示。所以U

28、=umn, umn代表m阶贝塞尔函数的根，n是根序号。可见 U的解除了和 m有关 之外也和贝塞尔函数的根序号n有关，从而光纤的场解（即模式）也和 m及n有关，不同（m,n）值给出不同的场解，即不同的模式。因此就将每一组（m,n）值对应的场解称为一个模式， 记作LPmn模，意思是线偏振模（Linearly Polarization），它表示弱导光纤中的模式基本上是线 偏振模。（m,n）也称为模指数。这里2m等于角向亮斑（极大值）数目，n等于径向亮斑数目。通常LPmn模的U值是在m-1阶贝塞尔函数的第 n个根和m阶贝塞尔函数的第 n个根之 间变化，这一点对于数值求解特征方程很重要，因为对于有多个根

29、的方程的求解，必须指定求解的区间，才能得到正确的根解。表2：贝塞尔函数的Jm根01234112. 4G53,8325,139& 3707.58825. 5207.016乱4177601L 0653B. 65410. 173It 62013.01514. 3134It 792IX 32114, 79616. 22017.6166. LPmn模的截止条件、归一化截止频率和单模传输条件当某一模式截止时，它就不能沿光纤传输了。可以用径向归一化衰减常数W来衡量某一模式是否截止。对于导波模（也称约束模，指约束在光纤中有效传输的模式），场在纤芯外是指数衰减的，保证了场约束在芯层，此时W2 >

30、0 ；当W2 < 0时，W为虚数，不再代表衰减常数，即场在纤芯外不再衰减了，也就不再约束在芯层中传播了，这种波叫做辐射波。而W2 =0正好是临界状态，以此作为导波模被截止的条件，此时的W记做Wc，相应的U和V值记做Uc和Vc，并将Vc称为归一化截止频率。由（3）式不难得到：Vc 二 U c （11）在截止条件（W=0）下，特征方程可以简化为Jm4（Uc）=0（12）即Uc是m-1阶贝塞尔函数的根。对于LP。!模（m=0 , n=1），截止时由（12）式得Vc=Uc=0,说明这种模式永不截止，是光纤中的最低阶模， 也称基模。第二个归一化截止频率较低的模式是LPii模式，称为二阶模，其归一化

31、截止频率Vc=Uc=2.40483。当光纤中只有基模工作而其他模式均被截止时，就称为单模光纤。要保证单模传输，必须要求光纤的归一化频率小于二 阶模LPii的截止频率 V，即让二阶模截止：V Vc(LP11H 2.40483(13)上式就是单模传输条件。根据(3)式，单模光纤的结构参数(a,二ni)就必须受到约束，比如纤芯半径a通常要较小。同样地，若要LPmn模式不被截止，则要求归一化频率V大于该模式的截止频率，即V Vc(LPmn)(14)7. LP mn模的功率密度以及功率平均坡印廷矢量Sav的大小就代表功率密度，使用S二！ Re( E H * )以及场解(5)(8)式,可以计算出光纤横截面

32、上的功率密度(Sav)近似为:Sz - A2 /(2Zo)cos2mv(15)n1Jm(Ur/a)/Jm(U)2 n2Km(Wr/a)/Km(W)2平均坡印廷矢量对横截面积分便得到光纤中的功率P：P =n2a2A2/(4Z°) (V2/)心 “(W) Km,W)/ K；(W) (16)令P=1 (归一化功率)，可算出以上各式的常数 A :21丄4Zo/(nr a2)2 (U/V) Km(W)/Km d(W) %1阿)2(17)8. 阶跃光纤中的模式数量以M表示光纤中能同时传播的模式数量，则有M -V2/n2k02a(18)可见M与归一化频率 V的平方成正比，因此，纤芯越大，折射率差越

33、大，光纤中传播的模 式数就越多。四、实验内容1. 计算横向电场Ey分布并将结果可视化显示分别计算单模光纤和多模光纤中LP34模的横向电场Ey分布并将结果可视化。给定实验参数：工作波长?0=1.55 m；纤芯折射率ni=1.5 ；相对折射率差A= 0.9 % ;纤芯半径a按(a)、(b)中有关要求输入，包层半径取为2a。要计算横向电场Ey,按式必须先计算归一化径向相位常数U和归一化径向衰减常数W(另一常数A由(17)式计算)，这就需要先求解特征方程(9)或者(10)式，可以先用(3)式将方程化为只含一个变量U。使用Matlab内建函数fzero可以求解特征方程：X=fzero(f, X0) 作用

34、是在 X0 附近搜索函数 f(x) 的零点 (根 )，返回值 X 为函数 f 的根。 f 是函数的句柄，可按照匿名函数的创建方法实现。对于单模光纤，基模LP01 模的 U 值介于0, 2.40483 之间，因此 X0 可取其间的某个值，比如1；对于多模光纤，当计算LP 34模时，根据第三部分的分析，其 U 值处于 2阶贝塞尔函数的第 4个根和 3阶贝塞尔函数的第 4 个 根之间，即区间 14.796, 16.22内，因此 X0 可取其间的某个值，如15。求出模式的U值后，使用(3)式算出 W，代入式计算Ey。(5)式中的m阶贝塞尔函数 Jm和m阶第二类修正的贝塞尔函数Km可以使用 Matlab

35、的内建函数besselj和besselk计算：J=besselj (NU, Z)计算Jm,输入参数NU指定贝塞尔函数的阶数，Z指定变量，J返回计算结果；K=besselk (NU, Z)计算Km,输入参数NU指定第二类贝塞尔函数的阶数，Z指定变量，K 返回计算结果；例如，besselj(0,U*r/a)计算 J°(Ur/a)的值。(a) 单模光纤单模光纤只有基模 LPoi(i.e.模指数m=0,n=1 )存在，应满足单模条件(13)式：V<2.40483， 再根据 (3)式,在给定实验参数的情况下,就要求芯径 a 有限制。输入参数 a 后,按上述方 法可计算出横向电场 Ey。按以下2种途径将数值结