偏光の計算 - その11

前回までやった数学を具体的に計算するためのMathematica関数をいくつか書いてパッケージにした。コードの簡単な説明から。

Mathematica は代数的な記述がやりやすいのでそのまま書き下せばよい。

偏光を扱うMathematica関数を書く

シンボル名の注意

MathematicaにはContextという名前空間のための構造があってしかもMathematicaがはじめから持っているシンボルのほとんどはProtected属性を持っているので間違って再定義してしまうということはまずない。でも苦労して書いた関数が名前の衝突で書き直すことになるのも煩わしい。ということで昔から僕はシンボルを小文字で始めることにしている。そうすれば少なくともMathematicaが定義しているシンボルと衝突することはない。これはMathematicaの流儀ではないけど面倒をちょっとだけ減らすことができる。

今回偏光計算用に定義したシンボルは全部（ふたつの例外を除いて）小文字で始めてある。大文字で始まっているシンボルはMathematicaが定義しているもので小文字のは今回の、として見てもらうと読みやすくなる。

偏光のMathematica 表現

まず偏光を表現するシンボルを作る。

 stokesParameter[s0, s1, s2, s3]
 poincareCoordinate[s1, s2, s3]
 ellipseRotated[psi, chi]
 jonesVector[ex, ey]
 ellipseCartesian[alpha, delta]

それぞれストークスパラメータ、ポアンカレ球上の座標値（直交座標）、楕円率と軸方向の組み、ジョーンズベクトル、振幅の比と位相の組み、を表すことにする。これらのシンボルは単なるタグで、このままでは評価されず、書き換えは起こらない。

数学を思い出せば明らかだけど、jonesVectorだけ複素数を引数にとる。それ以外は実数値だけになる。

偏光の表現間の変換

変換を実行するには、例えばジョーンズベクトルから振幅の比と位相の組みに変換するには

 ellipseCartesian[jonesVector[ex, ey]]

を評価する。例えば具体的に

 In[3]:= ellipseCartesian[jonesVector[1, I]]
 Out[3]= ellipseCartesian[Pi/4, Pi/2]
 In[4]:= ellipseCartesian[jonesVector[a1, a2]]
 Out[4]= ellipseCartesian[ArcTan[Abs[a1], Abs[a2]], -Arg[a1] + Arg[a2]]

などとなる。

ちょっとした扱いを楽にするために

 Attributes[jonesVector] = {Listable};
 jonesVector[jonesVector[ex_, ey_]] := jonesVector[ex, ey]

などとしておく。1行目はリストが渡されてきたときには自動的にリストの要素を変換することになる。2行目はジョーンズベクトルからさらにジョーンズベクトルに変換するのは恒等変換になるようにする。こうしておけば、変換済みかそうでないかをチェックする必要が無くなる。これはもちろん

  jonesVector[jv_jonesVector] := jv

などと書いても同じことになる。

変換規則は例えば

 ellipseCartesian[jonesVector[ex_, ey_]] := 
  ellipseCartesian[ArcTan[Abs[ex], Abs[ey]], Arg[ey] - Arg[ex]]
 jonesVector[ellipseCartesian[alpha_, delta_]] := 
  jonesVector[Cos[alpha], Sin[alpha] E^(I delta)]

などというようにそのまま書き下せばよい。割り算の分母が0になったりするのを気をつける必要がある。

これを必要な部分だけ書く。 先に書いたように、いきなり全部の変換規則を書き下すのは面倒だし、デバグが大変になるので図-6の矢印に対応するところだけの変換規則を定義する。

