(57)【要約】【課題】  電磁波が入射したときの母材中に散布された粒子による多重散乱波の強度を高速、高精度に計算すること。【解決手段】  粒子４と母材５とによって構成される材質２に電磁波６が入射したときの粒子４による多重散乱波の強度の計算方法において、電磁波６が材質２中では粒子４の屈折率と母材５の屈折率と粒子４の体積濃度から決まる有効電磁波として伝播するとして、有効電磁波の粒子単体による散乱波をＴ−ｍａｔｒｉｘを用いて計算し、全ての粒子４による散乱波の干渉を、粒子対確率密度分布関数を用いて統計的に足し合わせることにより観測方向の多重散乱波の強度を算出する。【選択図】  図１
2013218581.tif 000002
【特許請求の範囲】
【請求項１】
粒子と該粒子が散布されている母材とによって構成される材質に電磁波が入射したときの前記粒子による多重散乱波の強度の計算方法において、
前記粒子の半径、前記粒子の屈折率、前記母材の屈折率、前記材質の体積に対する前記粒子の体積濃度、前記電磁波の波長に基づいて、前記電磁波が前記材質中では前記粒子の屈折率と前記母材の屈折率と前記体積濃度から決まる有効電磁波として伝播するとして、前記有効電磁波の前記粒子単体による散乱波をＴ−ｍａｔｒｉｘを用いて計算し、全ての前記粒子による散乱波の干渉を、粒子対確率密度分布関数を用いて統計的に足し合わせることにより前記観測方向の多重散乱波の強度を算出することを特徴とする計算方法。
【請求項２】
粒子と該粒子が散布されている母材とによって構成される材質に電磁波が入射したときの前記粒子による多重散乱波の計算方法において、
前記粒子の半径をａ、前記粒子の屈折率をｎａ、前記母材の屈折率をｎｅ、前記材質の体積に対する前記粒子の体積濃度をｆ、前記母材中の波数をｋ、入射波の進行方向の前記材質の厚さをｌ、前記多重散乱波の観測方向を極座標表示したときの散乱角をθｓ，φｓとするとき、
前記多重散乱波の強度を以下の式を用いて算出することを特徴とする計算方法。
【数１】
2013218581.tif 000042
但し、ｎｅｆｆは前記粒子と前記母材によって決まる有効屈折率、ｋｅｆｆは前記粒子と前記母材によって決まる前記材質中の波数、Ｈ（ｐ）は粒子対確率密度分布関数の３次元フーリエ変換である。また、Ａ（θｓ，φｓ）は以下の式で表される。
【数２】
2013218581.tif 000043
【数３】
2013218581.tif 000044
【数４】
2013218581.tif 000045
【数５】
2013218581.tif 000046
【数６】
2013218581.tif 000047
【数７】
2013218581.tif 000048
【数８】
2013218581.tif 000049
【数９】
2013218581.tif 000050
【数１０】
2013218581.tif 000051
ここで、Ｔは粒子のＴ−ｍａｔｒｉｘでありａ、ｎａ、ｎｅ、ｋに依存し、Ｐｎｍ（ｃｏｓθｓ）はルジャンドル陪関数である。
【請求項３】
Ｈ（ｐ）が以下の式で表わされることを特徴とする請求項２に記載の計算方法。
【数１１】
2013218581.tif 000052
ただし、Ｃ（ｐ）は以下の式で表わされる。
【数１２】
2013218581.tif 000053
【数１３】
2013218581.tif 000054
【数１４】
2013218581.tif 000055
【数１５】
2013218581.tif 000056
【請求項４】
前記粒子と前記母材はそれぞれ粒子分散型の光学素子を構成する粒子と母材であり、前記粒子分散型の光学素子に入射する入射波を平面波として、前記多重散乱波の強度を計算することを特徴とする請求項２又は３に記載の計算方法。
【請求項５】
粒子と該粒子が散布されている母材とによって構成される材質に電磁波が入射したときの前記粒子による多重散乱波の計算方法において、
前記粒子の半径をａ、前記粒子の屈折率をｎａ、前記母材の屈折率をｎｅ、前記材質の体積に対する前記粒子の体積濃度をｆ、前記母材中の波数をｋ、入射波の進行方向の前記材質の厚さをｌ、前記多重散乱波の観測方向を極座標表示したときの散乱角をθｓ，φｓとするとき、
前記散乱角θｓとφｓの範囲内の多重散乱波の強度分布と散乱量を以下の式を用いて計算する計算方法。
【数１６】
2013218581.tif 000057
但し、ｎｅｆｆは前記粒子と前記母材によって決まる有効屈折率、ｋｅｆｆは前記粒子と前記母材によって決まる前記材質中の波数、Ｈ（ｐ）は粒子対確率密度分布関数の３次元フーリエ変換である。また、Ａ（θｓ，φｓ）は以下の式で表される。
【数１７】
2013218581.tif 000058
【数１８】
2013218581.tif 000059
【数１９】
2013218581.tif 000060
【数２０】
2013218581.tif 000061
【数２１】
2013218581.tif 000062
【数２２】
2013218581.tif 000063
【数２３】
2013218581.tif 000064
【数２４】
2013218581.tif 000065
【数２５】
2013218581.tif 000066
ここで、Ｔは粒子のＴ−ｍａｔｒｉｘでありａ、ｎａ、ｎｅ、ｋに依存し、Ｐｎｍ（ｃｏｓθｓ）はルジャンドル陪関数である。
【請求項６】
粒子と該粒子が散布されている母材とによって構成される材質に電磁波が入射したときの前記粒子による多重散乱波の計算方法において、
前記粒子の半径をａ、前記粒子の屈折率をｎａ、前記母材の屈折率をｎｅ、前記材質の体積に対する前記粒子の体積濃度をｆ、前記母材中の波数をｋ、入射波の進行方向の前記材質の厚さをｌ、前記多重散乱波の観測方向を極座標表示したときの散乱角をθｓ，φｓとするとき、
前記多重散乱波の強度を以下の式を用いて計算した計算結果と、粒子対確率密度分布関数を用いずに計算した計算結果との差分が許容範囲となるように、前記粒子対確率密度分布関数を用いない計算における入射波の進行方向に垂直な前記材質の幅ｗを決定することを特徴とする計算方法。
【数２６】
2013218581.tif 000067
但し、ｎｅｆｆは前記粒子と前記母材によって決まる有効屈折率、ｋｅｆｆは前記粒子と前記母材によって決まる前記材質中の波数、Ｈ（ｐ）は粒子対確率密度分布関数の３次元フーリエ変換である。また、Ａ（θｓ，φｓ）は以下の式で表される。
【数２７】
2013218581.tif 000068
【数２８】
2013218581.tif 000069
【数２９】
2013218581.tif 000070
【数３０】
2013218581.tif 000071
【数３１】
2013218581.tif 000072
【数３２】
2013218581.tif 000073
【数３３】
2013218581.tif 000074
【数３４】
2013218581.tif 000075
【数３５】
2013218581.tif 000076
ここで、Ｔは粒子のＴ−ｍａｔｒｉｘでありａ、ｎａ、ｎｅ、ｋに依存し、Ｐｎｍ（ｃｏｓθｓ）はルジャンドル陪関数である。
【請求項７】
前記粒子対確率密度分布関数を用いない計算は、Ｔ−ｍａｔｒｉｘ法であることを特徴とする請求項６に記載の計算方法。
【請求項８】
粒子と該粒子が散布されている母材とによって構成される材質で構成される粒子分散型の光学素子の設計方法において、
前記粒子の半径をａ、前記粒子の屈折率をｎａ、前記母材の屈折率をｎｅ、前記材質の体積に対する前記粒子の体積濃度をｆ、前記母材中の波数をｋ、入射波の進行方向の前記材質の厚さをｌ、材質に電磁波が入射したときの前記粒子による多重散乱波の観測方向を極座標表示したときの散乱角をθｓ，φｓとするとき、
以下の式を用いて計算した前記多重散乱波の散乱角θｓ，φｓの範囲内の散乱量が許容範囲となるように前記粒子の半径ａ、前記母材の屈折率ｎｅ、前記材質の体積に対する前記粒子の体積濃度ｆ、前記母材中の波数ｋ、前記材質の厚さｌのうち少なくとも１つの値を設定することを特徴とする粒子分散型の光学素子の設計方法。
【数３６】
2013218581.tif 000077
但し、ｎｅｆｆは前記粒子と前記母材によって決まる有効屈折率、ｋｅｆｆは前記粒子と前記母材によって決まる前記材質中の波数、Ｈ（ｐ）は粒子対確率密度分布関数の３次元フーリエ変換である。また、Ａ（θｓ，φｓ）は以下の式で表される。
【数３７】
2013218581.tif 000078
【数３８】
2013218581.tif 000079
【数３９】
2013218581.tif 000080
【数４０】
2013218581.tif 000081
【数４１】
2013218581.tif 000082
【数４２】
2013218581.tif 000083
【数４３】
2013218581.tif 000084
【数４４】
2013218581.tif 000085
【数４５】
2013218581.tif 000086
ここで、Ｔは粒子のＴ−ｍａｔｒｉｘでありａ、ｎａ、ｎｅ、ｋに依存し、Ｐｎｍ（ｃｏｓθｓ）はルジャンドル陪関数である。
【請求項９】
前記多重散乱波の散乱角θｓ，φｓの範囲内の散乱量が入射波のエネルギーの１％以下である請求項８に記載の粒子分散型の光学素子の設計方法。
【請求項１０】
粒子と該粒子が散布されている母材とによって構成される材質で構成される粒子分散型の光学素子の設計方法において、
前記粒子の半径をａ、前記粒子の屈折率をｎａ、前記母材の屈折率をｎｅ、前記材質の体積に対する前記粒子の体積濃度をｆ、前記母材中の波数をｋ、入射波の進行方向の前記材質の厚さをｌ、材質に電磁波が入射したときの前記粒子による多重散乱波の観測方向を極座標表示したときの散乱角をθｓ，φｓとするとき、
散乱角θｓ，φｓの範囲内の前記多重散乱波の以下の式を用いて計算した計算結果と、粒子対確率密度分布関数を用いずに計算した計算結果との差分が許容範囲となるように、前記粒子対確率密度分布関数を用いない計算における入射波の進行方向に垂直な前記材質の幅ｗを決定し、前記粒子対確率密度分布関数を用いない計算によって計算した前記多重散乱波の散乱角θｓ，φｓの範囲内の散乱量が許容範囲になるように前記粒子の半径ａ、前記母材の屈折率ｎｅ、前記材質の体積に対する前記粒子の体積濃度ｆ、前記母材中の波数ｋ、前記材質の厚さｌのうち少なくとも１つの値を設定することを特徴とする粒子分散型の光学素子の設計方法。
【数４６】
2013218581.tif 000087
但し、ｎｅｆｆは前記粒子と前記母材によって決まる有効屈折率、ｋｅｆｆは前記粒子と前記母材によって決まる前記材質中の波数、Ｈ（ｐ）は粒子対確率密度分布関数の３次元フーリエ変換である。また、Ａ（θｓ，φｓ）は以下の式で表される。
【数４７】
2013218581.tif 000088
【数４８】
2013218581.tif 000089
【数４９】
2013218581.tif 000090
【数５０】
2013218581.tif 000091
【数５１】
2013218581.tif 000092
【数５２】
2013218581.tif 000093
【数５３】
2013218581.tif 000094
【数５４】
2013218581.tif 000095
【数５５】
2013218581.tif 000096
ここで、Ｔは粒子のＴ−ｍａｔｒｉｘでありａ、ｎａ、ｎｅ、ｋに依存し、Ｐｎｍ（ｃｏｓθｓ）はルジャンドル陪関数である。
【請求項１１】
散乱波の散乱角θｓとφｓの範囲内の散乱量が入射波エネルギーの１％以下である請求項１０に記載の粒子分散型の光学素子の設計方法。
【請求項１２】
前記粒子対確率密度分布関数を用いない計算は、Ｔ−ｍａｔｒｉｘ法であることを特徴とする請求項１０又は１１に記載の粒子分散型の光学素子の設計方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、母材中に散布された複数の粒子からの多重散乱波を計算する方法及びその計算方法を利用した粒子分散型の光学素子の設計方法に関する。
【背景技術】
【０００２】
様々な光学分野において、波長未満の大きさを持つ散乱体からの散乱を解析、シミュレーションする技術が盛んになってきている。特に、複数の散乱体が存在する様な系では散乱体間の多重散乱を扱う必要がある。例えば粒子分散型光学素子の一つである屈折率分布型素子（Ｇｒａｄｉｅｎｔ  Ｉｎｄｅｘ）では母材と異なる屈折率を持つ粒子を母材内に散布させて屈折率を制御している。
【０００３】
ところが、粒子により散乱が起こってしまうため、その多重散乱量を見積もるために解析、シミュレーションが必要となる。散乱を解析、シミュレーションするための電磁場解析計算の技術には、ある条件の下で解析的に散乱波を解く手法や、数値的なシミュレーションにより散乱波を解く手法が考えられている。現実のモデルとの対応を取って解析、シミュレーションを実施する際には大規模領域に無数の散乱体が散布された状況での多重散乱の扱いが必要となる事もある。
【０００４】
散乱体による散乱計算を電磁場解析計算により行う場合、様々な手法が知られている。例えばＭｉｅ散乱による手法では孤立の球体粒子による散乱を高速、高精度に計算する事ができる。ＦＤＴＤ法（Ｆｉｎｉｔｅ−ｄｉｆｆｅｒｅｎｃｅ  ｔｉｍｅ−ｄｏｍａｉｎ  ｍｅｔｈｏｄ）や境界要素法に代表される電磁場解析手法では解析領域の全てまたは一部をメッシュ状に分割し、その格子点の全てに誘電率を定義することにより光の伝播計算を行う。
【０００５】
またＤＤＡ（ｄｉｓｃｒｅｔｅ−ｄｉｐｏｌｅ  ａｐｐｒｏｘｉｍａｔｉｏｎ）も注目されており、散乱体を双極子の集合体だとみなし分割する事により複数の散乱体による多重散乱を扱える（非特許文献１）。またＴ−ｍａｔｒｉｘ法においては、複数の球形粒子による多重散乱を厳密に扱える（非特許文献２）。さらにＥＦＡ（Ｅｆｆｅｃｔｉｖｅ  Ｆｉｅｌｄ  Ａｐｐｒｏｘｉｍａｔｉｏｎ）やＱＣＡ（Ｑｕａｓｉ−ｃｒｙｓｔａｌｌｉｎｅ  Ａｐｐｒｏｘｉｍａｔｉｏｎ）では粒子多重散乱を粒子配置の確率密度分布関数を用いる事により統計的に扱い、粒子が散布されている領域を進行する平均場を解析的に求める事ができる（非特許文献３）。
【０００６】
また複数の電磁場解析の手法を組み合わせて用いる事で各々の電磁場解析の長所を生かして多重散乱を解析、シミュレーションする手法も公開されている。例えば、Ｍｉｅ散乱の理論により孤立球体の散乱波の角度特性を求めておき、それを球体に電磁場が入射した時の光子進行方向の確率分布だとして定義し、光子を追跡する事により球体が複数存在する時の多重散乱波を算出するモンテカルロシミュレーション等がある。また孤立散乱体の散乱波をＦＤＴＤ法を用いて計算し、散乱体が複数ある時の多重散乱をＲＴＥ（Ｒａｄｉａｔｉｖｅ  Ｔｒａｎｓｆｅｒ  Ｅｑｕａｔｉｏｎ）を用いて計算する方法も開示されている（特許文献１）。
【先行技術文献】
【非特許文献】
【０００７】
【非特許文献１】Ｈ．ＤｅＶｏｅ，“Ｏｐｔｉｃａｌ  ｐｒｏｐｅｒｔｉｅｓ  ｏｆ  ｍｏｌｅｃｕｌａｒ  ａｇｇｒｅｇａｔｅｓ．Ｉ．Ｃｌａｓｓｉｃａｌ  ｍｏｄｅｌ  ｏｆ  ｅｌｅｃｔｒｏｎｉｃ  ａｂｓｏｒｐｔｉｏｎ  ａｎｄ  ｒｅｆｒａｃｔｉｏｎ，”Ｊ．Ｃｈｅｍ．Ｐｈｙｓ．４１，３９３−４００（１９６４）
【非特許文献２】Ｗａｔｅｒｍａｎ，Ｐ．Ｃ．，“Ｍａｔｒｉｘ  ｆｏｒｍｕｌａｔｉｏｎ  ｏｆ  ｅｌｅｃｔｒｏｍａｇｎｅｔｉｃ  ｓｃａｔｔｅｒｉｎｇ，”Ｐｒｏｃ．ＩＥＥＥ，５３，８０５−８１１（１９６５）
【非特許文献３】Ｌｅｕｎｇ  Ｔｓａｎｇ，Ｊｕｎ  Ａｕ  Ｋｏｎｇ，Ｋｕｎｇ−Ｈａｕ  Ｄｉｎｇ，ａｎｄ  Ｃｈｉ  Ｏｎ  Ａｏ，Ｓｃａｔｔｅｒｉｎｇ  ｏｆ  Ｅｌｅｃｔｒｏｍａｇｎｅｔｉｃ  Ｗａｖｅｓ  Ａｄｖａｎｃｅｄ  Ｔｏｐｉｃｓ，Ｗｉｌｅｙ，Ｎｅｗ  Ｙｏｒｋ，２０７−２１３
【特許文献】
【０００８】
【特許文献１】特開２０１１−２２８８９号公報
【発明の概要】
【発明が解決しようとする課題】
【０００９】
電磁場解析による散乱計算において従来の手法では、前提条件や原理から限られた条件に対しては高精度に計算できるが、実際の系に適用した場合には高精度に計算できないという問題があった。例えばＭｉｅ散乱による手法では孤立の球体しか扱えず、粒子が複数存在する場合の多重散乱を扱う事ができない。
【００１０】
また計算コストの問題から手法によってはモデルが計算機環境に大きく制限されるものがある。例えばＦＤＴＤ法や境界要素法に代表される電磁場解析計算においては解析領域の全てまたは一部をメッシュ状に分割しなければならない事から、計算時間やコンピュータのメモリ消費量が膨大になってしまう。また、ＤＤＡにおいては散乱体内部に微小な双極子を密に配置しなければならないため、やはり計算コストが膨大になってしまう。
【００１１】
これによりＦＤＴＤ法、境界要素法、ＤＤＡ等では設定できる粒子の数が制限されたり、粒子の構造を表現するために十分なメッシュ分割や双極子の配置ができなくなる事がある。また、これらの手法で粒子分散型光学素子等の散乱を解析しようとすると、解析対象のサイズが巨視的であるために有限の領域設定によるモデル化では高精度に計算できないことがある。
【００１２】
この問題に対処するために周期境界条件を課して疑似的に解析領域の大規模化を行う事がある。しかし、これによって散乱波にも周期性が現れてしまい、やはり高精度に計算ができなくなってしまう。また、解析対象のサイズが巨視的であると粒子配置のバリエーションを考慮しなければならず、粒子配置を何度も変更して繰り返し計算を行う必要があるため、計算時間が膨大となってしまう。
【００１３】
近年、粒子の多重散乱解析手法として注目されているＴ−ｍａｔｒｉｘ法では個々の粒子による散乱は厳密に解かれているものの、粒子数が大きくなってくるとやはり繰り返し計算が多数回必要となるため計算コストが膨大になってしまう。またモンテカルロシミュレーションにおいては粒子間の距離を平均自由光路長を用いて確率的に扱うため粒子配置を行う必要性がないものの、光線を無数に入射させる必要があるためやはり計算時間は膨大となってしまう。
【００１４】
本発明は、これらの従来の計算方法を鑑みなされたもので、電磁波が入射したときの粒子による多重散乱波を高速、高精度に計算することを目的とする。
【課題を解決するための手段】
【００１５】
本発明の１つの観点での計算方法は、粒子とその粒子が散布されている母材とによって構成される材質に電磁波が入射したときの粒子による多重散乱波の強度の計算方法において、
粒子の半径、粒子の屈折率、母材の屈折率、材質の体積に対する粒子の体積濃度、電磁波の波長に基づいて、電磁波が材質中では粒子の屈折率と母材の屈折率と体積濃度から決まる有効電磁波として伝播するとして、有効電磁波の粒子単体による散乱波をＴ−ｍａｔｒｉｘを用いて計算し、全ての粒子による散乱波の干渉を、粒子対確率密度分布関数を用いて統計的に足し合わせることにより観測方向の多重散乱波の強度を算出することを特徴としている。
【００１６】
本発明の別の観点での計算方法は、粒子とその粒子が散布されている母材とによって構成される材質に電磁波が入射したときの粒子による多重散乱波の計算方法において、
粒子の半径をａ、粒子の屈折率をｎａ、母材の屈折率をｎｅ、材質の体積に対する粒子の体積濃度をｆ、母材中の波数をｋ、入射波の進行方向の材質の厚さをｌ、多重散乱波の観測方向を極座標表示したときの散乱角をθｓ，φｓとするとき、
多重散乱波の強度を以下の式を用いて算出することを特徴としている。
【００１７】
【数１】
2013218581.tif 000003
【００１８】
但し、ｎｅｆｆは前記粒子と前記母材によって決まる有効屈折率、ｋｅｆｆは前記粒子と前記母材によって決まる前記材質中の波数、Ｈ（ｐ）は粒子対確率密度分布関数の３次元フーリエ変換である。また、Ａ（θｓ，φｓ）は以下の式で表される。
【００１９】
【数２】
2013218581.tif 000004
【００２０】
【数３】
2013218581.tif 000005
【００２１】
【数４】
2013218581.tif 000006
【００２２】
【数５】
2013218581.tif 000007
【００２３】
【数６】
2013218581.tif 000008
【００２４】
【数７】
2013218581.tif 000009
【００２５】
【数８】
2013218581.tif 000010
【００２６】
【数９】
2013218581.tif 000011
【００２７】
【数１０】
2013218581.tif 000012
【００２８】
ここで、Ｔは粒子のＴ−ｍａｔｒｉｘでありａ、ｎａ、ｎｅ、ｋに依存し、Ｐｎｍ（ｃｏｓθｓ）はルジャンドル陪関数である。
【発明の効果】
【００２９】
本発明によれば、複数の粒子による多重散乱波を高速、高精度に計算することができる。
【図面の簡単な説明】
【００３０】
【図１】解析モデルを表す図である。
【図２】極座標とデカルト座標の関係を表す図である。
【図３】散乱強度計算式を用いて特定の角度の散乱強度を算出する場合の手順を表すフローチャートである。
【図４】散乱強度計算式を用いて所定の散乱角範囲の散乱強度角度分布と散乱量を算出する場合の手順を表すフローチャートである。
【図５】数値的解法による電磁場計算を用いて散乱解析をするときに、材質幅を散乱強度計算式を用いて決定する際のユーザーの手順を表すフローチャートである。
【図６】散乱強度計算式を用いて粒子分散型の光学素子を設計する際のユーザーの手順を表すフローチャートである。
【図７】散乱強度計算式で球ベクトル関数展開項数を変化させた時の散乱強度の変化を表すグラフである。
【図８】散乱強度計算式により計算した散乱強度角度分布とＴ−ｍａｔｒｉｘ法を用いて材質の幅を変更して計算した散乱強度角度分布を表すグラフである。
【図９】散乱強度計算式により計算した散乱強度とＴ−ｍａｔｒｉｘ法を用いて材質幅を変更して計算した散乱強度を表すグラフである。
【図１０】散乱強度計算式を用いて計算した粒子半径を変更した時の散乱量の変化を表すグラフである。
【発明を実施するための形態】
【００３１】
以下に、本発明の好ましい実施形態を添付の図面に基づいて詳細に説明する。
【００３２】
図１、図２、図３、図４、図５、図６に基づいて、実施手順を説明する。図１に示しているのは、本発明の計算方法で対象としている解析モデルである。
【００３３】
図１において、１は設定した座標系を示している。２は母材に粒子が散布された領域を示しており、この粒子と粒子が散布されている母材とによって構成されるものを以下材質と呼ぶこととする。ｘ、ｙ方向（入射波の進行方向に垂直な方向）の長さを以下材質幅と呼ぶこととし、Ｔ−ｍａｔｒｉｘ法等の数値的解法による電磁場解析では有限の大きさを持つが、本発明に用いる散乱強度計算式では無限の大きさを持つとして定義してある。また、ｚ方向には３で示した境界で区切られた大きさを持ち、以下材質厚と呼ぶこととする。４は粒子、５は母材を示している。６は電磁波であるところの入射波を示している。入射波４は、本発明に用いる散乱強度計算式では振幅１、ｘ偏光、ｚ方向に進行する平面波である。７は複数の粒子で散乱した多重散乱波を示しており、本発明に用いる散乱強度計算式ではこの多重散乱波の強度を算出する。
【００３４】
図２は極座標とデカルト座標の関係を示している。８は位置（ｘ，ｙ，ｚ）に対してθを（ｘ，ｙ，ｚ）とｚ軸のなす角とし、φを（ｘ，ｙ，ｚ）をｘｙ平面に射影した時にｘ軸となす角である事を示している。
【００３５】
本実施形態においては、電磁波が材質２中では粒子４の屈折率と母材５の屈折率と、粒子４の体積濃度から決まる有効電磁波として伝播する想定している。そして、粒子単体による有効電磁波の散乱波をＴ−ｍａｔｒｉｘを用いて計算し、全ての粒子による散乱波の干渉を、粒子対確率密度分布関数を用いて統計的に足し合わせることにより観測方向の多重散乱波の強度を算出している。
【００３６】
具体的には、本実施形態において散乱強度計算式とは、
粒子の半径：ａ
粒子の屈折率：ｎａ
母材の屈折率：ｎｅ
粒子の分散されている領域の体積に対する粒子の体積濃度：ｆ
材質中の電磁波の波数：ｋ
入射波の進行方向の母材の厚さ：ｌ
散乱波の観測方向を極座標表示した時の散乱角：θｓ，φｓ
球ベクトル関数展開項数：ｎ
とするとき、以下の［数式１］を用いて、粒子とそれが散布された母材との集合体に光を含む電磁波が入射した場合の多重散乱波の強度Ｉｓ（θｓ，φｓ）を計算するものである。
【００３７】
【数１１】
2013218581.tif 000013
【００３８】
ただし、ｎｅｆｆは前記粒子と前記母材によって決まる有効屈折率であり、例として［数式２］、［数式３］、［数式４］、［数式５］の様な計算方法があるが他の計算方法を用いる事も可能である。本実施例では全て［数式４］を用いてｎｅｆｆを計算する。
【００３９】
【数１２】
2013218581.tif 000014
【００４０】
ただし、０＜β＜１である。
【００４１】
【数１３】
2013218581.tif 000015
【００４２】
【数１４】
2013218581.tif 000016
【００４３】
【数１５】
2013218581.tif 000017
【００４４】
ｋｅｆｆは前記粒子と前記母材によって決まる有効媒質中の波数であり、例として［数式６］の様な計算方法があるが他の計算方法を用いる事も可能である。本実施例では全て［数式６］を用いてｋｅｆｆを計算する。
【００４５】
【数１６】
2013218581.tif 000018
【００４６】
Ｈ（ｐ）は粒子対確率密度分布関数の３次元フーリエ変換であり、例として［数式７］の様な計算方法があるが他の計算方法を用いる事も可能である。本実施例では全て［数式７］を用いてＨ（ｐ）を計算する。
【００４７】
【数１７】
2013218581.tif 000019
【００４８】
ただし、Ｃ（ｐ）は以下の式で表わされる。
【００４９】
【数１８】
2013218581.tif 000020
【００５０】
【数１９】
2013218581.tif 000021
【００５１】
【数２０】
2013218581.tif 000022
【００５２】
【数２１】
2013218581.tif 000023
【００５３】
また、Ａ（θｓ，φｓ）は以下の式で表される。
【００５４】
【数２２】
2013218581.tif 000024
【００５５】
【数２３】
2013218581.tif 000025
【００５６】
【数２４】
2013218581.tif 000026
【００５７】
【数２５】
2013218581.tif 000027
【００５８】
【数２６】
2013218581.tif 000028
【００５９】
【数２７】
2013218581.tif 000029
【００６０】
【数２８】
2013218581.tif 000030
【００６１】
【数２９】
2013218581.tif 000031
【００６２】
【数３０】
2013218581.tif 000032
【００６３】
ここで、Ｔは粒子のＴ−ｍａｔｒｉｘでありａ、ｎａ、ｎｅ、ｋに依存し、Ｐｎｍ（ｃｏｓθｓ）はルジャンドル陪関数である。
【００６４】
〔第１実施形態〕
図３は散乱強度計算式［数式１］を用いて、球ベクトル関数展開項数ｎを決定し、特定の角度の散乱強度を算出する場合の手順を示すフローチャートである。
【００６５】
ステップ１ではユーザーが電磁波（入射波）の波長λ、粒子半径ａ、粒子体積濃度ｆ、母材屈折率ｎｅ、粒子屈折率ｎａ、材質厚ｌ、球ベクトル関数展開項数ｎ、散乱角θｓ、散乱角φｓ、散乱強度差分許容量δを入力する。ここで波長λと母材中の波数ｋはｋ＝２πｎｅ／λという関係がある。ステップ２では［数式１］に基づいて散乱強度を計算する。ステップ３ではステップ２で計算した散乱強度の値を保存する。ステップ４では球ベクトル関数展開次数ｎをｎ＋１に変更する。ステップ５では［数式１］を用いて散乱強度を計算する。ステップ６ではステップ５で計算した散乱強度の値を保存する。ステップ７ではステップ３で保存した球ベクトル関数展開項数ｎの時の散乱強度とステップ６で保存した球ベクトル関数展開項数ｎ＋１の時の散乱強度の差分がステップ１で設定した散乱強度差分許容量δより小さいか否かを判定する。差分が散乱強度差分許容量δより小さい場合には手順を終了し、小さくない場合にはステップ２に進む。散乱強度計算式［数式１］を用いる場合には全てこのフローにより決定された球ベクトル関数展開項数ｎを用いるのが望ましく、以下の実施例で用いるｎは全てこの手順で決定したものである。
【００６６】
（実施例１）
表１に示した値を用いて、図３で示したフローチャートに基づいて散乱強度を計算し、球ベクトル関数展開項数ｎを決定した。
【００６７】
【表１】
2013218581.tif 000033
【００６８】
図７は球ベクトル関数展開項数ｎを変化させた時の散乱強度の変化を示している。９は図３で示したフローチャートに基づいて散乱強度を計算した結果である。１０は９の計算結果を用いて球ベクトル関数展開項数ｎでの散乱強度から球ベクトル関数展開項数ｎ−１での散乱強度を引いた値である。これにより球ベクトル関数展開項数ｎを５以上にする事で散乱強度の差分を１．０×１０−５以下にすることができる。
【００６９】
（実施例２）
表２に示した値を用いて、散乱強度計算式［数式１］による計算と数値的解法である電磁場解析による計算で計算時間を比較した。数値的解法である電磁場解析はＴ−ｍａｔｒｉｘ法を採用し、以後の実施例でも全てその様にする。本発明の散乱強度計算式［数式１］で計算時間に寄与するのは球ベクトル関数展開項数ｎだけであり、Ｔ−ｍａｔｒｉｘ法では球ベクトル関数展開項数ｎと粒子数である。Ｔ−ｍａｔｒｉｘ法においては粒子濃度ｆを一定にしたまま粒子半径ａを小さくする、粒子半径ａを一定にしたまま粒子濃度ｆを大きくする、材質幅を大きくする、材質厚を大きくする事によって粒子数が大きくなる事で計算時間が増す。
【００７０】
【表２】
2013218581.tif 000034
【００７１】
表３、表４は計算時間の結果を表わしている。用いた計算機はどちらも同じものでプロセッサはＤｕａｌ−Ｃｏｒｅ  ＡＭＤ  ＯｐｔｅｒｏｎＴＭ  Ｐｒｏｃｅｓｓｏｒ  ２２２０  ２．８ＧＨｚである。
【００７２】
【表３】
2013218581.tif 000035
【００７３】
【表４】
2013218581.tif 000036
【００７４】
球ベクトル関数展開項数ｎは一定なので散乱強度計算式［数式１］ではいかなるパラメータを与えても計算を４秒程度で終える事ができる。しかし、Ｔ−ｍａｔｒｉｘ法では現実的な計算結果を得るために最低限必要な粒子数でも２時間程度の計算時間がかかり、高精度な計算結果を得ようとすると数週間の時間を要する。これにより散乱強度計算式［数式１］がＴ−ｍａｔｒｉｘ法よりも高速計算可能であることが分かる。
【００７５】
〔第２実施形態〕
図４は散乱強度計算式［数式１］を用いて所望の散乱角の範囲内の多重散乱波の強度分布と散乱量を算出する場合の手順を示したフローチャートである。
【００７６】
ステップ１ではユーザーが波長λ、粒子半径ａ、粒子体積濃度ｆ、母材屈折率ｎｅ、粒子屈折率ｎａ、材質厚ｌ、球ベクトル関数展開項数ｎ、散乱角θｓの最小値θｓｍｉｎ、散乱角θｓの最大値θｓｍａｘ、散乱角φｓの最小値φｓｍｉｎ、散乱角φｓの最大値φｓｍａｘ、散乱角θｓの分割数ｉｍａｘ、散乱角φｓの分割数ｊｍａｘを入力する。ここで波長λと母材中の波数ｋはｋ＝２πｎｅ／λという関係がある。ステップ２では散乱量の値をゼロに初期化する。ステップ３では散乱角θｓのループを開始する。ステップ４では散乱角θｓを計算する。ステップ５では散乱角φｓのループを開始する。ステップ６では散乱角φｓを計算する。ステップ７では散乱強度計算式［数式１］に基づいて散乱強度を計算する。ステップ８ではステップ７で計算した散乱強度を出力する。ステップ９では散乱量にステップ７で計算した散乱強度を加算する。ステップ１０では散乱角φｓのループを終了する。ステップ１１では散乱角θｓのループを終了する。ステップ１２では散乱量を出力する。ステップ１でθｓｍｉｎ＝０°、θｓｍａｘ＝１８０°、φｓｍｉｎ＝０°、φｓｍａｘ＝３６０°と入力した場合にはステップ９で算出する散乱量は全方位の総散乱量となる。また、θｓｍｉｎ＝θｓｍａｘ、φｓｍｉｎ＝φｓｍａｘである角度を入力した場合にはステップ９で算出する散乱量は散乱角θｓ＝θｓｍｉｎ、φｓ＝φｓｍｉｎの散乱強度となる。
【００７７】
（実施例３）
表２に示した値を用いて、図４で示したフローチャートに基づいて散乱角度分布を計算し、Ｔ−ｍａｔｉｒｘ法の計算結果と比較した。散乱角の範囲については表５に示した値を用いた。
【００７８】
【表５】
2013218581.tif 000037
【００７９】
Ｔ−ｍａｔｒｉｘ法では各パラメータは表２、表５と同じ値とし、材質幅については表６で示した値を用いた。
【００８０】
【表６】
2013218581.tif 000038
【００８１】
図８は散乱角度分布の計算結果を示している。１１は図４で示したフローチャートに基づいて散乱強度角度分布を計算した結果である。１２は材質幅が２７８ｎｍの時のＴ−ｍａｔｒｉｘ法の計算結果である。１３は材質幅が４８１ｎｍの時のＴ−ｍａｔｒｉｘ法の計算結果である。１４は材質幅が９９２ｎｍの時のＴ−ｍａｔｒｉｘ法の計算結果である。１５は材質幅が１９９９ｎｍの時のＴ−ｍａｔｒｉｘ法の計算結果である。Ｔ−ｍａｔｒｉｘ法で材質幅を大きくするにつれて散乱角度分布が散乱強度計算式［数式１］による計算結果に近づいているのが分かる。これにより材質幅が大きい時の散乱角度分布を計算する場合には散乱強度計算式［数式１］による計算の方が数値的解法である電磁場解析よりも高精度であることが分かる。
【００８２】
〔第３実施形態〕
図６は数値的解法による電磁場計算で散乱解析をする時に、散乱強度計算式［数式１］を用いて材質幅を決定する際のユーザーの手順を示しているフローチャートである。
【００８３】
ステップ１ではユーザーが波長λ、粒子半径ａ、粒子体積濃度ｆ、母材屈折率ｎｅ、粒子屈折率ｎａ、材質厚ｌ、球ベクトル関数展開項数ｎ、散乱角θｓ、散乱角φｓを入力する。ここで波長λと母材中の波数ｋはｋ＝２πｎｅ／λという関係がある。ステップ２では散乱強度計算式［数式１］に基づいて散乱強度を計算する。ステップ３ではステップ２で計算した散乱強度の値を保存する。ステップ４ではステップ１で入力した値の他に材質幅ｗ、数値的解法で粒子配置変更による散乱量平均値の変化量許容値δ１、数値的解法でレンズ幅変更した時の散乱量の散乱強度計算式［数式１］による結果との差分許容値δ２、材質幅変化量ｄｗを数値的解法による電磁場計算ソフトウェア等に入力する。ステップ５では０番目の粒子配置の散乱強度にゼロを代入する。ステップ６では粒子配置変更のループを開始する。ステップ７ではｉｒ番目の粒子配置に対して数値的解法により散乱強度を計算する。ステップ８ではｉｒ番目までの粒子配置のＥｎｓｅｍｂｌｅ平均により求めた散乱強度の値を保存する。ステップ９ではｉｒ番目までの粒子配置のＥｎｓｅｍｂｌｅ平均により求めた散乱強度とｉｒ−１番目までの粒子配置のＥｎｓｅｍｂｌｅ平均により求めた散乱強度の差が値δ１より小さい許容範囲であるか否かを判定し、δ１より小さいならステップ１０に進み粒子配置変更のループを終了する。ステップ１１ではステップ３で保存した散乱強度とｉｒ番目までの粒子配置のＥｎｓｅｍｂｌｅ平均により求めた散乱強度の差が値δ２より小さい許容範囲であるか否かを判定し、δ２より小さいなら処理を終了し、小さくないならステップ１２に進む。ステップ１２では材質幅をｄｗだけ変更してステップ６に進む。
【００８４】
（実施例４）
表２に示した値を用いて、図６で示したフローチャートに基づいて指定した散乱角での散乱強度が散乱強度計算式［数式１］による計算結果と数値的解法である電磁場解析の計算結果が所望の差分に収まる様に材質幅を決定した。散乱角は表７の値を用いた。
【００８５】
【表７】
2013218581.tif 000039
【００８６】
Ｔ−ｍａｔｒｉｘ法では、散乱角以外の各パラメータは表２と同じ値、散乱角は表７と同じ値とし、その他の値については表８で示した値を用いた。ただし、材質幅ｗは５００ｎｍを初期値として材質幅変化量ｄｗを５００ｎｍとしたが、実際に用いた値は粒子数をＮとした時に４πａ３Ｎ／（３ｌｗ２）＝ｆの関係を満たす様に材質幅ｗを調整した値となっている。
【００８７】
【表８】
2013218581.tif 000040
【００８８】
図９は散乱強度の計算結果を示している。横軸はＴ−ｍａｔｒｉｘ法による計算時の材質幅を表している。１６は散乱強度計算式［数式１］により散乱強度を計算した結果である。１７はＴ−ｍａｔｒｉｘ法の計算結果である。１８は散乱強度計算式［数式１］による散乱強度の計算結果とＴ−ｍａｔｒｉｘ法による散乱強度の計算結果の差分を散乱強度計算式［数式１］で割った値を百分率で表したものである。これによりＴ−ｍａｔｒｉｘ法による散乱量の計算結果は材質幅を大きくする程、散乱強度計算式［数式１］のそれに近づいている事が分かる。Ｔ−ｍａｔｒｉｘ法による散乱強度の計算結果と散乱強度計算式［数式１］による散乱強度の計算結果の差異は、材質幅を９回変更すると８．０Ｅ−８以下になり計算が終了する。この時の材質幅は２０００ｎｍであり、Ｔ−ｍａｔｒｉｘ法によって計算する時の材質幅が決定された。
【００８９】
〔第４実施形態〕
図７は散乱強度計算式［数式１］を用いて粒子分散型の光学素子を設計する際のユーザーの手順を示しているフローチャートである。
【００９０】
ステップ１ではユーザーが波長λ、粒子半径ａ、粒子体積濃度ｆ、母材屈折率ｎｅ、粒子屈折率ｎａ、材質厚ｌ、球ベクトル関数展開項数ｎ、散乱角θｓの最小値θｓｍｉｎ、散乱角θｓの最大値θｓｍａｘ、散乱角φｓの最小値φｓｍｉｎ、散乱角φｓの最大値φｓｍａｘ、散乱角θｓの分割数ｉｍａｘ、散乱角φｓの分割数ｊｍａｘ、指定した角度範囲の散乱量の許容量εを入力する。ここで波長λと母材中の波数ｋはｋ＝２πｎｅ／λという関係がある。ここで入力する球ベクトル関数展開項数ｎは図３のフローによって決定されたｎを用いるのが望ましい。また、指定した角度範囲の散乱量の許容量εは設計する粒子分散型素子を組み込んだ光学システムのフレア量の許容値から決定する事が可能である。ステップ２では指定した角度範囲の散乱量の値をゼロに初期化する。ステップ３では散乱角θｓのループを開始する。ステップ４では散乱角θｓを計算する。ステップ５では散乱角φｓのループを開始する。ステップ６では散乱角φｓを計算する。ステップ７では散乱強度計算式［数式１］を用いて散乱強度を計算する。ステップ８ではステップ７で計算した散乱強度の値を保存する。ステップ９では指定した角度範囲の散乱量の値を算出する。ステップ１０では散乱角φｓのループを終了する。ステップ１１では散乱角θｓのループを終了する。ステップ１２では指定した角度範囲の散乱量が値εより小さい許容範囲であるか否かを判定し、εより小さいならば終了し、小さくないならばステップ１３に進む。ステップ１３ではどのパラメータを変更するかを選択し、粒子半径を変更するならステップ１４に進み、粒子屈折率を変更するならステップ１５に進み、粒子体積濃度を変更するならステップ１６に進み、母材屈折率を変更するならステップ１７に進む。ステップ１４ではユーザーが粒子半径を小さくする量ｄａを指定し、変更した粒子半径を入力値としてステップ２に進む。ステップ１５ではユーザーが粒子屈折率を変化させる量ｄｎａを指定し、変更した粒子屈折率を入力値としてステップ２に進む。ステップ１６ではユーザーが粒子体積濃度を小さくする量ｄｆを指定し、変更した粒子体積濃度を入力値としてステップ２に進む。ステップ１７ではユーザーが母材屈折率を変化させる量ｄｎｅを指定し、変更した母材屈折率を入力値としてステップ２に進む。
【００９１】
（実施例５）
表９に示した値を用いて、図６で示したフローチャートに基づいて所望の散乱角θｓ、φｓの範囲の散乱量が入射波エネルギーの１％以下になる粒子半径を算出した。
【００９２】
【表９】
2013218581.tif 000041
【００９３】
図１０は粒子半径ａを５０ｎｍ、４０ｎｍ、３０ｎｍ、２０ｎｍ、１０ｎｍ、９ｎｍ、８ｎｍ、７ｎｍ、５ｎｍと変化させた時の散乱角θｓが０°〜１０°、φｓが０°〜３６０°の範囲の散乱量の計算結果を示している。１９は散乱強度計算式［数式１］により計算した散乱量の結果である。粒子半径ａが９ｎｍ以下になると散乱量が入射波のエネルギーの１％以下となることが分かる。つまり、粒子分散型光学素子において、同条件で設計する場合には粒子半径ａを９ｎｍ以下に設計する必要がある。
【図１】
2013218581.tif 000097
【図２】
2013218581.tif 000098
【図３】
2013218581.tif 000099
【図４】
2013218581.tif 000100
【図５】
2013218581.tif 000101
【図６】
2013218581.tif 000102
【図７】
2013218581.tif 000103
【図８】
2013218581.tif 000104
【図９】
2013218581.tif 000105
【図１０】
2013218581.tif 000106