照度分布計算その4 - 実装の方針

さくっと1週間ぐらいで終わらせたかったのに、他のことやってるせいでブツ切れになって、前やったことを忘れてしまう回折計算による照度分布の計算の実装。まず、FFTを回折計算に使うときの注意。

パーセバルの定理を調整する

fftwなんかのFFTでの「パーセバルの定理」は1/Nの係数分だけ違っている。物理の分野ではFourier変換の定義を変更して振幅の絶対値の2乗の和（総エネルギーに対応する）が変化しないようにするのが普通なので、FFT、逆FFTとも変換後にすべての要素を√Nで割る必要がある。今回は2回Fourier変換を行うのでそのままでは数字が大きくなってしまって、計算上の問題は無いけど、表示やデバグには不便である。ということで、FFTを行った後、毎回全要素を√Nで割って総エネルギーが一定になるようにする。
このへんはSAmplitudePlaneのクラスに埋め込んでしまう。

光軸位置の調整とSCoordinateConverterクラス

それにfftwを使うと（ってFFTだとみんな同じだけど）1次元でも多次元でも配列の最初の要素が原点になる。添字が負の方向i<0に対してはN+iに配置される（というか、サンプリング定理から同じものである）。従って回折計算に使おうとした場合、光軸中心がデータの端っこになる。瞳面か焦点面かで並べ直したりしなかったりするのは不便なので、瞳面でも光軸中心を配列の原点に持ってくることにしておく。

SCoordinateConverterはこの変換もコミでやる。つまり

• 正の実座標にはそれに対応するインデクス（<N/2）を返す
• 負の実座標にはN/2 < i < Nを返す
• 逆にN/2より小さいインデクス値には正の実座標を
• N/2 < i < Nの範囲にあるインデクス値iには負の実座標を

それぞれ返すことにする。

つまりこのクラスは

として実座標からインデクス値を、この逆関数の

でインデクス値から実座標を返す。式-40のround()は四捨五入して整数に丸める関数。

また、配列から値を取り出すとき、実座標を指定してインデクス値を得る必要があるが、このときバイリニア補完できるように小数のインデクス値f_xも返すようにする。つまり

これには横シフトは関係ない（配列自身を取り出すための座標なので）ので付加されていない。実際には式-40は必要ない。単にデバグのために作っておく。

さてそのSCoordinateConverterのヘッダは

@interface SCoordinateConverter : NSObject {
    int     iRange;
    double  indexUnitv;
    double  indexUnith;
    double  shiftv;
    double  shifth;
}

- (id)initWithVerticalUnit:(double)vUnit
            horizontalUnit:(double)hUnit
             verticalShift:(double)vShift
           horizontalShift:(double)hShift
             andIndexRange:(int)indexRange;


- (double)horizontalFromIndex:(int)index;
- (double)verticalFromIndex:(int)index;
- (int)horizontalIndexFrom:(double)hor;
- (int)verticalIndexFrom:(double)ver;
- (double)horizontalFractionalIndexFrom:(double)hor;
- (double)verticalFractionalIndexFrom:(double)ver;
@end

というようなもの。単に変換の定数を保持して、それぞれ変換、あるいは逆変換するだけ。indexRangeはNのことで、*FractionalIndexFrom:は式-42用のもの。

これにあと、自分でかってにSParametersのインスタンスを読みに行って自分でかってに設定する

- (id)initForLaser;
- (id)initForAperture;
- (id)initForPinhole;
- (id)initForSubstrate;

というのを作っとけばいいだろう。これはSParametersのインスタンスに自分の必要なパラメータの値を問い合わせて、それをもとに自分のinitWithVerticalUnit:horizontalUnit:verticalShift:horizontalShift:andIndexRange:を呼ぶ。計算はこのそれぞれのinit*に任せる。こないだのSAmplitudeFunctionも同じように

@interface SGaussian (SParameters)
- (id)initForLaser;
@end
@interface SCircularAperture  (SParameters)
- (id)initForPinhole;
@end

なんかを作ってそれぞれ自分でSParametersのインスタンスを読みに行くようにしよう。計算が局所的になってデバグが楽になってよろしい。

ピンホールのあとの照度計算

「その2」でやった式-22の計算を引き受けるクラスSGeometricIssuesを定義する。

@interface SGeometricIssues : NSObject {
....
}
- (id)initWithDistance:(double)distance
        substrateAngle:(double)substrateAngle
  andSubstrateDecenter:(double)decenter;

- (void)setCurrentPositionOnSubstrateX:(double)x andY:(double)y;

- (double)illuminationCoefficient;
- (double)lengthFromPinhole;
- (double)horizontalDirection;
- (double)verticalDirection;

これはinit*で定数を設定し、基板上のどの位置を計算するか、をsetCurrentPositionOnSubstrateX: andY:で設定する。そのあと、illuminationCoefficientは式-22の値、lengthFromPinholeはピンホールからの実距離、などを返す。ちょっとダサいけどまあええでしょ。 ついでに、他のクラスと同じようにSParametersからかってに読むための

- (id)initForGeometry;

も作っておく。

あとはいわゆるコントローラクラスと、入出力のクラスを書けばいい。それは適当にやればいいので、このあとはこれまでヘッダを書いたクラスの実装ファイルを書こう。