第 2 章 チュートリアル 33
2.4 伝搬回折計算
2.4.1 角スペクトル法による円形開口からの回折像計算
次のソースコード例は，直径1mmの円形開口からの回折像を回折距離20mmで求めるものである．
バイナリの円形開口
まず，開口パターンを2値で与えた円形開口からの計算例を示す．
Example 円形開口からの回折像計算(1)
1 # i n c l u d e < wfl . h >
2 u s i n g n a m e s p a c e wfl ;
3
4 v o i d m a i n ()
5 {
6 S t a r t ();
7 W a v e F i e l d a (256 , 256 , 2 e - 6 ) ;
8
9 // = = = = = 円 形 開 口 をaに 生 成 す る
10 int i , j ;
11 for ( i = 0; i < a . Nx (); i ++)
12 {
13 for ( j = 0; j < a . Ny (); j ++)
14 {
15 d o u b l e x = a . i t o x ( i );
16 d o u b l e y = a . j t o y ( j );
17 d o u b l e r = s q r t ( x * x + y * y );
18 if ( r < 0.1 e -3) //半径0.1 m mの 円 形 開 口
19 a . S e t P i x e l ( i , j , C o m p l e x (1.0 , 0 ) ) ; //開口内は1.0
20 e l s e
21 a . S e t P i x e l ( i , j , C o m p l e x (0 , 0 ) ) ; //開口外は0.0
22
23 }
24 }
25 a . S a v e A s B M P ( " A p e r t u r e . bmp " , A M P L I T U D E ); //開口形状の保存 26
27 // = = = = = 角 ス ペ ク ト ル 法 に よ る 回 折 伝 搬 計 算
28 a . A s m P r o p (1 e - 3 ) ; // 1 m mの 伝 搬 回 折 計 算 を 実 行
29 a . N o r m a l i z e ();
30 a . S a v e A s B m p ( " D i f f r a c t i o n . bmp " , A M P L I T U D E ); //回折像(振 幅 像)の 保 存 31 }
ここで，10∼24行はオブジェクトa内に円形開口を設定している．ここでは，単純にオブジェクトa の全サンプリング点をスキャンして，物理座標の原点から半径0.1mm以内ではサンプリング値を1，そ れ以外では0に設定している．作成した円形開口を図 2.6(a)に示す．
( a ) A p e r t u r e . b m p ( b ) D i f f r a c t i o n . b m p 図2.6 バイナリの円形開口とその回折像
回折計算は，28行でWaveFieldクラスの AsmProp()メンバー関数を呼び出して行っている．オブ
ジェクトaの内容はAsmProp() メンバー関数の実行により，回折計算の結果に置き換わる．つまり，
WaveFieldオブジェクトを一種のスクリーンと考え，スクリーンが指定の距離だけ平行移動していると考
えれば良い．求めた回折像(振幅像)を図2.6(b)に示す．
ここでは，回折計算の手法として角スペクトル法 [1]を用いてる．日本の光学の参考書にはあまり記述 がないが，角スペクトル法は平面波展開法とも呼ばれ，フレネル近似条件を用いずに波動方程式の解を求 める手法の一つである．フレネル近似を用いないため，伝播距離などの制限がなく，少なくとも解析的に はヘルムホルツ方程式の完全な解となっている．実際ここでも，わずか1mmと短距離の伝搬計算を行っ ている．
ジャギーの影響を軽減した円形開口
よく見ると，図2.6(b)の回折像には本来光がないはずの周辺部にも放射状の光が見られる．これは，円 形開口が透過度0と1の2値パターンであり，開口の縁で，ビットマップ画像におけるジャギーと同様の 効果が発生しているためであると考えられる．ジャギーを軽減するために，高次ガウス関数を用いアンチ エイリアシングと同様の効果を持たせたのが次のコード例である．
( a ) A p e r t u r e . b m p ( b ) D i f f r a c t i o n . b m p
図2.7 アンチエイリアス効果を持たせた円形開口とその回折像
Example 円形開口からの回折像計算(2)
1 # i n c l u d e < wfl . h >
2 u s i n g n a m e s p a c e wfl ;
3
4 v o i d m a i n ()
5 {
6 S t a r t ();
7 W a v e F i e l d a (256 , 256 , 2 e - 6 ) ;
8
9 // = = = = = 円 形 開 口 をaに 生 成 す る
10 a . S e t G a u s s i a n ( 0 . 1 e -3 , 5 0 ) ; //半径0.1 m mの 円 形 開 口
11 a . S a v e A s B M P ( " A p e r t u r e . bmp " , A M P L I T U D E ); //開口形状の保存 12
13 // = = = = = 角 ス ペ ク ト ル 法 に よ る 伝 搬 計 算
14 a . A s m P r o p (1 e - 3 ) ; // 1 m mの 伝 搬 回 折 計 算
15 a . N o r m a l i z e ();
16 a . S a v e A s B m p ( " D i f f r a c t i o n . bmp " , A M P L I T U D E ); //回折像(振 幅 像)の 保 存 17 }
この例では，10行目でWaveFieldクラスのSetGaussian()メンバー関数を用いて，円形開口を作成 している．ガウス関数は本来指数部が2次の関数であるが，SetGaussian()メンバー関数は第2引数で
44 第2章 チュートリアル 2次を超える高次(ここでは50次)を指定*1できる．この様にすると，非常に急峻な立ち上がりを有する がバイナリではない円形開口ができる．
開口パターンとその回折像を図 2.6に示す．この回折像では，図 2.6(b)に見られたような放射状のノ イズが消えていることがわかる．
伝搬距離が長い場合
角スペクトル法による伝搬回折計算は，FFTを用いた離散高速畳み込み演算により行われる．しかし，
WaveFieldクラスのAsmProp()メンバー関数は単純な円状畳み込みを行うため，伝搬距離が伸び，光が サンプリング領域の境界近くまで広がってくると誤差が生じる．AsmProp()メンバー関数を用い，伝播 距離を変えて図 2.6(a)の開口からの伝搬計算を行った結果を図 2.8に示す．用いたソースコードは以下 のとおりである．
Example 伝搬距離を変えた回折像計算
1 # i n c l u d e < wfl . h >
2 u s i n g n a m e s p a c e wfl ;
3
4 v o i d m a i n ()
5 {
6 S t a r t ();
7 W a v e F i e l d a (256 , 256 , 2 e - 6 ) ;
8
9 a . S e t G a u s s i a n ( 0 . 1 e -3 , 5 0 ) ;
10 a . S a v e A s B M P ( " A p e r t u r e . bmp " , A M P L I T U D E ); //開口形状の保存 11
12 // = = = = = 角 ス ペ ク ト ル 法 に よ る 伝 搬 計 算
13 for ( int i = 1; i <= 3; i ++)
14 {
15 W a v e F i e l d b = a ; //オブジェクトbに 開 口aを コ ピ ー す る
16 d o u b l e d = 10 e -3 * i ; // dは 伝 搬 距 離
17 b . A s m P r o p ( d ); // bを 距 離dだ け 伝 搬 す る
18 c h a r f n a m e [ 8 0 ] ; //伝搬距離を入れたファイル名を作成
19 s p r i n t f ( fname , " D i f f r a c t i o n -% d . bmp " , ( int )( d /1 e - 3 ) ) ;
20 b . N o r m a l i z e ();
21 b . S a v e A s B m p ( fname , A M P L I T U D E ); //回折像(振 幅 像)の 保 存
22 }
23 }
( a ) - 1 0 . b m p ( b ) D i f f r a c t i o n - 2 0 . b m p ( b ) D i f f r a c t i o n - 3 0 . b m p
図2.8 AsmProp()で伝搬距離を増加した場合の回折像．伝搬距離を，(a)10mm，(b)20mm，(c)30mm とした結果．
*1このような関数はしばしばスーパーガウス関数と呼ばれる．
15行目では二つ目のWaveFieldオブジェクトとしてbを定義し，開口データの入っているaをコピー しているが，これはAsmProp()メンバー関数の実行によりそのオブジェクトの内容が回折計算結果に置 き換わり，元の開口データが破壊されるためである．
図 2.8(c)では，回折パターンが円形でなくなってきていることがわかる．円形開口からの回折像とし
てはこれは明らかにおかしい．これは，円状畳み込み演算のいわゆる“端効果”により，サンプリング領 域の端で畳み込みが正しく行われていないためである．
この様な端効果を防ぐためには，回折計算前にサンプリング数を4倍(縦横2倍)に拡張して周辺を0 で埋め，回折計算後に中心の4分の1を切り取って元に戻せばよい．その処理にはWaveFieldクラスの Embed()メンバー関数とExtract()メンバー関数を用いる．すなわち，次のソースのとおり，これらで AsmProp()メンバー関数を挟めばよい．
Example より正確な回折像計算(ソースの一部のみ記載)
1 // = = = = = 角 ス ペ ク ト ル 法 に よ る 伝 搬 計 算
2 for ( int i = 1; i <= 3; i ++)
3 {
4 W a v e F i e l d b = a ;
5 d o u b l e d = 10 e -3 * i ;
6 b . E m b e d ( 1 ) ; //サンプリング領域を4倍 拡 張
7 b . A s m P r o p ( d ); //伝搬する
8 b . E x t r a c t ( 1 ) ; //サンプリング領域を4分 の1に 縮 小
9 c h a r f n a m e [ 8 0 ] ;
10 s p r i n t f ( fname , " D i f f r a c t i o n -% d . bmp " , ( int )( d /1 e - 3 ) ) ;
11 b . N o r m a l i z e ();
12 b . S a v e A s B m p ( fname , A M P L I T U D E );
13 }
( a ) - 1 0 . b m p ( b ) D i f f r a c t i o n - 2 0 . b m p ( b ) D i f f r a c t i o n - 3 0 . b m p
図2.9 4倍拡張して伝搬距離を増加した場合の回折像．伝搬距離を，(a)10mm，(b)20mm，(c)30mm とした結果．
この場合の結果を図 2.9に示す．正しく円形の回折パターンになっていることがわかる．
Note
• 上に示したとおり，4倍拡張・縮小を行うとより正確な計算が行えるが，その代償として一時的に オブジェクトの4倍のメモリを消費する．また，計算時間も4倍以上遅くなることに注意しなけれ ばならない．
• 4倍拡張と縮小を一つの関数で行うのが，ExactAsmProp()メンバー関数である．
• 畳み込み演算におけるFFTを分離した関数としてAsmPropFs()メンバー関数がある．フーリエ
46 第2章 チュートリアル 空間で畳み込み以外の処理が必要な場合はこちらが便利である．