それ以外の変換規則はpoincareCoordinateを経由するように書く。具体的にはたとえば

 ellipseCartesian[hd_[args__]] := 
  Check[ellipseCartesian[
    checkConvert[poincareCoordinate[hd[args]], 1] /. 
     checkConvert -> (#1 &)], $Failed]

とする。 これは変換したい任意の形式を一旦checkConvertに渡す。 Check[]という組み込み関数はその引数を評価して、Messageが発行されたら第2引数の方をCheckの結果として返す。Messageがなかったら第1引数の評価結果を返す、というもの。

こうしておいてcheckConvertは変換ができなかった場合にMessageを発行すればいい。

 poincareCoordinate::cannotconvert = 
   "Cannot convert `1` to generic poincareCoordinate.";

 checkConvert[poincareCoordinate[hd_[args__]], n_] := 
  If[n >= 2, Message[poincareCoordinate::cannotconvert, hd], 
   checkConvert[poincareCoordinate[hd[args]], n + 1]]

とする。これはcheckConvertが2回入れ子に呼ばれると変換できなかったと見なしてMessageを出力する。 Mathematica は無限評価系（評価結果が変化しなくまで評価を繰り返す）なので、例えば単純な

 ellipseCartesian[unknown_] := 
  ellipseCartesian[poincareCoordinate[unknown]]

のようなやりかただと、変換できなかった場合はpoincareCoordinateが無限に入れ子になる（実際には$RecursionLimit繰り返されてエラーで止まる）。

こうしておけば

  In[84]:= ellipseCartesian[ellipseRotated[Pi/4, 0]]
 Out[84]= ellipseCartesian[Pi/4, Pi/2]
  In[85]:= ellipseCartesian[g[0, 0]]
  poincareCoordinate::cannotconvert: Cannot convert g to generic poincareCoordinate.
 Out[85]= $Failed

というふうになる。

checkConvertはもう少し簡単な書き方ができると思うけど、こういう変換の記述はMathematica らしくて慣れるとシンプルでわかりやすい。

偏光変換素子を表すシンボル

素子のシンボルとしてとりあえず、

 retarder[phse,dir]
 polarizer[dir]

を書いた。これもそのままでは未評価のまま残るが、jonesMatrixに変換できる。

 In[95]:= jonesMatrix[retarder[Pi/2, Pi/4]]
 Out[95]= jonesMatrix[1/2 -I/2, 1/2 + I2,  1/2 + I/2, 1/2 - I/2]

このjonesMatrixは2×2のマトリクスの要素をただ並べただけ。

ある偏光に対していくつかの偏光変換素子を通した後の偏光を計算するapplyOpticsを書いた。

 In[96]:= applyOptics[jonesVector[1, 0],
     {retarder[Pi/2, Pi/4], retarder[Pi/2, Pi/2]}]
 Out[96]= jonesVector[-(1/2) - I/2, 1/2 + I/2]

のように、引数のリストにある素子を順番に作用させてその結果のジョーンズベクトルを返す。

またintensityは

 In[98]:= intensity[
     applyOptics[jonesVector[1, 0], {retarder[Pi/2, Pi/4],
      retarder[Pi/2, Pi/2], polarizer[Pi/2]}]]
 Out[98]= 1/2

のように光強度を計算する。これは素子のリストの最後がpolarizerでないと意味がない。

applyOpticsにはverctorListというオプションがあってTrueをセットすると

 In[99]:= applyOptics[
      jonesVector[1, 0], {retarder[Pi/2, Pi/4], retarder[Pi/2, Pi/2]}, 
       vectorList -> True]
 Out[99]= {jonesVector[1, 0], jonesVector[1/2 - I/2, 1/2 + I/2], 
      jonesVector[-(1/2) - I/2, 1/2 + I/2]}

のようにそれぞれの素子の途中の結果も出力する。

Mathematica の関数と組み合わせて例えばグラフの表示も当然できる。

 Plot[Evaluate[
   intensity[
    applyOptics[
     jonesVector[1, 0], {retarder[Pi*550/lambda, Pi/4], 
      polarizer[Pi/2]}]]], {lambda, 400, 700}, PlotRange -> {0, 1}]

とすれば図-7のように描ける。