第3章 射线追踪理论及其MATLAB实现 当物体的尺寸远远大于电磁波的波长时，可以采用几何光学的方法对电磁问题进行近似分析。在这种情况下，用光线来表示电磁波的传输。传统的透镜、成像等都是这方面的典型应用。当材料特性比较复杂时，如非均匀材料、各向异性材料等，其中光线的传播路径不是直线，通常要使用射线追踪理论研究光线在这些媒质中的传播特性。本章重点介绍媒质中的色散方程、光学哈密顿函数、哈密顿正则方程等，并介绍利用MATLAB数值分析光线路径的方法。作为实例，介绍了电磁隐形衣、光学“黑洞”以及三种经典透镜中光线的追踪方法和结果。 3.1从费马原理到哈密顿原理 媒质中光线的传播满足费马原理，基于费马原理即可得到光射线方程。这个原理与分析力学中的哈密顿原理有相似之处，因此，可以利用类比的方法得到光学哈密顿函数和哈密顿正则方程。 3.1.1泛函及其极值 在介绍费马原理之前，首先介绍一下与其密切相关的泛函的概念、泛函取极值的条件、欧拉拉格朗日方程以及哈密顿正则方程等。 1. 泛函的定义 设对于某一函数集合内的任意一个函数y(x)，按照一定的规则，都有另一个数J［y］与之对应，则称J［y］为y(x)的泛函。这里的函数集合，即泛函的定义域，通常要求y(x)满足一定的边界条件并且具有连续的二阶导数。这样的y(x)称为可取函数，通常也称为变量函数。泛函的形式可以多种多样，用积分定义的泛函应用最为广泛，即 J［y］=∫x1x0F(x，y，y′)dx(31) 其中： F是它的宗量的已知函数； y′表示对自变量x求导数。 从上述定义可以看出，只有当可取函数y(x)作为一个整体全部给定时(而不是仅已知个别点的函数值)，泛函才有确定的值，这与复合函数是完全不同的。 如果泛函的变量函数是一个二元函数，如u(x，y)，则可以定义泛函为 J［u］=SF(x，y，u，ux，uy)dxdy(32) 也可以定义多个函数所构成的泛函，如 J［x，y，z］=∫t1t0L(t，x，y，z，x·，y·，z·)dt(33) 其中： L为已知函数(之所以使用L，是为了与分析力学中的表达形式相一致)； x、y、z都是变量t的函数，且x·=dxdt，其他加点变量含义相同。 泛函(33)也可以简写为 J［r］=∫t1t0L(t，r，r·)dt(34) 如果t表示时间，x、y、z表示空间坐标，那么加点变量就是质点的速度矢量。在很多条件下，光线就可以看作一个质点以一定速度运动形成的轨迹。因此，掌握式(33)或式(34)所表示的泛函，对于理解光线和费马原理尤为重要。 2. 泛函的极值 与一元函数可以取极小值类似，泛函也可以在定义域的某个特定函数上取极小值。“当变量函数为y(x)时，泛函J［y］取极小值”的含义就是： 对于极值函数y(x)及其“附近”的变量函数y(x)+δy(x)，恒有J［y+δy］≥J［y］。 所谓函数y(x)+δy(x)在另一个函数y(x)的“附近”，指的是|δy(x)|<ε，ε是一个任意小的正数。有时还要求|δy′(x)|<ε。这里的δy(x)称为函数y(x)的变分。 3. 泛函取极值的必要条件 众所周知，一元函数取极值的条件是函数在该点对应的一阶导数为0，该极值点称为函数的驻点。泛函也有类似的性质。考虑泛函的全部变量函数都通过两个定点，即y(x0)=a，y(x1)=b，则可以证明，泛函取极值的必要条件为δJ=0，即泛函的一阶变分等于0，具体可以表示为 Fy-ddxFy′=0(35) 也称泛函在极值函数处是“驻定”的。式(35)称为欧拉拉格朗日方程，一般来讲，这是一个二阶常微分方程。 当泛函的变量函数为二元函数时，对应的必要条件为 Fu-xFux-yFuy=0(36) 对于多个函数构成的泛函∫t1t0L(t，x，y，z，x·，y·，z·)dt，其所对应的必要条件为 Lx-ddtLx·=0； Ly-ddtLy·=0； Lz-ddtLz·=0(37) 或者简写为 Lr-ddtLr·=0(38) 4. 哈密顿函数和正则方程 泛函理论具有重要的应用价值。在分析力学中，哈密顿最小作用原理就是基于泛函理论表示的。该原理表明： 在一个力学体系内，质点所有可能的运动中，真实运动使得哈密顿作用量取极值。如果定义体系的作用量泛函为 I=∫t1t0L(t，q，q·)dt=∫t1t0L(t，q1，q2，…，qn，q·1，q·2，…，q·n)dt 那么，取该泛函的变分为0，即可得到质点的真实运动轨迹q(t)，即 δI=δ∫t1t0L(t，q，q·)dt=0 其中： L=T-V是拉格朗日函数，表示体系动能T和势能V的差值； q表示广义坐标； q·表示广义速度。 类比式(37)和式(38)，即可得到其对应的欧拉拉格朗日方程如下： Lqi-ddtLq·i=0i=1，2，…，n(39) 这是关于广义坐标的二阶常微分方程组。 在分析力学中，还可以采用哈密顿正则方程来表示体系的运动。此时，对式(39)采用勒让德变换，定义体系的广义动量，即 pi=Lq·i i=1,2,…,n(310) 并定义系统所对应的哈密顿函数为 H(t，qi，pi)=∑iq·ipi-L i=1,2,…,n(311) 则下面的方程组成立： dqidt=Hpi dpidt=-Hqii=1,2,…,n(312) 方程组(312)称为哈密顿正则方程组。与方程组(39)相比，方程的个数增加了一倍，但是阶数降为一阶。求解式(312)即可得到质点的运动轨迹。 3.1.2费马原理和光线方程 众所周知，光线在真空中是沿直线传播的。光线为什么会沿直线传播？光在普通媒质中，尤其是非均匀介质或各向异性材料中是如何传播的？这些问题都可以由费马原理来解释。 费马原理表明： 在连接空间任意两点A、B的所有曲线中，光线实际经过的路线所对应的光程取极值。如果定义光程为如下的泛函，即 J［r］=∫BAn(r)ds(313) 费马原理意味着光程对应的变分为0，即 δJ=δ∫BAn(r)ds=0 其中： n是折射率分布； ds是沿曲线的弧长； n(r)ds表示光程的微元。 式(313)实际上表示的是一个三元函数的泛函。该泛函可取函数的集合为连接A、B的任意曲线，其取极值时的函数(曲线)正好是实际的光线。在真空或均匀介质中，由于折射率是1或常数，且两点之间线段最短，所以光线只能沿连接A、B的直线传播。当介质不均匀分布或者为各向异性材料时，必须利用泛函求极值的方法确定光线在其中的传播路径，这就是射线追踪的理论基础。 可以将上述泛函的极值曲线也就是光线用弧长表示为参数方程，则式(313)还可以写作 J［r］=∫BAn(r(s))drds·drdsds(314) 其中 drds=1 根据式(37)或式(38)，式(314)所对应的欧拉拉格朗日方程为 ddsndrds=n(315) 式(315)就是非均匀各向同性介质中的光线方程。 3.1.3光射线的哈密顿函数和正则方程 类比分析力学的方法，可以采用式(310)和式(311)得到光线对应的哈密顿函数。但简单的推演发现，哈密顿函数H≡0，这意味着直接套用分析力学的结论得不到相应的正则方程。通常情况下，这是由于光线参数选择不当造成的。如果不用弧长作为参数描述光线，例如采用空间坐标z作为光线的参数，就可以修正上述问题。此时光程的表达式为 J=∫BAn(x，y，z)1+x·2+y·2dz(316) 其中： x·=dxdz，y·=dydz。 上式中被积函数可以看作系统的拉格朗日函数L，其所对应的欧拉拉格朗日方程为 xn(x，y，z)=11+x·2+y·2ddzn(x，y，z)x·1+x·2+y·2 yn(x，y，z)=11+x·2+y·2ddzn(x，y，z)y·1+x·2+y·2 在这种情况下，继续应用勒让德变换，取 px=Lx·，py=Ly·(317) 则系统所对应的哈密顿函数可以表示为 H(x，y，px，py，z)= x·px+y·py-L(318) 这也就是 H(x，y，px，py，z)=-n2-p2x-p2y(319) 且有下面的常微分方程组成立 dxdz=Hpx dydz=Hpy dpxdz=-Hx dpydz=-Hy (320) 3.1.4其他形式的哈密顿函数和正则方程 实际应用中，考虑到3.1.3节中的分析方式对z坐标有特殊的对待，没有注意到空间坐标的平等性，因此，人们往往选择其他形式的哈密顿函数。例如，对弧长作为参数的光程进行参数变换，不妨取dt=nds，并考虑 drds=1，则式(313)变为 J［r］=∫BAn(r(s))drds·drdsds=∫BAn(r(s))drdtdtds·drdtdtdsds=∫BAn2(r(t))drdt·drdtdt 事实上，令上述泛函取极值(费马原理)，并采用勒让德变换，就可以得到 H(t，x，y，z，px，py，pz)=p2x+p2y+p2z4n2(x，y，z)(321) 同样，有如下的哈密顿正则方程成立： dxdt=Hpx dydt=Hpy dzdt=Hpz dpxdt=-Hx dpydt=-Hy dpzdt=-Hz (322) 由上面的推导过程可知，在光学里面，哈密顿函数不是唯一的。在后面的章节中，还会给出哈密顿函数的其他形式及其具体获得过程，并给出它们之间的联系。 3.2色散方程及其求解 前面给出了哈密顿函数的几种表达形式及其由来。事实上，在多数情况下，光学哈密顿函数还可以通过媒质中的色散方程获得。本节首先推导一般材料中的色散方程，然后利用特征线方程对其进行求解。在多数情况下，利用色散方程获得哈密顿函数，进而得到哈密顿正则方程的做法，更容易理解，也是常用的方法。 3.2.1非均匀材料的色散方程 考虑一种一般类型的各向异性材料，且材料分布不均匀，电磁参数随空间位置的不同而有一些缓慢的变化。假设电磁场在该材料中传输时具有如下类似均匀平面波的形式： E(r)=E0(r)· exp［-jk0(r)］ H(r)=H0(r)·exp［-jk0(r)］ (323) 其中： -k0(r)表示电磁波的相位； k0=2π/λ0是真空中的波数； λ0是真空中的波长； (r)是一个标量函数，表示波程。 对于位置矢量r，有 r=xex+yey+zez 显然，在各向同性材料中，波程可以表示为 (r)=∫n(r)·ds 对式(323)计算旋度，有 ×E(r)=×{E0(r)·exp［-jk0(r)］}(324) 考虑到 ×(fA)=f×A+f×A 则式(324)可以写为 ×E(r)={exp［-jk0(r)］}×E0(r)+［×E0(r)］exp［-jk0(r)］ =［-jk0(r)×E0(r)］exp［-jk0(r)］+［×E0(r)］exp［-jk0(r)］ 如果考虑波长很短，即λ0→0，则k0=2π/λ0→+∞； 且材料是空间位置的缓变函数，因此可以忽略上式中的第二项，而仅仅考虑第一项，于是有 ×E(r)=［-jk0(r)×E0(r)］exp［-jk0(r)］ 将结果代入麦克斯韦方程组，并考虑无源情况，即 ×H=jωD ×E=-jωB 则有 k0(r)×E0(r)=ωμ·H0(r)(325) 同理，有 k0(r)×H0(r)=-ωε·E0(r)(326) 式(325)和式(326)可以表示为 ［k0(r)］×-ωμ· +ωε·［k0(r)］×E0(r) H0(r)=0(327) 其中 ［k0(r)］×=k00- zy z0-x -yx0 假设材料的电磁参数形式为 μ=μ0μ100 0μ20 00μ3，ε=ε0ε100 0ε20 00ε3 则式(327)可以进一步详细表示为 0-k0zk0y-ωμ0μ100 k0z0-k0x0-ωμ0μ20 -k0yk0x000-ωμ0μ3 ωε0ε1000-k0zk0y 0ωε0ε20k0z0-k0x 00ωε0ε3-k0yk0x0E0x E0y E0z H0x H0y H0z=0 或 0-zy-Z0μ100 z0-x0-Z0μ20 -yx000-Z0μ3 1Z0ε1000-zy 01Z0ε20z0-x 001Z0ε3-yx0E0x E0y E0z H0x H0y H0z=0 此方程要有非零解，则其系数行列式必须为0，于是可以得到如下的方程： ε1μ14x+ε2μ24y+ε3μ34z+(ε1μ2+ε2μ1)2x2y+(ε1μ3+ε3μ1)2x2z+ (ε2μ3+ε3μ2)2y2z- ε1μ1(ε2μ3+ε3μ2)2x-ε2μ2(ε1μ3+ε3μ1)2y-ε3μ3(ε1μ2+ε2μ1)2z+ ε1ε2ε3μ1μ2μ3=0(328) 式(328)就是一般各向异性材料中的色散方程。 如果考虑非磁性材料，则上式中的磁导率均为1，于是有 ε14x+ε24y+ε34z+(ε1+ε2)2x2y+(ε1+ε3)2x2z+(ε2+ε3)2y2z- ε1(ε2+ε3)2x-ε2(ε1+ε3)2y-ε3(ε1+ε2)2z+ε1ε2ε3=0 (329) 对于单轴晶体，ε2=ε3，则式(329)变为 (2x+2y+2z-ε2) (ε12x+ε22y+ε22z-ε1ε2 )=0(330) 当材料为各向同性材料时，ε1=ε2=ε3，于是有 2x+2y+2z-ε=0(331) 上式也就是 ||2=n2(332) 式(332)就是大家所熟知的程函方程，其中n是各向同性材料的折射率。 上述各方程都可以用下面的函数形式来描述，即 H(x，y，z，，x，y，z)=0(333) 在3.2.3节将会对式(333)进行特征线分析。 关于式(328)的推导，也可以使用MATLAB进行符号计算。具体代码如下： syms phix phiy phiz Z0 mu1 mu2 mu3 eps1 eps2 eps3%定义符号变量 A=［0 -phiz phiy -Z0*mu1 0 0;%定义系数矩阵 phiz 0 -phix 0 -Z0*mu2 0; -phiy phix 0 0 0 -Z0*mu3; 1/Z0*eps1 0 0 0 -phiz phiy 0 1/Z0*eps2 0 phiz 0 -phix; 0 0 1/Z0*eps3 -phiy phix 0］; det(A)%计算系数行列式的值 MATLAB计算结果如下： eps1*mu1*phix^4+eps2*mu2*phiy^4+eps3*mu3*phiz^4+eps1*mu2*phix^2*phiy^2+eps2*mu1*phix^2*phiy^2+eps1*mu3*phix^2*phiz^2+eps3*mu1*phix^2*phiz^2+eps2*mu3*phiy^2*phiz^2+eps3*mu2*phiy^2*phiz^2-eps1*eps2*mu1*mu3*phix^2-eps1*eps3*mu1*mu2*phix^2-eps1*eps2*mu2*mu3*phiy^2-eps2*eps3*mu1*mu2*phiy^2-eps1*eps3*mu2*mu3*phiz^2-eps2*eps3*mu1*mu3*phiz^2+eps1*eps2*eps3*mu1*mu2*mu3 对比上面的结果和式(328)，可以发现二者完全一致。 3.2.2均匀材料的色散方程 需要说明的是，3.2.1节的推导过程是电磁场在几何光学情况下的近似； 但是当材料为均匀分布时，由于材料中有严格的均匀平面波存在，故上述推导严格成立，且有 k0(r)=k·r=kxx+kyy+kzz k0(r)=k0xex+k0yey+k0zez=k 此时，对于各向异性的非磁性材料，将方程(329)两边同乘以k60，则有 k21k4x+k22k4y+k23k4z+(k21+k22)k2xk2y+(k21+k23)k2xk2z+(k22+k23)k2yk2z- k21(k22+k23)k2x-k22(k21+k23)k2y-k23(k21+k22)k2z+k2 1k22k23=0 (334) 其中 k2i=k20εii=1，2，3 对于单轴晶体，ε2=ε3，则上式变为 k2x+k2y+k2z-k22k21k2x+k22k2y+k22k2z-k21k22=0(335) 所以 k2x+k2y+k2z=k22(336) 或者 k21k2x+k22k2y+k 22k2z=k21k22 也就是 k2xk22+k2y+k2zk21=1(337) 式(336)和式(337)对应的分别是寻常光(o光)和非寻常光(e光)的色散方程。 3.2.3色散方程的特征线法分析 在数学上，形如式(333)的微分方程可以用特征线的方法求解。将其重写为 H(x，y，z，，x，y，z)=0 可以简写为 H(r，，p)=0(338) 其中： p==(x，y，z)，称为光学方向余弦，并将其视为自由变量。 p的方向与等相位面垂直，也就是波矢的方向。由于k0p=k0=k，所以p也可以看作相对于真空中的波数k0做归一化的波矢量。 可以设想构造一条所谓的特征线，其可以用一个含参的曲线来表述，即r(τ)，且对于特征线上的各点，方程(338)恒成立。所以有 dHdτ=Hr·drdτ+Hr·drdτ+Hp·dpdτ=0(339) 也即 dHdτ=Hr+Hr·drdτ+Hp·dpdτ=0(340) 为了保证上式成立，可以考虑 drdτ=Hp(341) dpdτ=-Hr+Hr(342) 当然有 ddτ=r·drdτ=p·Hp(343) 于是 =∫p·Hpdτ(344) 由3.2.1节的推证过程可知，一般情况下，H中不显含，因此，上面的方程组可以简化为 drdτ=Hp dpdτ=-Hr (345) 通过对上面常微分方程组的求解，即可得到特征线，在几何光学里面，也就是光线。方程组(345)就是哈密顿正则方程。所以色散方程(338)的左侧函数即可看作哈密顿函数。 需要指出的是，色散方程(333)不是唯一的。例如，可以考虑 H′=f(H)=0(346) 其中： f为已知函数。 则式(333)和式(346)对应的是同一色散方程。可以证明，此时，对上述方程进行特征线分析，所得特征线不变，只不过特征线的参量进行了变换而已。 drdτ=H′p=dfdHHp dpdτ=-H′r=-dfdHHr (347) 可知，只要定义考虑dτ′=dτdfdH，则可以将原参量变换为新参量，而特征线的方程不变。 同样的道理，如果选择H′=f(r，p)H(r，，p)=0，其中f为任意一个已知的二元函数，则 drdτ=H′p=fpH+fHp=fHp dpdτ=-H′r=-frH-fHr=-fHr (348) 类似地，可以选择dτ′=fdτ，从而改变曲线参数的选择，但并不改变相应的特征线。上式中用到了H(r，，p)=0的条件。 例如，对于程函方程，即式(331)或式(332) 2x+2y+2z=n2 显然，有 2x+2y+2z4n2-14=0 那么根据特征线法的推证，可以选择哈密顿函数为等号左侧的函数，即 H=2x+2y+2z4n2-14 观察上面的结果，其与式(321)除了一个常数外，完全相同。换句话说，用色散方程的方法也可以得到哈密顿函数的表达式。事实上，当两个哈密顿函数相差一个常数时，二者对应的正则方程，即式（345）完全相同。后面章节中，也会采用类似手段舍去此类常数。 3.2.4媒质分界面的处理 在真空中光线沿一条直线传播，很容易进行射线追踪。但在两种媒质的分界面处，由于材料的特性发生了变化，光线一般会发生方向偏转。因此，在分界面处必须确定界面两侧p的大小和方向关系。在这种情况下，可以假设媒质1中的波矢p1(相对于真空中的波数归一化处理)已知，需要得到光线在媒质2中的波矢 p2。从3.2.1节中推导色散方程的过程来看，射线追踪时，一般是把电磁波看作准平面波，所以界面两侧必须满足切向波矢的连续性，也就是 (p1-p2)×n=0(349) 且在媒质2中，色散方程必然成立，即 H(p2)=0(350) 此外，考虑光线是从媒质1进入媒质2，还应该有 Hp·n>0(351) 图31光线在两种媒质分 界面的情况 其中： n是两种媒质分界面的单位法矢量，n=(nx，ny，nz)，其方向是从媒质1指向媒质2。图31给出了各个物理量的相互关系。 一般情况下，波矢p1已知时，通过式(349)、式(350)和式(351)，可以唯一确定p2。当光线进入媒质2后，就可以将入射点作为起始点，将p2作为起始波矢，从而继续进行射线追踪的过程。 3.3隐形衣中的射线追踪 电磁隐形衣是由特殊设计的电磁材料构成的，一般是非均匀各向异性的磁性材料，它可以包裹在物体周围，实现对目标的隐形。在电磁隐形衣内部，材料可以导引电磁波绕过被隐形的物体； 而在外部，隐形衣不对外界电磁场产生任何扰动。本节首先推导隐形衣中的哈密顿函数，由此得到光线的哈密顿正则方程，并利用数值方法对射线方程进行求解。射线追踪的结果表明了电磁隐形装置的有效性。 3.3.1隐形衣中的哈密顿函数 在电磁领域，可以使用变换光学理论来设计各种电磁隐形装置。变换光学本质上就是广义的坐标变换或者空间变换，它涉及两个变换空间，一个是虚拟空间，另一个是物理空间。虚拟空间也就是所希望模拟的电磁空间，而物理空间就是变换以后电磁器件实际上存在的空间。因为麦克斯韦方程组的协变性，可以通过在物理空间中填充某种经过变换后所得的具有特殊参数的材料来获得虚拟空间中的场分布效果，从而设计出种类繁多的特殊电磁器件，如隐形衣。 下面以图32所示的任意隐形斗篷作为示例简单说明变换光学的原理。 图32(a)所示是虚拟空间，此处以真空为例； 图32(b)所示是实际的物理空间，为填充有物理材料的实际空间。在虚拟空间，光线显然是沿着直线传播的。如果将虚拟空间中的任意点P(x，y，z)通过变换函数 X′=X′(x，y，z)映射为物理空间中的点P′(x′，y′，z′)，从而使图32(a)所示的一个闭合单连通区域从中间撕裂开，形成图32(b)所示的复连通区域； 同时，保持两个空间的外边界一致，从而保证边界的连续性，不会因为变换 图32变换光学原理示意图 的引入而使得边界处产生反射。在这种情况下，如果根据一定的变换规则在图32(b)中的套层区域填充特殊材料，则就可以实现封闭区域外场的无扰动传播，从而达到隐形的目的。 当虚拟空间和物理空间都以直角坐标系进行表示时，根据麦克斯韦方程组的协变性，变换后的电磁材料参数有如下关系式 ε′(x′，y′，z′)=A·ε(x，y，z)·ATdetA，μ′(x′，y′，z′)=A·μ(x，y，z)·ATdetA(352) 式(352)就是虚拟空间与物理空间材料分布的变换公式。其中，A是雅可比矩阵，A=X′/X=(x′，y′，z′)/(x，y，z)，体现了两个空间的变换关系。从式(352)可以看出，隐形材料对应的介电常数和磁导率都是对称的张量形式，二者都是空间坐标的函数且大小相同(相对值)。 接下来考虑光线在这种所谓的“变换媒质”中的传播特性。依然采用类似于3.2.1节的方法，假设光以准平面波的形式传播，并考虑缓变条件，则由麦克斯韦方程组可得 ［k0(r)］×-ωμ· +ωε·［k0(r)］×E0(r) H0(r)=0(353) 与3.2.1节不同的是，其中 μ=μ0N=μ0 n11n12n13 n12n22n23 n13n23n33， ε=ε0N=ε0 n11n12n13 n12n22n23 n13n23n33(354) 因此，在电磁隐形衣内部，相对介电常数张量和磁导率张量相同，用矩阵N表示，且 ［k0(r)］×=k00-zy z0-x -yx0(355) 如果用分量表示，式(353)也就是 0-k0zk0y-ωμ0n11-ωμ0n12-ωμ0n13 k0z 0 -k0x-ωμ0n12-ωμ0n22-ωμ0n23 -k0yk0x0-ωμ0n13-ωμ0n23-ωμ0n33 ωε0n11ωε0n12ωε0n130-k0zk0y ωε0n12ωε0n22ωε0n23k0z 0 -k0x ωε0n13ωε0n23ωε0n33-k0yk0x 0E0x E0y E0z H0x H0y H0z=0(356) 即 0 -z y-Z0n11-Z0n12-Z0n13 z 0 -x-Z0n12-Z0n22-Z0n23 -y x 0-Z0n13-Z0n23-Z0n33 1Z0n111Z0n121Z0n130-zy 1Z0n121Z0n221Z0n23z0-x 1Z0n131Z0n231Z0n33-yx0E0x E0y E0z H0x H0y H0z=0(357) 式(357)有非零解的条件是系数行列式的值为0，于是得到 |A|=(pNpT-|N|)2=0(358) 为简单计算，并根据式(358)，不妨考虑取 H=pNpT-|N|=0(359) 式(359)就是隐形衣等器件中光射线所满足的哈密顿函数。利用哈密顿正则方程，则容易得到射线的轨迹。 上述推证相对比较复杂，可以使用MATLAB符号工具箱进行辅助推演。代码如下： syms phix phiy phiz Z0 n11 n12 n13 n22 n23 n33%定义符号变量 A=［ 0 -phiz phiy -Z0*n11 -Z0*n12 -Z0*n13;%定义系数矩阵 phiz 0 -phix -Z0*n12 -Z0*n22 -Z0*n23; -phiy phix 0 -Z0*n13 -Z0*n23 -Z0*n33; 1/Z0*n11 1/Z0*n12 1/Z0*n13 0-phizphiy 1/Z0*n12 1/Z0*n22 1/Z0*n23 phiz0-phix; 1/Z0*n13 1/Z0*n23 1/Z0*n33 -phiyphix0］; n=［n11 n12 n13;%定义材料参数矩阵N n12 n22 n23; n13 n23 n33;］ factor(det(A))%计算行列式的值并做因式分解 B=［phix phiy phiz］*(n)*［phix phiy phiz］.'-det(n)%计算pNp^T-det(N) factor(det(A)-(B^2))%对比二者差值是否为0 MATLAB给出的系数矩阵的行列式的值如下： (-n33*n12^2+2*n12*n13*n23-2*n12*phix*phiy-n22*n13^2-2*n13*phix*phiz-n11*n23^2-2*n23*phiy*phiz-n11*phix^2-n22*phiy^2-n33*phiz^2+n11*n22*n33)^2 MATLAB最终计算得到的差值为0。因此，式(359)是严格成立的。 3.3.2球形隐形衣的射线追踪 1. 射线方程和起始条件 设计一个球形隐形装置，其内外半径分别为a和b，并考虑如下形式的变换： r′=b-abr+a，θ′=θ，φ′=φ(360) 则利用变换光学的理论，可以得到 N=bb-aI-2ar-a2r4xx(361) |N|=bb-a3r-ar2(362) 其中： I表示单位阵。 选择如下的哈密顿函数： H=12b-ab(pN pT-|N|)(363) 并将式(361)代入，那么球形隐形衣中的哈密顿函数为 H=12p·p-122ar-a2r4(x·p)2-12 b(r-a)r(b-a)2(364) 于是有 Hp=p-2ar-a2r4(x·p)x(365) Hx=-2ar-a2r4(x·p)p+3ar-2a2r6(x·p)2x-bb-aar-a2r4x(366) 其中： x为(x，y，z)； p为(px，py，pz)。 那么上面的两个方程就写成了6个微分方程。由式(365)和式(366)，得 dxdt=px-2ar-a2b4(xpx+ypy+zpz)x dydt=py-2ar-a2b4(xpx+ypy+zpz)y dzdt=pz-2ar-a2b4(xpx+ypy+zpz)z dpxdt=2ar-a2r4(xpx+ypy+zpz)px-3ar-2a2r6(xpx+ypy+zpz)2x+bb-aar-a2r4x dpydt=2ar-a2r4(xpx+ypy+zpz)py-3ar-2a2r6(xpx+ypy+zpz)2y+bb-aar-a2r4y dpzdt=2ar-a2r4(xpx+ypy+zpz)pz-3ar-2a2r6(xpx+ypy+zpz)2z+bb-aar-a2r4z 图33球形隐形衣的射线追踪示意 解这6个微分方程，并考虑起始条件，就会得到光线在介质中的传播路径。 假设光线沿x方向入射，入射点分布在球形隐形衣外表面一个半径为r0的圆上，即p1=(1，0，0)，如图33所示。球形隐形装置表面的单位法向矢量满足下面关系： nx=-cosθ ny=-sinθcosφ nz=-sinθsinφ (367) 其中的角度定义如图33所示。θ为矢径与x轴正方向的夹角，φ为矢径在yOz平面上的投影与y轴正方向的夹角。假设p2=(px，py，pz)，并考虑在外表面处r=b，则有 pzny=pynz pxnz-pznx=nz p2x+p2y+p2z-2ab-a2b4(xpx+ypy+zpz)2-1=0 Hp·n>0 (368) 解这个方程组就会得到px、py和pz的值，即 px0=-g+g2-4hu2h py0=ny0nx0(px0-1) pz0=nz0nx0(px0-1) 其中 h=1+ny0nx02+nz0nx02-mx20-my20ny0nx02-mz20 nz0nx02-2mx0z0nz0nx0- 2my0z0ny0nz0n2x0-2mx0y0ny0nx0 g=-2ny0nx02-2nz0nx02+2my20ny0nx02+2mz20nz0nx02+4my0z0ny0nz0n2x0+ 2mx0z0nz0nx0+2mx0y0ny0nx0 u=ny0nx02+nz0nx02-my20ny0nx02-mz20nz0nx02-2my0z0ny0nz0n2x0-1 把t=0时x、y、z、px、py和pz的值(初始条件)代入6个微分方程组就会解出x、y、z的解，这就是射线在介质中的传播轨迹，用MATLAB画出其轨迹就能得到射线在球形隐形介质中的传播路径。 2. MATLAB代码 基于上面的分析，可以编制如下的MATLAB程序，实现球形隐形衣内部的射线追踪。从图33可以看出，一束光线从x=-d处入射到x=-c对应的球面上，然后进入隐形衣内部。在绘制过程中使用均分圆周的方法确定30个入射点位置，如图33所示。经过隐形衣的导引，光线从x=c处出射。由于光线在真空中沿直线传播，所以在隐形衣外部采用直接绘制直线的方式绘制光线； 在隐形衣表面处，通过进行波矢匹配计算得到进入隐形衣内部的波矢起始值； 在隐形衣内部，采用哈密顿正则方程来描述光线路径，基于ode45函数进行数值求解； 光线出射后，依然采用绘制直线的方法绘制光线。具体代码如下： tspan=［0 20］;%设置光线对应的参数区间 no_lines=30;%设置光线条数 a=2;%隐形衣内径 b=4;%隐形衣外径 ［x,y,z］=sphere(100);%计算单位球面所对应的空间坐标 surf(a*x,a*y,a*z);hold on;%绘制内部的球面，并设置叠加绘图模式 surf(b*x,b*abs(y),b*z);%绘制外部球面。由于使用了绝对值，仅绘制半球面 shading interp; light('Position',［-1,0,0］);%设置光照位置 lighting gouraud;%设置光照模式 view(-40,17);%设置视角 r0=1;%入射点所在圆的半径为1 d=6; c=sqrt(b^2-r0^2);%入射点所在圆面到原点的距离 phi0=linspace(0,2*pi,no_lines);%在圆面上均匀选择30个点，作为入射点 theta0=pi-asin(r0/b);%计算入射点对应的theta角度 Y0=r0*cos(phi0);%入射点的y坐标 Z0=r0*sin(phi0);%入射点的z坐标 X0=linspace(-c,-c,no_lines);%入射点的x坐标 xx=linspace(-d,-d,no_lines); zz=Z0;yy=Y0;%光线起始点的坐标 for l=1:no_lines line(［X0(l),xx(l)］,［Y0(l),yy(l)］,［Z0(l),zz(l)］,'color','g','linewidth',1.5); %绘制每条光线 end for lp=1:no_lines x0=X0(lp);y0=Y0(lp); z0=Z0(lp);%每条光线的入射点坐标 ny0=-sin(theta0)*cos(phi0(lp)); nx0=-cos(theta0);nz0=-sin(theta0)*sin(phi0(lp)); %方向分量 m=(2*a*b-a^2)/b^4;%定义常数 h=1+ny0^2/nx0^2+nz0^2/nx0^2-m*y0^2*ny0^2/nx0^2-m*z0^2*(nz0^2/nx0^2)-2*m*x0*z0*(nz0/nx0)-2*m*y0*z0*(ny0*nz0/nx0^2)-m*x0^2-2*m*x0*y0*ny0/nx0; %定义h g=-2*ny0^2/nx0^2-2*nz0^2/nx0^2+2*m*y0^2*ny0^2/nx0^2+2*m*z0^2*nz0^2/nx0^2+4*m*y0*z0*nz0*ny0/nx0^2+2*m*x0*y0*ny0/nx0+2*m*x0*z0*nz0/nx0; %定义g u=ny0^2/nx0^2+nz0^2/nx0^2-m*y0^2*ny0^2/nx0^2-m*z0^2*(nz0^2/nx0^2)-2*m*y0*z0*(ny0*nz0/nx0^2)-1; %定义u kx0=(-g+sqrt(g^2-4*h*u))/(2*h);%计算波矢的起始x分量 ky0=(ny0/nx0)*(kx0-1);%计算波矢的起始y分量 kz0=(nz0/nx0)*(kx0-1);%计算波矢的起始z分量 Yi=［x0 y0 z0 kx0 ky0 kz0］;%起始矢量 options=odeset('RelTol',1e-10,'AbsTol',［1e-10］);%ode45函数的附加设置项 ［T,Y］=ode45(@rayham,tspan,Yi,options);%利用ode45函数进行计算 plot3(Y(:,1),Y(:,2),Y(:,3),'color','g','linewidth',1.5);%绘制光线 end X0=linspace(c,c,no_lines);%以下用于绘制出射的光线 xx=linspace(d,d,no_lines);zz=Z0;yy=Y0; for l=1:no_lines line(［X0(l),xx(l)］,［Y0(l),yy(l)］,［Z0(l),zz(l)］,'color','g','linewidth',1.5); end axis equal;%设置三个方向等比例显示 上面代码执行后所得到的结果如图34所示。 图34球形隐形衣射线追踪结果 另外，在ode45函数求解中，用到了描述哈密顿正则方程的函数rayham。其定义如下。读者可以结合哈密顿方程的表达式，进一步加深理解。 function dy=rayham(t,y); dy =zeros(6,1); a=2;b=4;%内外半径数值 r=(y(1).^2+y(2).^2+y(3).^2).^(1/2);%球坐标系下矢径的大小 s=y(1)*y(4)+y(2)*y(5)+y(3)*y(6);%定义波矢与矢径的点积 dy(1)=y(4)-((2*a*r-a^2)/r^4)*s*y(1);%哈密顿方程组的前三个方程 dy(2)=y(5)-((2*a*r-a^2)/r^4)*s*y(2); dy(3)=y(6)-((2*a*r-a^2)/r^4)*s*y(3); %哈密顿方程组的后三个方程 dy(4)= -((2*a*r-a^2)/r^4)*s*y(4)+((3*a*r-2*a^2)/r^6)*(s^2)*y(1)-((b/(b-a))^2)*((a*r-a^2)/r^4)*y(1); dy(4)=-dy(4); dy(5)=-((2*a*r-a^2)/r^4)*s*y(5)+((3*a*r-2*a^2)/r^6)*(s^2)*y(2)-((b/(b-a))^2)*((a*r-a^2)/r^4)*y(2); dy(5)=-dy(5); dy(6)=-((2*a*r-a^2)/r^4)*s*y(6)+((3*a*r-2*a^2)/r^6)*(s^2)*y(3)-((b/(b-a))^2)*((a*r-a^2)/r^4)*y(3); dy(6)=-dy(6); 3.3.3柱状隐形衣的射线追踪 1. 哈密顿正则方程 在柱坐标系下，选择柱状隐形衣的哈密顿函数为 H=12ρ-aρ(pNpT-|N|) 如果考虑采用如下的变换函数： ρ′=ρb-ab+a，φ′=φ，z′=z(369) 同样，利用变换光学理论，可以得到 N=ρ-aρT-2aρ-a2ρ3(ρ-a)ρρ+bb-a2ρ-aρZ 其中 T=100 010 001，Z=000 000 001，ρ=xex+yey 且有 |N|= bb-a2ρ-aρ 那么哈密顿方程可以写成下面的形式： H=12pTpT-122aρ-a2ρ4(ρ·p)2+12 b(ρ-a)ρ(b-a)(pZpT-1)(370) 对H分别求偏导，得 Hp=TpT-2aρ-a2ρ4(ρ·p)ρ+ b(ρ-a)ρ(b-a)2ZpT Hx=3aρ-2a2ρ6(ρ·p)2ρ-2aρ-a2ρ4(ρ·p)TpT+bb-a2aρ-a2ρ4(pZpT-1)ρ 写成分量形式，有 dxdt=px-2aρ-a2ρ4(xpx+ypy)x(371) dydt=py-2aρ-a2ρ4(xpx+ypy)y(372) dzdt= b(ρ-a)ρ(b-a)2pz(373) dpxdt=-3aρ-2a2ρ6(xpx+ypy)2x+2aρ-a2ρ4(xpx+ypy)px- bb-a2aρ-a2ρ4( p2z-1)x(374) dpydt=-3aρ-2a2ρ6(xpx+ypy)2y+2aρ-a2ρ4(xpx+ypy)py- bb-a2aρ-a2ρ4(p2z-1)y(375) dpzdt=0(376) 考虑初始条件，数值求解这6个方程就会得到射线在柱形隐形介质中的传播路径。 假设n是两种媒质分界面单位法矢量，则在圆柱表面，对应的内法线矢量可以表示为 nx=-cosφ ny=-sinφ nz=0 (377) 假设入射光方向p1=22，0，22，把n和p1代入式(349)~式(351)这三个方程来确定射线在介质中的方向余弦p，得到 nxpy-nypx=-22ny pz=22 p2x+p2y+p2z-2ab-a2b4(xpx+ypy)2-1=0 Hp·n>0 (378) 解这个方程组就能得到px、py和pz的初始值，即 px0=-g+g2-4hu2h py0=ny0nx0px0-22 pz0=22 其中 h=1+ny0nx02-mx20-my20ny0nx02-2mx0y0ny0nx0 g=-2ny0nx02+2my20ny0nx02+2mx0y0ny0nx0 u=12ny0nx02-my202ny0nx02-12 m=2ab-a2b4 2. MATLAB代码 同球形隐形衣一样，把t=0时x、y、z、px、py和pz的值(初始条件)代入6个微分方程组就会解出x、y、z的方程，这就是光线在柱形 隐形衣中的传播方程，用MATLAB画出其轨迹就能得到光线在柱形隐形衣中的传播路径。 图35柱形隐形衣中z=0平面 上的入射点分布示意图 为了方便演示，仅考虑光线在隐形衣内部的传输情况。图35给出了z=0的平面上隐形衣外表面上分布的入射点和入射光线的示意图。如前所述，光线从真空中入射到隐形衣表面，入射光方向为p1=22，0，22。下面的代码给出了射线追踪的全过程。 tspan=［0 150］;%定义光线参数区间 no_lines=20;%定义光线根数 a=2; b=4;%内外半径 phi=linspace(0,2*pi,no_lines);%角度分布，用于判断点是否在圆内，见plot3语句 ［x,y,z］=cylinder(1);%高度为1，半径为1的柱面坐标 surf(a*x,a*y,8*z);%绘制隐形衣内部表面，半径为2，高度为8 hold on;%设置为叠加绘制模式 surf(b*x,b*abs(y),8*z);%绘制隐形衣外表面，使用绝对值仅绘制曲面一半 shading interp; light('Position',［-1,0,0］);%光源位置 lighting gouraud;%设置光照模式为gouraud colormap('Spring');%设置颜色图为Spring view(-14,17);%设置视角 phi00=linspace(pi/2,3*pi/2,no_lines);%入射点对应的角度分布 y00=b*sin(phi00);z00=linspace(0,0,no_lines); x00=b*cos(phi00);%入射点坐标 for lp=4:length(x00)-2 x0=x00(lp);y0=y00(lp); phi0=phi00(lp);z0=z00(lp);%选择入射点追踪 nx0=-cos(phi0); ny0=-sin(phi0);nz0=0;%法线方向 c=sqrt(2);%以下定义常数 m=(2*a*b-a^2)/b^4; h=1+ny0^2/nx0^2-m*x0^2-m*y0^2*ny0^2/nx0^2-2*m*x0*y0*ny0/nx0; g=-c*(ny0^2/nx0^2)+c*m*y0^2*(ny0^2/nx0^2)+c*m*x0*y0*ny0/nx0; u=ny0^2/(2*nx0^2)-(1/2)*m*y0^2*(ny0^2/nx0^2)-1/2; kx0=(-g+sqrt(g^2-4*h*u))/(2*h); %计算波矢初值 ky0=(ny0/nx0)*(kx0-c/2); kz0=c/2; Y0=［x0 y0 z0 kx0 ky0 kz0］;%哈密顿方程初值 options=odeset('RelTol',1e-10,'AbsTol',［1e-10］);%求解参数设置 ［T,Y］=ode45(@cylinderham,tspan,Y0,options);%数值求解 IN=inpolygon(Y(:,1),Y(:,2) ,［b*cos(phi)］,［b*sin(phi)］);%判断光线是否在隐形衣内部 tmpx=Y(:,1); tmpy=Y(:,2);tmpz=Y(:,3);%光线坐标 plot3(tmpx(IN),tmpy(IN),tmpz(IN),'b','linewidth',2); hold on;%绘制内部光线 end axis equal;%同比例显示 图36(a)给出了柱状隐形衣的光线追踪的结果，图36(b)给出了俯视观察的情景。 图36柱状隐形衣中的射线追踪结果 上面代码中所涉及的哈密顿方程组是通过使用函数cylinderham来具体描述的。其对应的代码如下： function dy=cylinderham(t,y); a=2;b=4;%隐形衣内外半径 r=(y(1)^2+y(2)^2)^(1/2);%柱坐标系下的变量定义 s=y(1)*y(4)+y(2)*y(5);%定义表达式 dy =zeros(6,1); dy(1)=y(4)-((2*a*r-a^2)/r^4)*s*y(1)+0;%式(371) dy(2)=y(5)-((2*a*r-a^2)/r^4)*s*y(2)+0;%式(372) dy(3)=0-((2*a*r-a^2)/r^4)*s*0+(b*(r-a)/(r*(b-a)))^2*y(6);%式(373) dy(4)=-((2*a*r-a^2)/r^4)*s*y(4)+((3*a*r-2*a^2)/r^6)*(s^2)* y(1)+((b/(b-a))^2)*((a*r-a^2)/r^4)*y(1)*(y(6)^2-1); dy(4)=-dy(4);%式(374) dy(5)=-((2*a*r-a^2)/r^4)*s*y(5)+((3*a*r-2*a^2)/r^6)*(s^2)* y(2)+((b/(b-a))^2)*((a*r-a^2)/r^4)*y(2)*(y(6)^2-1); dy(5)=-dy(5);%式(375) dy(6)=-((2*a*r-a^2)/r^4)*s*0+((3*a*r-2*a^2)/r^6)*(s^2)*0+((b/(b-a))^2)*((a*r-a^2)/r^4)*0*(y(6)^2-1); dy(6)=-dy(6);%式(376) 3. 通过射线追踪来深入理解变换光学理论 变换光学理论通过变换函数实现虚拟空间与物理空间的映射，同时也把虚拟空间中的电磁场映射为物理空间内的电磁场。依旧考虑柱状隐形装置。假设虚拟空间为真空，则光线沿直线传播，如图37(a)所示。在实现柱状隐形衣映射的过程中，如式(369)所示，这些光线上的点也同样被映射到物理空间，形成物理空间的新的光线。图37(b)给出了光线映射后的状况。可以看出，光线在隐形装置内部发生弯曲，光线绕过中间圆形区域； 在隐形衣外部，光线未发生变化。 但正如前面的分析一样，物理空间中的电磁隐形装置可以通过射线追踪的方式加以理解。为此，选择隐形衣内部的8条光线进行射线追踪，如图38所示。注意到光线从真空入射到隐形衣表面上时，真空中的光线不发生变化。因此，可以利用隐形衣外表面上的入射点，获得虚拟空间中的光线。与此同时，还可以采用坐标变换的方法，即式(369)，映射出物理空间的同一组光线。为了便于区别，映射过来的点都用空心圆点表示，如图38所示。从图中可以明确看出，采用变换光学理论和射线追踪方法 与采用空间点映射的方法得到了相同的光射线。这个结论对于深入理解变换光学理论具有重要意义。 图37变换光学中光线的变化 图38光线的坐标映射与射线追踪的对比 3.4人工“黑洞”中的射线追踪 所谓人工“黑洞”，是指具有特定电磁参数和分布的特殊材料，当光线入射到这种材料表面时，反射非常小，绝大部分能量进入材料内部并且被吸收。这个特性与天文学中的“黑洞”相类似，因此被称为人工“黑洞”。基于人工“黑洞”理论，可以设计性能优异的吸波材料等。 3.4.1“黑洞”中的哈密顿函数和正则方程 由3.2.1节可知，在各向同性的非磁性材料中，程函方程可以表示为 ||2=n2=ε(379) 如果考虑二维情况，并考虑圆柱坐标系，有 =err+eφrφ=erpr+eφpφr(380) 则程函方程可表示为 p2r+pφr2=ε(381) 也就是 p2r2ε+p2φ2r2ε=12(382) 因此，可以考虑二维的哈密顿函数为 H(r，p)=p2r2ε+p2φ2r2ε(383) 人工光学“黑洞”的介电常数ε满足下面的关系式： ε(r)=Rrn，0