GRINレンズの光線追跡

[Maple Plot]

GRINレンズはGRaded INdexレンズから作られた言葉で屈折率が均一でない媒質(Inhomogenious Media)のことです。屈折率分布型レンズとも呼ばれます。昔から有名な例としてMaxwellレンズやLunebergレンズがあります。

屈折率が不均一の場合、光線は一般的に直進しないので光線追跡は簡単ではありません。

ここではMapleを使って不均一媒質内の光線追跡について紹介します。

>    restart:

不均一媒質の光線追跡

不均一媒質内の光線の軌跡は次の方程式の解によって与えられます。

>    Diff(n(r(s))*Diff(r(s),s),s)=grad(n(r));

Diff(n(r(s))*Diff(r(s),s),s) = grad(n(r))

ここで r = (x, y, z) は光線上の点の位置ベクトルです。n(r)は位置 r での屈折率で ds は光線にそった線分です。

時間に相当するパラメータtを導入して ds を次のように置き換えます。

>    dt=ds/n;

dt = ds/n

光線方程式は次のようになります。

>    diff(r(t),t$2)=n*grad(n);

diff(r(t),`$`(t,2)) = n*grad(n)

maple で光線追跡はこの微分方程式を解きプロットします。視覚的には DEplot3d という関数が適しています。

光線追跡例

では具体的な例で光線追跡を行ってみましょう。

微分方程式のツールパッケージを呼び出しておきます。

>    with(DEtools):

屈折率を与える関数nを定義します。z軸方向には変化しないケースです。

>    n:=(x,y)->2*exp(-(x^2+y^2));

n := proc (x, y) options operator, arrow; 2*exp(-x^2-y^2) end proc

z 軸から離れるにしたがって屈折率が下がっています。屈折率分布をプロットしてみましょう。

>    plot3d(n(x,y),x=-3..3,y=-3..3,shading=zhue);

[Maple Plot]

n(x, y, z) の勾配 (grad) を与える関数 nx, ny, nz を定義します。z 軸方向には変化しないので nz は 0 です。

>    nx:=unapply(diff(n(x,y),x),x,y);ny:=unapply(diff(n(x,y),y),x,y);nz:=0;

nx := proc (x, y) options operator, arrow; -4*x*exp(-x^2-y^2) end proc

ny := proc (x, y) options operator, arrow; -4*y*exp(-x^2-y^2) end proc

nz := 0

plots パッケージを呼びだしておきます。

>    with(plots):

Warning, the name changecoords has been redefined

屈折率分布を色で表現してプロットしてみましょう。

>    c1:=plot3d([0,x,y],x=-2..2,y=-2..2,color=n(x,y)):

>    c2:=plot3d([1,x,y],x=-2..2,y=-2..2,color=n(x,y)):

>    c3:=plot3d([2,x,y],x=-2..2,y=-2..2,color=n(x,y)):

>    c4:=plot3d([3,x,y],x=-2..2,y=-2..2,color=n(x,y)):

>    display({c1,c2,c3,c4},style=patchnogrid);

[Maple Plot]

同じ屈折率断面が z 軸方向につながっています。

上で示したように微分方程式は次のようになります。

>    deq:={diff(x(t),t$2)=n(x(t),y(t))*nx(x(t),y(t)),diff(y(t),t$2)=n(x(t),y(t))*ny(x(t),y(t)),diff(z(t),t)=1};

deq := {diff(z(t),t) = 1, diff(x(t),`$`(t,2)) = -8*exp(-x(t)^2-y(t)^2)^2*x(t), diff(y(t),`$`(t,2)) = -8*exp(-x(t)^2-y(t)^2)^2*y(t)}

いくつかの光線を表示するために初期値を複数定義します。

>    init:=NULL:

>    init:=init,[x(0)=0.1,y(0)=0.1,z(0)=0,D(x)(0)=0,D(y)(0)=0];

init := [x(0) = .1, y(0) = .1, z(0) = 0, D(x)(0) = 0, D(y)(0) = 0]

>    init:=init,[x(0)=0.1,y(0)=0.2,z(0)=0,D(x)(0)=0,D(y)(0)=0]:

