照度分布計算その3（ちょっと修正と実装開始）

さくっと終わらせたい回折計算による照度分布の計算の実装。
昨日の最後のところ、エイリアジング比を定義したが、これはN.A.が確定している場合には便利だけど、今回の場合、N.A.はアパチャの大きさで変わってしまうので便利さが無くなってしまう。そこでN.A.が変わっても他にその影響が及びにくいようなパラメータの定義を選ぶことにする。このへんは顕微鏡対物や光ディスクの計算と違うところやね。まあ、たいした差ではないけど。

Fourier変換は解析的には実行できないのでFFTを使う。そのため有限の格子点上での値だけにサンプリングされ、また座標の単位が格子の区切り方に依存するということになる。

座標の単位は瞳面の格子点の数をN、格子点間隔の実距離をqとするとFFTを行ったあとの格子点の間隔pはfraunhofer回折の理論から

なので

ここでfは焦点面を作るレンズの焦点距離であり

である。エイリアジング比の代わりにこっちのQの値を使うことにする。

実装

今回のような、さくっと書いてさくっと計算してさくっと結果を得たい場合には、あまり効率がどうのとか、サイズがこうのとか考えない方がいい。それよりもデバグのやりやすさ、それだけにつきる。だからなるべくミニマルに書くということを心がける。ようするにあるひとつのクラスだけ見て、他のクラスまで同時に参照しなくてもデバグできるように書く。もちろん、Objective-Cで書く。Mathematicaで書いてもいいけど計算を繰り返して最適な条件を探したいので計算はさくっと終わってくれる方がいいので。

まず、とにかく計算できて結果を出すことを前提に、Mathematicaをfront endに使い、ヘッドレスの状態で実装する。MathematicaからはRunThrough[]で起動することを前提に実装する。その後、必要があれば（さらに今後何度も計算する必要が出るなら）ユーザインターフェイスを作ることにする。

パラメータ

計算するときに必要なパラメータは、昨日の算数のまとめを思い出しながら

光源波長	lambda
光源FWHM(X-Y)	fwhmv, fwhmh
リレー倍率	magnificationOfRelay
リレーシフト(X-Y)	relayShiftv, relayShiftd
アパチャサイズ	apertureDiameter
アパチャシフト(X-Y)	apertureShiftv, apertureShiftd
絞り込みレンズ焦点距離	focalLengtdOfConversion
ピンホール径	pinholeDiameter
ピンホールシフト(X-Y)	pinholeShiftv, pinholeShiftd
ピンホールデフォーカス	pinholeDefocus
距離	distance
入射角	incidentAngle
基板サイズ	substrateDiameter
ミラー反射率	reflactanceOfLloyd
ミラー基板ギャップ	substrateGap
横シフト	lateralShift
瞳面サイズ	pupilAreaSize
FFTサイズ	fftSize

ぐらいかな。これらをユーザから受けて計算を行う。数値計算のソフトとしてはそれほどパラメータの数が多い訳ではないので、構造化せずフラットに保持することにする。ユーザからパラメータ設定を受けて、それぞれのオブジェクトに配分するuserDefaultsのようなオブジェクトをひとつ用意して、それはパラメータの名前とその値をNSDictionaryで持つことにする。

ほかのオブジェクトは自分の必要なパラメータの値をそのパラメータ設定オブジェクトに問い合わせることで得る。複雑になるとこういうやり方は訳わからなくなるけど、この程度の計算なら問題ないでしょ。

具体的には

@interface SParameters : NSObject {
    NSMutableDictionary	*parameters;
    NSDictionary	*defaults;
}
+ (id)currentParameters;
- (void)setValue:(id)value for:(NSString *)parameter;
- (id)querryFor:(NSString *)parameter;

こんなの。シングルトンで設定と照会を受け付ける。設定されていないパラメータはデフォルトの値を返すようにする。デフォルトの値はこの際、面倒なのでハードコードする。

SAmplitudePlaneクラス

計算の中心はSAplitudePlaneと言うクラスに実装する。これは

typedef struct {
    double	re;
    double	im;
} complex;

@interface SAmplitudePlane : NSObject {
    complex	*plane;
    int		sSize;	//	square shape
}

のようなもの。中身はパックした2次元のcomplexの正方配列で、大きさはsSize×sSizeである。complexはfftwの複素数型fftw_complex型とコンパチであって、読みやすさのために使っているだけである。

このクラスのインスタンスは

- (id)initWithAmplitudeFunction:(SAmplitudeFunction *)aFunc
            coordinateConverter:(SCoodinateConverter *)converter
                         ofSize:(int)squareSize;

で初期化される。SAmplitudeFunctionクラスは振幅分布を初期化するためのクラスで、ついでにアパチャなんかの透過率分布も同じクラスで表現する。たとえばgaussianの分布は

@interface SGaussian : SApmlitudeFunction {
    double	av;
    double	ah;
}
- (id)initGaussianFWHMVertical:(double)fwhmVertical
    andHorizontal:(double)fwhmHorizontal;
- (complex)amplitudeAtVertical:(double)v andHorizontal:(double)h;

みたいな感じ。

SCoodinateConverterクラスは、SAmplitudeFunctionが実数（double）の座標を受けるようになっているけど、それを2次元格子の整数座標にしたり、その逆を計算する変換用のオブジェクト。

SAmplitudePlaneのインスタンスは自分の配列を初期化するときに、配列の添字をSCoodinateConverterのインスタンスに渡して実数値を取り出し、その値をまたSAmplitudeFunctionに渡して配列のある場所の値を決める、と言う手順を踏むことにする。そのあとはこのふたつのクラスのインスタンスは不要なので、init...を出た後は参照も持たない。SAmplitudePlaneはその他に

- (id)addByPlane:(SAmplitudePlane *)othrePlane;
- (id)multiplyByPlane:(SAmplitudePlane *)otherPlane;
- (id)scaleBy:(double)factor;
- (id)rotateHalf;
- (id)fourier;
- (id)inverseFourier;
- (id)intensityPlane;

などのインスタンスメソッドを持つ。これはどれも、足したり掛たりした結果を保持しているSAmplitudePlaneを新しく作ってidとして返す。自分の中身を書き換える方が手っ取り早くてメモリ領域も無駄にならないけど、計算する前も残しておいた方がデバグは楽になるのでこうする。

基本的には回折計算の部分はこのクラスでやってしまう。FFTを計算するのにはfftwのライブラリを使うつもりなのでそのスタイルに合わせるけど、このクラスのインターフェイスにはそれは見えなくなる。べつに見える見えないはどっちでもいいが、本当は効率を考えると同じ大きさの2次元FFTばかりを繰り返すことになるので、演算は自分の中身を書き換えて、しかもfftw用の計算設定構造体

fftw_plan	p;

をインスタンス変数に持っておく方が断然いい。でも今回は効率は二の次なのでfourier変換メソッドが呼ばれるたびにfftw_planを作ってFourier変換して終わったら毎回fftw_planは破棄する、と言う方針で行くことにする。僕の経験からはこういう一見無駄なことが数値計算のデバグのしやすさになる。効率は必要になった時点で上げることを考えれば良い。特に、自分は頭が良くない、とか自分は人よりかなり記憶力が劣っている、とか思っている人にはおすすめ。僕は自分自身に強く勧めたい。

ああ、もう遅くなった。続きは明日にしよ。