Mathematica/ミラー光学系 の履歴(No.3)
更新概要†
ミラー光学系の収差を Mathematica を使って計算してみる。
前提†
まずは収差のない場合のおさらいから。
焦点距離と曲率†
曲率が で表わされる凹面レンズの の位置に垂直に平行光が入射するとき、
反射面の傾きは なので、小角近似では
入射光は だけ曲がって反射されるので、
すなわち、焦点距離は となる。
曲率と曲率半径†
原点を通る半径 の円
を について の小さなところで解くと、
したがって、曲率半径 に対して、 となり、焦点距離は
となる。
点光源からの光を集光する†
凹面鏡の中心軸上、 の距離に置かれた点光源から、 凹面鏡上の の位置に到達する光が法線から だけずれているとすると、
の関係がなりたつ。平面鏡であればこの光は反射角 で反射されるが、 凹面鏡の傾きが であれば、実際の反射角は となる。
このようにして中心軸上に集められる光の焦点位置 は、
となるから、
の関係が得られる。
斜め入射の場合†
中心線に対して の角度で入射する平行光に対する焦点位置を考える。
中心位置で反射される光は反射角 で反射される
位置 で反射される光は反射角 で反射される
これらの光の為す角は だから、色を付けた三角形を考えることで、
となって、 のファクターだけ焦点距離が短くなる。
一方で入射面に垂直な方向の焦点距離は変わらないから、 入射面内での光の収束位置と入射面に垂直な方向への光の収束位置とがずれてしまうことになる。 すなわち凹面鏡への斜め入射は大きな収差(非点収差)を生むことになる。
この収差を無くすため、斜め入射で用いるための凹面鏡は入射面内の曲率半径と、 垂直方向の曲率半径を だけ違えて設計され、トロイダル凹面鏡などと呼ばれる。
トロイダル凹面鏡の面は数学的にはドーナツ面(トーラス面)の側面として表せる。
数値計算で確かめてみる†
以下は Wolfram Mathematica 9.0 による結果です。
まずはトーラス面(ドーナツ面、トロイダル面)を作成してみます。
パラメータ表示
&math( \left\{\begin{array}{l} x=o_x+r_1\cos s+r_2\cos t\cos s\\ y=o_y+r_1\sin s+r_2\cos t\sin s\\ z=o_z+\phantom{r_1\sin s+}\ \,r_2\sin t\\ \end{array}\right . );
と、
方程式表示
(* ドーナツ面のパラメータ表示 *) (* r1=0 で球面を表わせる *) torusPar[m_,s_,t_]:= Module[{o,r1,r2,ss,tt}, {o,r1,r2,srange,trange}=m; ss = s srange[[1]] + (1-s) srange[[2]]; tt = t trange[[1]] + (1-t) trange[[2]]; o + { r1 Cos[ss]+r2 Cos[tt] Cos[ss], r1 Sin[ss]+r2 Cos[tt] Sin[ss], r2 Sin[tt] } ] (* ドーナツ面の方程式表示 *) (* r1=0 で球面を表わせる *) torusEq[m_] := Module[{o, r1, r2, sr, tr, d}, {o, r1, r2, sr, tr} = m; (Sqrt[(x - o[[1]])^2 + (y - o[[2]])^2] - r1)^2 + (z - o[[3]])^2 == r2^2 ] (* どちらも引数 m には、 - o: 中心位置 - r1: 半径1 ( = 0 なら球になる) - r2: 半径2 - sr: パラメータ s の範囲 - tr: パラメータ t の範囲 をリストにした物を与える *)
試しに表示してみると、
With[ {m = {{0, 0, 0}, 3, 1, {0, 2 Pi}, {0, 2 Pi}}}, { ParametricPlot3D[torusPar[m, s, t], {s, 0, 1}, {t, 0, 1}], ContourPlot3D[Evaluate[torusEq[m]], {x, -4, 4}, {y, -4, 4}, {z, -4, 4}] }]
うまくいっている。
範囲を指定することでミラーの領域を限ることもできる
With[{m = {{-4, 0, 0}, 0, 4, {-Pi/8, Pi/8}, {-Pi/8, Pi/8}}}, ParametricPlot3D[ torusPar[m, s, t], {s, 0, 1}, {t, 0, 1}, PlotRange -> {{-2, 2}, {-2, 2}, {-2, 2}}] ]
次に直線のパラメータ表示と方程式表示。
(* 直線のパラメータ表示 *) linePar[l_, t_] := Module[{p, v}, {p, v} = l; p + t*v ] (* 直線の方程式表示 *) lineEq[l_] := Module[{p, v}, {p, v} = l; If[v[[1]] == 0, If[v[[2]] == 0, If[v[[3]] == 0, {x, y, z} == p, {x, y} == {p[[1]], p[[2]]} ], If[v[[3]] == 0, {x, z} == {p[[1]], p[[3]]}, x == p[[1]] && (y - p[[2]])/v[[2]] == (z - p[[3]])/v[[3]] ] ], If[v[[2]] == 0, If[v[[3]] == 0, {y, z} == {p[[2]], p[[3]]}, y == p[[2]] && (x - p[[1]])/v[[1]] == (z - p[[3]])/v[[3]] ], If[v[[3]] == 0, z == p[[3]] && (x - p[[1]])/v[[1]] == (y - p[[2]])/v[[2]], (x - p[[1]])/v[[1]] == (y - p[[2]])/v[[2]] == (z - p[[3]])/v[[3]] ] ] ] ] (* どちらも引数 l には、 - p: 開始位置 - v: 方向ベクトル をリストにした物を与える *)
方程式表示はちょっと面倒だ。。。
これらのトーラス面と直線を使って計算していく。
(以下暫定版で、まだ動きません・・・)
(* 直線とトーラス面との交点を求めてトーラス面のパラメータで返す *) intersection[l_, m_] := Module[{points, params, o, r1, r2, srange, trange, delta}, delta = 10^(-5); {o, r1, r2, srange, trange} = m; (* 交点を求め *) points = NSolve[torusEq[m] && lineEq[l], {x, y, z}]; (* パラメータ表示に直す *) params = Map[ Function[ point, Minimize[{ Norm[torusPar[m, s, t] - {x, y, z}] /. point (*, 0-delta < s < 1+delta && 0-delta < t < 1+delta*) }, {s, t} ][[2]] ], points]; (* 範囲外の物を除く *) Select[ params, 0 <= s < 1 && 0 <= t < 1 /. #1 & ] ] (* トーラス面上の s, t で表わされる点における法線ベクトルを求める *) normalVector[m_, s_, t_] := Normalize[ Cross[ D[torusPar[m, ss, tt], ss], D[torusPar[m, ss, tt], tt] ] /. {ss -> s, tt -> t} ] (* 入射方向ベクトルと法線ベクトルから反射方向ベクトルを求める *) reflection[v_, n_] := Module[ { vv = Normalize[v], nn = Normalize[If[v . n < 0, -n, n]] }, Normalize[vv - 2*vv . nn*nn] ] (* 直線 l がトーラス面 m で反射する様子を追跡する *) (* 反射点と反射方向ベクトルを返す *) trace[l_, m_] := Module[{p, v, st, q, n, r}, {p, v} = l; list = intersection[l, m]; Map[ Function[ st, q = torusPar[m, s, t] /. st; n = normalVector[m, s, t] /. st; r = reflection[v, n]; {q, r} ], list ] ] (* ベクトル v を v x z 方向に a, v x (v x z) 方向に b ずらす *) offsetVector[v_, a_, b_] := Module[{ez = {0, 0, 1}, ea, eb}, ea = Normalize[Cross[v, ez]]; eb = Normalize[Cross[v, ea]]; v + a*ea + b*eb ] (* 点 a と点 b を結ぶ線分を表示する *) showSegment[a_, b_] := ParametricPlot3D[ t * a + (1 - t) * b, {t, 0, 1}]
を用意しておき、例えば
With[{ n = 3, d = 0.1, m1 = {{8, 0, 0}, 0, 8, {-Pi/16, Pi/16}, {Pi(1-1/16), Pi(1+1/16)}} }, Show[ ParametricPlot3D[torusPar[m1, s, t], {s, 0, 1}, {t, 0, 1}], Function[ij, Module[{i, j, p, v, q, r, qrs}, {i, j} = ij; p = offsetVector[{10, -2, 0}, i d, j d]; v = Normalize[{-10, 2, 0}]; Map[ Function[qr, Module[{q, r}, {q, r} = qr; {showSegment[p, q], showSegment[q, q + 5*r]} ] ], trace[{p, v}, m1] ] ] ] /@ Flatten[ Array[ {#1 - n - 1, #2 - n - 1} &, {2 n + 1, 2 n + 1} ], 1 ], PlotRange -> {{-10, 10}, {-2, 2}, {-2, 2}} ] ]
半径 8 の球面の焦点距離が予想通り 4 程度になっているのが見て取れる。