>    init:=init,[x(0)=0.1,y(0)=0.3,z(0)=0,D(x)(0)=0,D(y)(0)=0]:

>    init:=init,[x(0)=0.1,y(0)=0.4,z(0)=0,D(x)(0)=0,D(y)(0)=0]:

>    init:=init,[x(0)=0.1,y(0)=0.5,z(0)=0,D(x)(0)=0,D(y)(0)=0]:

>    init:=init,[x(0)=0.1,y(0)=0.6,z(0)=0,D(x)(0)=0,D(y)(0)=0]:

>    init:=init,[x(0)=0.1,y(0)=0.7,z(0)=0,D(x)(0)=0,D(y)(0)=0]:

>    init:=[init];

init := [[x(0) = .1, y(0) = .1, z(0) = 0, D(x)(0) = 0, D(y)(0) = 0], [x(0) = .1, y(0) = .2, z(0) = 0, D(x)(0) = 0, D(y)(0) = 0], [x(0) = .1, y(0) = .3, z(0) = 0, D(x)(0) = 0, D(y)(0) = 0], [x(0) = .1, ...
init := [[x(0) = .1, y(0) = .1, z(0) = 0, D(x)(0) = 0, D(y)(0) = 0], [x(0) = .1, y(0) = .2, z(0) = 0, D(x)(0) = 0, D(y)(0) = 0], [x(0) = .1, y(0) = .3, z(0) = 0, D(x)(0) = 0, D(y)(0) = 0], [x(0) = .1, ...
init := [[x(0) = .1, y(0) = .1, z(0) = 0, D(x)(0) = 0, D(y)(0) = 0], [x(0) = .1, y(0) = .2, z(0) = 0, D(x)(0) = 0, D(y)(0) = 0], [x(0) = .1, y(0) = .3, z(0) = 0, D(x)(0) = 0, D(y)(0) = 0], [x(0) = .1, ...
init := [[x(0) = .1, y(0) = .1, z(0) = 0, D(x)(0) = 0, D(y)(0) = 0], [x(0) = .1, y(0) = .2, z(0) = 0, D(x)(0) = 0, D(y)(0) = 0], [x(0) = .1, y(0) = .3, z(0) = 0, D(x)(0) = 0, D(y)(0) = 0], [x(0) = .1, ...

2階の微分方程式なので微分係数(角度)の初期値も与えなければなりません。Dは微分演算子でこれを使って微係数の初期値を与えています。

準備ができましたのでいよいよ DEplot3d で光線の軌跡(解)をプロットします。

>    DEplot3d(deq,{x(t),y(t),z(t)},t=0..1,init,scene=[x(t),z(t),y(t)],
stepsize=.05,linecolor=t,orientation=[20,60],scaling=constrained,axes=frame);

[Maple Plot]

色はt(時間)の変化を表しています。

今度は2次元的な初期値のセットを用意してみます。

>    n:=8;

n := 8

>    init:=NULL:

>    for j to n do
  for k to n do
    init:=init,[x(0)=evalf(j/n),y(0)=evalf(k/n),z(0)=0,D(x)(0)=0,D(y)(0)=0]:
  end do:
end do:

>    init:=[init]:

>    DEplot3d(deq,{x(t),y(t),z(t)},t=0..1,init,scene=[x(t),z(t),y(t)],
stepsize=.05,linecolor=t,orientation=[-20,60]);

[Maple Plot]

今度は初期値を円状に配置してみましょう。

>    n:=20;

n := 20

>    init:=NULL:

>    for j to n do
  init:=init,[x(0)=evalf(cos(2*Pi*j/n)),y(0)=evalf(sin(2*Pi*j/n)),z(0)=0,D(x)(0)=0,D(y)(0)=0]:
end do:

>    init:=[init]:

>    DEplot3d(deq,{x(t),y(t),z(t)},t=0..4,init,scene=[x(t),z(t),y(t)],
stepsize=.05,linecolor=t,orientation=[-20,60]);

[Maple Plot]

今度は 1 点から光線が出ますが角度 (微係数) を配置してみます。

>    n:=20;

n := 20

>    init:=NULL:

>    for j to n do
  init:=init,[x(0)=0.1,y(0)=0,z(0)=0,D(x)(0)=0,D(y)(0)=evalf((2*j-n)/n)]:
end do:

>    init:=[init]:

>    DEplot3d(deq,{x(t),y(t),z(t)},t=0..2,init,scene=[x(t),z(t),y(t)],
stepsize=.05,linecolor=t,orientation=[-10,80]);

[Maple Plot]

ルーネベルグ・レンズ

もうひとつの例として有名なルーネベルグ・レンズをとりあげてみましょう。

屈折率は半径方向に次のように定義されます。

>    n:=r->sqrt(2-(r/R)^2);

n := proc (r) options operator, arrow; sqrt(2-r^2/R^2) end proc

R = 3 のケースを考えます。まずプロットしてみましょう。

>    R:=3:

>    plot(n(r),r=0..4);

[Maple Plot]

2 次元的な関数で定義してプロットしてみましょう。

>    n:=(x,y)->sqrt(2-(x^2+y^2)/R^2);

n := proc (x, y) options operator, arrow; sqrt(2-(x^2+y^2)/R^2) end proc

>    plot3d(n(x,y),x=-4..4,y=-4..4,shading=zhue,grid=[40,40]);

[Maple Plot]

色で屈折率の分布を表現してみましょう。

>    c1:=plot3d([x,y,0],x=-2..2,y=-2..2,color=n(x,y)):

>    c2:=plot3d([x,0,y],x=-2..2,y=-2..2,color=n(x,y)):

>    c3:=plot3d([0,x,y],x=-2..2,y=-2..2,color=n(x,y)):

>    display({c1,c2,c3},style=patchnogrid,axes=normal);

[Maple Plot]

実際の解のために x, y, z の関数として定義します。

>    n:=(x,y,z)->sqrt(2-(x^2+y^2+z^2)/R^2);

n := proc (x, y, z) options operator, arrow; sqrt(2-(x^2+y^2+z^2)/R^2) end proc

勾配 (grad) を計算します。unapply は関数化するたいへん便利な関数です。

>    nx:=unapply(diff(n(x,y,z),x),x,y,z);ny:=unapply(diff(n(x,y,z),y),x,y,z);nz:=unapply(diff(n(x,y,z),z),x,y,z);

nx := proc (x, y, z) options operator, arrow; -1/3*1/(18-x^2-y^2-z^2)^(1/2)*x end proc

ny := proc (x, y, z) options operator, arrow; -1/3*1/(18-x^2-y^2-z^2)^(1/2)*y end proc

nz := proc (x, y, z) options operator, arrow; -1/3*1/(18-x^2-y^2-z^2)^(1/2)*z end proc

微分方程式は次のようになります。

>    deq:={diff(x(t),t$2)=n(x(t),y(t),z(t))*nx(x(t),y(t),z(t)),diff(y(t),t$2)=n(x(t),y(t),z(t))*ny(x(t),y(t),z(t)),diff(z(t),t$2)=n(x(t),y(t),z(t))*nz(x(t),y(t),z(t))};

deq := {diff(x(t),`$`(t,2)) = -1/9*x(t), diff(z(t),`$`(t,2)) = -1/9*z(t), diff(y(t),`$`(t,2)) = -1/9*y(t)}

初期値を用意しましょう。

>    n:=10:
init:=NULL:

>    for k to n do
  init:=init,[x(0)=0.1,y(0)=evalf(2*k/n),z(0)=-2,D(x)(0)=0,D(y)(0)=0,D(z)(0)=1];
end do:

>    init:=[init];

init := [[x(0) = .1, y(0) = .2000000000, z(0) = -2, D(x)(0) = 0, D(y)(0) = 0, D(z)(0) = 1], [x(0) = .1, y(0) = .4000000000, z(0) = -2, D(x)(0) = 0, D(y)(0) = 0, D(z)(0) = 1], [x(0) = .1, y(0) = .600000...
init := [[x(0) = .1, y(0) = .2000000000, z(0) = -2, D(x)(0) = 0, D(y)(0) = 0, D(z)(0) = 1], [x(0) = .1, y(0) = .4000000000, z(0) = -2, D(x)(0) = 0, D(y)(0) = 0, D(z)(0) = 1], [x(0) = .1, y(0) = .600000...
init := [[x(0) = .1, y(0) = .2000000000, z(0) = -2, D(x)(0) = 0, D(y)(0) = 0, D(z)(0) = 1], [x(0) = .1, y(0) = .4000000000, z(0) = -2, D(x)(0) = 0, D(y)(0) = 0, D(z)(0) = 1], [x(0) = .1, y(0) = .600000...
init := [[x(0) = .1, y(0) = .2000000000, z(0) = -2, D(x)(0) = 0, D(y)(0) = 0, D(z)(0) = 1], [x(0) = .1, y(0) = .4000000000, z(0) = -2, D(x)(0) = 0, D(y)(0) = 0, D(z)(0) = 1], [x(0) = .1, y(0) = .600000...
init := [[x(0) = .1, y(0) = .2000000000, z(0) = -2, D(x)(0) = 0, D(y)(0) = 0, D(z)(0) = 1], [x(0) = .1, y(0) = .4000000000, z(0) = -2, D(x)(0) = 0, D(y)(0) = 0, D(z)(0) = 1], [x(0) = .1, y(0) = .600000...
init := [[x(0) = .1, y(0) = .2000000000, z(0) = -2, D(x)(0) = 0, D(y)(0) = 0, D(z)(0) = 1], [x(0) = .1, y(0) = .4000000000, z(0) = -2, D(x)(0) = 0, D(y)(0) = 0, D(z)(0) = 1], [x(0) = .1, y(0) = .600000...
init := [[x(0) = .1, y(0) = .2000000000, z(0) = -2, D(x)(0) = 0, D(y)(0) = 0, D(z)(0) = 1], [x(0) = .1, y(0) = .4000000000, z(0) = -2, D(x)(0) = 0, D(y)(0) = 0, D(z)(0) = 1], [x(0) = .1, y(0) = .600000...
init := [[x(0) = .1, y(0) = .2000000000, z(0) = -2, D(x)(0) = 0, D(y)(0) = 0, D(z)(0) = 1], [x(0) = .1, y(0) = .4000000000, z(0) = -2, D(x)(0) = 0, D(y)(0) = 0, D(z)(0) = 1], [x(0) = .1, y(0) = .600000...
init := [[x(0) = .1, y(0) = .2000000000, z(0) = -2, D(x)(0) = 0, D(y)(0) = 0, D(z)(0) = 1], [x(0) = .1, y(0) = .4000000000, z(0) = -2, D(x)(0) = 0, D(y)(0) = 0, D(z)(0) = 1], [x(0) = .1, y(0) = .600000...
init := [[x(0) = .1, y(0) = .2000000000, z(0) = -2, D(x)(0) = 0, D(y)(0) = 0, D(z)(0) = 1], [x(0) = .1, y(0) = .4000000000, z(0) = -2, D(x)(0) = 0, D(y)(0) = 0, D(z)(0) = 1], [x(0) = .1, y(0) = .600000...

さっそく光線をプロットしてみましょう。

>    DEplot3d(deq,{x(t),y(t),z(t)},t=0..5,init,scene=[x(t),z(t),y(t)],
stepsize=.05,linecolor=t,orientation=[20,60],scaling=constrained,axes=frame);

[Maple Plot]

前回と同様 2 次元的な初期値のセットを用意しましょう。

>    n:=8:init:=NULL:

>    for j to n do
  for k to n do
    init:=init,[x(0)=evalf(2*j/n),y(0)=evalf(2*k/n),z(0)=-2,D(x)(0)=0,D(y)(0)=0,D(z)(0)=1];
  end do:
end do:

>    init:=[init]:

これで光線を表示します。

>    DEplot3d(deq,{x(t),y(t),z(t)},t=0..5,init,scene=[x(t),z(t),y(t)],
stepsize=.05,linecolor=t,orientation=[-20,60]);

[Maple Plot]