(57)【要約】【課題】媒質の誘電率が空間的に不連続に変化する直線状の境界を有する２次元平面において、その境界が、２次元平面を格子状に離散化したときの互いに直交する格子軸のいずれに対しても平行でない場合、その離散化した格子状の平面を伝搬する電磁界成分の高精度なシミュレーションを可能にする電磁界シミュレーション方法、電磁界シミュレーションプログラム、およびプログラム記録媒体を提供する。【解決手段】互いに割り当てられる電磁界成分が異なる３種類の格子点をそれぞれ等間隔に配置して離散化した平面（Ｙｅｅ格子）を構成し、このうち電界成分が割り当てられる格子点のうち境界近傍に位置する格子点に対して所定の関係式を満たす誘電率をそれぞれ割り当てる。【選択図】      図４
2004341644.tif 000002
【特許請求の範囲】
【請求項１】
一つの平面内に誘電率ε１  、ε２  をそれぞれ有する二つの媒質が一つの直線を境界として接し、前記平面内に、当該平面に平行な磁界成分の値がゼロである電磁界が存在するときに、
前記平面に対応する離散化された仮想平面を、互いに割り当てられる電磁界が異なるそれぞれ複数の第１種、第２種、および第３種の格子点を用いて、同種の格子点がそれぞれ格子間隔Δの正方形状をなすとともに、第２種格子点が最近接する二つの第１種格子点の中点に位置し、なおかつ第３種格子点が最近接する二つの第１種格子点の中点のうち第２種格子点が位置する中点とは異なる中点に位置するように配置した上で、第１種格子点に前記平面と垂直な磁界成分を割り当てる一方で、第２種格子点には当該第２種格子点と最近接する第１種格子点を結ぶ線分に直交し、なおかつ前記平面に平行な電界成分を割り当て、さらに、第３種格子点には当該第３種格子点と最近接する第１種格子点を結ぶ線分に直交し、なおかつ前記平面に平行な電界成分を割り当てることによって構成し、
この構成した仮想平面における電磁界の時間変化を、コンピュータを用いた数値計算によって求める電磁界シミュレーション方法であって、
第１種格子点と第３種格子点とが交互に配置される列のいずれかをｘ軸とする一方、第１種格子点と第２種格子点とが交互に配置される列のいずれかをｙ軸とし、このｙ軸の負方向に誘電率ε１  を有する媒質が位置するようにした座標系を用いて前記境界をなす点のｙ座標をｘ座標の１次式で表すとき、当該直線の傾きの絶対値が２であり、なおかつ前記直線のｙ切片が、−１．５より大きく１．５より小さい範囲の任意の実数値をとる定数ｓおよび前記格子間隔Δを用いてｓΔと表される場合、
前記境界を表す直線と平行であり、なおかつｙ切片が−Δである直線上に配置される第３種格子点の誘電率をａ、
前記境界を表す直線と平行であり、なおかつｙ切片が−Δ／２である直線上に配置される第２種格子点の誘電率をｂ、
前記境界を表す直線と平行であり、なおかつｙ切片が０である直線上に配置される第３種格子点の誘電率をｃ、
前記境界を表す直線と平行であり、なおかつｙ切片がΔ／２である直線上に配置される第２種格子点の誘電率をｄ、
前記境界を表す直線と平行であり、なおかつｙ切片がΔである直線上に配置される第３種格子点の誘電率をｅ、
前記ａ，ｂ，ｃ，ｄ，ｅのいずれかが割り当てられていない第２種または第３種格子点のうち、前記仮想平面に対応する前記平面内で誘電率ε１の媒質に対応する格子点の誘電率をε１  、
前記ａ，ｂ，ｃ，ｄ，ｅのいずれかが割り当てられていない第２種または第３種格子点のうち、前記仮想平面に対応する前記平面内で誘電率ε２の媒質に対応する格子点の誘電率をε２として、
前記コンピュータが、
関係式
【数１】
2004341644.tif 000023
（ここで、係数ｆｉｊ（ｉ，ｊ＝０，１，２）は、
【数２】
2004341644.tif 000024
係数ｇｉｊｋ（ｉ＝０，１；ｊ，ｋ＝０，１，２）は、
【数３】
2004341644.tif 000025
係数ｈｉｊ（ｉ，ｊ＝０，１）は、
【数４】
2004341644.tif 000026
と定義される）を満たすように誘電率ａ，ｂ，ｃ，ｄ，ｅ，ε１  ，  ε２  を設定するステップと、
このステップで設定した誘電率ａ，ｂ，ｃ，ｄ，ｅ，ε１  ，  ε２  を用いることにより、各格子点に割り当てられた電磁界を、離散化された時間の変化に応じて順次計算してメモリに格納するステップと、
このステップで計算した電磁界を前記メモリから読み出して当該電磁界に関連する情報を計算するステップと
を実行することを特徴とする電磁界シミュレーション方法。
【請求項２】
一つの平面内に誘電率ε１  、ε２  をそれぞれ有する二つの媒質が一つの直線を境界として接し、前記平面内に、当該平面に平行な磁界成分の値がゼロである電磁界が存在するときに、
前記平面に対応する離散化された仮想平面を、互いに割り当てられる電磁界が異なるそれぞれ複数の第１種、第２種、および第３種の格子点を用いて、同種の格子点がそれぞれ格子間隔Δの正方形状をなすとともに、第２種格子点が最近接する二つの第１種格子点の中点に位置し、なおかつ第３種格子点が最近接する二つの第１種格子点の中点のうち第２種格子点が位置する中点とは異なる中点に位置するように配置した上で、第１種格子点に前記平面と垂直な磁界成分を割り当てる一方で、第２種格子点には当該第２種格子点と最近接する第１種格子点を結ぶ線分に直交し、なおかつ前記平面に平行な電界成分を割り当て、さらに、第３種格子点には当該第３種格子点と最近接する第１種格子点を結ぶ線分に直交し、なおかつ前記平面に平行な電界成分を割り当てることによって構成し、
この構成した仮想平面における電磁界の時間変化を、コンピュータを用いた数値計算によって求める電磁界シミュレーション方法であって、
第１種格子点と第３種格子点とが交互に配置される列のいずれかをｘ軸とする一方、第１種格子点と第２種格子点とが交互に配置される列のいずれかをｙ軸とし、このｙ軸の負方向に誘電率ε１  を有する媒質が位置するようにした座標系を用いて前記境界をなす点のｙ座標をｘ座標の１次式で表すとき、当該直線の傾きの絶対値が２であり、なおかつ前記直線のｙ切片が、−１より大きく１より小さい範囲の任意の実数値をとる定数δおよび前記格子間隔Δを用いて（１／２＋δ）Δと表される場合、
前記境界を表す直線と平行であり、なおかつｙ切片が０である直線上に配置される第３種格子点の誘電率をｕ、
前記境界を表す直線と平行であり、なおかつｙ切片がΔ／２である直線上に配置される第２種格子点の誘電率をｖ、
前記境界を表す直線と平行であり、なおかつｙ切片がΔである直線上に配置される第３種格子点の誘電率をｗ、
前記ｕ，ｖ，ｗのいずれかが割り当てられていない第２種または第３種格子点のうち、前記仮想平面に対応する前記平面内で誘電率ε１の媒質に対応する格子点の誘電率をε１  、
前記ｕ，ｖ，ｗのいずれかが割り当てられていない第２種または第３種格子点のうち、前記仮想平面に対応する前記平面内で誘電率ε２の媒質に対応する格子点の誘電率をε２として、
前記コンピュータが、
関係式
【数５】
2004341644.tif 000027
を満たすように設定される誘電率ｕ，ｖ，ｗ，ε１  ，  ε２  を設定するステップと、
このステップで設定した誘電率ｕ，ｖ，ｗ，ε１  ，  ε２  を用いることにより、各格子点に割り当てられた電磁界を、離散化された時間の変化に応じて順次計算してメモリに格納するステップと、
このステップで計算した電磁界を前記メモリから読み出して当該電磁界に関連する情報を計算するステップと
を実行することを特徴とする電磁界シミュレーション方法。
【請求項３】
２次精度を有する時間差分式、および２次精度を有する空間差分式を用いた時間領域差分法によって離散化された時間の変化に応じた電磁界を順次求めることを特徴とする請求項１または２記載の電磁界シミュレーション方法。
【請求項４】
一つの平面内に誘電率ε１  、ε２  をそれぞれ有する二つの媒質が一つの直線を境界として接し、前記平面内に、当該平面に平行な磁界成分の値がゼロである電磁界が存在するときに、
前記平面に対応する離散化された仮想平面を、互いに割り当てられる電磁界が異なるそれぞれ複数の第１種、第２種、および第３種の格子点を用いて、同種の格子点がそれぞれ格子間隔Δの正方形状をなすとともに、第２種格子点が最近接する二つの第１種格子点の中点に位置し、なおかつ第３種格子点が最近接する二つの第１種格子点の中点のうち第２種格子点が位置する中点とは異なる中点に位置するように配置した上で、第１種格子点に前記平面と垂直な磁界成分を割り当てる一方で、第２種格子点には当該第２種格子点と最近接する第１種格子点を結ぶ線分に直交し、なおかつ前記平面に平行な電界成分を割り当て、さらに、第３種格子点には当該第３種格子点と最近接する第１種格子点を結ぶ線分に直交し、なおかつ前記平面に平行な電界成分を割り当てることによって構成し、
この構成した仮想平面における電磁界の時間変化を、コンピュータを用いた数値計算によって求めるために、
第１種格子点と第３種格子点とが交互に配置される列のいずれかをｘ軸とする一方、第１種格子点と第２種格子点とが交互に配置される列のいずれかをｙ軸とし、このｙ軸の負方向に誘電率ε１  を有する媒質が位置するようにした座標系を用いて前記境界をなす点のｙ座標をｘ座標の１次式で表すとき、当該直線の傾きの絶対値が２であり、なおかつ前記直線のｙ切片が、−１．５より大きく１．５より小さい範囲の任意の実数値をとる定数ｓおよび前記格子間隔Δを用いてｓΔと表される場合、
前記境界を表す直線と平行であり、なおかつｙ切片が−Δである直線上に配置される第３種格子点の誘電率をａ、
前記境界を表す直線と平行であり、なおかつｙ切片が−Δ／２である直線上に配置される第２種格子点の誘電率をｂ、
前記境界を表す直線と平行であり、なおかつｙ切片が０である直線上に配置される第３種格子点の誘電率をｃ、
前記境界を表す直線と平行であり、なおかつｙ切片がΔ／２である直線上に配置される第２種格子点の誘電率をｄ、
前記境界を表す直線と平行であり、なおかつｙ切片がΔである直線上に配置される第３種格子点の誘電率をｅ、
前記ａ，ｂ，ｃ，ｄ，ｅのいずれかが割り当てられていない第２種または第３種格子点のうち、前記仮想平面に対応する前記平面内で誘電率ε１の媒質に対応する格子点の誘電率をε１  、
前記ａ，ｂ，ｃ，ｄ，ｅのいずれかが割り当てられていない第２種または第３種格子点のうち、前記仮想平面に対応する前記平面内で誘電率ε２の媒質に対応する格子点の誘電率をε２として、
前記コンピュータに、
関係式
【数６】
2004341644.tif 000028
（ここで、係数ｆｉｊ（ｉ，ｊ＝０，１，２）は、
【数７】
2004341644.tif 000029
係数ｇｉｊｋ（ｉ＝０，１；ｊ，ｋ＝０，１，２）は、
【数８】
2004341644.tif 000030
係数ｈｉｊ（ｉ，ｊ＝０，１）は、
【数９】
2004341644.tif 000031
と定義される）を満たすように誘電率ａ，ｂ，ｃ，ｄ，ｅ，ε１  ，  ε２  を設定するステップと、
このステップで設定した誘電率ａ，ｂ，ｃ，ｄ，ｅ，ε１  ，  ε２  を用いることにより、各格子点に割り当てられた電磁界を、離散化された時間の変化に応じて順次計算してメモリに格納するステップと、
このステップで計算した電磁界を前記メモリから読み出して当該電磁界に関連する情報を計算するステップと
を実行させることを特徴とする電磁界シミュレーションプログラム。
【請求項５】
一つの平面内に誘電率ε１  、ε２  をそれぞれ有する二つの媒質が一つの直線を境界として接し、前記平面内に、当該平面に平行な磁界成分の値がゼロである電磁界が存在するときに、
前記平面に対応する離散化された仮想平面を、互いに割り当てられる電磁界が異なるそれぞれ複数の第１種、第２種、および第３種の格子点を用いて、同種の格子点がそれぞれ格子間隔Δの正方形状をなすとともに、第２種格子点が最近接する二つの第１種格子点の中点に位置し、なおかつ第３種格子点が最近接する二つの第１種格子点の中点のうち第２種格子点が位置する中点とは異なる中点に位置するように配置した上で、第１種格子点に前記平面と垂直な磁界成分を割り当てる一方で、第２種格子点には当該第２種格子点と最近接する第１種格子点を結ぶ線分に直交し、なおかつ前記平面に平行な電界成分を割り当て、さらに、第３種格子点には当該第３種格子点と最近接する第１種格子点を結ぶ線分に直交し、なおかつ前記平面に平行な電界成分を割り当てることによって構成し、
この構成した仮想平面における電磁界の時間変化を、コンピュータを用いた数値計算によって求めるために、
第１種格子点と第３種格子点とが交互に配置される列のいずれかをｘ軸とする一方、第１種格子点と第２種格子点とが交互に配置される列のいずれかをｙ軸とし、このｙ軸の負方向に誘電率ε１  を有する媒質が位置するようにした座標系を用いて前記境界をなす点のｙ座標をｘ座標の１次式で表すとき、当該直線の傾きの絶対値が２であり、なおかつ前記直線のｙ切片が、−１より大きく１より小さい範囲の任意の実数値をとる定数δおよび前記格子間隔Δを用いて（１／２＋δ）Δと表される場合、
前記境界を表す直線と平行であり、なおかつｙ切片が０である直線上に配置される第３種格子点の誘電率をｕ、
前記境界を表す直線と平行であり、なおかつｙ切片がΔ／２である直線上に配置される第２種格子点の誘電率をｖ、
前記境界を表す直線と平行であり、なおかつｙ切片がΔである直線上に配置される第３種格子点の誘電率をｗ、
前記ｕ，ｖ，ｗのいずれかが割り当てられていない第２種または第３種格子点のうち、前記仮想平面に対応する前記平面内で誘電率ε１の媒質に対応する格子点の誘電率をε１  、
前記ｕ，ｖ，ｗのいずれかが割り当てられていない第２種または第３種格子点のうち、前記仮想平面に対応する前記平面内で誘電率ε２の媒質に対応する格子点の誘電率をε２として、
前記コンピュータに、
関係式
【数１０】
2004341644.tif 000032
を満たすように設定される誘電率ｕ，ｖ，ｗ，ε１  ，  ε２  を設定するステップと、
このステップで設定した誘電率ｕ，ｖ，ｗ，ε１  ，  ε２  を用いることにより、各格子点に割り当てられた電磁界を、離散化された時間の変化に応じて順次計算してメモリに格納するステップと、
このステップで計算した電磁界を前記メモリから読み出して当該電磁界に関連する情報を計算するステップと
を実行させることを特徴とする電磁界シミュレーションプログラム。
【請求項６】
２次精度を有する時間差分式、および２次精度を有する空間差分式を用いた時間領域差分法によって離散化された時間の変化に応じた電磁界を順次求めることを特徴とする請求項４または５記載の電磁界シミュレーションプログラム。
【請求項７】
請求項４乃至６のいずれか１項に記載した電磁界シミュレーションプログラムを記録したことを特徴とするプログラム記録媒体。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、コンピュータを用いて電磁界を解析する電磁界シミュレーション方法、電磁界シミュレーションプログラム、およびプログラム記録媒体に関する。
【０００２】
【従来の技術】
近年、“社会の情報化”が著しく進展している。その中で、情報のキャリアとしての電磁波は、ますます高度化された形で情報の伝達に利用されている。例えば、光導波路やアンテナ等は、それぞれの目的に合わせ、より精密に設計することが必要とされてきており、このような状況においては、電磁界の時間変化を詳細にシミュレーションすることが求められる。
【０００３】
上述した目的のために使用される電磁界シミュレーション方法としては、時間領域差分法（ＦＤＴＤ：  Ｆｉｎｉｔｅ−Ｄｉｆｆｅｒｅｎｃｅ  Ｔｉｍｅ−Ｄｏｍａｉｎ  ｍｅｔｈｏｄ，  例えば，  Ｙｅｅ，  Ｋ．  Ｓ．，  “Ｎｕｍｅｒｉｃａｌ  ｓｏｌｕｔｉｏｎ  ｏｆ  ｉｎｉｔｉａｌ  ｂｏｕｎｄａｒｙ  ｖａｌｕｅ  ｐｒｏｂｌｅｍｓ  ｉｎｖｏｌｖｉｎｇ  Ｍａｘｗｅｌｌ’ｓ  ｅｑｕａｔｉｏｎｓ  ｉｎ  ｉｓｏｔｒｏｐｉｃ  ｍｅｄｉａ”，  ＩＥＥＥ  Ｔｒａｎｓａｃｔｉｏｎｓ  ｏｎ  Ａｎｔｅｎｎａｓ  ａｎｄ  Ｐｒｏｐａｇａｔｉｏｎ，  Ｖｏｌ．Ａｐ−１４，  ｐｐ．  ３０２−３０７，  １９６６．  を参照）がある。このシミュレーション方法では、電磁界の時間変化を記述する基本方程式であるＭａｘｗｅｌｌの方程式を、Ｙｅｅ格子と呼ばれる２次元格子を用いて空間的に差分化し、コンピュータを用いて電磁界の時間変化を計算する。
【０００４】
時間領域差分法は、電磁界の時間変化を追跡するために必要な同一時刻の電磁界の空間偏微分を、離散化した空間差分式を用いて近似し、空間差分式の結果に基づいて、電磁界の時間変化を時間差分式を用いて順次計算するものである。
【０００５】
この時間領域差分法は、計算に用いる差分式の精度次数によって分類される。差分式の精度次数は、差分間隔（時間差分式の場合は時間ステップ、空間差分式の場合はセル幅）をτとし、τを０に近づけたときに、誤差がτｎ  に比例して減少する場合、ｎ次精度と呼ばれる。空間差分の場合は、誘電率と透磁率が空間的に一定な状況や、誘電率と透磁率の空間的変化が連続的な状況等において、この識別を適用する。
【０００６】
最も一般的に用いられている時間領域差分法は、時間２次精度、空間２次精度の時間領域差分法である。従来、シミュレーション領域内での誘電率の空間的な変化が連続的でない場合、すなわち、誘電率が不連続な境界が存在する場合に、２次精度を達成するためのさまざまな手法が知られており、例えば特許文献１では、境界が離散化された格子状の平面における格子軸のいずれかと平行な場合に、境界付近の電界成分の時間変化を計算するための誘電率を工夫して割り当てることにより、２次精度を達成する手法が開示されている。
【０００７】
【特許文献１】
特開２００１−２４３２１５号公報
【０００８】
【発明が解決しようとする課題】
しかしながら、上述した従来の時間２次精度、空間２次精度の時間領域差分法において、シミュレーション領域に異なる誘電率を有する媒質の境界があり、この境界が離散化した格子状の平面で互いに直交する格子軸のいずれに対しても平行でない場合、境界付近の誘電率をどのように割り当てるかが明らかではなかった。
【０００９】
本発明はこのような事情に鑑みてなされたものであり、その目的は、媒質の誘電率が空間的に不連続に変化する直線状の境界を有する２次元平面において、その境界が、２次元平面を格子状に離散化したときの互いに直交する格子軸のいずれに対しても平行でない場合、その離散化した格子状の平面を伝搬する電磁界成分の高精度なシミュレーションを可能にする電磁界シミュレーション方法、電磁界シミュレーションプログラム、およびプログラム記録媒体を提供することにある。
【００１０】
【課題を解決するための手段】
上記目的を達成するために、請求項１記載の発明は、一つの平面内に誘電率ε１  、ε２  をそれぞれ有する二つの媒質が一つの直線を境界として接し、前記平面内に、当該平面に平行な磁界成分の値がゼロである電磁界が存在するときに、前記平面に対応する離散化された仮想平面を、互いに割り当てられる電磁界が異なるそれぞれ複数の第１種、第２種、および第３種の格子点を用いて、同種の格子点がそれぞれ格子間隔Δの正方形状をなすとともに、第２種格子点が最近接する二つの第１種格子点の中点に位置し、なおかつ第３種格子点が最近接する二つの第１種格子点の中点のうち第２種格子点が位置する中点とは異なる中点に位置するように配置した上で、第１種格子点に前記平面と垂直な磁界成分を割り当てる一方で、第２種格子点には当該第２種格子点と最近接する第１種格子点を結ぶ線分に直交し、なおかつ前記平面に平行な電界成分を割り当て、さらに、第３種格子点には当該第３種格子点と最近接する第１種格子点を結ぶ線分に直交し、なおかつ前記平面に平行な電界成分を割り当てることによって構成し、この構成した仮想平面における電磁界の時間変化を、コンピュータを用いた数値計算によって求める電磁界シミュレーション方法であって、第１種格子点と第３種格子点とが交互に配置される列のいずれかをｘ軸とする一方、第１種格子点と第２種格子点とが交互に配置される列のいずれかをｙ軸とし、このｙ軸の負方向に誘電率ε１  を有する媒質が位置するようにした座標系を用いて前記境界をなす点のｙ座標をｘ座標の１次式で表すとき、当該直線の傾きの絶対値が２であり、なおかつ前記直線のｙ切片が、−１．５より大きく１．５より小さい範囲の任意の実数値をとる定数ｓおよび前記格子間隔Δを用いてｓΔと表される場合、前記境界を表す直線と平行であり、なおかつｙ切片が−Δである直線上に配置される第３種格子点の誘電率をａ、前記境界を表す直線と平行であり、なおかつｙ切片が−Δ／２である直線上に配置される第２種格子点の誘電率をｂ、前記境界を表す直線と平行であり、なおかつｙ切片が０である直線上に配置される第３種格子点の誘電率をｃ、前記境界を表す直線と平行であり、なおかつｙ切片がΔ／２である直線上に配置される第２種格子点の誘電率をｄ、前記境界を表す直線と平行であり、なおかつｙ切片がΔである直線上に配置される第３種格子点の誘電率をｅ、前記ａ，ｂ，ｃ，ｄ，ｅのいずれかが割り当てられていない第２種または第３種格子点のうち、前記仮想平面に対応する前記平面内で誘電率ε１の媒質に対応する格子点の誘電率をε１  、前記ａ，ｂ，ｃ，ｄ，ｅのいずれかが割り当てられていない第２種または第３種格子点のうち、前記仮想平面に対応する前記平面内で誘電率ε２の媒質に対応する格子点の誘電率をε２として、前記コンピュータが、関係式
【数１１】
2004341644.tif 000003
（ここで、係数ｆｉｊ（ｉ，ｊ＝０，１，２）は、
【数１２】
2004341644.tif 000004
係数ｇｉｊｋ（ｉ＝０，１；ｊ，ｋ＝０，１，２）は、
【数１３】
2004341644.tif 000005
係数ｈｉｊ（ｉ，ｊ＝０，１）は、
【数１４】
2004341644.tif 000006
と定義される）を満たすように誘電率ａ，ｂ，ｃ，ｄ，ｅ，ε１  ，  ε２  を設定するステップと、このステップで設定した誘電率ａ，ｂ，ｃ，ｄ，ｅ，ε１  ，  ε２  を用いることにより、各格子点に割り当てられた電磁界を、離散化された時間の変化に応じて順次計算してメモリに格納するステップと、このステップで計算した電磁界を前記メモリから読み出して当該電磁界に関連する情報を計算するステップとを実行することを要旨とする。
【００１１】
請求項２記載の発明は、一つの平面内に誘電率ε１  、ε２  をそれぞれ有する二つの媒質が一つの直線を境界として接し、前記平面内に、当該平面に平行な磁界成分の値がゼロである電磁界が存在するときに、前記平面に対応する離散化された仮想平面を、互いに割り当てられる電磁界が異なるそれぞれ複数の第１種、第２種、および第３種の格子点を用いて、同種の格子点がそれぞれ格子間隔Δの正方形状をなすとともに、第２種格子点が最近接する二つの第１種格子点の中点に位置し、なおかつ第３種格子点が最近接する二つの第１種格子点の中点のうち第２種格子点が位置する中点とは異なる中点に位置するように配置した上で、第１種格子点に前記平面と垂直な磁界成分を割り当てる一方で、第２種格子点には当該第２種格子点と最近接する第１種格子点を結ぶ線分に直交し、なおかつ前記平面に平行な電界成分を割り当て、さらに、第３種格子点には当該第３種格子点と最近接する第１種格子点を結ぶ線分に直交し、なおかつ前記平面に平行な電界成分を割り当てることによって構成し、この構成した仮想平面における電磁界の時間変化を、コンピュータを用いた数値計算によって求める電磁界シミュレーション方法であって、第１種格子点と第３種格子点とが交互に配置される列のいずれかをｘ軸とする一方、第１種格子点と第２種格子点とが交互に配置される列のいずれかをｙ軸とし、このｙ軸の負方向に誘電率ε１  を有する媒質が位置するようにした座標系を用いて前記境界をなす点のｙ座標をｘ座標の１次式で表すとき、当該直線の傾きの絶対値が２であり、なおかつ前記直線のｙ切片が、−１より大きく１より小さい範囲の任意の実数値をとる定数δおよび前記格子間隔Δを用いて（１／２＋δ）Δと表される場合、前記境界を表す直線と平行であり、なおかつｙ切片が０である直線上に配置される第３種格子点の誘電率をｕ、前記境界を表す直線と平行であり、なおかつｙ切片がΔ／２である直線上に配置される第２種格子点の誘電率をｖ、前記境界を表す直線と平行であり、なおかつｙ切片がΔである直線上に配置される第３種格子点の誘電率をｗ、前記ｕ，ｖ，ｗのいずれかが割り当てられていない第２種または第３種格子点のうち、前記仮想平面に対応する前記平面内で誘電率ε１の媒質に対応する格子点の誘電率をε１  、前記ｕ，ｖ，ｗのいずれかが割り当てられていない第２種または第３種格子点のうち、前記仮想平面に対応する前記平面内で誘電率ε２の媒質に対応する格子点の誘電率をε２として、前記コンピュータが、関係式
【数１５】
2004341644.tif 000007
を満たすように設定される誘電率ｕ，ｖ，ｗ，ε１  ，  ε２  を設定するステップと、このステップで設定した誘電率ｕ，ｖ，ｗ，ε１  ，  ε２  を用いることにより、各格子点に割り当てられた電磁界を、離散化された時間の変化に応じて順次計算してメモリに格納するステップと、このステップで計算した電磁界を前記メモリから読み出して当該電磁界に関連する情報を計算するステップとを実行することを要旨とする。
【００１２】
請求項３記載の発明は、２次精度を有する時間差分式、および２次精度を有する空間差分式を用いた時間領域差分法によって離散化された時間の変化に応じた電磁界を順次求めることを要旨とする。
【００１３】
請求項４記載の発明は、一つの平面内に誘電率ε１  、ε２  をそれぞれ有する二つの媒質が一つの直線を境界として接し、前記平面内に、当該平面に平行な磁界成分の値がゼロである電磁界が存在するときに、前記平面に対応する離散化された仮想平面を、互いに割り当てられる電磁界が異なるそれぞれ複数の第１種、第２種、および第３種の格子点を用いて、同種の格子点がそれぞれ格子間隔Δの正方形状をなすとともに、第２種格子点が最近接する二つの第１種格子点の中点に位置し、なおかつ第３種格子点が最近接する二つの第１種格子点の中点のうち第２種格子点が位置する中点とは異なる中点に位置するように配置した上で、第１種格子点に前記平面と垂直な磁界成分を割り当てる一方で、第２種格子点には当該第２種格子点と最近接する第１種格子点を結ぶ線分に直交し、なおかつ前記平面に平行な電界成分を割り当て、さらに、第３種格子点には当該第３種格子点と最近接する第１種格子点を結ぶ線分に直交し、なおかつ前記平面に平行な電界成分を割り当てることによって構成し、この構成した仮想平面における電磁界の時間変化を、コンピュータを用いた数値計算によって求めるために、第１種格子点と第３種格子点とが交互に配置される列のいずれかをｘ軸とする一方、第１種格子点と第２種格子点とが交互に配置される列のいずれかをｙ軸とし、このｙ軸の負方向に誘電率ε１  を有する媒質が位置するようにした座標系を用いて前記境界をなす点のｙ座標をｘ座標の１次式で表すとき、当該直線の傾きの絶対値が２であり、なおかつ前記直線のｙ切片が、−１．５より大きく１．５より小さい範囲の任意の実数値をとる定数ｓおよび前記格子間隔Δを用いてｓΔと表される場合、前記境界を表す直線と平行であり、なおかつｙ切片が−Δである直線上に配置される第３種格子点の誘電率をａ、前記境界を表す直線と平行であり、なおかつｙ切片が−Δ／２である直線上に配置される第２種格子点の誘電率をｂ、前記境界を表す直線と平行であり、なおかつｙ切片が０である直線上に配置される第３種格子点の誘電率をｃ、前記境界を表す直線と平行であり、なおかつｙ切片がΔ／２である直線上に配置される第２種格子点の誘電率をｄ、前記境界を表す直線と平行であり、なおかつｙ切片がΔである直線上に配置される第３種格子点の誘電率をｅ、前記ａ，ｂ，ｃ，ｄ，ｅのいずれかが割り当てられていない第２種または第３種格子点のうち、前記仮想平面に対応する前記平面内で誘電率ε１の媒質に対応する格子点の誘電率をε１  、前記ａ，ｂ，ｃ，ｄ，ｅのいずれかが割り当てられていない第２種または第３種格子点のうち、前記仮想平面に対応する前記平面内で誘電率ε２の媒質に対応する格子点の誘電率をε２として、前記コンピュータに、関係式
【数１６】
2004341644.tif 000008
（ここで、係数ｆｉｊ（ｉ，ｊ＝０，１，２）は、
【数１７】
2004341644.tif 000009
係数ｇｉｊｋ（ｉ＝０，１；ｊ，ｋ＝０，１，２）は、
【数１８】
2004341644.tif 000010
係数ｈｉｊ（ｉ，ｊ＝０，１）は、
【数１９】
2004341644.tif 000011
と定義される）を満たすように誘電率ａ，ｂ，ｃ，ｄ，ｅ，ε１  ，  ε２  を設定するステップと、このステップで設定した誘電率ａ，ｂ，ｃ，ｄ，ｅ，ε１  ，  ε２  を用いることにより、各格子点に割り当てられた電磁界を、離散化された時間の変化に応じて順次計算してメモリに格納するステップと、このステップで計算した電磁界を前記メモリから読み出して当該電磁界に関連する情報を計算するステップとを実行させることを要旨とする。
【００１４】
請求項５記載の発明は、一つの平面内に誘電率ε１  、ε２  をそれぞれ有する二つの媒質が一つの直線を境界として接し、前記平面内に、当該平面に平行な磁界成分の値がゼロである電磁界が存在するときに、前記平面に対応する離散化された仮想平面を、互いに割り当てられる電磁界が異なるそれぞれ複数の第１種、第２種、および第３種の格子点を用いて、同種の格子点がそれぞれ格子間隔Δの正方形状をなすとともに、第２種格子点が最近接する二つの第１種格子点の中点に位置し、なおかつ第３種格子点が最近接する二つの第１種格子点のうち第２種格子点が位置する中点とは異なる中点に位置するように配置した上で、第１種格子点に前記平面と垂直な磁界成分を割り当てる一方で、第２種格子点には当該第２種格子点と最近接する第１種格子点を結ぶ線分に直交し、なおかつ前記平面に平行な電界成分を割り当て、さらに、第３種格子点には当該第３種格子点と最近接する第１種格子点を結ぶ線分に直交し、なおかつ前記平面に平行な電界成分を割り当てることによって構成し、この構成した仮想平面における電磁界の時間変化を、コンピュータを用いた数値計算によって求めるために、第１種格子点と第３種格子点とが交互に配置される列のいずれかをｘ軸とする一方、第１種格子点と第２種格子点とが交互に配置される列のいずれかをｙ軸とし、このｙ軸の負方向に誘電率ε１  を有する媒質が位置するようにした座標系を用いて前記境界をなす点のｙ座標をｘ座標の１次式で表すとき、当該直線の傾きの絶対値が２であり、なおかつ前記直線のｙ切片が、−１より大きく１より小さい範囲の任意の実数値をとる定数δおよび前記格子間隔Δを用いて（１／２＋δ）Δと表される場合、前記境界を表す直線と平行であり、なおかつｙ切片が０である直線上に配置される第３種格子点の誘電率をｕ、前記境界を表す直線と平行であり、なおかつｙ切片がΔ／２である直線上に配置される第２種格子点の誘電率をｖ、前記境界を表す直線と平行であり、なおかつｙ切片がΔである直線上に配置される第３種格子点の誘電率をｗ、前記ｕ，ｖ，ｗのいずれかが割り当てられていない第２種または第３種格子点のうち、前記仮想平面に対応する前記平面内で誘電率ε１の媒質に対応する格子点の誘電率をε１  、前記ｕ，ｖ，ｗのいずれかが割り当てられていない第２種または第３種格子点のうち、前記仮想平面に対応する前記平面内で誘電率ε２の媒質に対応する格子点の誘電率をε２として、前記コンピュータに、関係式
【数２０】
2004341644.tif 000012
を満たすように設定される誘電率ｕ，ｖ，ｗ，ε１  ，  ε２  を設定するステップと、このステップで設定した誘電率ｕ，ｖ，ｗ，ε１  ，  ε２  を用いることにより、各格子点に割り当てられた電磁界を、離散化された時間の変化に応じて順次計算してメモリに格納するステップと、このステップで計算した電磁界を前記メモリから読み出して当該電磁界に関連する情報を計算するステップとを実行させることを要旨とする。
【００１５】
請求項６記載の発明は、２次精度を有する時間差分式、および２次精度を有する空間差分式を用いた時間領域差分法によって離散化された時間の変化に応じた電磁界を順次求めることを要旨とする。
【００１６】
請求項７記載の発明は、請求項４乃至６のいずれか１項に記載した電磁界シミュレーションプログラムを記録したことを要旨とする。
【００１７】
この発明におけるプログラム記録媒体とは、ハードディスク、フレキシブルディスク、ＣＤ−ＲＯＭ（Ｃｏｍｐａｃｔ  Ｄｉｓｃ  Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）、ＤＶＤ（Ｄｉｇｉｔａｌ  Ｖｅｒｓａｔｉｌｅ  Ｄｉｓｋ）、光磁気ディスク、ＰＣカード等のコンピュータ読み取り可能な記録媒体を意味する。
【００１８】
【発明の実施の形態】
次に、図面を参照して本発明の実施の形態を説明する。
【００１９】
図１は、本発明の一実施形態に係る電磁界シミュレーション方法を実行する電子的な装置である電磁界シミュレーション装置の構成を示す機能ブロック図である。同図に示す電磁界シミュレーション装置１は、シミュレーション用基礎データと各種パラメタを入力するためのキーボード、マウス等の入力装置から成る入力部１１、入力された基礎データや各種パラメタを用いて電磁界シミュレーションを行うために、演算機能および制御機能を備えた中央処理装置（ＣＰＵ：Ｃｅｎｔｒａｌ  Ｐｒｏｃｅｓｓｉｎｇ  Ｕｎｉｔ）から構成される演算部１２、入力部１１からの入力情報や演算部１２からの演算結果等を記憶する記憶部１３、記憶部１３で記憶した情報をシミュレーション結果として出力するためのディスプレイ装置等から成る出力部１６を少なくとも備えている。
【００２０】
記憶部１３は、ＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）等から構成される主記憶部１４と、ハードディスクドライブ、フレキシブルディスクドライブ、ＣＤ−ＲＯＭ（Ｃｏｍｐａｃｔ  Ｄｉｓｃ  Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）ドライブ、ＤＶＤ（Ｄｉｇｉｔａｌ  Ｖｅｒｓａｔｉｌｅ  Ｄｉｓｋ）ドライブ、光磁気ディスクドライブ、ＰＣカードドライブ等から構成される補助記憶部１５とを備えており、後述する計算結果を随時格納するために必要なメモリ領域が確保されている。
【００２１】
以後に説明する本実施形態の各種処理は、一つの電子的な装置が実行する場合だけでなく、各ステップの実行を適宜分割して二つ以上の電子的な装置から構築されたシステムが全体で実行する場合も含む。そこで、本実施形態においては、一つまたは複数の電子的な装置から成る装置（コンピュータ）を「電磁界シミュレーション装置」と呼ぶ。この意味で、本実施形態には、電磁界シミュレーション処理を行うためのプログラムの指示にしたがって、例えば電磁界シミュレーション装置１上で稼動するＯＳ（Ｏｐｅｒａｔｉｎｇ  Ｓｙｓｔｅｍ）が実際の処理を行う場合等が含まれることはいうまでもない。この点は、本発明の全ての実施の形態に共通する事項である。
【００２２】
図２は、本実施形態の電磁界シミュレーションを実行する際に、誘電率が不連続的に異なる二つの領域（誘電率をε１  、ε２  とする）が直線状の境界で隣接している連続媒質の平面を離散化して得られる格子状の仮想平面としての２次元Ｙｅｅ格子（以後、単にＹｅｅ格子と記載する）の概略構成を示す説明図である。同図に示すＹｅｅ格子２０は、３種類の異なるタイプの格子点２１（第１種格子点）、格子点２２（第２種格子点）、および格子点２３（第３種格子点）が交互に２次元的に配置され、各々が正方格子を単位格子とする２次元格子（同種の格子点間の格子間隔である格子定数をΔとする。ここでΔは正の定数）を構成することによって離散化された平面（仮想平面）を形成している。より具体的には、格子点を用いて、同種の格子点がそれぞれ格子間隔Δの正方形状をなすとともに、格子点２２が最近接する二つの格子点２１の中点に位置し、なおかつ格子点２３が最近接する二つの格子点２１の中点のうち格子点２２が位置する中点とは異なる中点に位置するように配置した上で、格子点２１に前記平面と垂直な磁界成分を割り当てることによってその離散化された平面を構成している。したがって、このＹｅｅ格子２０において、最近接する異種格子点間の距離はΔ／２である。
【００２３】
図２では、格子点２１と格子点２３とが交互に配置される列に平行な方向（図の水平方向）をｘ軸方向とし、この列のうちの一つをｘ軸としている。他方、格子点２１と格子点２２とが交互に配置される列に平行な方向（図の鉛直方向）をｙ軸とし、この列のうちの一つをｙ軸としている。このような座標軸の原点には格子点２１のいずれかが配置される。図２では、以上の設定を満たす座標系としてΣ０を取った場合を図示しているが、この座標系の取り方が任意であることはいうまでもない。
【００２４】
図２に示す座標系Σ０を用いると、任意の整数ｉおよびｊを用いることにより、格子点２１の座標を（ｉΔ，ｊΔ）、格子点２２の座標を（ｉΔ，（ｊ＋１／２）Δ）、格子点２３の座標を（（ｉ＋１／２）Δ，ｊΔ）とそれぞれ表すことができる。以後の説明においては、記述を簡単にするために、Δを明記しなくともその意味が十分理解できるとみなしうる場合には、各格子点のｘ、ｙ座標をそれぞれΔで割り、単なる整数ｉ、ｊの組み合わせとして、例えば上記格子点２１、２２、２３の座標をそれぞれ（ｉ，ｊ）、（ｉ，ｊ＋１／２，ｊ）、（ｉ＋１／２，ｊ）と表すことにする。
【００２５】
本実施形態においては、格子点２１には図２の平面に垂直な方向（ｚ方向）の磁界成分Ｈｚ（ｉ，ｊ）を割り当て、格子点２２には、ｘ方向の電界成分Ｅｘ（ｉ，ｊ＋１／２）を、格子点２３には、ｙ方向の電界成分Ｅｙ（ｉ＋１／２，ｊ）をそれぞれ割り当てる。これにより、Ｙｅｅ格子面に垂直な方向への振動を磁界のみが有する２次元ＴＭ（Ｔｒａｎｓｖｅｒｓｅ  Ｍａｇｎｅｔｉｃ）モードの電磁界をモデル化できる。ちなみに、この２次元ＴＭモードでは、ｚ方向の電界成分Ｅｚ  、ｘ方向の磁界成分Ｈｘ  、およびｙ方向の磁界成分Ｈｙは、時間によらず常に０（一定）である。
【００２６】
二つの媒質の境界２４は、図２にも示すように、格子の配列方向であるｘ軸およびｙ軸方向に対して所定の傾斜角をなす方向を指向する。より具体的には、ｙをｘの１次関数として境界２４を表現する場合、その直線の傾きの大きさが２である。そこで、以後の説明の便宜上、この境界をなす直線の傾きを−２とし、Δおよびｓ（−１．５＜ｓ＜１．５）を用いて
ｙ＝−２ｘ＋ｓΔ            （１）
とするが、本実施形態の電磁界シミュレーション方法は、傾きが＋２であってもよいし、ｙ切片が上記以外の値をとってもかまわない。これらは、座標系の取り方の任意性にも依存するものである。しかしながら、以上のような限定を加えても、本実施形態に係る電磁界シミュレーション方法の一般性が失われないことはいうまでもない。
【００２７】
式（１）で表される境界２４の近傍に配置される格子点２２および２３（ともにゼロでない電界成分を有する格子点）の誘電率は、次のように設定する：
（Ｉ−１）直線ｙ＝−２ｘ−Δ上に配置される格子点２３の誘電率をａとする。
（Ｉ−２）直線ｙ＝−２ｘ−Δ／２上に配置される格子点２２の誘電率をｂとする。
（Ｉ−３）直線ｙ＝−２ｘ上に配置される格子点２３の誘電率をｃとする。
（Ｉ−４）直線ｙ＝−２ｘ＋Δ／２上に配置される格子点２２の誘電率をｄとする。
（Ｉ−５）直線ｙ＝−２ｘ＋Δ上に配置される格子点２３の誘電率をｅとする。
【００２８】
これら以外の格子点２２および２３の誘電率は、次のように設定する：
（Ｉ−６）ｙ＜−２ｘ＋ｓΔの領域（この領域を領域Ｄ１  とする）に配置される上記以外の格子点２２および２３の誘電率をε１  とする。
（Ｉ−７）ｙ＞−２ｘ＋ｓΔの領域（この領域を領域Ｄ２  とする）に配置される上記以外の格子点２２および２３の誘電率をε２  とする。
【００２９】
なお、透磁率は、媒質に関わらず一定であるとする。
【００３０】
図３は、以上の誘電率の割り当てを一つの表３０にまとめた図である。同図に示す表３０において、例えば「ｙ＝−２ｘ−Δ」と記載されている列は、この直線上に配置される各格子点に割り当てられる誘電率を意味している。すなわち、上記（Ｉ−１）より、格子点２３には誘電率ａが割り当てられる。他方、直線ｙ＝−２ｘ−Δ上には格子点２２が存在しないため、表３０では格子点２２の誘電率を記載する欄に横線（−）を引いている。他の欄の記載についても同様である。
【００３１】
図４は、境界２４付近の各格子点に割り当てられる誘電率を概念的に示す説明図である。なお、図２と同じ記号によって記載されており、割り当てられている誘電率が明記されていない格子点２２および２３については、上記（Ｉ−６）および（Ｉ−７）にしたがって、適宜ε１またはε２が割り当てられているものとする。
【００３２】
領域Ｄ１側から領域Ｄ２側へ平面波が入射した時、磁界振幅に基づく境界での反射率Ｒおよび透過率Ｔは、格子定数Δの級数として、次のように表すことができる。
【数２１】
2004341644.tif 000013
ここで、θｉは入射角、θｔは反射角をそれぞれ表している。境界２４を誘電率ε１の領域と誘電率ε２の領域の境界とすると、式（２）および（３）の右辺第１項はともに真値である。
【００３３】
本実施形態では、１次の誤差項の係数Ｒ１およびＴ１を入射角θｉの値に関わらずゼロとする２次精度のシミュレーションを行う。そこで、式（２）および（３）で、Δの１次の係数が入射角θｉの値に関わらずにゼロとなる条件を求めると、次の式が得られる。
【数２２】
2004341644.tif 000014
ここで、式（４）の係数ｆｉｊ（ｉ，ｊ＝０，１，２）、式（５）の係数ｇｉｊｋ（ｉ＝０，１；ｊ，ｋ＝０，１，２）、および式（６）の係数ｈｉｊ（ｉ，ｊ＝０，１）は、上記（Ｉ−１）〜（Ｉ−５）に基づいて設定される誘電率ａ，ｂ，ｃ，ｄ，ｅを用いてそれぞれ次のように表される。
【００３４】
【数２３】
2004341644.tif 000015
【００３５】
【数２４】
2004341644.tif 000016
【００３６】
【数２５】
2004341644.tif 000017
【００３７】
ちなみに、領域Ｄ２から領域Ｄ１へ平面波が入射した場合に２次精度を与える条件式を求めても、式（４）〜（６）と同じ式が得られる。
【００３８】
以上説明した誘電率の設定法を、「誘電率設定法Ｉ」と呼ぶ。
【００３９】
ところで、この誘電率設定法Ｉにおいて、さらに、
【数２６】
2004341644.tif 000018
とおき、これらをｕ，ｖ，ｗについて解くことにより、次の式（７）〜（９）を得る。
【数２７】
2004341644.tif 000019
【００４０】
これらの式（７）〜（９）を用いる場合には、境界２４を表す直線を、
【数２８】
2004341644.tif 000020
とし、格子点２２および２３の誘電率を、次のように割り当てる：
まず、上記（Ｉ−３）〜（Ｉ−５）のｃ，ｄ，ｅをそれぞれｕ，ｖ，ｗとし、
（ＩＩ−１）直線ｙ＝−２ｘ上に配置される格子点２３の誘電率をｕとする。
（ＩＩ−２）直線ｙ＝−２ｘ＋Δ／２上に配置される格子点２２の誘電率をｖとする。
（ＩＩ−３）直線ｙ＝−２ｘ＋Δ上に配置される格子点２３の誘電率をｗとする。
さらに、
（ＩＩ−４）上記３つの条件に含まれない格子点２２および２３のうち、ｙ＜−２ｘ＋（１／２＋δ）Δの領域（この領域を領域Ｄ１とする）に配置される格子点２２および２３の誘電率をε１  とする。
（ＩＩ−５）上記４つの条件に含まれない格子点２２および２３の誘電率をε２  とする。この誘電率ε２  が割り当てられる領域は、ｙ＞−２ｘ＋（１／２＋δ）Δの領域（この領域を領域Ｄ２とする）内に含まれる。
【００４１】
条件（ＩＩ−１）〜（ＩＩ−５）および式（７）〜（９）に基づく誘電率設定法を「誘電率設定法ＩＩ」と呼ぶことにする。この誘電率設定法ＩＩが、誘電率設定法Ｉの特別な場合であることは、以上の説明からも明らかである。
【００４２】
図５は、上記各条件を一つの表４０にまとめて記載した図であり、表４０の見方は図３と同様である。
【００４３】
また、図６の説明図には、境界２４付近の各格子点に割り当てられる誘電率を概念的に示した。同図においても、図２と同じ記号で記載され、割り当てられている誘電率が明記されていない格子点２２および２３には、上記（ＩＩ−４）および（ＩＩ−５）にしたがって、適宜ε１またはε２が割り当てられているものとする。
【００４４】
電磁界の各成分は、以上説明したように、空間的に異なった位置に割り当てられるが、時間的な意味においても、Ｙｅｅ格子２０の平面に時間の次元を加えて成る３次元時空間の中の異なった位置に割り当てられる。例えば、磁界が、離散化された時刻ｔ１，ｔ２，…，ｔｎ−１，ｔｎ，ｔｎ＋１，…（ｎは正の整数）に割り当てられるとすると、電界はｔ１／２，ｔ１＋１／２，ｔ２＋１／２，…，ｔｎ−１／２，ｔｎ＋１／２，ｔｎ＋３／２，…というように、磁界成分が割り当てられた時刻の中間の時刻に割り当てられる。このように時刻を設定することによって、磁界と電界を交互に計算することが可能になる。なお、電界と磁界に割り当てられる時刻を逆に定義することも勿論可能である。
【００４５】
次に、以上の設定に基づいた電磁界シミュレーション方法について詳細に説明する。
【００４６】
本実施形態で適用する時間２次精度、空間２次精度の時間領域差分法では、電磁界シミュレーション装置１の記憶部１３内のメモリに記憶された電磁界成分から、時系列的な意味で次の時刻の電磁界が以下に示す演算によって順次求められる。
【００４７】
【数２９】
2004341644.tif 000021
ここで、Δｔは時間ステップ間隔、Ｅｘｎ＋１／２（ｉ，ｊ＋１／２）は、格子点２２（ｉ，ｊ＋１／２）の時刻ｔｎ＋１／２  における電界のｘ方向成分、Ｅｙｎ＋１／２（ｉ＋１／２．ｊ）は格子点２３の時刻ｔｎ＋１／２  における電界のｙ方向成分、Ｈｚｎ＋１（ｉ，ｊ）は格子点２１の時刻ｔｎ＋１  における磁界のｚ方向成分である。また、ε（ｉ，ｊ＋１／２）は座標（ｉ，ｊ＋１／２）の格子点２２に付与される誘電率、ε（ｉ＋１／２，ｊ）は座標（ｉ＋１／２，ｊ）の格子点２３に付与される誘電率である。さらに、μは透磁率であり、前述したように媒質に関わらず一定値をとる。
【００４８】
なお、各式の右辺にある［  ］は、その［  ］内に記載された偏微分を近似する空間２次精度の差分式を意味しており、より具体的には以下の差分式で定義される：
【数３０】
2004341644.tif 000022
【００４９】
媒質の誘電率が一様な場合のように、各電磁界成分が少なくとも２回以上偏微分可能としてよい場合は、これらの差分式が空間２次精度を与える事を各電磁界成分の空間座標（ｘ，ｙ）についてのテイラー展開を用いて示す事ができる。一方、境界２４の近傍では、各電磁界成分が境界２４と垂直な方向（ｘ方向）に対して一般に２回偏微分可能であるとは限らない。特に、境界２４に垂直な方向の電界成分は不連続であるため、これまで高精度の計算が可能であるかどうかは自明ではなかった。
【００５０】
これに対して、本実施形態においては、境界２４が、座標系Σ０において傾き２または−２の直線で表される場合、境界２４近傍の格子点２２、２３に割り当てる誘電率を条件式（４）〜（６）または式（７）〜（９）を満たすように設定し、通常の２次精度の空間差分式（式（１４）〜（１７）を参照）を用いることにより、境界２４に入射する電磁波の入射角に関わらず境界２４での反射率と透過率を２次精度で近似可能な２次元ＴＭモードの電磁界シミュレーションを実現することができる。
【００５１】
図７は、本実施形態に係る電磁界シミュレーション方法として、時間２次精度、空間２次精度の時間領域差分法を用いる際の具体的な処理手順を示すフローチャート図である。同図に示す処理を行うにあたっては、シミュレーション領域の格子点への誘電率の設定等、シミュレーション用基礎データおよび各種パラメタが予め入力部１１から入力され、記憶部１３に記憶されており、これらの入力内容に基づいて、誘電率の具体的な設定が演算部１２によってなされているものとする。
【００５２】
図７に示す時間領域差分法においては、まず、電界と磁界の初期値を設定する（ステップＳ１０１）。ここで、初期磁界の時刻をｔ＝−１／２、初期電界の時刻をｔ＝０とする。
【００５３】
次に、初期時刻（時間ステップ）をｎ＝０と設定する（ステップＳ１０２）。いうまでもなく、図７の記号「：＝」は、「ｎを０と定義する」ことを意味している。以後の時間ステップは、上述したように、電界がｔ１／２，ｔ１＋１／２，ｔ２＋１／２，…，ｔｎ−１／２，ｔｎ＋１／２，ｔｎ＋３／２，…、磁界がｔ１，ｔ２，…，ｔｎ−１，ｔｎ，ｔｎ＋１，…と推移して時間が経過していく。
【００５４】
この後、時間の経過にしたがって電磁界成分を演算部１２において順次計算する。そこで、以後のステップは、一般のｎ（正の整数）で表される任意の時間ステップについて説明する。
【００５５】
最初に、式（１１）、（１２）を用いて時間ステップｔ＝ｎ＋１／２の電界の各成分Ｅｘｎ＋１／２（ｉ，ｊ＋１／２）、Ｅｙｎ＋１／２（ｉ＋１／２，ｊ）を計算し、その計算結果を格子点と時間ステップの組（３次元時空間の座標に相当）に対応づけて記憶部１３内のメモリに格納する（ステップＳ１０３）。
【００５６】
式（１３）を用いて時間ステップｔ＝ｎ＋１の磁界成分Ｈｚｎ＋１（ｉ，ｊ）を計算し、その計算結果をステップＳ１０３と同様にメモリに格納する（ステップＳ１０４）。
【００５７】
電磁界成分の計算結果は、順次記憶部１３内のメモリに格納されて記憶されることになるが、前の時刻（時間ステップ）の電界または磁界を、次の時刻（時間ステップ）の電界または磁界に置換しても、以後の電磁界の計算に支障は無い。なぜならば、電界と磁界の計算は、式（１１）〜（１３）からも明らかなように、時間ステップごとに交互に実行されるからである。このようにして電磁界成分を順次更新することにより、計算時に用いるメモリ領域を削減して、負荷を軽減することが可能になる。
【００５８】
続いて、ステップＳ１０３およびＳ１０４でそれぞれ求めた電界成分Ｅｘｎ＋１／２（ｉ，ｊ＋１／２）、Ｅｙｎ＋１／２（ｉ＋１／２，ｊ）、および磁界成分Ｈｚｎ＋１（ｉ，ｊ）から、それらに関するデータを抽出してメモリへ格納する（ステップＳ１０５）。ちなみに、この抽出された電磁界成分に関するデータは、後述するステップにおいて電磁界の各モードの強度や位相等のデータを計算する際に使用することもできる。この意味で、これらのデータを格納するメモリ領域は、電磁界成分自体を格納するメモリ領域とは異なる領域であることが好ましい。なお、実行するシミュレーションの内容によっては、このステップＳ１０５が必須ではなく、ステップＳ１０４終了後、以下に説明するステップＳ１０６に推移するようにしてもよい。
【００５９】
この後、ｎの値を１増やして時間ステップを進め（ステップＳ１０６）、ステップＳ１０７において、ｎが予め設定した閾値ｎｍａｘに達していない（ｎ＜ｎｍａｘ）と判定した場合には、ステップＳ１０３に戻って上述したステップＳ１０３〜Ｓ１０６の処理を繰り返す。なお、ステップＳ１０７の判定が演算部１２で行われることはいうまでもない。
【００６０】
他方、ステップＳ１０７でｎがｎｍａｘに達したと判定した場合には、電磁界分布を出力後、電磁界の振幅や位相、反射率等を計算して計算結果を出力部１６から出力する（ステップＳ１０８）。
【００６１】
以上の処理によって得られるシミュレーション結果、およびシミュレーションのために入力される基礎データおよび各種パラメタ等は、必要に応じて記憶部１３から読み出して出力部１６から出力させることができる。
【００６２】
以上説明したように、本発明の一実施形態によれば、境界が、所定の座標系を用いてある傾きを有する直線で表される場合、境界近傍の格子点に対して条件式（４）〜（６）または式（７）〜（９）を満たすような誘電率を設定し、時間２次精度、空間２次精度の時間領域差分法を用いることにより、境界に入射する電磁波の入射角に関わらず境界での反射率と透過率を２次精度で近似可能な２次元ＴＭモードの電磁界シミュレーションを実現することができる。
【００６３】
なお、本実施形態の電磁界シミュレーション処理を実現するための電磁界シミュレーションプログラムを記録したコンピュータ読み取り可能なプログラム記録媒体をコンピュータに装着し、そのプログラム記録媒体に格納されているプログラムを読み出すことによって、コンピュータが上述した処理を実行するようにしてもよい。ここで、「コンピュータ読み取り可能な」プログラム記録媒体としては、ハードディスク、フレキシブルディスク、ＣＤ−ＲＯＭ、ＤＶＤ、光磁気ディスク、ＰＣカード等を用いることができる。このようなプログラム記録媒体を提供することによって、本実施形態の電磁界シミュレーションプログラムを広く流通させることができるようになる。
【００６４】
（その他の実施形態）
ところで、本発明は、上述した一実施形態においてのみ特有の効果を奏するものと理解されるべきではない。
【００６５】
例えば、境界近傍の誘電率を設定するに際して、より広範囲の境界近傍の格子点に対して、もとの連続媒質における誘電率（ε１  、ε２）とは異なる誘電率を割り当てることができる。具体的には、境界近傍において、その境界をなす直線と平行でｙ切片が異なるもののうち、境界をなす直線との距離が短い順にｍ個（ｍは２以上の任意の整数）の直線上を通過する格子点に対し、その直線ごとに誘電率の設定を行うことができる（誘電率設定法Ｉがｍ＝５の場合であり、誘電率設定法ＩＩがｍ＝３の場合に相当）。
【００６６】
また、境界をなす直線のＹｅｅ格子の軸方向に対する傾きは、上述した場合に限定されるわけではなく、軸方向に対して他の傾きを有する境界に対しても、同様に電磁界シミュレーションを実行することが可能である。
【００６７】
さらに、正方格子ではないＹｅｅ格子を用いて電磁界シミュレーションを行うことも勿論可能である。
【００６８】
このように、本発明は、上記一実施形態と同様に、時間２次精度、空間２次精度の時間領域差分法を用いることによって高精度の電磁界シミュレーションを実現することのできる様々な実施の形態等を含みうるものである。
【００６９】
【実施例】
以下、本発明の一実施形態に係る電磁界シミュレーション方法を用いた電磁界シミュレーションの実施例について説明する。
【００７０】
＜シミュレーション領域＞
図８は、本実施例において適用されるシミュレーション領域の概要を示す説明図である。なお、離散化の対象となる連続媒質の２次元平面は、誘電率がε２＝１０．２４の領域である半導体スラブ導波路５１と誘電率がε１＝１．０の領域であるバックグラウンド５２が、図８に示すように配置されたものである。
【００７１】
図８に示すシミュレーション領域５０は長方形をなしており、半導体スラブ導波路５１のバックグラウンド５２との境界では、境界をなす直線が、その境界付近に設けられる適当な局所座標系において、式（１０）を満たすようになっている。すなわち、境界５３では、図中鉛直下向きをｘ軸、水平方向左向きをｙ軸、図面に垂直で図の裏側に向かう方向をｚ軸とする局所座標系Σ１を取る一方で、境界５４では、図中鉛直上向きをｘ軸、水平方向右向きをｙ軸、図面に垂直で図の裏側に向かう方向をｚ軸とする局所座標系Σ２を取ることにより、各境界５３および５４が式（１０）の直線として表される。
【００７２】
シミュレーション領域５０内には図２と同様のＹｅｅ  格子２０が形成されており、局所座標系Σ１もしくはΣ２に基づいて各点の座標を定めたとき、磁界成分Ｈｚ（ｉ，ｊ）を割り当てる格子点２１（ｉ，ｊ）、電界成分Ｅｘ（ｉ，ｊ＋１／２）を割り当てる格子点２２（ｉ，ｊ＋１／２）、電界成分Ｅｙ（ｉ＋１／２，ｊ）を割り当てる格子点２３（ｉ＋１／２，ｊ）が交互に配列され、各格子点はそれぞれ同じ格子定数Δを有する正方格子をなしている。
【００７３】
＜シミュレーション用基礎データおよび各種パラメタの設定＞
シミュレーション用基礎データおよび各種パラメタについて説明する。本実施例では、各格子点に割り当てる誘電率の設定を誘電率設定法ＩＩに基づいて設定するものとする。
【００７４】
境界５３の局所座標系Σ１における直線の式（１０）のδの値を０とし、シミュレーション領域５０の対角線の長さを、伝搬シミュレーションを行う固有モードの整数倍、すなわち（２π／［伝搬シミュレーションを行うモードの伝搬定数］）の整数倍とした。また、シミュレーション領域５０の周囲には周期境界条件を課し、半導体スラブ導波路５１の幅を１．５μｍ（１．５×１０−６  ｍ）とした。
【００７５】
これに対して、境界５４の局所座標系Σ２における直線の式（１０）のδの値は、格子定数Δの値に応じて若干変化するが、その変化率がたかだか２％程度であるようにδの値を定めることができる。このδの具体的な値の例については後述する。
【００７６】
なお、以上の設定とは逆にして、境界５４の局所座標系Σ２における直線の式（１０）のδの値を０と先に定めた後、境界５３の局所座標系Σ１における直線の式（１０）のδの値を定めても勿論かまわない。
【００７７】
次に、伝搬シミュレーションを行う電磁波の設定について説明する。
【００７８】
まず、伝搬する電磁波の自由空間波長を１．５５μｍとした。この場合、磁界成分Ｈｚ  、および二つの方向の電界成分Ｅｘ  ，Ｅｙ  によって表される電磁波の伝搬モードは、６つのＴＭモード（０次〜５次）を保持することができる。各モードは、半導体スラブ導波路５１を伝搬する波の導波路境界への入射角が異なる。本実施例では、入射角が６１．４度である２次モードと入射角が２２．１度である５次モードについて伝搬シミュレーションを行った。
【００７９】
各格子点に割り当てられる電磁界成分の初期値については、解析的に得られる結果を初期値として各格子点に割り当てた。
【００８０】
図９は、半導体スラブ導波路５１とバックグラウンド５２の境界５３または５４近傍の格子点２２または２３に付与される誘電率の設定例を示す図である。同図に示す表６０においては、δ＝０．０となる境界（境界５３、５４のうち少なくとも１つの境界ではこの値をとる）において、ｕ＝０．８２６３６，  ｖ＝１．９３９３，  ｗ＝１６．０６１とした。
【００８１】
他方、δ≠０となる境界において、半導体スラブ導波路５１内の波長あたりのセル数（以後、導波路内波長あたりのセル数）が２３．１９５でδ＝０．４の場合、ｕ＝０．８１７９３，  ｖ＝０．９６９３１，  ｗ＝７．１７２２となる。表６０では、δ＝０．４のところに、括弧書きで０．３９９４７（＝０．４−０．００５３）という値が記載されているが、これは、δ≠０の場合には、格子定数Δ（導波路内波長あたりのセル数は格子の細かさの度合いを表す指標である）の値に応じて若干変化するためである。しかしながら、この変化は上述したようにたかだか２％程度なので、これが本実施例の電磁界シミュレーション結果に対して特段の影響を及ぼすことはない。
【００８２】
また、半導体スラブ導波路５１を構成する導波路内波長あたりのセル数が１７．５０３でδ＝−０．２の場合、ｕ＝０．９６９７０，  ｖ＝２．７４７１７，  ｗ＝２３．４４１６となる。この場合に表６０に括弧書きで記載されている値−０．１９８２８（＝−０．２＋０．００１７２）は、導波路内波長あたりのセル数が１７．５０３のときにシミュレーションに用いられる正確な値を示している。
【００８３】
ｕ，ｖ，ｗ以外の格子点の誘電率については、バックグラウンド５２を真空としてその誘電率ε１を１．０とする一方、上記以外の半導体スラブ導波路５１の誘電率ε２は１０．２４とした。
【００８４】
なお、本実施例のシミュレーションにおいては、Ｙｅｅ格子による離散化を行う際に、空間のスライス幅に対して時間のスライス幅をどの程度にするかの指標を与えるＣＦＬ数（Ｃｏｕｒａｎｔ−Ｆｒｉｅｄｒｉｃｈｓ−Ｌｅｖｙ  ｎｕｍｂｅｒ）を０．６とした。
【００８５】
＜シミュレーションの内容＞
以上説明したシミュレーション用基礎データおよび各種パラメタの設定に基づいて、２次元ＴＭモードのうちの２次モードと５次モードにつき、さまざまな格子定数Δを有するＹｅｅ格子に対してシミュレーションを行った。
【００８６】
反射率Ｒおよび透過率Ｔの精度は、電磁波の位相速度の精度により評価できるので、本実施例では、この位相速度の誤差を調べた。
【００８７】
シミュレーション結果としての位相速度は、伝搬シミュレーションの電磁界分布から位相情報を抽出することによって求められる（図７のステップＳ１０８）。
【００８８】
他方、半導体スラブ導波路５１を伝搬する電磁界の伝搬モード（導波路モード）は、前述したように解析的に計算することができる。
【００８９】
そこで、これら二つの結果を用いることにより、位相速度誤差を厳密に評価した。
【００９０】
＜評価結果＞
図１０は、半導体スラブ導波路５１の導波路内波長あたりのセル数（横軸）と位相速度誤差（縦軸）の関係を示す図である。同図に示すグラフ７０では、半導体スラブ導波路５１内での導波路内波長当たりのセル数に対する位相速度誤差を百分率で表示した。
【００９１】
位相速度誤差は、２次モードが直線２１０（図で□の点を通過する実線で表示）、５次モードが直線５１０（図で○の点を通過する実線で表示）でそれぞれ示されている。これらの直線は、セル数を２倍にすると、位相速度誤差が１／４になっており、その傾きから２次精度を実現していることが分かる。
【００９２】
これに対し、曲線２２０（図で×の点を通過する破線で表示）および５２０（図で△の点を通過する破線で表示）は、従来法による電磁界シミュレーションに基づいて位相速度誤差を評価した結果を与えるものである。
【００９３】
ここで、従来法の概要を説明する。この方法は、使用する格子点の種類と位置ならびに電磁界の初期値、電磁界の時間変化を計算する差分式は本発明と同一である一方で、境界近傍の格子点への誘電率の割り当て方のみが異なっており、誘電率を割り当てるべき格子点の近傍の誘電率を算術平均し、その平均値を割り当てるものである。より具体的には、ｘ軸もしくはｙ軸に平行な長さΔの辺から成る正方形において、その正方形内部の誘電率を面積で重み付けした平均値を、正方形の対角線の交点に設けられる格子点に割り当てる。
【００９４】
この従来法は、２次元ＴＥ（Ｔｒａｎｓｖｅｒｓｅ  Ｅｌｅｃｔｒｉｃ）モードについて適用可能であることが知られている（例えば、Ａｌｌｅｎ  Ｔａｆｌｏｖｅ，  ”Ｃｏｍｐｕｔａｔｉｏｎａｌ  Ｅｌｅｃｔｒｏｄｙｎａｍｉｃｓ，  Ｔｈｅ  Ｆｉｎｉｔｅ−Ｄｉｆｆｅｒｅｎｃｅ  Ｔｉｍｅ−Ｄｏｍａｉｎ  Ｍｅｔｈｏｄ”，  ｐ．  ３７４，  Ａｒｔｅｃｈ  Ｈｏｕｓｅ，  １９９５を参照）が、ここでは本実施例と比較するために２次元ＴＭモードに適用し、半導体スラブ導波路５１（を構成するＹｅｅ格子上）における電磁波の伝搬シミュレーションを実行した。その際、導波路内波長当たりのセル数とδの値は、本実施例のシミュレーションと同じ値を採用した。
【００９５】
図１０からも明らかなように、従来法によって得られた曲線２２０および曲線５２０は、直線２１０や直線５１０のような直線ではないことに加え、その傾きの絶対値も直線２１０や５１０に比べて明らかに小さい。この結果から、本発明の電磁界シミュレーション方法が、従来法よりも高次の精度（２次精度）を実現していることが明らかになった。
【００９６】
以上説明した本発明の一実施例によれば、格子点２２および格子点２３のうち、境界近傍の格子点に対して式（７）〜（９）の関係式を満たす誘電率を付与し、時間差分式として式（１１）〜（１３）で定義される２次精度の差分式を用いることにより、誘電率が不連続に変化する二つの媒質の境界が、シミュレーションを行うＹｅｅ格子軸に対して所定角傾斜している場合にも、２次元ＴＭモードに対する精度の高い電磁界シミュレーションを実現できることが分かる。
【００９７】
【発明の効果】
以上の説明からも明らかなように、本発明によれば、媒質の誘電率が空間的に不連続に変化する直線状の境界を有する２次元平面において、その境界が、２次元平面を格子状に離散化したときの互いに直交する格子軸のいずれに対しても平行でない場合、その離散化した格子状の平面を伝搬する電磁界成分の高精度なシミュレーションを可能にする電磁界シミュレーション方法、電磁界シミュレーションプログラム、およびプログラム記録媒体を提供することができる。
【図面の簡単な説明】
【図１】本発明の一実施形態に係る電磁界シミュレーション装置の機能ブロック図である。
【図２】本発明の一実施形態に係る電磁界シミュレーション方法で用いられるＹｅｅ格子の概略構成を表す説明図である。
【図３】図２のＹｅｅ格子の第２種および第３種格子点に割り当てられる誘電率の設定例（誘電率設定法Ｉ）を示す説明図である。
【図４】誘電率設定法Ｉに従って割り当てられたＹｅｅ格子の境界付近の格子点の誘電率を例示する図である。
【図５】図２のＹｅｅ格子の第２種および第３種格子点に割り当てられる誘電率の設定例（誘電率設定法ＩＩ）を示す説明図である。
【図６】誘電率設定法Ｉに従って割り当てられたＹｅｅ格子の境界付近の格子点の誘電率を例示する図である。
【図７】本発明の一実施形態に係る電磁界シミュレーション方法の概要を示すフローチャート図である。
【図８】本発明の一実施例として用いられるシミュレーション領域を示す説明図である。
【図９】図８のシミュレーション領域の半導体スラブ導波路とバックグラウンドの境界近傍の第２種および第３種格子点に付与される誘電率を示す図である。
【図１０】本発明の一実施形態に係る電磁界シミュレーション方法によって位相速度誤差を評価した結果を示す図である。
【符号の説明】
１  電磁界シミュレーション装置
１１  入力部
１２  演算部
１３  記憶部
１４  主記憶部
１５  補助記憶部
１６  出力部
２０  Ｙｅｅ  格子
２１、２２、２３  格子点
２４、５３、５４  境界
５０  シミュレーション領域
５１  半導体スラブ導波路
５２  バックグラウンド
２１０、５１０  直線
２２０、５２０  曲線
【図１】
2004341644.tif 000033
【図２】
2004341644.tif 000034
【図３】
2004341644.tif 000035
【図４】
2004341644.tif 000036
【図５】
2004341644.tif 000037
【図６】
2004341644.tif 000038
【図７】
2004341644.tif 000039
【図８】
2004341644.tif 000040
【図９】
2004341644.tif 000041
【図１０】
2004341644.tif 000042