照度分布計算その12 - 計算制御

回折による照度計算の続き。あともうちょっと。今回は計算の流れを制御するクラス。その2でやった、数学に従って順に計算を進めるSPropagateToSubstrateというクラス。

材料はだいたいそろっているので、計算手順に従って進めるだけ。 計算手順をおさらいすると

1. レーザの出力の振幅分布を表す複素2次元配列を作る
2. アパチャの透過形状を合成する
3. 焦点位置の振幅分布をFFTで求める
4. ピンホールの透過形状を合成する
5. ピンホール射出後の放射分布をFFTで求める
6. 強度に変換する
7. 基板の特定位置での照度を計算する

ということ順番にやればいい。

たとえば、一番最初のレーザの振幅分布を作るには

  SGaussian            *lfunc = [[[SGaussian alloc] initForLaser]
                                                    autorelease];
  SCoordinateConverter *lconv = [[[SCoordinateConverter alloc]
                                        initForLaser] autorelease];
  SAmplitudePlane      *laser = [[[SAmplitudePlane alloc]
                                initWithAmplitudeFunction:lfunc
                                      coordinateConverter:lconv
                                                   ofSize:sSize]
                                                    autorelease];

とすればいい。(1)振幅分布関数を用意して、(2)座標変換のインスタンスを作って、(3)複素配列を初期化する。決まりきった書き方だけど、長くてうるさいので、書き換え用に

#define	gInit(gmethod)	[[[SGaussian alloc] gmethod] autorelease]
#define	dInit(dmethod)	[[[SCircularApertureWithDefocus alloc]\
                                    dmethod] autorelease]
#define	aInit(amethod)	[[[SCircularAperture alloc] amethod]\
                                    autorelease]
#define	cInit(cmethod)	[[[SCoordinateConverter alloc] cmethod]\
                                    autorelease]
#define	pInit(f, c)     [[[SAmplitudePlane alloc]\
                        initWithAmplitudeFunction:f\
                              coordinateConverter:c\
                                           ofSize:sSize]\
                                    autorelease]

とした。いまさらかっこわるいけど、しょうがない、読みにくいよりまし。これを使って

- (void)propagate
{
    SParameters	*par = [SParameters currentParameters];
    int         sSize = [[par querryFor:kFftSize] intValue];

    SGaussian               *lfunc = gInit(initForLaser);
    SCoordinateConverter    *lconv = cInit(initForLaser);
    SAmplitudePlane         *laser = pInit(lfunc, lconv);
    [self addToResult:laser IfChecked:kOutputLaserPlane];
	
    SCircularApertureWithDefocus    *afunc = dInit(initForAperture);
    SCoordinateConverter    *aconv = cInit(initForAperture);
    SAmplitudePlane         *aper = pInit(afunc, aconv);
	
    SAmplitudePlane         *afterAper = [laser multiplyByPlane:aper];
    [self addToResult:afterAper IfChecked:kOutputAperturePlane];

    SAmplitudePlane         *prepin = [afterAper fourier];
	[self addToResult:prepin IfChecked:kOutputPrePinholePlane];
	
    SCircularAperture       *pfunc = aInit(initForPinhole);
    SCoordinateConverter    *pconv = cInit(initForPinhole);
    SAmplitudePlane         *pin = pInit(pfunc, pconv);
    [self addToResult:pin IfChecked:kOutputPinholePlane];
	
    SAmplitudePlane         *postpin = [prepin multiplyByPlane:pin];
    [self addToResult:postpin IfChecked:kOutputPostPinholePlane];
    SAmplitudePlane         *presubst = [postpin fourier];
	
    intensity = [[presubst intensityMap] retain];
    [self addToResult:intensity IfChecked:kOutputIntensityPlane];
    gissue = [[SGeometricIssues alloc] initForGeometry];
    converter = [[SCoordinateConverter alloc] initForSubstrate];

    [self scanSubstrateBy];
}

とした。これでも長いけど、まあ、いっこずつ見れば順番通りにやってるだけ。ここにでてくる-addToResult:IfChecked:メソッドはouput*パラメータが設定されていればNSMutableDictionaryのresultインスタンス変数にそのパラメータの名前をキーにしてセットするような

- (void)addToResult:(id)obj IfChecked:(NSString *)name
{
    SParameters	*par = [SParameters currentParameters];
    if (![[par querryFor:name] isEqualToNumber:nFalse])
        [result setObject:obj forKey:name];
}

なんてもの。ここにあるnFalseはstaticに確保された変数でinitの中で

    if (nFalse == nil)
        nFalse = [[NSNumber alloc] initWithInt:0];

としている。ようするにoutput*のフラグをチェックするための定数。

それと最後の-scanSubstrateByメソッドはメッシュに切った基板上での値を計算して配列に入れる

- (void)scanSubstrateBy
{
    NSMutableArray  *direct = nil;
    SParameters     *par = [SParameters currentParameters];
    double          sradius = [[par querryFor:kSubstrateDiameter]
                               doubleValue] * 0.5;
    double          rad2 = sradius * sradius;
    double          step = [[par querryFor:kOutputMeshStep] doubleValue];
    double          h, v;
    int             capacityHint = (int)(rad2 / step / step * 3.14);
    int             startp = -((int)(sradius / step)) * step;
	
    if (![[par querryFor:kOutputDirectIllumination] isEqualToNumber:nFalse]) {
        direct = [NSMutableArray arrayWithCapacity:capacityHint]; 
        [self addToResult:direct IfChecked:kOutputDirectIllumination];
    }
    for (v = startp ; v <= sradius ; v += step) {
        double	v2 = v * v;
        for (h = startp ; h <= sradius ; h += step) {
            if (v2 + h * h <= rad2) 
                [self appendTo:direct horizontal:h andVertical:v];
        }
    }
}

みたいなの。ああ、汚い。やってることはループを回すことだけ。最後の-appendTo:horizontal:andVertical:メソッドは

- (void)appendTo:(NSMutableArray *)ret horizontal:(double)h andVertical:(double)v
{
    [gissue setCurrentPositionOnSubstrateX:h andY:v];
    double  x = [converter horizontalFractionalIndexFrom:
                                  [gissue horizontalDirection]];
    double  y = [converter verticalFractionalIndexFrom:
                                  [gissue verticalDirection]];
    double  illum = [gissue illuminationCoefficient];
    double  intat = [intensity bilinearedRealValueAtIndexX:x andY:y];
    SDistribution   *dist = [[SDistribution alloc]
                             initX:h Y:v andItsValue:illum * intat];
    [ret addObject:dist];
}

というもので、具体的な計算をやってる。gissue、converter 、intensityはインスタンス変数でそれぞれSGeometricIssues、SCoordinateConverter、SAmplitudePlaneのインスタンス。どんどん醜くなるな。かっこわるう。自己弁護する訳ではないが、影響を受けるパラメータの数も多いので、しょうがない面もある。

最後のSDistributionというのは

@interface SDistribution : NSObject  {
    double  x;
    double  y;
    double  val;
}
- (id)initX:(double)xi Y:(double)yi andItsValue:(double)value;
@end

なんてもので、計算結果を保持するためだけのクラス。SMathematicaStyleDescriptionプロトコルに従ってMathematica用に{x,y,val}の形に出力できるようにしている。この実装は、まあここに出すまでもないな。

さて、あとはmain()の中身を書いて実行デバグするだけ。よし、行こう。