第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. 泛函取极值的必要条件
众所周知,一元函数取极值的条件是函数在该点对应的一阶导数为零,该极值点称为函数的驻点。泛函也有类似的性质。考虑泛函的全部变量函数都通过两个定点,即y(x0)=a,y(x1)=b,则可以证明,泛函取极值的必要条件为δJ=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


那么,取该泛函的变分为零,即可得到质点的真实运动轨迹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)


费马原理意味着光程对应的变分为零,即


δ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


此方程要有非零解,则其系数行列式必须为零,于是可以得到如下的方程


ε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(b)中的套层区域填充特殊材料,则就可以实现封闭区域外场的无扰动传播,从而达到隐形的目的。


图32变换光学原理示意图


当虚拟空间和物理空间都以直角坐标系进行表示时,根据麦克斯韦方程组的协变性,变换后的电磁材料参数有如下关系式


ε′(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)有非零解的条件是系数行列式的值为零,于是得到


|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))%对比二者差值是否为零

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)。
那么上面的两个方程就写成了六个微分方程。由式(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


解这六个微分方程,并考虑起始条件,就会得到光线在介质中的传播路径。
假设光线沿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)




图33球形隐形衣的射线追踪示意


解这个方程组就会得到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的值(初始条件)代入六个微分方程组就会解出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)


考虑初始条件,数值求解这六个方程就会得到射线在柱形隐形介质中的传播路径。
假设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的值(初始条件)代入六个微分方程组就会解出x,y,z的方程,这就是光线在柱形隐形衣中的传播方程,用MATLAB画出其轨迹就能得到光线在柱形隐形衣中的传播路径。
为了方便演示,仅考虑光线在隐形衣内部的传输情况。图35给出了z=0的平面上隐形衣外表面上分布的入射点和入射光线的示意图。如前所述,光线从真空中入射到隐形衣表面,入射光方向为p1=22,0,22。下面的代码给出了射线追踪的全过程。



图35柱形隐形装置中z=0平面上的入射点分布示意图



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)给出了光线映射后的状况。可以看出,光线在隐形装置内部发生弯曲,光线绕过中间圆形区域; 在隐形衣外部,光线未发生变化。


图37变换光学中光线的变化



但正如前面分析一样,物理空间中的电磁隐形装置可以通过射线追踪的方式加以理解。为此,选择隐形衣内部的8条光线进行射线追踪,如图38所示。注意到光线从真空入射到隐形衣表面上时,真空中的光线不发生变化。因此,可以利用隐形衣外表面上的入射点,获得虚拟空间中的光线。与此同时,还可以采用坐标变换的方法,即式(369),映射出物理空间的同一组光线。为了便于区别,映射过来的点都用空心圆点表示,如图38所示。从图中可以明确看出,采用变换光学理论和射线追踪方法
与采用空间点映射的方法,得到了相同的光射线。这个结论对于深入理解变换光学理论具有重要意义。


图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<r<R(384)


其中: R是光学黑洞的半径; n是常数,n取不同值,其对光线的吸收情况不同。
为了方便研究黑洞内部射线轨迹,选取简单的具有轴对称分布的圆形介质进行研究,则相应的哈密顿函数为


H=p2r2ε(r)+P2φ2ε(r)r2(385)


由哈密顿原理,可以得到如下四个微分方程


drdt=rnRnpr(386)

dφdt=rnr2Rnpφ(387)

dprdt=-nrn-1p2r2Rn-n-2rn-3p2φ2Rn(388)

dpφdt=0(389)


利用方程(386)~(389),并结合t=0时的初始条件,就可以求解得到光射线在“黑洞”中的轨迹方程。大多数情况下,可以使用MATLAB下的ode45等函数,利用数值计算得到射线轨迹。
3.4.2MATLAB代码
在下面的程序中,平行光线从x=-3的位置发出,照射在半径为2的人工“黑洞”结构上,进入“黑洞”


图39初始条件的确定示意图

