(57)【要約】【課題】波動の伝播の初期段階の、入射波が全領域に行き渡らない間の計算量を低減し、計算を高速化することができる波動場解析方法、装置および波動場解析方法をコンピュータに実行させるプログラムを記録したコンピュータに読み取り可能な記録媒体を提供する。【解決手段】事前に解析領域内における波動が到達している到達領域を予測し、所定の計算ステップにおける波動場の時間変化の計算を、当該計算ステップに解析領域内の波動が到達していると予測された到達領域を当該計算ステップの空間計算範囲として行うことにより、上記課題を解決する。【選択図】図６
2008250956.tif 000002
【特許請求の範囲】
【請求項１】
波動が伝播する解析領域内の波動場の時間変化の計算を初期状態から各計算ステップ毎に連続して順次行い、前記解析領域内を前記波動が伝播する伝播状態を解析する波動場解析方法であって、
事前に前記解析領域内における前記波動が到達している到達領域を予測し、
所定の計算ステップにおける前記波動場の時間変化の計算を、当該計算ステップに前記解析領域内の前記波動が到達していると予測された前記到達領域を当該計算ステップの空間計算範囲として行うことを特徴とする波動場解析方法。
【請求項２】
前記波動は、前記解析領域の一端から略前記解析領域の所定の軸方向へ伝播するものである請求項１に記載の波動場解析方法。
【請求項３】
前記波動の前記到達領域の予測は、前記解析領域を伝播する前記波動の伝播速度を予測することにより行われるものである請求項１または２に記載の波動場解析方法。
【請求項４】
前記波動の前記到達領域の予測は、前記解析領域内の屈折率分布の情報に基づいて行われるものである請求項１〜３のいずれかに記載の波動場解析方法。
【請求項５】
前記解析領域内の屈折率分布は、前記到達領域内における前記波動の伝播方向と直交する面の屈折率の代表値を、この前記波動の伝播方向と直交する面内において最も小さい屈折率とすることにより得られるものである請求項４に記載の波動場解析方法。
【請求項６】
前記波動の前記到達領域の予測は、前記解析領域内の屈折率を１であると仮定して行われるものである請求項１〜４のいずれかに記載の波動場解析方法。
【請求項７】
前記波動の前記到達領域の予測は、前記解析領域の前記一端から前記波動が伝播する伝播長を計算することにより行われるものである請求項１〜６のいずれかに記載の波動場解析方法。
【請求項８】
前記波動は、電磁波であり、前記波動場は、電磁場である請求項１〜７のいずれかに記載の波動場解析方法。
【請求項９】
波動が伝播する解析領域内の波動場の時間変化の計算を初期状態から各計算ステップ毎に連続して順次行い、前記解析領域内を前記波動が伝播する伝播状態を解析する波動場解析装置であって、
前記波動が伝播する前記解析領域内の波動場の時間変化を初期状態から連続して計算する波動場計算手段と、
事前に、前記解析領域内における前記波動が到達している到達領域を予測する予測手段と、
前記波動場計算手段は、所定の計算ステップにおける前記波動場の時間変化の計算を、前記所定の計算ステップに前記解析領域内の前記波動が到達していると予測された前記到達領域を当該計算ステップの空間計算範囲として行うことを特徴とする波動場解析装置。
【請求項１０】
前記波動は、前記解析領域の一端から略前記解析領域の所定の軸方向へ伝播するものであり、前記波動場計算手段は、前記所定の軸方向に沿って前記波動場の時間変化を計算するものである請求項９に記載の波動場解析装置。
【請求項１１】
前記予測手段は、前記波動の前記到達領域の予測を、前記解析領域を伝播する前記波動の伝播速度を予測することにより行う請求項９または１０に記載の波動場解析装置。
【請求項１２】
前記予測手段は、前記波動の前記到達領域の予測を、前記解析領域内の屈折率分布の情報に基づいて行う請求項９〜１１のいずれかに記載の波動場解析装置。
【請求項１３】
前記予測手段は、前記到達領域内における前記波動の伝播方向と直交する面の屈折率の代表値を、この前記波動の伝播方向と直交する面内において最も小さい屈折率として、前記解析領域内の屈折率分布を得、得られた前記解析領域内の屈折率分布に基づいて前記波動の前記到達領域の予測を行う請求項９〜１２のいずれかに記載の波動場解析装置。
【請求項１４】
前記予測手段は、前記波動の前記到達領域の予測を、前記解析領域内の屈折率を１であると仮定して行う請求項９〜１２のいずれかに記載の波動場解析装置。
【請求項１５】
前記予測手段は、前記波動の前記到達領域の予測を、前記解析領域の前記一端から前記波動が伝播する伝播長を計算することにより行う請求項９〜１４のいずれかに記載の波動場解析装置。
【請求項１６】
前記波動は、電磁波であり、前記波動場は、電磁場である請求項９〜１５のいずれかに記載の波動場解析装置。
【請求項１７】
請求項１〜８のいずれかに記載の波動場解析方法をコンピュータに実行させるプログラムを記録したことを特徴とするコンピュータに読み取り可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、波動場解析方法および装置ならびに波動場解析方法をコンピュータに実行させるプログラムを記録したコンピュータに読み取り可能な記録媒体に関し、詳しくは、ＦＤＴＤ法（Finite Difference Time Domain method（時間領域差分法））を用いて、電磁波などの波動の伝播、すなわち電磁場（電磁界）などの波動場の時間的変化を数値解析する波動場解析方法および波動場解析装置ならびに波動場解析方法をコンピュータに実行させるプログラムを記録したコンピュータ読み取り可能な記録媒体に関する。
【背景技術】
【０００２】
近年、光ディスク、レンズ、特にマイクロレンズなどの種々の光学部品、光コネクタ、光導波路、光ファイバなどの光通信用部品、アンテナ、レーザ装置、顕微鏡、半導体集積回路（ＩＣ）、圧電素子、構造物等々の設計のために、流体内や固体内の波動解析シミュレーション法として、例えば、光やＸ線などを含む電磁波、音波や、弾性波なども含め、波動の伝播の時間的な変化を解析するシミュレーション法として、ＦＤＴＤ法が利用されてきている。
【０００３】
このＦＤＴＤ法は、特に、マックスウェル（Maxwell）方程式に従う電磁波の伝播の時間的な変化の解析ツールとして利用されている。このＦＤＴＤ法を用いて、対象とする伝播系、すなわち、解析領域内において、マックスウェル方程式をコンピュータで数値的に解くことにより、この解析領域内において、例えば空気中を電磁場が伝播し、対象とする解析領域内に存在する物質表面で反射したり、屈折して透過したりする伝播状態を理論的に予測することができる。
このような電磁波の伝播状態の論理的な予測は、上述の光学部品、光通信用部品、アンテナを含む種々の部品等の設計に極めて有用な情報や必要な情報を与える。
【０００４】
このようなＦＤＴＤ法を用いた電磁場解析装置および方法が、特許文献１に提案されている。この特許文献１には、解析領域の全領域（全空間）の電磁場がゼロ「０」の状態を初期状態として、解析領域の一部分に、入射光を発生する励震部を設定し、励震部から連続波、パルス波などを発生させ、ＦＤＴＤ法によって、解析領域内の電磁場の時間変化を計算することが開示されている。また、特許文献１では、入射波に関する情報を入力し、入力された情報に基づいて入射波を計算することにより入射波の誤差を正確に再現し、再現された誤差を含んだ値を計算された入射波とし、この計算された入射波に基づいて電磁場を計算することによりＦＤＴＤ法による誤差を相殺し、ＦＤＴＤ法による誤差のない計算結果を出力することができることを開示している。
【０００５】
この特許文献１に開示されるようなＦＤＴＤ法を用いた従来の電磁場解析方法のフローチャートを図１２に示す。なお、この例では、従来の電磁場解析方法で用いられているＦＤＴＤ法計算プログラムの計算モデルが設定されている。
まず、ステップＳ２００において、この計算モデルの解析領域内の全空間に対して屈折率の空間分布を定義する。
また、ステップＳ２０２において、この計算モデルの空間計算範囲を解析領域の全空間（全空間離散化数：Ｓ０）に設定するなど、解析領域の空間離散化数（１〜Ｓ０）（空間離散化方法および条件）、境界条件などの解析領域の（状態や条件の）設定や、入射波（入力電磁波）の（状態や条件の）設定およびその他、全計算ステップ数（Ｎ０）電磁場計算に必要となる条件の設定をするために、計算条件を入力する。
ここで、屈折率の空間分布定義ステップＳ２００と計算条件の入力ステップＳ２０２とは、いずれを先に行っても良いし、同時に行っても良い。こうして、まず、電磁場計算を開始する前に、この計算モデルにおける電磁場計算の対象となる解析領域の屈折率の空間分布の定義および解析領域の全空間を空間計算範囲に設定することなどの必要な計算条件の設定がなされる。
【０００６】
その後、ステップＳ２０４およびＳ２０６において、解析領域の屈折率の空間分布の定義および必要な計算条件の設定に基づいて、ＦＤＴＤ法計算プログラムを用いて、電磁場が存在しない状態を初期状態とし、解析領域の空間上に定めた励震部からの電磁波を入射させ（ｊ＝１）、電磁波が解析空間の全空間（ｉ＝１〜Ｓ０）を伝播する時間（ｊ＝Ｎ０）まで、マックスウェル方程式をコンピュータで数値的に解く電磁場の計算を行い、電磁場の時間的、空間的変化を計算する。すなわち、各時間計算ステップ（ｊ）について、前の計算ステップ（ｊ−１）の解析領域の全空間（１〜Ｓ０）の電磁場の分布に基づいて、ＦＤＴＤ法計算プログラムを用い、固定化された空間計算範囲（１〜Ｓ０）である解析領域の全空間（１〜Ｓ０）について電磁場の時間的な変化を計算して、当該計算ステップ（ｊ）の解析領域の全空間（１〜Ｓ０）の電磁場の分布を求めることを、電磁波が解析空間の全空間（１〜Ｓ０）を伝播するまでの全時間ステップ（ｊ＝１〜Ｎ０）について繰り返し、解析領域の全空間（１〜Ｓ０）の電磁場の時間的、空間的変化を計算する。
【０００７】
具体的には、各時間計算ステップ（ｊ）について、ステップＳ２０４において、この時間ステップ（ｊ）の前の計算ステップ（ｊ−１）の解析領域の全空間の電磁場の分布に基づいて、ＦＤＴＤ法計算プログラムを用い、解析領域の全空間（１〜Ｓ０）について電磁場の時間的な変化を計算して、当該計算ステップ（ｊ）の解析領域の全空間（１〜Ｓ０）の電磁場の分布を求める。
なお、ｊ−１が０（ｊ−１＝０）の場合には、この前の計算ステップ（ｊ−１）における解析領域の全空間の電磁場の分布は、初期状態、例えば、０に設定する。
次に、ステップＳ２０６において、全時間計算ステップの計算が終了したか否かを判断する。図示例では、計算ステップが全時間計算ステップ数Ｎ０（ｊ＝Ｎ０）であるか、Ｎ０より大（ｊ＞Ｎ０）であるか否かを判断する。全計算ステップの計算が終了していれば（図示例では、ｊがＮ０以上（ｊ≧Ｎ０）であれば）、Ｙｅｓとして、ステップＳ２０８に抜け出し、計算後処理を行なう。一方、全計算ステップの計算が終了していなければ（図示例では、ｊがＮ０未満（ｊ＜Ｎ０）であれば）、Ｓ０：ｊ＝ｊ＋１として、ステップＳ２０４に戻る。
【０００８】
続いて、計算ステップ（ｊ＋１）について、ステップＳ２０４において、前の計算ステップ（ｊ）の解析領域の全空間の電磁場の分布に基づいて、同様に、ＦＤＴＤ法計算プログラムを用い、解析領域の全空間（１〜Ｓ０）について電磁場の時間的な変化を計算して、当該計算ステップ（ｊ＋１）の解析領域の全空間（１〜Ｓ０）の電磁場の分布を求め、次に、ステップＳ２０６において、全計算ステップの計算が終了した（ｊ＋１≧Ｎ０）か否かを判断し、終了していなければ（ｊ＋１＜Ｎ０）ステップＳ２０４に戻り、全計算ステップの計算が終了するまで、上記ステップＳ２０４〜Ｓ２０６を繰り返し、終了していれば（ｊ＋１≧Ｎ０）ステップＳ２０８に抜け出して、計算後処理を行なう。
【０００９】
続いて、ステップＳ２０８において、計算後処理として、こうして計算された電磁場の時間、空間分布に基づいて、解析目的に応じた物理量を計算する。また、計算結果を、表示デバイスのディスプレー上にグラフィック表示したり、プリンタにハードコピー出力する。
この後、ステップＳ２１０（ＥＮＤ）において、電磁場解析を終了する。
こうして、この計算モデルの解析領域の全空間を伝播した電磁波についての電磁場計算が終了する。
このような従来の電磁場解析技術では、１つの計算モデルについて、解析領域の全空間（１〜Ｓ０）おける電磁場の計算を、単純に各計算ステップ（ｊ）において、すなわち、全ての時間計算ステップ（ｊ＝１〜Ｎ０）において、逐次実施している。
【００１０】
また、ＦＤＴＤ法、伝送線路行列法、空間回路網法および有限要素時間領域法などの、マックスウエル方程式に従って有限の速度で伝搬していく電磁波などの波動の変化を、逐次的、かつ過渡的に計算する波動データの計算方法を用いた従来の波動解析方法および装置が、特許文献２に開示されている。
この特許文献２においては、解析領域を複数のブロックに分割し、各ブロック中の計算対象となる離散点の内の所定の判定離散点の波動データが所定値よりも大きいか否かを各ブロック毎に判定して、判定離散点の波動データが所定値よりも大きいブロックについてのみ、当該ブロックに含まれる全ての離散点の波動データの空間的および時間的変化を計算することにより、波動データの計算時間を短縮できることを開示している。
また、特許文献２では、計算対象ブロックの波動データの計算に必要なデータ、すなわち当該ブロックに含まれる全ての離散点の波動データおよびその近傍の所要の波動データを主記憶手段に記憶し、他の全てのブロックに含まれる離散点の波動データを補助記憶手段に記憶しておき、計算対象ブロックの計算が終了した時には、計算結果を補助記憶手段に記憶させ、次の計算対象となる他のブロックの波動データの計算に必要なデータを補助記憶手段から読み出し、主記憶手段に記憶させることにより、計算に必要なコンピュータの主記憶容量を少なくすることができることを開示している。
【００１１】
【特許文献１】特開２０００−２２７４５０号公報
【特許文献２】特開平９−１６５４８号公報
【発明の開示】
【発明が解決しようとする課題】
【００１２】
ところで、特許文献１に開示されたような技術では、上述したように、所定の解析領域の電磁波の伝播モデルについて、ＦＤＴＤ法によって電磁波の伝播の時間変化を時間計算ステップ毎に逐次計算する場合、解析領域の全空間について計算している。
このため、いずれの計算ステップにおいても、すなわち計算ステップに関わらず、解析領域の全空間の電磁場を計算するため、計算初期の電磁波が到達しない領域（空間）も計算を行うという問題がある。
また、特許文献２に開示されたような技術では、上述したように、計算対象の解析領域を複数のブロックに分割し、所定点の波動データが規定値よりも大きなブロックだけを計算対象とすることで、計算量を低減することができる。しかしながら、この場合、計算ステップの度に、各ブロック毎にその中の所定点の波動データと規定値との大小を調べる必要があり、その分の計算負荷を生じてしまうという問題がある。
【００１３】
本発明の課題は、上記従来技術の問題点を解消し、ＦＤＴＤ法などの波動データの計算方法を用いて、１つの波動解析モデルの解析領域を伝播する波動場の時間的な変化を計算する場合に、波動の伝播の初期段階、すなわち波動場の計算の初期段階における、入射波が全領域に行き渡らない間の計算量を低減し、計算を高速化することができる波動場解析方法および装置ならびに波動場解析方法をコンピュータに実行させるプログラムを記録したコンピュータに読み取り可能な記録媒体を提供することにある。
【課題を解決するための手段】
【００１４】
上記課題を解決するために、本発明の第１の態様は、波動が伝播する解析領域内の波動場の時間変化の計算を初期状態から各計算ステップ毎に連続して順次行い、前記解析領域内を前記波動が伝播する伝播状態を解析する波動場解析方法であって、事前に前記解析領域内における前記波動が到達している到達領域を予測し、所定の計算ステップにおける前記波動場の時間変化の計算を、当該計算ステップに前記解析領域内の前記波動が到達していると予測された前記到達領域を当該計算ステップの空間計算範囲として行うことを特徴とする波動場解析方法を提供するものである。
【００１５】
ここで、前記波動は、前記解析領域の一端から略前記解析領域の所定の軸方向へ伝播するものであるのが好ましい。
また、前記波動の前記到達領域の予測は、前記解析領域を伝播する前記波動の伝播速度を予測することにより行われるものであるのが好ましい。
また、前記波動の前記到達領域の予測は、前記解析領域内の屈折率分布の情報に基づいて行われるものであるのが好ましい。
また、前記解析領域内の屈折率分布は、前記到達領域内における前記波動の伝播方向と直交する面の屈折率の代表値を、この前記波動の伝播方向と直交する面内において最も小さい屈折率とすることにより得られるものであるのが好ましい。
また、前記波動の前記到達領域の予測は、前記解析領域内の屈折率を１であると仮定して行われるものであるのが好ましい。
また、前記波動の前記到達領域の予測は、前記解析領域の前記一端から前記波動が伝播する伝播長を計算することにより行われるものであるのが好ましい。
また、前記波動は、電磁波であり、前記波動場は、電磁場であるのが好ましい。
【００１６】
また、上記課題を解決するために、本発明の第２の態様は、波動が伝播する解析領域内の波動場の時間変化の計算を初期状態から各計算ステップ毎に連続して順次行い、前記解析領域内を前記波動が伝播する伝播状態を解析する波動場解析装置であって、前記波動が伝播する前記解析領域内の波動場の時間変化を初期状態から連続して計算する波動場計算手段と、事前に、前記解析領域内における前記波動が到達している到達領域を予測する予測手段と、前記波動場計算手段は、所定の計算ステップにおける前記波動場の時間変化の計算を、前記所定の計算ステップに前記解析領域内の前記波動が到達していると予測された前記到達領域を当該計算ステップの空間計算範囲として行うことを特徴とする波動場解析装置を提供するものである。
【００１７】
ここで、前記波動は、前記解析領域の一端から略前記解析領域の所定の軸方向へ伝播するものであり、前記波動場計算手段は、前記所定の軸方向に沿って前記波動場の時間変化を計算するものであるのが好ましい。
また、前記予測手段は、前記波動の前記到達領域の予測を、前記解析領域を伝播する前記波動の伝播速度を予測することにより行うのが好ましい。
また、前記予測手段は、前記波動の前記到達領域の予測を、前記解析領域内の屈折率分布の情報に基づいて行うのが好ましい。
また、前記予測手段は、前記到達領域内における前記波動の伝播方向と直交する面の屈折率の代表値を、この前記波動の伝播方向と直交する面内において最も小さい屈折率として、前記解析領域内の屈折率分布を得、得られた前記解析領域内の屈折率分布に基づいて前記波動の前記到達領域の予測を行うのが好ましい。
また、前記予測手段は、前記波動の前記到達領域の予測を、前記解析領域内の屈折率を１であると仮定して行うのが好ましい。
また、前記予測手段は、前記波動の前記到達領域の予測を、前記解析領域の前記一端から前記波動が伝播する伝播長を計算することにより行うのが好ましい。
また、前記波動は、電磁波であり、前記波動場は、電磁場であるのが好ましい。
【００１８】
また、上記課題を解決するために、本発明の第３の態様は、上記第１の態様の波動場解析方法をコンピュータに実行させるプログラムを記録したことを特徴とするコンピュータに読み取り可能な記録媒体を提供するものである。
【発明の効果】
【００１９】
本発明によれば、上記構成により、ＦＤＴＤ法などの波動データの計算方法を用いて、１つの波動解析モデルの解析領域を伝播する波動場の時間的な変化を各時間計算ステップ毎に計算する場合に、ある計算ステップにおいて、入射波が解析領域内の到達可能な到達領域を予め計算しておき、計算ステップ毎に波動場の計算対象となる到達領域である空間計算範囲を変化させ、波動場の反復計算を行なう空間計算範囲を、当該時間ステップにおいて入射波が到達可能な到達領域に限定することにより、解析領域内の、入射波が到達していない未到達領域の波動場の計算を省略できるので、計算量を低減し、計算を高速化することができる。
また、本発明によれば、入射波が到達可能な到達領域の予測や入射波が到達していない未到達領域の予測は、電磁場などの波動場の計算を行なう前に、簡易な計算で済むため、時間計算ステップ毎の計算領域を決める上での負荷は小さい。
その結果、本発明によれば解析領域内の、入射波が到達していない未到達領域の波動場の計算を省略による計算量の低減効果および計算の高速化効果を享受できる。
【発明を実施するための最良の形態】
【００２０】
以下に、本発明に係る波動場解析方法および装置ならびに波動場解析方法をコンピュータに実行させるプログラムを記録したコンピュータに読み取り可能な記録媒体を、添付の図面に示す好適実施形態に基づいて詳細に説明する。
【００２１】
本発明の波動場解析方法および装置は、マックスウェル方程式などの波動方程式に従う解析領域内の電磁波等の波動の伝播状態の１つの計算モデルを、ＦＤＴＤ法などの波動データの計算方法を用いて各時間計算ステップ毎に逐次計算する場合に、まず、最初に、ある計算ステップにおいて入射波が入射した初期状態から解析領域内を伝播する波動が到達可能な到達領域を予め計算して、全ての計算ステップにおける到達領域を求めておき、計算ステップ毎に、電磁場等の波動場の計算対象となる到達領域である空間計算範囲を変化させ、波動場の反復計算を行なう空間計算範囲を、当該時間ステップにおいて入射波が到達可能な到達領域に限定して、解析領域内の、入射波が到達していない未到達領域の波動場の計算を省略することにより、計算量を低減し、計算を高速化するものである。
【００２２】
以下では、波動の代表例として、光や電子線やＸ線などを含む電磁波を用い、波動場および波動方程式の代表例として、電磁場およびマックスウェル方程式を用い、ＦＤＴＤ法によってマックスウェル方程式をコンピュータで数値的に解き、対象とする解析領域の電磁波の空間および時間分布を計算し、電磁場の時間的な変化を求めて電磁場の解析を行う例を代表例として説明するが、本発明はこれに限定されるわけではない。本発明は、波動方程式をコンピュータで数値的に解き、対象とする解析領域の波動の空間および時間分布を計算し、波動場の時間的な変化を求めて波動場の解析を行うものであれば、例えば、音波や、弾性波などの電磁波以外の波動にも適用可能であることはもちろんである。また、本発明は、マックスウェル方程式などの波動方程式を解く数値計算方法として、ＦＤＴＤ法に限定されず、特許文献２に開示されている伝送線路行列法、空間回路網法および有限要素時間領域法などの他の波動データの数値計算方法を用いることができることはもちろんである。
【００２３】
すなわち、本発明は、電磁波の空間、時間分布を、マックスウェル方程式に基き、屈折率（誘電率）の空間分布を境界条件として、コンピュータで数値的に解くためのＦＤＴＤ法などの数値計算プログラムを用いたシミュレーション方法および手段を提供するもので、主に、光ディスク、レンズ、特に、マイクロレンズ等の種々の光学部品、光コネクタ、光導波路、光ファイバ等の光通信用部品、アンテナ、レーザ装置、顕微鏡、半導体集積回路（ＩＣ）、圧電素子、構造物等々の解析、微小領域の光学現象解析などに用いられ、開始状態の電磁場を０とし、外部から入射された電磁波の時間変化を計算する手法を提供するものである。
【００２４】
まず、本発明の一実施形態に係る波動場解析方法を実施する波動場解析装置が適用される電磁場解析装置の機能的構成について説明する。
図１は、本発明の波動場解析装置が適用される電磁場解析装置の機能的構成の一実施例を示すブロック図である。
図１に示す本発明の電磁場解析装置１０は、パソコンやワークステーションなどの計算機（コンピュータ）環境で動作するプログラムを含み、コンピュータによって構成することができるものである。
【００２５】
図１に示すように、電磁場解析装置１０は、複数の計算モデルに対する屈折率の空間分布を定義する屈折率分布定義部１２と、解析領域を含む電磁場計算に必要となる条件（全時間計算ステップ数および全空間離散化層数）を設定するために計算条件を入力する計算条件入力部１４と、屈折率分布定義部１２によって定義された計算モデルの屈折率分布および計算条件入力部１４によって入力された計算条件に基づいて、予め、解析領域内の入射電磁波が所定の計算ステップにおいて到達する到達領域を予測する到達領域予測部１６と、上述の計算モデルの屈折率分布および上述の計算条件に基づいて、ＦＤＴＤ法を用いてマックスウエル方程式を各計算ステップ毎に逐次数値的に解き、伝播する電磁場を計算する電磁場計算部１８と、到達領域予測部１６において予測された入射電磁波の予測到達領域と計算ステップの関係を記憶保存する記憶部２０と、到達領域予測部１６からの計算ステップ毎の予測到達領域の記憶部２０への書き込みおよび記憶部２０から電磁場計算部１８への所定の計算ステップにおける予測到達領域の読み出しを制御して行う書込／読出制御部２２と、電磁場計算部１８によって計算された解析領域の全空間を伝播した電磁場の計算結果を出力する計算結果出力部２４とを有する。
【００２６】
屈折率分布定義部１２は、解析対象となる解析領域の形状、材料に基づいて、計算モデルとなる屈折率の空間分布を定義するものである。
計算条件入力部１４は、空間離散化方法（例えば、電磁波が面状に解析領域の一端面から他端面まで解析領域の所定の軸方向へ伝播する場合、離散化層：１〜Ｓ０）、計算範囲（例えば、全空間離散化層数Ｓ０）、境界条件などの解析領域の入力や設定、空間分布、時間分布、偏光条件などの入力電磁波（入射波）の入力や設定、および、その他の、全時間計算ステップ数（Ｎ０）、時間離散化条件（計算ステップ：１〜Ｎ０）などの計算に必要となる条件の入力や設定を行うものである。
【００２７】
これらの屈折率分布定義部１２および計算条件入力部１４は、電磁場計算部１８によって解析領域内の電磁場計算を開始する前に、１つの計算モデルについて、電磁場計算の対象となる解析領域の屈折率の空間分布の定義および必要な計算条件の入力・設定を行う。
なお、屈折率分布定義部１２による屈折率の空間分布定義と、計算条件入力部１４による計算条件の入力および／または設定とは、互いに対応させながら同時に行っても良いし、いずれかを先に行っても良いし、また、独立して行っても良い。
【００２８】
到達領域予測部１６は、本発明の特徴とする部分であって、屈折率分布定義部１２によって定義された計算モデルの屈折率分布、およびこの屈折率分布および計算条件入力部１４によって設定された計算条件に基づいて、解析領域内に入射した電磁波が、解析領域内を伝播し、入射してからある時間後に到達する、すなわち、計算開始からある計算ステップ（ｊ）において到達する到達領域を予測するものである。
なお、到達領域予測部１６では、図２に実線で示すように、計算ステップ（ｊ）（ｊ＝１〜Ｎ０）と予測到達領域との関係を解析領域全体に亘って求めておく。
【００２９】
ここで、上述したように、例えば、電磁波が面状に解析領域の一端面から他端面まで解析領域の所定の軸方向へ伝播するとし、解析領域が、離散化層が一端面から他端面まで、すなわち１層〜Ｓ０層（全空間離散化層数Ｓ０）まで、空間的に離散化され、全時間計算ステップ数（Ｎ０）の方が、全空間離散化層数（Ｓ０）よりも大きい場合には、到達領域予測部１６による予測到達領域が、複数の計算ステップにおいて、同一の離散化層に該当する場合がある。このため、図２に点線で示すように、予測到達領域が、離散化層（ｉ）（ｉ＝１〜Ｓ０）に到達する時間に該当する計算ステップ（ｊ）をｔ（ｉ）として、予測到達領域（離散化層：ｉ＝１〜Ｓ０）と計算ステップｔ（ｉ）との関係を全離散化層に亘って求めておくのが好ましい。ここで、予測到達領域（離散化層：ｉ）と計算ステップｔ（ｉ）とは一対一に対応する。
【００３０】
なお、到達領域予測部１６による所定の計算ステップ（ｊ）における解析領域内の到達領域の予測は、電磁場計算部１８による解析領域内の電磁場計算の開始前に行われる。
こうして、電磁場計算の開始前に求められた予測到達領域（例えば、離散化層：ｉ＝１〜Ｓ０）と計算ステップ（ｔ（ｉ））との関係は、ＬＵＴ（ルックアップテーブル）や、関係式として、書込／読出制御部２２によって記憶部２０に記憶保存される。
上述した例では、到達領域予測部１６は、予測到達領域と計算ステップとの関係を解析領域の全空間に亘って求め、この関係を記憶部２０に記憶保存し、後述するように、電磁場計算部１８が、この関係式を読み出して、電磁場計算を行うべき空間計算範囲を設定しているが、本発明はこれに限定されず、到達領域予測部１６において、計算ステップ（ｊ）に対し、計算ステップ（ｊ）において計算すべき空間計算範囲を設定しておき、計算ステップ（ｊ）と空間計算範囲との関係をＬＵＴなどとして記憶部２０に記憶しておき、電磁場計算部１８が、記憶部２０から、直接、計算ステップ（ｊ）において計算すべき空間計算範囲を読み出すようにしても良い。
【００３１】
なお、本発明においては、到達領域予測部１６は、到達領域の予測を、解析領域を伝播する電磁波の伝播速度を予測することにより行うのが好ましい。なお、電磁波の伝播速度（ｖ）は、電磁波が伝播する媒質の屈折率ｎに反比例する（ｖ＝ｃ０／ｎ：ここで、ｃ０は、真空中の光速（ｃ０＝２．９９８×１０８ｍ／ｓ）である）ので、到達領域予測部１６は、到達領域の予測を、解析領域内の屈折率分布の情報に基づいて行うのが良い。また、屈折率ｎは、誘電率εおよび透磁率μの関数である（ｎ＝√（εμ／ε０μ０）：ε０およびμ０は、それぞれ真空中の誘電率（10７／4πc０２）および透磁率（4π×10−７）である）ので、到達領域予測部１６は、到達領域の予測を、解析領域内の媒質の誘電率εおよび透磁率μの情報に基づいて行っても良い。
【００３２】
具体的には、到達領域予測部１６は、解析領域内の到達領域における電磁波の伝播方向と直交する面を代表する屈折率の代表値として、この電磁波の伝播方向と直交する面内において最も小さい屈折率を採用し、この解析領域内の電磁波の伝播方向における屈折率分布を得、得られた屈折率分布に基づいて電磁波の到達領域の予測を行うのが好ましい。
その理由は、上述したように、電磁波の伝播速度（ｖ）は屈折率ｎに反比例するので、電磁波の伝播方向と直交する面内における屈折率の内の最も小さい屈折率を採用することにより、予測される伝播速度は、実際の伝播速度より少し速くなり、予測到達領域は、実際の到達領域より少し大きい、すなわち、少しより遠くまで到達した領域となり、効率よく、かつ、確実に、実際の到達領域を含む領域を予測できるからである。
【００３３】
また、到達領域予測部１６は、このような解析領域内の電磁波の伝播方向における屈折率分布の代わりに、解析領域内の屈折率を１であるとして、電磁波の到達領域の予測を行うこともできる。その理由は、上述したように、定義により真空の屈折率は1であり、1より小さい屈折率を持つ物質、すなわち媒質は存在しないし、また、空気の屈折率は１に近いので、解析領域内に真空や空気を含む場合には、特に、より簡単に、かつ、確実に、実際の到達領域を含む領域を予測できるからである。
なお、本発明においては、電磁波は、解析領域の一端から、ほぼ、解析領域の所定の軸方向へ伝播するものであるとするのが好ましい。もちろん、この場合には、電磁場計算部１８は、この所定の軸方向に沿って電磁場の時間変化を計算する。
また、到達領域予測部１６は、電磁波の到達領域の予測を、解析領域の一端から電磁波が伝播する伝播長を計算することにより行うのが好ましい。
なお、到達領域の予測の詳細および具体例については、後述する。
【００３４】
電磁場計算部１８は、屈折率分布定義部１２によって定義された計算モデルの屈折率分布および計算条件入力部１４によって設定された計算条件に基づいて、ＦＤＴＤ法計算プログラムによって、マックスウエル方程式をパソコンやワークステーションなどのコンピュータで数値的に解くものであり、解析領域（解析空間）内の電磁場の時間的、空間的変化を計算する。
このとき、本発明においては、電磁場計算部１８は、所定の計算ステップの電磁場計算を開始する前に、記憶部２０に記憶されている予測到達領域と計算ステップとの関係から当該計算ステップにおいて計算すべき解析領域内の空間計算範囲を設定し、設定された空間計算範囲内において当該計算ステップの電磁場計算を実施する。または、電磁場計算部１８は、所定の計算ステップの電磁場計算を開始する前に、記憶部２０に記憶されている計算ステップと空間計算範囲との関係から、直接、当該計算ステップで計算すべき空間計算範囲を読み出し、読み出された空間計算範囲内において当該計算ステップの電磁場計算を実施しても良い。
【００３５】
すなわち、本発明においては、電磁場計算部１８は、計算ステップ（ｊ＝Ｎ）における電磁場の計算を開始する前に、記憶部２０に記憶されている予測到達領域（例えば、離散化層：ｉ＝１〜Ｓ０）と計算ステップ（ｔ（ｉ））との関係を書込／読出制御部２２を介して読み出し、読み出した関係を用いて、計算ステップ（Ｎ）における予測到達領域（離散化層（ｉ））を求める。具体的には、電磁場計算部１８は、図２に示すような関係式またはＬＵＴなどを読み出し、計算ステップ（Ｎ）より大きい計算ステップｔ（ｉ）の中の最小の計算ステップｔ（Ｓ）、すなわち最小の離散化層Ｓを求め、離散化層１〜Ｓを計算ステップ（Ｎ）における空間計算範囲Ｓとして設定する。
なお、本発明において、電磁波の予測到達領域のみを空間計算範囲Ｓとして設定することができるのは、電磁場計算開始直後は、入射波が入射した直後であり、入射波が到達しない領域の電磁場は初期状態と同じ０のままであり、計算の必要がないため、省略できるからである。
この後、電磁場計算部１８は、計算ステップ（Ｎ）における電磁場の計算を、離散化層１〜Ｓの空間計算範囲に制限して行う。その結果、本発明においては、計算ステップ（Ｎ）において、残りの離散化層Ｓ＋１〜Ｓ０についての電磁場の計算を省略することができるので、電磁場計算部１８による電磁場の計算量を低減し、この計算を高速化することができる。
【００３６】
なお、電磁場計算部１８が実施するＦＤＴＤ法計算プログラムは、電磁場の物理現象を表すマックスウエル方程式をコンピュータで数値的に解くための１つの手法であり、電磁場が存在しない状態を初期状態とし、解析空間上に定めた励震部（波源）からの電磁波を入射させ、電磁場の時間変化を逐次更新しながら計算するシミュレーション手法である。（なお、ＦＤＴＤ法やＦＤＴＤ法計算プログラムについては、「ＦＤＴＤ法による電磁界およびアンテナ解析」宇野  亨、コロナ社発行の文献を参照することができる。）
【００３７】
ＦＤＴＤ法やＦＤＴＤ法計算プログラムについては、上記文献を始めとして、種々の従来公知の文献に開示された技術を用いれば良いので、詳細な説明は省略するが、以下に、本発明において適用されるＦＤＴＤ法の一例について簡単に説明する。
ＦＤＴＤ法は、下記式（１）および（２）に示すマックスウェル方程式を、時間、空間方向に差分化した電磁場の数値計算手法である。
電磁場計算では、下記式（１）および（２）に示すマックスウェル方程式を、時間および空間を離散化して差分方程式、すなわち、下記式（３）および（４）に示す電磁場計算差分式(一例として電磁場の１成分のみが示されている)に変換する。
【００３８】
【数１】
2008250956.tif 000003
【００３９】
【数２】
2008250956.tif 000004
【００４０】
本発明に適用されるＦＤＴＤ法においては、図３に示すように、入射光（電場Ｅ、磁場Ｈの電磁波）が入射する解析空間（領域）をメッシュ状に分割して離散化し、解析空間の形状や材料に応じて、各格子（微小セル）点に、屈折率ｎ（誘電率ε）を割り当てる。図３において、解析空間から１つの格子（単位セル）を抜き出して示すように、Ｙｅｅセルと呼ばれる１つの格子上の点に、電場および磁場のｘ，ｙ，ｚ成分、すなわち、Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚを定義する。
【００４１】
ここで、上記式（３）から明らかなように、磁場（磁界）の変化（差分）は、周りの電場（電界）の和および差の定数倍、すなわち、図４に示されるように、周りの電場が形成する渦の定数倍で表され、上記式（４）から明らかなように、電場の変化（差分）は、周りの磁場の和および差の定数倍、すなわち、図４に示されるように、周りの磁場が形成する渦の定数倍で表される。
このため、ＦＤＴＤ法では、各Ｙｅｅセルにおいて、ある時間の電場から磁場を上記式（３）で計算し、これで求まった磁場を上記式（４）に代入して、次の時間の電場を計算し、さらに、これで求まった電場を再び上記式（３）に代入して、磁場を計算する、あるいは、ある時間の磁場から電場を上記式（４）で計算し、これで求まった電場を上記式（３）に代入して次の時間の磁場を計算し、さらに、これで求まった磁場を再び上記式（４）に代入して、電場を計算する。
【００４２】
ＦＤＴＤ法においては、一般的には、解析領域の全空間、本発明では、計算ステップ（ｊ）における電磁波の予測到達領域に該当する解析領域の空間計算範囲の全Ｙｅｅセルについて計算を繰り返し行うことにより、電磁場の時間変化を計算することができる。すなわち、まず、ある時間の、解析空間における電場の空間分布から、半ステップ後の磁場の空間分布を計算する。引続き、この解析空間における磁場分布から、半ステップ後の電場を計算することで、1ステップ後の電場に更新される。このように、電場、磁場を交互に更新しながら、電磁場の時間変化を計算する。
したがって、ＦＤＴＤ法では、このように、電場と磁場とを交互に計算する計算手順（リープフロッグ（蛙飛び）アルゴリズム）を繰り返すことにより、電磁場の計算を時系列に沿って行うことができる。
【００４３】
以上から、本発明では、ＦＤＴＤ法において、解析領域内に励震源（波源）を定め、入射波を与え、解析領域内の全空間の電磁場を０とした状態を初期（開始）状態として、上述した式（３）および（４）による計算を繰り返し行うことにより、入射波が、計算モデル内を伝搬する様子を計算することができる。
また、本発明では、入射波として、連続波を与え、変化が収束するまで計算を続ける、もしくは、パルス波を与えることもできる。また、本発明において、境界としては、吸収境界条件、周期境界条件などを設定することができる。
本発明においては、上述したように、ＦＤＴＤ法やＦＤＴＤ法計算プログラムを適用することができる。
【００４４】
記憶部２０は、到達領域予測部１６によって計算ステップ（ｊ）（ｊ＝１〜Ｎ０）または（ｔ（ｉ））と予測到達領域（例えば、離散化層：ｉ＝１〜Ｓ０）との関係などを、解析領域全体に亘る関係式や一次元ＬＵＴなどとして記憶保存するものである。あるいは、記憶部２０は、計算ステップ（ｊ）と空間計算範囲との関係などを、解析領域全体に亘る関係式や一次元ＬＵＴなどとして記憶保存しても良い。
書込／読出制御部２２は、このような、計算ステップと予測到達領域または空間計算範囲などとの関係などを記憶部２０に記憶保存させるために、そのような関係を示す関係式やＬＵＴとして記憶部２０に書き込むと共に、その書き込み制御を行い、かつ、電磁場計算部１８によって記憶部２０に記憶保存された予測された計算ステップと予測到達領域または空間計算範囲との関係を記憶部２０から読み出すと共に、その読み出し制御を行うものである。
【００４５】
計算結果出力部２４は、電磁場計算部１８による１つの計算モデルの電磁場計算終了後の後処理として、電磁場計算部１８によって行われた解析領域内の磁場の内時間および空間分布の計算結果に基づいて、反射率、透過率、吸収率および遠方放射分布などの解析目的に応じた物理量を計算する。また、計算結果出力部２４は、こうして計算された物理量や、電磁場の時間および空間分布などの計算結果を、表示デバイスのディスプレー上にグラフィック表示する機能などを備えていても良いし、プリンタでハードコピーを出力する機能などを備えていても良い。
【００４６】
次に、図１に示す電磁場解析装置のハードウエア構成について説明する。
図５は、図１に示す電磁場解析装置のハードウエア構成の一実施例を示すブロック図である。
図５に示す電磁場解析装置３０は、図１に示す電磁場解析装置１０の機能的構成をハードウエア構成として持つもので、上述したように、図１に示す本発明の電磁場解析装置１０を、パソコンやワークステーションなどのコンピュータなどで構成した場合の一実施例である。
【００４７】
図５に示す電磁場解析装置３０は、解析領域や、屈折率や、計算条件などのデータの入力を行う入力装置３２と、本発明の電磁場解析方法をコンピュータなどで実行するプログラム、ＦＤＴＤ法計算プログラムなどのプログラムや、解析領域、屈折率、計算条件などのデータを記憶した着脱可能なデータ記録媒体からデータを読み取り、またデータ記録媒体にデータを書き込む読取／書込装置３４と、本発明の特徴とする解析領域内を伝播する電磁波の到達領域を予測し、時間または時間計算ステップと予測到達領域または空間計算範囲との関係を求め、または計算ステップと予測到達領域との関係から空間計算範囲を求めるなどの演算、ＦＤＴＤ法計算プログラムを実行して電磁場の時間的および空間的な変化を計算するなどの演算等の処理を行う、また、時間または時間計算ステップと予測到達領域または空間計算範囲との関係の書き込み／読み出しの制御を初めとして、全体を制御する演算処理装置（例えば、ＣＰＵ）３６と、本発明の特徴とする時間または時間計算ステップと予測到達領域または空間計算範囲との関係、ＦＤＴＤ法計算プログラムや、本発明の電磁場解析方法の実行プログラムなどのプログラムや、必要な計算条件などのデータを記憶するメモリ３８と、計算結果等を表示する表示装置４０と、計算結果等を記録する記録装置（プリンタ）４２と、これらの各装置を接続するためのバス４４とを有する。
【００４８】
ここで、入力装置３２は、コンピュータ等のキーボードやマウスなどを含み、また、読取／書込装置３４は、各種のデータ記録媒体（メディア）の読み書きをそれぞれ行うメディアドライブなどを含む。また、入力装置３２および読取／書込装置３４は、図１に示す電磁場解析装置１０の複数の計算モデルに対する屈折率の空間分布を定義する屈折率分布定義部１２および解析領域を含む電磁場計算に必要となる条件を設定するために計算条件（例えば、全空間離散化層数（Ｓ０）や、全時間計算ステップ数（Ｎ０））を入力する計算条件入力部１４としての機能を含む。
また、演算処理装置３６は、例えば、コンピュータ等のＣＰＵ（中央処理ユニット）などで構成され、時間または時間計算ステップと解析領域内の電磁波の予測到達領域または空間計算範囲との関係を予測して求める到達領域予測部１６およびＦＤＴＤ法プログラムを実行してマックスウエル方程式を数値的に解く電磁場計算部１８、ならびに上記関係の記憶部２０への書き込および記憶部２０からの上記関係の読み出しを制御して行う書込／読出制御部２２としての機能を含む。なお、演算処理装置３６は、電磁場計算部１８による電磁場計算後の目的に応じた物理量の算出などの計算結果出力部２４の機能を備えていても良い。
【００４９】
メモリ３８は、コンピュータ等のハードディスクや、各種の内蔵または外付けのメモリを含み、時間または時間計算ステップと解析領域内の電磁波の予測到達領域または空間計算範囲との関係や、計算ステップの前のステップの電磁場のデータ（時間、空間分布）などを記憶保存する記憶部２０としての機能を含む。
また、表示装置４０は、計算結果などを表示画面に表示するためのもので、コンピュータなどの液晶表示デバイスやＣＲＴなどのディスプレーやモニタを含み、また、記録装置４２は、計算結果などを紙などに記録したハードコピーを出力するためのもので、いわゆるプリンタを含む。表示装置４０および記録装置（プリンタ）４２は、解析領域の一端面から初期状態で開始され、他端面で終了した全空間を伝播した電磁場の計算結果を出力する計算結果出力部２４としての機能を持つ。
【００５０】
本発明の波動場解析装置が適用される電磁場解析装置は、基本的に以上のような機能的構成およびハードウエア構成を有するものであるが、以下に、上述した電磁場解析装置の作用および電磁場解析の実施例を説明することにより、本発明の波動場解析方法が適用される電磁場解析方法を説明する。
図６は、本発明の波動場解析方法が適用される電磁場解析方法の一例を示すフローチャートである。なお、図６に示す電磁場解析方法は、図１に示す電磁場解析装置によって実施されるものである。
【実施例】
【００５１】
以下では、図７に示すような、ガラス基板５２上に配置されるシリンドリカルレンズアレー５４が空気５６中に存在する場合に、入射光が、空気層５６に入射し、空気層５６、シリンドリカルレンズアレー５４およびガラス基板５２を順次伝播する計算モデル５０を実施例とする。なお、この計算モデル５０においては、空気層５６の厚みは、屈折率ｎは１（ｎ＝１）であり、シリンドリカルレンズアレー５４の各シリンドリカルレンズ５４ａの頂点において１μｍであり、その基部（隣接するシリンドリカルレンズ５４ａの接続部）において２μｍである。すなわち、この計算モデル５０においては、各シリンドリカルレンズ５４ａは、屈折率ｎが１．５（ｎ＝１．５）ある、半径１μｍの半円筒形状の部材であり、屈折率ｎが１．５（ｎ＝１．５）である、厚さ３μｍのガラス基板上に接続されている。なお、この計算モデル５０では、紙面に垂直な方向には、同じ構造が無限に連続すると仮定する。
【００５２】
本実施例では、本計算モデル５０の解析領域５０ａとして、図８に示すように、シリンドリカルレンズアレー５４を構成する１つのシリンドリカルレンズ５４ａを含み、その上部に空気層５６ａが存在し、シリンドリカルレンズ５４ａの下部にガラス基板５２ａが存在する２μｍ×５μｍの空間を設定する。この解析領域５０ａにおいては、紙面に垂直な方向には、同じ構造が無限に連続すると仮定した周期境界条件を設定し、上限方向の境界には吸収境界条件を設定する。
なお、以下では、図７および８に示すような、ガラス基板５２に支持され、空気中に存在するシリンドリカルレンズアレー５４に空気層５６から入射光が入射して伝播する計算モデル５０および解析領域５０ａを代表例として説明を行なうが、本発明はこれに限定されないことはいうまでもないことである。
【００５３】
まず、図６に示す電磁場解析方法においては、ＦＤＴＤ法計算プログラムを用いて電磁場計算を行うものであるが、そのステップＳ１００において、図１に示す装置１０の屈折率分布定義部１２によって、２μｍ×５μｍの解析領域（空間）５０ａに対して、その各領域が異なる屈折率を持つ屈折率の空間分布が定義される。
例えば、図８に示すように、２μｍ×５μｍの解析領域（空間）５０ａでは、上部の最小厚み１μｍの空気層５６ａでは、屈折率ｎが１．０となり、その下の半径１μｍの半円筒形状のシリンドリカルレンズ５４ａで、屈折率ｎが１．５となり、シリンドリカルレンズ５４ａを支持する下部の厚さ３μｍのガラス基板では、屈折率ｎが１．５となる屈折率分布が定義される。
【００５４】
また、ステップＳ１０２において、計算条件入力部１４によって、解析領域５０ａの計算範囲として２μｍ×５μｍの解析領域（空間）５０ａ、空間離散化方法として、空間グリッドΔｄが２０ｎｍに設定され、したがって、Ｙｅｅセルとして２０ｎｍ×２０ｎｍ×２０ｎｍの単位セル（格子）が設定される。
また、入射光（電磁波）条件として、解析領域５０ａの一端面（図７および図８中、計算モデル５０および解析領域５０ａの空気層５６ａの上側端面）を励震部５８とし、励震部５８から波長５００ｎｍの平面波が連続波として入射され、解析領域５０ａの空気層５６ａの上側端面からガラス基板５２ａの下側端面まで伝播することが設定される。
したがって、電磁波が伝播する解析領域５０ａの上側端面から下側端面に向かう進行方向と直交する方向においては、電磁波は一様であるので、解析領域５０ａでは、電磁波の到達を予測する予測到達領域として、進行方向において２０ｎｍの空間グリッドΔｄを持つ空間離散化層を２５０層（全空間離散化層数Ｓ０＝２５０）に設定することができる。
また、その他の計算条件として、紙面に垂直な方向には、上述した周期境界条件を設定し、上限方向の境界には上述した吸収境界条件が設定される。
【００５５】
また、ステップＳ１０２において、初期状態として、解析領域５０ａの全空間の電磁場が０に設定され、本発明の電磁場計算部１８によるの計算開始と同時に励震部５８による励震を開始するように設定される。
さらに、本発明の電磁場計算において、ＦＤＴＤ法計算プログラムで電磁場を更新する１計算ステップ（ｊ）当たりの時間をΔｔとする。Δｔは、クーラン条件と呼ばれる上限値以下に設定される定数であり、本実施例では、例えば、Δｔ＝０．０２９７８［ｆｓ／ｓｔｅｐ］に設定され、全計算ステップ数Ｎ０は、１０００［ｓｔｅｐ］に設定される。
こうして、ステップＳ１０２においては、必要な計算条件が入力される。
ここで、ステップＳ１００とステップＳ１０２とは、いずれを先に行っても良いし、同時に行っても良い。これらのステップＳ１００およびＳ１０２によって、電磁場計算を開始する前に、まず、各計算モデルについて、電磁場計算の対象となる解析領域の屈折率の空間分布の定義および必要な計算条件の設定がなされる。
【００５６】
次に、電磁場計算部１８によって電磁場計算を開始する前に、ステップＳ１０４において、到達領域予測部部１６によって時間計算ステップ（ｊ）毎の電磁波の到達領域を予測する。
具体的には、到達領域予測部部１６では、まず、図８に示す解析領域５０ａにおいて、入射波の進行方向を横断する格子の面（ｉ）毎に、屈折率の最小値 ｎ（ｉ）を求める。その結果、励震部５８から１μｍまでの１μｍの厚みの領域、すなわちｉ＝１〜５０［ｇｒｉｄ］（空間離散化層として、ｉ＝０〜４９層）には、空気層５６ａのみが存在するので、屈折率の最小値 ｎ（ｉ）は１．０となる。次に、励震部５８から１μｍ〜２μｍまでの１μｍの厚みの領域（ｉ＝５０〜９９層）には、屈折率ｎ＝１．５のシリンドリカルレンズ５４ａと屈折率ｎ＝１．０の空気層５６ａとが混在するので、この領域の屈折率の最小値 ｎ（ｉ）も１．０となる。次に、励震部５８から２μｍ〜５μｍまでの３μｍの厚みの領域（ｉ＝１００〜２５０層）には、屈折率ｎ＝１．５のガラス基板のみが存在するので、この領域の屈折率の最小値 ｎ（ｉ）は１．５となる。
すなわち、解析領域５０ａにおいて、ｉ＝０〜９９では、ｎ（ｉ）＝１．０であり、ｉ＝１００〜２５０では、ｎ（ｉ）＝１．５となる。
【００５７】
次に、各面（各空間離散化層）（ｉ）までの最速到達計算ステップｔ（ｉ）を下記式（５）に従って計算する。
【数３】
2008250956.tif 000005
ここで、ｎ（ｋ）およびｄ（ｋ）は、上述の空間離散化層（ｋ）の屈折率の最小値および空間離散化層（ｋ）のメッシュの距離（グリッド間隔）（ｎ（０）〜ｎ（９９）＝１．０、ｄ（０）〜ｄ（２４９）＝２０ｎｍ）であり、ｃ０は、真空中の光速（ｃ０＝２．９９８×１０８ｍ／ｓ）であり、Δｔ は、上述したように、ＦＤＴＤ法計算プログラムで電磁場を更新する１計算ステップ（ｊ）当たりの時間であり、例えば、Δｔ＝０．０２９７８［ｆｓ／ｓｔｅｐ］である。
【００５８】
こうして得られた空間離散化層数（ｉ）と、時間計算ステップｔ（ｉ）とをプロットすると、概略、図２に点線で示す関係が得られる。ここで、Ｓ０＝２５０、Ｎ０＝１０００である。
ステップＳ１０６において、このようにステップＳ１０４で得られた図２に点線で示す時間計算ステップと空間離散化層数との関係は、ＬＵＴや関係式などとして、書込／読出制御部２２によって記憶部２０に記憶保存される。
【００５９】
次に、ステップＳ１０８において、ｊ＝１として、電磁場計算部１８は、計算ステップ（ｊ＝１）における電磁場の計算を開始する前に、まず、記憶部２０に記憶保存されている図２に点線で示す計算ステップと空間離散化層数との関係を読み出し、計算ステップ（ｊ＝１）における空間計算範囲Ｓを設定する。
なお、本発明において、電磁波の予測到達領域のみを空間計算範囲Ｓとして設定することができるのは、電磁場計算の初期段階では、入射波は十分に伝播しておらず、入射波が到達しない領域の電磁場は、初期状態と同じ０のままであり、計算の必要がないため、省略可能だからである。その結果、本発明においては、以下の電磁場計算部１８による電磁場の計算量を低減し、この計算を高速化することができる。
【００６０】
なお、ここでは、一般的に、計算ステップ（ｊ＝Ｎ）における空間計算範囲Ｓの設定について具体的に説明する。
すなわち、本発明においては、電磁場計算部１８は、計算ステップ（ｊ＝Ｎ）における電磁場の計算を開始する前に、図２に示す空間離散化層数（予測到達領域：ｉ＝１〜Ｓ０）と計算ステップ（ｔ（ｉ））との関係を記憶部２０から読み出す。
図２において、Ｎｔ＜（ｉ）を満たす最小の計算ステップｔ（Ｓ）求め、離散化層１〜Ｓを計算ステップ（Ｎ）における空間計算範囲（Ｓ）として設定する。
こうして設定された解析領域５０ａ内の空間計算範囲（Ｓ）６０を図９に示す。
なお、空間計算範囲（Ｓ）６０の外側のｉ＝Ｓ＋１層目からＳ０層目までは、光速の限界で、入射波が到達し得ないため、全ての電磁場が０であり、計算の必要はない。
なお、Ｎ≧ｔ（Ｓ０）となった場合には、解析領域５０ａ内の全領域を空間計算範囲（Ｓ）として設定し、解析領域５０ａ内の全領域を計算する。
【００６１】
次に、ステップＳ１１０において、ｉ＝１として、電磁場計算部１８によって、図９に示すステップＳ１０８で設定された空間計算範囲（Ｓ）６０の１層目（ｉ＝１）の電磁場を計算する。
なお、ここでは、一般的に、設定空間計算範囲（Ｓ）６０のｉ層目（ｉ）の電磁場の計算について具体的に説明する。
すなわち、電磁場計算部１８によって、図９に示すステップＳ１０８で設定された空間計算範囲（Ｓ）６０について、まず、解析領域５０ａのｉ層目（ｉ）の計算ステップ（Ｎ−１）の電磁場に基づいて計算ステップ（Ｎ）の電磁場の変化を計算し、計算ステップ（Ｎ）における電磁場を計算する。
【００６２】
電磁場計算部１８によって、図８および図９に示す解析領域５０ａについて、屈折率の空間分布の定義および必要な計算条件に基づいて、ＦＤＴＤ法計算プログラムを用いて、電磁場が存在しない状態を初期状態とし、解析領域の上端面に定めた励震部５４からの入射光を入射させ、マックスウェル方程式をコンピュータで数値的に解き、設定空間計算範囲（Ｓ）の電磁場の計算を行い、設定空間計算範囲（Ｓ）の電磁場の時間的、空間的変化を計算する。
【００６３】
なお、電磁場計算部１８による電磁場の計算では、解析領域５０ａの全空間の電磁場が０の状態を開始状態として計算を進める。電磁場計算部１８では、励震部５８で与えた入射光（電磁波）が、徐々に解析空間５０ａ内を伝播して行くので、シリンドリカルレンズ５４ａに当たっても、シリンドリカルレンズ５４ａによる収束を計算できる。電磁場計算部１８では、設定空間計算範囲（Ｓ）や全時間計算ステップ（Ｎ０）などのように、一定の計算ステップ、もしくは、計算の収束状態の判定条件など終了条件を満たした段階で計算を終了させる。
【００６４】
次に、ステップＳ１１２において、全空間計算範囲（Ｓ）の計算が終了したか否か、すなわち、ｉ≧Ｓを満足するか否かを判断する。全空間計算範囲（Ｓ）の計算が終了していなければ（図示例では、ｉ＜Ｓであれば）、Ｎｏ：ｉ＝ｉ＋１として、ステップＳ１１０に戻り、同様にして、空間計算範囲（Ｓ）６０の（ｉ＋１）層目の電磁場を計算することを、全空間計算範囲（Ｓ）の計算が終了するまで繰り返す。
一方、全空間計算範囲（Ｓ）の計算が終了していれば（図示例では、ｉ≧Ｓを満足すれば）、ＹｅｓとしてステップＳ１１４に抜け出す。
【００６５】
次に、ステップＳ１１４において、全計算ステップ（Ｎ０）の計算が終了したか否か、すなわち、ｊ≧Ｎ０を満足するか否かを判断する。全計算ステップ（Ｎ０）の計算が終了していなければ（図示例では、ｊ＜Ｎ０であれば）、Ｎｏ：ｊ＝ｊ＋１として、ステップＳ１０８に戻り、同様にして、計算ステップ（ｊ＋１）における空間計算範囲（Ｓ）を新たに設定し、ステップＳ１１０およびＳ１１２において、空間計算範囲Ｓについて、ｉ＝１層目〜Ｓ層目まで、電磁場計算部１８によって電磁場を計算することを全空間計算範囲（Ｓ）の計算が終了するまで繰り返し、全空間計算範囲（Ｓ）の計算が終了した時点で、再びステップＳ１１４に抜け出す。こうして、全計算ステップ（Ｎ０）の各空間計算範囲における電磁場の計算が終了するまで、ステップＳ１０８、ステップＳ１１０、Ｓ１１２およびＳ１１４を繰り返す。
【００６６】
なお、ある時間計算ステップ（Ｎ）以降では、Ｎ≧ｔ（Ｓ０）となるので、解析領域５０ａ内の全領域を空間計算範囲（Ｓ）として設定し、電磁場計算部１８によって解析領域５０ａ内の全領域が計算される。
一方、全計算ステップ（Ｎ０）の計算が終了していれば（図示例では、ｊ≧Ｎ０を満足すれば）、ＹｅｓとしてステップＳ１１６に抜け出す。
【００６７】
次に、ステップＳ１１６において、計算後処理として、こうして計算された電磁場の時間、空間分布に基づいて、解析目的に応じた物理量を計算する。また、計算結果出力部２４において、計算結果を、表示デバイスのディスプレー上にグラフィック表示したり、プリンタにハードコピー出力する。
次に、ステップＳ１１８（ＥＮＤ）において、１つ計算モデル５０の解析領域５０ａについての電磁場計算が終了する。
【００６８】
このようにして、本実施例において、１つの計算モデル５０の解析領域５０ａについて行われた計算結果を、伝播ステップ数が５０（Ｎ＝５０ｓｔｅｐ）毎に、図１０（ａ），（ｂ）および（ｃ）に示す。なお、図１０（ｄ）には、計算モデル５０の解析領域５０ａのサイズを具体的に示す。
また、本実施例において求めた時間計算ステップ数（ｔ（ｉ），Ｎ）［ｓｔｅｐ］と、空間計算範囲（計算領域）［ｇｒｉｄ］との関係を示すグラフを図１１に示す。この図１１に示すグラフは、図２に示すグラフと同様であるが、図１１においては、グラフの上側の面積が計算量の削減割合を示している点で異なっている。
【００６９】
これらの図から、初期段階では（０計算ステップから約７３０ステップ目までの間では）、解析領域５０ａ内に、入射光が到達していない領域が存在し、その電磁場は０であるので、それらの計算ステップにおいては、入射光が到達していない領域の電磁場の計算を省略できるので、図１０（ａ）０計算ステップから７００ステップ目までの実線で囲まれた計算省略部分の電磁場の計算を省略できる。
その結果、図１１に示す例においては、削減された計算量を示すグラフの上側の面積が全体の計算量を示す全面積に対する割合が約３０％を示しており、約３０％の計算量削減率を達成することができる。
【００７０】
上述した計算例によれば、屈折率が低く、入射光の伝搬速度の速い領域を通過した光に対して、入射光の到達ステップを計算しているため、計算初期の各ステップにおいて、必要最小限の計算領域を設定することができ、未到達領域の電磁場計算を省略することで、計算量を低減し、高速化することができる。
なお、上述した実施例では、解析領域内の屈折率分布から入射波の進行方向と直交する面内の最小屈折率を求めて、入射波の到達領域の予測計算や、未到達領域の予測計算を行っているが、本発明はこれに限定されず、単純に、全領域を空気と仮定して、上述した入射波の到達領域の予測計算や、未到達領域の予測計算を行なうことも可能である。通常の媒質の屈折率ｎは、１以上（ｎ≧１）であり、空気と仮定しておけば、屈折率条件は十分満足する。高屈折率が存在する場合に、必要以上の計算を要するが、処理をより簡便にできる。
【００７１】
上述したように、本発明は、入射波が所定方向に進行する計算モデルにおける入射波未到達の領域の計算量の削減による計算の効率化についての技術であり、上述した光学、通信、半導体、検出素子、構造物などの技術分野に適用可能なことはもちろんである。
【００７２】
なお、上述した波動場解析方法は、予め用意されたプログラムをパーソナルコンピュータやワークステーションなどのコンピュータで実行することによっても実現されるものである。このプログラムは、波動場解析方法をコンピュータに実行させるプログラムであって、ハードディスク、フロッピーディスク、ＣＤ−ＲＯＭ、ＭＯ、ＤＶＤ等のコンピュータで読み取り可能な記録媒体に記録され、コンピュータによって記録媒体から読み出されることによって実行される。またこのプログラムは、上記記録媒体を介して、インターネット等のネットワークを介して配布することができる。
【００７３】
以上、本発明に係る波動場解析方法および装置ならびに波動場解析方法をコンピュータに実行させるプログラムを記録したコンピュータに読み取り可能な記録媒体について詳細に説明したが、本発明は、以上の実施形態に限定されるものではなく、本発明の要旨を逸脱しない範囲において、各種の改良や変更を行ってもよい。
【図面の簡単な説明】
【００７４】
【図１】本発明の波動場解析装置が適用される電磁場解析装置の機能的構成の一実施例を示すブロック図である。
【図２】図１に示す波動場解析装置の到達領域予測部で得られる予測到達領域と計算ステップとの関係の一例を示すグラフである。
【図３】本発明に適用されるＦＤＴＤ法を説明する説明図である。
【図４】本発明に適用されるＦＤＴＤ法よる磁場および電場の計算を説明する説明図である。
【図５】図１に示す電磁場解析装置のハードウエア構成の一実施例を示すブロック図である。
【図６】本発明の波動場解析方法が適用される電磁場解析方法の一例を示すフローチャートである。
【図７】図６に示す電磁場解析方法が適用される計算モデルの一実施例を示す模式図である。
【図８】図７に示す計算モデルの解析領域とその屈折率分布および最小屈折率とを説明するための模式図である。
【図９】図８に示す解析領域における空間計算範囲の一例を示す模式図である。
【図１０】（ａ）は、図７に示す計算モデルの図８に示す解析領域について行われた計算結果を、所定の時間計算ステップ数毎に示す図であり、（ｂ）は、図８に示す解析領域の具体的な寸法を説明する説明図である。
【図１１】図８に示す解析領域における空間計算範囲と計算ステップとの関係の一例を示すグラフおよび計算量削減率を示す説明図である。
【図１２】従来の電磁場解析方法を示すフローチャートである。
【符号の説明】
【００７５】
１０，３０  電磁場解析装置
１２  屈折率分布定義部
１４  計算条件入力部
１６  到達領域予測部
１８  電磁場計算部
２０  記憶部
２２  書込／読出（制御）部
２４  計算結果出力部
３２  入力装置
３４  読取／書込装置
３６  演算処理装置（ＣＰＵ）
３８  メモリ
４０  表示装置
４２  記録装置（プリンタ）
４４  バス
５０  計算モデル
５０ａ  解析領域
５２，５２ａ  ガラス基板  
５４  シリンドリカルレンズアレー
５４ａ  シリンドリカルレンズ
５６，５６ａ  空気層
５８  励震部（波源）
６０  空間計算範囲
【図１】
2008250956.tif 000006
【図２】
2008250956.tif 000007
【図３】
2008250956.tif 000008
【図４】
2008250956.tif 000009
【図５】
2008250956.tif 000010
【図６】
2008250956.tif 000011
【図７】
2008250956.tif 000012
【図８】
2008250956.tif 000013
【図９】
2008250956.tif 000014
【図１１】
2008250956.tif 000015
【図１２】
2008250956.tif 000016
【図１０】
2008250956.tif 000017