Mathematicaお役立ちPackage - その1

Mathematicaを6から10にすることができたので、7以降に定義された新しい機能を探索（Unitパッケージはobsoleteになったのね。便利だったのでけっこう慣れてたんだけど）しながら、普段使っている自分で書いた古いパッケージを10に変換するという作業をしている。5から6の時は結構苦労したんだけど、今回はファイルフォーマットを変更するだけで、ほとんど修正の必要はなさそうである。この際なのでできるところはちゃんと整備して（できないのはほったらかしで）、汎用性の高そうなの（と言ってもごく一部の人にとってだけだろうけど）はここで公開することにする。また、普通のMathematicaユーザはあまりパッケージにまとめることはしないらしいので、パッケージプログラムの例としてみてもらってもいいかもしれない。第1回の今日はガウシアンビーム伝搬計算用パッケージGaussianBeam。普段使っているパッケージの中では比較的新しいもの．．．．

僕のMathematicaパッケージの書き方

まず最初に、僕のMathematicaパッケージの書き方をまとめておく。

1. パッケージは大分類の直下に配置する。例えば今回のGaussianBeamパッケージのコンテクストは「Optics`GaussianBeam`」である
2. パッケージは名前、バージョン、ヒストリ、簡単な使い方の文字列を必ず定義する。今回のGaussianBeamパッケージでは名前がGaussianBeamで、そのコンテクストにGaussianBeamVersion、GaussianBeamVersionHistory、GaussianBeamUsageが定義されている。このパターンはどのパッケージでも同じにする
3. 上記のパターンのシンボルは大文字で始める。それ以外のパッケージ内で定義されたシンボルはすべて小文字で始める。例えばcreateGaussianBeamなどである。これはMathematica本体で定義されたシンボルではなくパッケージ内のシンボルであることがすぐわかるようにするためである
4. usage messageやコメントに日本語は使わず、英語表記のみにする
5. パッケージのファイルは「*.m」でこれだけあれば使えるが、「*.nb」ファイルも同時に配置する。「*.nb」にはパッケージ定義と、簡単なExampleを書く。Exampleセルは初期化セルではなく、「*.nb」ファイルから「*.m」を改めて生成しても影響しない

などである。これは僕が20年来守ってきたやり方である。

日本語を使わないのはCやObjective-Cのソースも同じで、僕はソースにはタブと改行以外はASCIIで0x20〜0x7Eの範囲だけにしている。最近の人にはわからないだろうけど、僕は若いころシフトJISエンコードの日本語をコメントに埋めたばっかりにつまらないバグに悩まされた（日本語コメントのせいだとわかった後に日本語をソースから全部削除してもバグがおさまらないという不思議な現象に悩まされ続けた）。

今ではUTF-8とかを使うし、コンパイラもUTFは通るようになっているのが普通なので問題ないのだろうけど、そういう昔からの癖と、僕の嫌いな

    i ++;   // iをインクリメントする

なんてコメントを書く気力をなくさせてほんとに必要なコメントだけを書くようにして、コードそのものに語らせる意識を失わないようにするためである。それと、英語は技術屋のデファクトの共通言語なので、それを尊重する姿勢を表す意味もある。下手くそな英語のせいで何言ってるのかわからなくなることも多いけど。コメントは未来の自分自身のためのものである、ということを身にしみてわかっているので内容は厳選したい。

Gaussian Beamのパラメータ計算

ガウシアンビームはレーザ光の伝搬を考えるときによく使う波動方程式の近似解である。僕が今いる会社は半導体レーザを扱っていて、それを利用した簡単な応用製品を作っている。モジュールの先10mのところに1$\mu$mのスポットを作るようなレーザモジュールを作れ、とごむたいなことをときどき言われるので、その場でサクッと計算できるようにパッケージを作った。

ダウンロードはここから

このzipを解凍して、$PathにあるパスにOpticsディレクトリを作って、その中に入れると呼び出せるようになる。OS Xの場合、自分のホームディレクトリのLibrary/Mathematica/Applicationsの下に作るのがオススメのはずである。例えば僕の場合

/Users/decafish/Library/Mathematica/Applications/Optics/GaussianBeam.m

となっている。

何ができるかというと、波長とビームウェスト半径を指定すると、他のパラメータが計算できる、というもの。まあそれだけ。具体的には、まずパッケージを読み込む。

Needs["Optics`GaussianBeam`"];

成功すると、usage messageが表示できるはずである。

?Optics`GaussianBeam`*