内部,并被“吸入”黑洞中心。由于黑洞边缘对应的折射率为1,与空气匹配,因此在分界面处光线不发生偏折。但由于采用了柱坐标系,因此,需要对初始条件做深入考虑。图39给出了初始条件确定过程中所涉及的几何关系。假设黑洞半径为R,对应入射点的极角为φ,则初始位置为


r=R(390)

φ=π+arctanyx(391)


初始光学方向余弦为


pr=cosφ(392)

pφ=-Rsinφ(393)


光线入射到“黑洞”的初始位置式(390)、式(391)是容易确定的。现简单论述一下初始光学方向余弦的确定。
用柱坐标系下的单位矢量与式(380)两端做点积运算,并考虑梯度的性质,可知



er·=r=pr

reφ·=φ=pφ



由程函方程式(379)可知,||=n,且该梯度的方向就是光射线的方向,同时考虑在“黑洞”的边沿折射率为1,半径为R,于是


pr=er·=1·n·cosφ=cosφ

pφ=reφ·=R·1·n·cosφ+π2=-Rsinφ




基于ode45函数所编制的具体代码如下: 

clc; clear;%清屏,清内存

tspan=[0 25];%设置光线对应的参数所在的区间

no_lines=5;%设置光线数量

R=2;%黑洞半径

phi=linspace(0,2*pi,100);

x=R*cos(phi); y=R*sin(phi);

plot(x,y,'r','linewidth',3);hold on;%绘制黑洞对应的圆

theta=linspace(pi/2+0.2,pi-0.1,no_lines);%光线入射点在圆周上的分布角度

for lp=1:length(theta)

line([R*cos(theta(lp)),-3],[R*sin(theta(lp)),R*sin(theta(lp))],'linewidth',2);

%绘制空气中的光线

p1=cos(theta(lp));%径向的光学方向余弦

p2=-sin(theta(lp))*R;%角向的光学方向余弦

Y0=[R theta(lp) p1 p2];%初始位置和光学方向余弦

options=odeset('RelTol',1e-10,'AbsTol',[1e-10]);%求解参数设置

[T,Y]=ode45(@blackholeham,tspan,Y0,options);%利用ode45函数进行数值求解

plot(Y(:,1).*cos(Y(:,2)),Y(:,1).*sin(Y(:,2)),'linewidth',2);%绘制光线

end

axis equal;%设置x、y方向等比例绘制

在上述程序中,ode45函数引用了描述哈密顿正则方程组的函数blackholeham。该函数的实现代码具体如下: 

function dy=blackholeham(t,y)

dy=zeros(4,1);

R=2;

m=4;%折射率分布中的m值

dy(1)=(y(1)^m/R^m)*y(3);%式(3-86)

dy(2)=y(1)^(m-2)*y(4)/R^m;%式(3-87)

dy(3)=-m*y(3)^2*y(1)^(m-1)/2/R^m-(m-2)*y(4)^2*y(1)^(m-3)/2/R^m;%式(3-88)

dy(4)=0;%式(3-89)

改变上述函数中m的取值,可以实现不同的折射率分布,进而得到不同情况下的射线追踪情景。具体情况如图310所示。从图310(a)中可以看出,光线进入到“黑洞”内部之后,以螺旋形轨迹绕中心旋转并进入中心; 当m增大时,如图310(b)所示,螺旋的圈数减少,光线沿弧线射入圆心。因此,在圆心位置处布置吸波材料,即可将光线导入“黑洞”内部并吸收。


图310人工“黑洞”内部的射线追踪



3.5龙伯透镜、麦克斯韦鱼眼透镜和伊顿透镜的射线追踪实现
龙伯透镜、麦克斯韦鱼眼透镜和伊顿透镜是光学器件中最具代表性的三种透镜,并具有重要的实际应用价值。通过射线追踪分析光线在这些透镜内部的传输轨迹,可以准确了解透镜的特性。本节针对二维情况,分析三种透镜的哈密顿函数和正则方程,
采用MATLAB进行数值求解、射线追踪并绘制相应的射线轨迹。
3.5.1透镜中的哈密顿函数和正则方程
对程函方程,考虑二维情况,并采用直角坐标系,则有


