Mathematica/ミラー光学系 の履歴(No.5)
更新概要†
ミラー光学系の収差を Mathematica を使って計算してみたい。
Mathematica によるプログラミングはあまり慣れていないので、 いろいろとまどろっこしいことをしていると思いますが、 ご容赦ください。
前提†
まずは収差のない場合のおさらいから。
焦点距離と曲率†
曲率が で表わされる凹面レンズの の位置に垂直に平行光が入射するとき、
反射面の傾きは なので、小角近似では
入射光は だけ曲がって反射されるので、
すなわち、焦点距離は となる。
曲率と曲率半径†
原点を通る半径 の円
を について の小さなところで解くと、
したがって、曲率半径 に対して、 となり、焦点距離は
となる。
点光源からの光を集光する†
凹面鏡の中心軸上、 の距離に置かれた点光源から、 凹面鏡上の の位置に到達する光が法線から だけずれているとすると、
の関係がなりたつ。平面鏡であればこの光は反射角 で反射されるが、 凹面鏡の傾きが であれば、実際の反射角は となる。
このようにして中心軸上に集められる光の焦点位置 は、
となるから、
の関係が得られる。
斜め入射の場合†
中心線に対して の角度で入射する平行光に対する焦点位置を考える。
中心位置で反射される光は反射角 で反射される
位置 で反射される光は反射角 で反射される
これらの光の為す角は だから、色を付けた三角形を考えることで、
&math( f_\parallel=\frac{x\cos\phi}{2\theta}=\frac{x\cos\phi}{4ax}=\frac{\cos\phi}{4a}=f\cos\phi );
となって、 のファクターだけ焦点距離が短くなる。
一方で入射面に垂直な方向には、入射面に垂直な面へ投影した際の焦点距離が 元の焦点距離と等しくなることから、斜めに進む分だけ焦点距離が長くなる。
&math( f_\perp=f/\cos\phi );
したがって入射面内での光の収束位置と入射面に垂直な方向への光の収束位置とがずれてしまうことになる。 すなわち凹面鏡への斜め入射は大きな収差(非点収差)を生むことになる。
この収差を無くすため、斜め入射で用いるための凹面鏡は入射面内の曲率半径と、 垂直方向の曲率半径を だけ違えて設計され、トロイダル凹面鏡などと呼ばれる。
トロイダル凹面鏡の面は数学的にはドーナツ面(トーラス面)の側面として表せる。
数値計算で確かめてみる†
以下は 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 程度になっているのが見て取れる。
この時の様子を入射面方向と入射面に垂直方向とで比べると、 入射面に平行に見た焦点距離が計算通り 4 なのに対して、 入射面に垂直に見た焦点距離が明らかに 4 より短いことも分かる。
perpendicularPlanePar[l_, d_, a_, s_, t_] := Module[{p, v, ea, eb}, {p,v} = l; eb = Normalize[Cross[v, a]]; ea = Normalize[Cross[eb, v]]; p + d v + s ea + t eb ] planeSection[l0_, d_, a_, l_] := Module[{p, v, stuList}, stu = NSolve[ perpendicularPlanePar[l0, d, a, s, t] == linePar[l, u], {s, t, u} ]; {s, t} /. stu ]
を用意して、反射光線
lines = Module[{ n = 4, d = 0.1, m1 = {{8, 0, 0}, 0, 8, {-Pi/16, Pi/16}, {Pi(1-1/16), Pi(1+1/16)}}, indice }, indice = DeleteDuplicates[ Flatten[ Array[ {(#1 - n - 1), (#2 - n - 1)} &, {2 n + 1, 2 n+1} ], 1 ]]; Map[ 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}]; trace[{p, v}, m1] ] ], Sort[indice, Norm[#1]<Norm[#2]&] ] // Flatten[#, 1]& ]
焦点付近でミラーからの距離を変えながら 光路に直交する平面との交点をプロットする。
Animate[ ListPlot[ Map[ Function[line, planeSection[lines[[1]], d, {0,0,1}, line] ], lines ] // Flatten[#, 1] &, PlotRange->{{-0.06,0.06},{-0.06,0.06}} ], {d,3.6,4.4} ]
アニメーションにすると変化が分かりやすい。
SetDirectory["C:\\Users\\osamu\\Desktop"] Table[ ListPlot[ Map[ Function[line, planeSection[lines[[1]], d, {0,0,1}, line] ], lines ] // Flatten[#, 1] &, PlotRange->{{-0.06,0.06},{-0.06,0.06}} ], {d,3.6,4.4,0.01} ]; Export["reflection.gif", %, "GIF"]
アニメーションを Web に載せるには、このようにして動画 gif を作ればいい。
Table[ ListPlot[ Map[ Function[line, planeSection[lines[[1]], d, {0,0,1}, line] ], lines ] // Flatten[#, 1] &, PlotRange->{{-0.06,0.06},{-0.06,0.06}} ], {d,3.87,4.15,0.02} ]
縦と横の焦点距離が違うことが明らかに分かるのと同時に、 縦も横も、球面収差などのせいで完璧な直線には収束しないことが分かる。
Table[ HistogramList[ Map[ Function[line, planeSection[lines[[1]], d, {0,0,1}, line][[1]][[1]] ], lines ], {-0.0525,0.0525,0.005} ][[2]], {d,3.6,4.4,0.05} ] // ArrayPlot[#, FrameTicks -> { Transpose[{Range@17, Table[x, {x, 3.6, 4.4, 0.05}]}], Transpose[{Range@21, Table[If[Abs[Round[100 x] - 100 x] < 0.001, x, ""], {x, -0.05, 0.05, 0.005}]}] }] & Table[ HistogramList[ Map[ Function[line, planeSection[lines[[1]], d, {0,0,1}, line][[1]][[2]] ], lines ], {-0.0525,0.0525,0.005} ][[2]], {d,3.6,4.4,0.05} ] // ArrayPlot[#, FrameTicks -> { Transpose[{Range@17, Table[x, {x, 3.6, 4.4, 0.05}]}], Transpose[{Range@21, Table[If[Abs[Round[100 x] - 100 x] < 0.001, x, ""], {x, -0.05, 0.05, 0.005}]}] }] &
などとして、ミラーからの距離に対する縦方向の幅、横方向の幅をプロットしてみると、 各方向に対する焦点距離を比較的正確に見ることができる。
それぞれ、およそ (4 ± 0.075) 程度と見積もれる。
元々の入射角が だから、 であり、
見積り通りと言える。
トラス面にすることで非点収差を補正してみる。
With[{ n = 3, d = 0.1, m1 = {{8, 0, 0}, 8 ArcTan[2/10]^2/2, 8 (1-ArcTan[2/10]^2/2), {-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, 0, -2}, i d, j d]; v = Normalize[{-10, 0, 2}]; 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 -> {{-0.1, 10}, {-2, 2}, {-2, 2}} ] ]