僕がときどきやるパターンで、ひとまとまりのパラメータをひとつのMathematicaシンボルに結びつけるやり方を使っている。オブジェクト指向風にすることで複数の組み合わせを独立に扱いやすいのでこうしている。このGaussianBeamパッケージはまず、適当なシンボルを用意して専用にする。

createGaussianBeam[gbSymbol, 0.000633, 0.05]

これはgbSymbolというシンボルをGaussianBeam用に使い、波長が633nm、ウエスト半径が0.05mmのガウシアンビームである、という宣言である（もちろん次元と単位の一貫性はユーザが保証する）。もしこの前にgbSymbolになにか定義されていたら、それはクリアされる。

また、

createGaussianBeam[gbSymbol, 0.000633, 10,specifiedArgument->rayleighLength]

とオプションをつければ最後の引数をRayleigh長さとみなして設定する。

こうやってシンボルを設定しておいて、知りたいパラメータの名前を引数にしてシンボルを呼べばその値が返ってくる。例えばこのビームの広がり角は

gbSymbol[divergenceAngle]

を評価すると値を返す。光軸に沿ったビーム半径は

gbSymbol[beamRadius]

が引数をひとつとる純関数を返すのでこれをプロットすればいい。

Plot[Evaluate[gbSymbol[beamRadius]][z], {z, 0, 200}]

他も同じなので、usage messageとExampleを見て欲しい。

これだけではつまらないので、後ひとつ工夫がある。僕の場合、レンズの出口から何mmのところに一番小さなスポットを作りたい、という要求が多い。そのときのウェスト半径や広がり角を知りたい、という。それには

createGaussianBeam[epup, 0.000633, {0.5, 100}]

という呼び出し形式を使う。これは波長633nmで、レンズの射出瞳のところでビーム半径が0.5mmのとき、100mm先でウェストが来るようなガウシアンビームパラメータを設定する。そのときのパラメータは

Map[epup,{beamWaist,divergenceAngle,rayleighLength}]

とすれば、ウェスト半径、広がり角、Rayleigh長さがわかる、というようなものである。射出瞳位置でのビーム半径を決めたとき、ある距離以上にウェストを作ることはできないので、その場合には最大値はこれだけという値を返して、シンボルを変更しないで返る。

また、一番最初の呼び出し形式で作ったシンボルを

setFocalDistance[gbSymbol,5]

で、ウェスト位置を変更することもできる。最初の呼び出し形式ではウェストと瞳位置が一致していて原点にある、とみなしていたのを、場所として分離するという作業をする。このとき、setFocalDistanceに

preserveWaist->False

というオプションをつけると、瞳位置でのビーム半径を保持して、ウェストはその先にする、従ってより絞り込まれる、ということにできる。

複数のレンズでリレーするような光学系の場合は近軸マトリクスを使ったやり方が簡単なんだけど、僕の場合あまりそういう光学系を要求されることは少なく、これだけで十分な場合がほとんどである。

実はウェスト位置と射出瞳位置でのビーム半径を決めると解はふたつできる。太い方の解が必要になることが今までなかったので、その一方のウェストの細い方だけ設定するようになっている。汎用性を考えるとちょっと寂しい対応ではある。

そのためにもうひとつオプションを用意してある。

createGaussianBeam[epup, 0.000633, {0.5, 0.01},specifiedArgument -> curvatureValue]

とすると、引数のリストのふたつめ（例では0.01）は射出瞳での波面の曲率とみなす（本来符号があるけど、ここでは射出瞳を出てからウェストが来るように勝手に符号を選ぶ）。こっちの形式を使えば太い方の解も得られることになる。波面の曲率に関しては解は連続にひとつだけ決まるからである。なぜそうなのかは計算してもらえばいいけど、以前ガウシアンビームの解をまったりと導こうと思って始めたのが途中になっている。再開のタイミングを失っていただけなので、この際なので再開して、この波面の曲率だと解がひとつに決まることなどを示したい。

たかだかこれだけのパッケージで、最初のバージョンはたぶん2時間ぐらいで作った。でもこれは僕が日頃一番たくさん呼び出すパッケージになっているので、一番最初に10化した。と言っても10で読み込んで、「更新されたスタイルで続行」というボタンが現れるので、それを押して保存するだけでよかった。Mathematicaの関数は抽象性が高くて、汎用的にできているのでバージョン間の互換性に悩むことは少ないのでありがたい。逆に言えば、10の機能を使っていなくてバージョンアップは無駄、ということでもあるんだけど。

ということで、このあとも他の人に役立ちそうだと思えたパッケージはここで公開することにする。でもまあ、ほとんど光学系の技術屋以外には無用の長物だろうけど。