x2+y2=p2x+p2y=n2(394)


因此,考虑采用如下形式的哈密顿函数


 H(x,y,px,py)=(p2x+p2y-n2(x,y))/2(395)


其中: n代表折射率; p代表光学方向余弦。
在射线追踪过程中,n是已知的关于x,y的函数。
由哈密顿函数,可得哈密顿正则方程为


dxidt=Hpii=x,y(396)

dpidt=-Hxii=x,y(397)


方程式(396)、式(397)并结合起始条件,即可求得光线的踪迹,完成射线追踪的过程。
3.5.2龙伯透镜的射线追踪分析
在二维情况下,龙伯透镜为一单位圆,对应的折射率分布为


n=2-x2-y2(398)


可以看出,在单位圆上,折射率为1,与周围的环境相匹配; 在圆心处,折射率取最大值2。对于放置在单位圆上的点光源,龙伯透镜能够将其发出并进入透镜的光线准直,并且沿光源和圆心连线的方向平行出射; 改变点光源的位置,平行光线的方向随之改变; 由于光路可逆,当平行光线入射到龙伯透镜时,该透镜能够把平行光线汇聚在圆周上,且圆心和焦点的连线与光线相平行。正是由于这些特性,龙伯透镜在卫星通信中具有广泛的用途。
把折射率n的表达式代入式(395)中,得哈密顿函数的表达式为


H(x,y,px,py)=(p2x+p2y+x2+y2-2)/2(399)


于是,对应的哈密顿正则方程为


dxdt=Hpx=px(3100)

dydt=Hpy=py(3101)

dpxdt=-Hx=-x(3102)

dpydt=-Hy=-y(3103)


对上述常微分方程组进行数值求解,即可得到龙伯透镜内的光线踪迹。下面的程序中,10条光线从x=2的位置发出,沿x轴负方向照射到龙伯透镜,然后经过透镜汇聚到(-1,0)的焦点位置。由于龙伯透镜与空气是匹配的,因此在交界面处波矢的方向不发生变化,光学方向余弦也保持不变,其数值为(-1,0)。

clc; clear;%清屏,清理内存

tspan=[0 1.6];%设置光线对应的参数区间

no_lines=10;%设置光线根数为10

a=1;%设置半径为1

y0=linspace(-a,a,no_lines);x0=sqrt(1-y0.^2);%设置光线照射到透镜上的位置坐标

x1=ones(1,length(x0))*2*a; y1=y0;%设置空气中光线的初始位置坐标

line([x0;x1],[y0;y1],'color','g'); hold on;%绘制空气中的光线,绿色

phi=linspace(0,2*pi,100);%以下语句用于绘制龙伯透镜

x=a*cos(phi); y=a*sin(phi);

plot(x,y,'r','linewidth',3);%绘制单位圆

for lp=1:length(x0)

Y0=[x0(lp) y0(lp) -1 0];%定义光线对应的起始位置和光学方向余弦

options=odeset('RelTol',1e-5,'AbsTol',[1e-6]);%设置求解参数

[T,Y]=ode45(@ham,tspan,Y0,options);%利用ode45函数进行数值求解

plot(Y(:,1),Y(:,2),'b'); hold on;%绘制对应的光线,蓝色

end

axis equal;%设置x、y方向等比例显示

射线追踪的结果如图311所示。从图311中可以看出,水平入射的光线进入龙伯透镜后发生弯曲,并被聚焦到龙伯透镜表面一点。


图311龙伯透镜对平行光线实现汇聚


在上述程序中,ode45函数引用了描述哈密顿正则方程的函数ham。该函数的实现代码如下: 

function dy=ham(t,y)

dy=zeros(4,1);%dy为一个列向量

dy(1)=y(3);%式(3-100)

dy(2)=y(4);%式(3-101)

dy(3)=-y(1);%式(3-102)

dy(4)=-y(2);%式(3-103)


3.5.3麦克斯韦鱼眼透镜的射线追踪分析
在二维情况下,麦克斯韦鱼眼透镜也是一个单位圆,对应的折射率分布为


n=2x2+y2+1(3104)


可以看出: 在单位圆上,折射率为1,与周围的空气相匹配; 在圆心处,折射率取最大值2。对于放置在单位圆上的点光源,麦克斯韦鱼眼透镜能够将光线汇聚到与光源相对应的单位圆上的另外一点,且二者连线经过圆心。改变点光源的位置,焦点位置随之改变。
把折射率n的表达式代入式(395)中,得哈密顿函数的表达式为



H(x,y,px,py)=12p2x+p2y-4(x2+y2+1)2(3105)


于是,哈密顿正则方程为


dxdt=Hpx=px(3106)

dydt=Hpy=py(3107)

dpxdt=-Hx=-8x(x2+y2+1)3(3108)

dpydt=-Hy=-8y(x2+y2+1)3(3109)


对常微分方程式(3106)~式(3109)进行数值求解,即可得到麦克斯韦鱼眼透镜内的光线踪迹。在下面的程序中,15根光线从单位圆上一点(1,0)发出,分别沿不同的方向射入麦克斯韦鱼眼透镜,经过透镜汇聚到(-1,0)的焦点位置。由于鱼眼透镜与空气是匹配的,因此在交界面处波矢的方向不发生变化,光学方向余弦也保持不变。

clc; clear;%清屏,清内存

tspan=[0 4];%设置光线对应的参数范围

no_lines=15;%15根光线

a=1;%透镜半径为1

y0=0;x0=1;%光源位置(1,0)

phi0=linspace(pi/2,3*pi/2,no_lines);%光源发出的光线的角度

phi=linspace(0,2*pi,100);%绘制透镜

x=a*cos(phi); y=a*sin(phi);

plot(x,y,'r-','linewidth',3); hold on;%绘制单位圆

for lp=1:no_lines

Y0=[x0 y0 cos(phi0(lp)) sin(phi0(lp))];%起始位置和光学方向余弦

options=odeset('RelTol',1e-7,'AbsTol',[1e-8]);%设置求解参数

[T,Y]=ode45(@ham,tspan,Y0,options);%数值求解光线轨迹

IN=inpolygon(Y(:,1),Y(:,2) ,[a*cos(phi)],[a*sin(phi)]);%判断光线是否在透镜内部

tmpx=Y(:,1); tmpy=Y(:,2);%光线坐标

plot(tmpx(IN),tmpy(IN),'b','linewidth',2); hold on;%绘制内部光线

end

axis equal;%x、y方向等比例绘制

射线追踪的结果如图312所示。从图312中可以看出,位于圆周上的光源,其发出的光线,经过麦克斯韦鱼眼透镜之后,汇聚到圆周的另外一点,二者的连线正好是圆的一条直径。


图312麦克斯韦鱼眼透镜射线追踪的结果


在上述程序中,ode45函数引用了描述哈密顿正则方程的函数ham。该函数的实现代码如下: 

function dy=ham(t,y)

dy=zeros(4,1);%dy为一个列向量

dy(1)=y(3);%式(3106)

dy(2)=y(4);%式(3107)

dy(3)=-8*y(1)./(y(1).^2+y(2).^2+1).^3;%式(3108)

dy(4)=-8*y(2)./(y(1).^2+y(2).^2+1).^3;%式(3109)

3.5.4伊顿透镜的射线追踪分析
在二维情况下,伊顿透镜为一单位圆,对应的折射率分布为


n=2x2+y2-1(3110)


由式(3110)可见: 在单位圆上,透镜折射率为1,与周围的环境相匹配; 在圆心处,折射率为无穷大,具有奇异性。对于入射的光线,伊顿透镜能够在内部将光线绕圆心弯曲,并沿与入射方向平行但相反的方向出射。这个特性与三个两两垂直的镜面所构成的角反射器非常相似。
把折射率n的表达式代入式(395)中,得伊顿透镜的哈密顿函数的表达式为


H(x,y,px,py)=12p2x+p2y-2x2+y2+1(3111)


于是,对应的哈密顿正则方程为


dxdt=Hpx=px(3112)

dydt=Hpy=py(3113)

dpxdt=-Hx=-x(x2+y2)3/2(3114)

dpydt=-Hy=-y(x2+y2)3/2(3115)


对上述常微分方程组进行数值求解,即可得到伊顿透镜内的光线踪迹。在下面的程序中,10根光线从x=2的位置发出,沿x轴负方向照射到伊顿透镜,然后经过透镜弯曲后返回。由于伊顿透镜与空气相匹配,因此在交界面处波矢的方向不发生变化,光学方向余弦也保持不变。

clc; clear;%清屏,清内存

tspan=[0 3];%设置光线对应的参数区间

no_lines=10;%设置10根光线

a=1;%透镜半径为1

y0=linspace(0,a*0.95,no_lines);x0=sqrt(1-y0.^2);%设置光线入射到透镜的位置

x1=ones(1,length(x0))*2*a; y1=y0;%设置空气中相应光线的起始位置

line([x0;x1],[y0;y1],'color','g','linewidth',2); hold on;%绘制空气中的光射线

phi=linspace(0,2*pi,50);%绘制透镜

x=a*cos(phi); y=a*sin(phi);

plot(x,y,'r','linewidth',3);%绘制透镜对应的单位圆

for lp=1:no_lines

Y0=[x0(lp) y0(lp) -1 0];%设置起始位置和相应的光学方向余弦

options=odeset('RelTol',1e-7,'AbsTol',[1e-8]);%设置求解参数

[T,Y]=ode45(@ham,tspan,Y0,options);%利用ode45函数进行数值求解

IN=inpolygon(Y(:,1),Y(:,2) ,[a*cos(phi)],[a*sin(phi)]);

%判断光线是否在透镜内部

tmpx=Y(:,1); tmpy=Y(:,2);%光线坐标

plot(tmpx(IN),tmpy(IN),'b','linewidth',2); hold on;
%绘制透镜内部光线

end

axis equal;%设置x、y方向等比例绘制

y0=-y0;x0=x0;%将入射点位置关于x轴对称

x1=ones(1,length(x0))*2*a; y1=y0;%设置空气对应光线的终止位置

line([x0;x1],[y0;y1],'color','g','linewidth',2); hold on;%绘制空气中反射回去的光射线

具体射线追踪的效果如图313所示。从图313中可以看出,水平光线从右侧(绿色光线,圆外上面部分)入射到伊顿透镜后,沿透镜中心弯曲并实现180°转向; 光线从透镜出射,依旧呈平行状态并沿原方向返回(红色光线,圆外下面部分)。


图313伊顿透镜对应的射线追踪结果


在上述程序中,ode45函数引用了描述哈密顿正则方程组的函数ham。该函数的实现代码具体如下: 

function dy=ham(t,y)

dy=zeros(4,1);%dy为一个列向量

dy(1)=y(3);%式(3112)

dy(2)=y(4);%式(3113)

dy(3)=-y(1)./(y(1).^2+y(2).^2).^(3/2);%式(3114)

dy(4)=-y(2)./(y(1).^2+y(2).^2).^(3/2);%式(3115)

3.6小结
本章从一般材料中的麦克斯韦方程组出发,推导了电磁波的色散方程; 以此为基础,构造光学哈密顿函数; 利用特征线方法,得到了光线的哈密顿正则方程组。此方程组可以使用MATLAB中的ode45函数进行数值求解,进而得到材料中的光线轨迹。基于光学中的费马原理,并类比分析力学的哈密顿原理,也可以得到类似的结论。此外,对边界分界面上的波矢匹配等给出了结果。射线追踪是分析许多光学元器件的基础,具有重要的理论和应用价值。