(57)【特許請求の範囲】
【請求項１】
１つの平面内に透磁率μ・誘電率ε１を有する第１の媒質と、透磁率μ・誘電率ε２を有する第２の媒質とが１つの直線を境界として接し、前記平面内に、前記平面に平行な交流電界成分を持たない電磁波が存在する場合に、前記平面と等価な仮想平面を前記平面と同一視し、
前記仮想平面内に、タイプ１の格子点を、前記境界に垂直なｘ方向及び平行なｙ方向に、それぞれ、間隔△ｘ及び△ｙで二次元配置し、タイプ２の格子点を、前記タイプ１の格子点の位置からｘ方向に△ｘ/２だけ移動した位置に配置し、タイプ３の格子点を、前記タイプ１の格子点の位置からｙ方向に△ｙ/２だけ移動した位置に配置し、
前記タイプ１の格子点には、該格子点における前記平面と垂直なｚ方向の電界成分と、該電界成分の時間変化を計算するための誘電率とを割り当て、前記タイプ２の格子点には該格子点におけるｙ方向の磁界成分を割り当て、前記タイプ３の格子点には該格子点におけるｘ方向の磁界成分を割り当て、前記電界成分及び磁界成分の時間変化を時間領域差分法によって求める二次元電磁界シミュレーション方法であって、
前記境界に最も近い位置にある前記タイプ１の格子点を最近接格子点とし、該最近接格子点から前記境界までの距離を｜δ｜△ｘとし、前記最近接格子点が前記境界上または誘電率ε１の領域内にある場合に０≦δとし、前記最近接格子点が誘電率ε２の領域内にある場合にδ＜０とし、
誘電率ε１の領域にあって前記最近接格子点に最も近い位置にある前記タイプ１の格子点を第１近接格子点とし、誘電率ε２の領域にあって前記最近接格子点に最も近い位置にある前記タイプ１の格子点を第２近接格子点とし、
式、
(εｉ１+εｉ２+εｉ３)/3 = [(1+ 2δ/3)ε１ + (1- 2δ/3)ε２]/2、（１）
及び式、
12(εｉ３-εｉ１) = (13- 6δ２)(ε２-ε１)                   （２）
を満足するεｉ２、εｉ１及びεｉ３を、それぞれ、前記最近接格子点、前記第１近接格子点及び前記第２近接格子点に割り当て、前記最近接格子点と前記第１近接格子点と前記第２近接格子点とを除く前記タイプ１の格子点に、該格子点が誘電率ε１の領域にある場合にはε１を、該格子点が誘電率ε２の領域にある場合にはε２を、それぞれ割り当てる事を特徴とする二次元電磁界シミュレーション方法。
【請求項２】
１つの平面内に透磁率μ・誘電率ε１を有する第１の媒質と、透磁率μ・誘電率ε２を有する第２の媒質とが１つの直線を境界として接し、前記平面内に、前記平面に平行な交流電界成分を持たない電磁波が存在する場合に、前記平面と等価な仮想平面を前記平面と同一視し、
前記仮想平面内に、前記境界に平行なｙ軸と、前記境界に垂直なｘ軸と、前記ｘ軸とｙ軸との交点を通り前記ｘ軸とｙ軸とに垂直なｚ軸とを設け、ｉ、ｊを整数とし、△ｘ、△ｙを正の定数とし、前記仮想平面内に、ｘ座標がｉ△ｘであり、ｙ座標がｊ△ｙであるタイプ１の格子点(ｉ，ｊ）、ｘ座標が(ｉ＋１/２）△ｘであり、ｙ座標がｊ△ｙであるタイプ２の格子点（ｉ＋１/２，ｊ）、ｘ座標がｉ△ｘであり、ｙ座標が（ｊ＋１/２）△ｙであるタイプ３の格子点（ｉ，ｊ＋１/２）を配置し、
前記タイプ１の格子点（ｉ，ｊ）に該格子点におけるｚ方向の電界成分Ｅｚ(ｉ，ｊ）と、該電界成分Ｅｚ(ｉ，ｊ）の時間変化を計算する為の誘電率ε(ｉ，ｊ）とを割り当て、前記タイプ２の格子点（ｉ＋１/２，ｊ）に該格子点におけるｙ方向の磁界成分Ｈｙ(ｉ＋１/２，ｊ）を割り当て、前記タイプ３の格子点（ｉ，ｊ＋１/２）に該格子点におけるｘ方向の磁界成分Ｈｘ(ｉ，ｊ＋１/２）を割り当て、前記電界成分及び磁界成分の時間変化を時間領域差分法によって求める二次元電磁界シミュレーション方法であって、
前記境界に最も近い位置にある前記タイプ１の格子点を最近接格子点（ｉｍ ，ｊ）とし、前記最近接格子点（ｉｍ ，ｊ）から前記境界までの距離を｜δ｜△ｘとし、前記最近接格子点（ｉｍ ，ｊ）が前記境界上または誘電率ε１の領域内にある場合に０≦δとし、前記最近接格子点（ｉｍ ，ｊ）が誘電率ε２の領域内にある場合にδ＜０とし、
誘電率ε１の領域にあって前記最近接格子点（ｉｍ ，ｊ）に最も近い位置にある前記タイプ１の格子点を第１近接格子点（ｉ１ ，ｊ）とし、誘電率ε２の領域にあって前記最近接格子点（ｉｍ ，ｊ）に最も近い位置にある前記タイプ１の格子点を第２近接格子点（ｉ２，ｊ）とし、
式、
(εｉ１+εｉ２+εｉ３)/3 = [(1+ 2δ/3)ε１ + (1- 2δ/3)ε２]/2、（３）
及び式、
12(εｉ３-εｉ１) = (13- 6δ２)(ε２-ε１)                   （４）
を満足するεｉ１、εｉ２及びεｉ３が、それぞれ、ε(ｉ１ ，ｊ）、ε(ｉｍ ，ｊ）及びε(ｉ２，ｊ）に相等しいとし、前記最近接格子点（ｉｍ ，ｊ）と前記第１近接格子点（ｉ１ ，ｊ）と前記第２近接格子点（ｉ２，ｊ）とを除く前記タイプ１の格子点（ｉ，ｊ）に対して、該格子点（ｉ，ｊ）が誘電率ε１の領域にある場合にはε(ｉ，ｊ）＝ε１とし、該格子点（ｉ，ｊ）が誘電率ε２の領域にある場合にはε(ｉ，ｊ）＝ε２とする事を特徴とする二次元電磁界シミュレーション方法。
【請求項３】
前記電界成分と、該電界成分の時間変化を計算する為の誘電率と、前記ｙ方向の磁界成分と、前記ｘ方向の磁界成分とをコンピュータハードウェア資源であるメモリに数値として蓄え、該数値を引き出して計算に用いる事を特徴とする請求項１または請求項２に記載の二次元電磁界シミュレーション方法。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は電磁波の電磁界を計算に依って求める二次元電磁界シミュレーション方法に関し、特に、二次元ＴＥモード電磁界シミュレーション方法に関する。
【０００２】
【従来の技術】
近年、“社会の情報化”が著しく進展している。その中で、情報のキャリアとしての電磁波の利用は、ますます、高度化している。具体的には、光導波路やアンテナをそれぞれの目的に合わせ、精密に設計することが必要とされるようになってきている。そこで、複雑な環境の中で、電磁界の時間変化を詳細にシミュレーションすることが求められている。そのような目的のために使用される手法として、時間領域差分法（K. S. Yee,“Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media,”IEEE Transactions on Antennas and Propagation, vol. AP-14, pp. 302-307, 1966.）がある。この手法では、電磁界の時間変化を記述する基本方程式である Maxwell の方程式を、Yee 格子を用いて空間的に差分化し、電磁界の時間変化を計算する。この計算にはコンピュータが用いられる。
【０００３】
図１はコンピュータにより電磁界シミュレーションを行う装置主要部の機能ブロック図である。図においては、１はシミュレーション用基礎データと各種パラメタを入力する為の入力装置であり、３は電磁界シミュレーションを行う演算装置であり、２は入力装置１からの入力情報や演算装置３からの演算結果等を記憶する記憶装置であり、４は記憶装置２の記憶情報をシミュレーション結果として出力する為の出力装置である。図中、太い矢印はデータの流れを示している。シミュレーションする領域の誘電率分布、初期電磁界成分や計算に必要な時間ステップ等のパラメタが入力装置１から入力され、記憶装置２に記憶される。演算装置３において、記憶装置２に記憶されている電磁界成分のデータと誘電率分布のデータから次の時間ステップでの電磁界成分データが求められる。求められた電磁界成分データは再び記憶装置２に記憶され、更に次の時間ステップ電磁界成分を求める為に使用される。記憶装置２に記憶された電磁界成分データは、必要に応じて、出力装置４から出力させる事ができる。
【０００４】
図２に電磁界シミュレーションの手順の概要をフローチャートに纏めた。入力手順では、シミュレーション用基礎データと各種パラメタを入力装置から入力する。記憶・演算手順では、入力されたデータを記憶装置に記憶させ、記憶させたデータをもとに演算を行い電磁界を計算し、その結果を再び記憶装置に記憶させる。また、この記憶・演算手順では、計算した電磁界をもとに電磁界のモード強度や位相等のパラメタを計算する。また、総演算回数の制御を行う。出力手順では、記憶装置に記憶されたデータを出力する。出力手順は、入力手順や記憶・演算手順の開始後、随時、行っても構わない。
【０００５】
時間領域差分法では、電磁界の時間変化を追跡する為に必要な同一時刻の電磁界の空間偏微分を空間差分式を用いて近似する。空間差分式の結果をもとに、電磁界の時間変化を時間差分式を用いて計算する。時間領域差分法は、用いる差分式の精度次数により分類される。差分式の精度次数は、差分間隔（時間差分の場合は時間ステップ、空間差分の場合はセル幅、ここに、セルとは差分演算のための格子における単位格子のことである）をτとし、τを０に近づけた時に、誤差がτｎに比例して減少する場合、ｎ次精度と呼ばれる。空間差分の場合は、誘電率・透磁率が空間的に一定な状況、または、誘電率・透磁率の空間的変化が連続的な状況に於いて、この識別を適用する。
【０００６】
最も一般的に用いられている時間領域差分法は時間２次精度・空間２次精度の時間領域差分法である。この手法では、空間２次精度の差分式を用いる為、セル幅を半分にしても、空間差分近似の誤差は１/４にしかならない。
【０００７】
そこで、高精度の計算を行う場合、セル幅を大層小さくする必要が有り、コンピュータのメモリ使用量が膨大となった。従って、このメモリ使用量の増大を抑えるためには、この差分法における精度の次数を高める必要がある。
【０００８】
実際、誘電率が一様な場合においては、セル幅が同一の条件の下で、時間２次精度・空間４次精度の時間領域差分法は時間２次精度・空間２次精度の時間領域差分法より高精度であることが示されている（K. L. Shlager, J. G. Moloney, S. L. Ray, and A. F. Peterson, "Relative accuracy of several finite-difference time-domain method in two and three dimensions," IEEE Transactions on Antennas and Propagation, vol. 41, pp. 1732-1737, Dec. 1993.）。また、誘電率が一様な場合、セル幅が同一な条件のもとで、時間４次精度・空間４次精度の時間領域差分法は時間２次精度・空間２次精度の時間領域差分法の時間領域差分法より高精度であることが示されている（特開２０００−１０５２５９号公報「電磁界シミュレーション方法および当該シミュレーション方法の制御プログラムを記憶した記憶媒体」）。
【０００９】
【発明が解決しようとする課題】
しかしながら、空間４次精度の時間領域差分法においては、シミュレーション領域に異なる誘電率の媒質の境界が有る場合、即ち、誘電率が不連続に変化する境界が有る場合、境界付近ではどのような差分式を用いるべきか指針が得られていなかった。
【００１０】
本発明はこの点に鑑みてなされたものであり、誘電率が異なる２つの領域が１つの直線を境界として接している場合の二次元ＴＥモード電磁界を高い精度で計算することを可能とする二次元電磁界シミュレーション方法を提供する事を目的としている。
【００１１】
【課題を解決するための手段】
上記課題を解決するために、本発明においては、請求項１に記載のように、
１つの平面内に透磁率μ・誘電率ε１を有する第１の媒質と、透磁率μ・誘電率ε２を有する第２の媒質とが１つの直線を境界として接し、前記平面内に、前記平面に平行な交流電界成分を持たない電磁波が存在する場合に、前記平面と等価な仮想平面を前記平面と同一視し、
前記仮想平面内に、タイプ１の格子点を、前記境界に垂直なｘ方向及び平行なｙ方向に、それぞれ、間隔△ｘ及び△ｙで二次元配置し、タイプ２の格子点を、前記タイプ１の格子点の位置からｘ方向に△ｘ/２だけ移動した位置に配置し、タイプ３の格子点を、前記タイプ１の格子点の位置からｙ方向に△ｙ/２だけ移動した位置に配置し、
前記タイプ１の格子点には、該格子点における前記平面と垂直なｚ方向の電界成分と、該電界成分の時間変化を計算するための誘電率とを割り当て、前記タイプ２の格子点には該格子点におけるｙ方向の磁界成分を割り当て、前記タイプ３の格子点には該格子点におけるｘ方向の磁界成分を割り当て、前記電界成分及び磁界成分の時間変化を時間領域差分法によって求める二次元電磁界シミュレーション方法であって、
前記境界に最も近い位置にある前記タイプ１の格子点を最近接格子点とし、該最近接格子点から前記境界までの距離を｜δ｜△ｘとし、前記最近接格子点が前記境界上または誘電率ε１の領域内にある場合に０≦δとし、前記最近接格子点が誘電率ε２の領域内にある場合にδ＜０とし、
誘電率ε１の領域にあって前記最近接格子点に最も近い位置にある前記タイプ１の格子点を第１近接格子点とし、誘電率ε２の領域にあって前記最近接格子点に最も近い位置にある前記タイプ１の格子点を第２近接格子点とし、
式、
(εｉ１+εｉ２+εｉ３)/3 = [(1+ 2δ/3)ε１ + (1- 2δ/3)ε２]/2、（１）
及び式、
12(εｉ３-εｉ１) = (13- 6δ２)(ε２-ε１)                   （２）
を満足するεｉ２、εｉ１及びεｉ３を、それぞれ、前記最近接格子点、前記第１近接格子点及び前記第２近接格子点に割り当て、前記最近接格子点と前記第１近接格子点と前記第２近接格子点とを除く前記タイプ１の格子点に、該格子点が誘電率ε１の領域にある場合にはε１を、該格子点が誘電率ε２の領域にある場合にはε２を、それぞれ割り当てる事を特徴とする二次元電磁界シミュレーション方法を構成する。
【００１２】
また、本発明においては、請求項２に記載のように、
１つの平面内に透磁率μ・誘電率ε１を有する第１の媒質と、透磁率μ・誘電率ε２を有する第２の媒質とが１つの直線を境界として接し、前記平面内に、前記平面に平行な交流電界成分を持たない電磁波が存在する場合に、前記平面と等価な仮想平面を前記平面と同一視し、
前記仮想平面内に、前記境界に平行なｙ軸と、前記境界に垂直なｘ軸と、前記ｘ軸とｙ軸との交点を通り前記ｘ軸とｙ軸とに垂直なｚ軸とを設け、ｉ、ｊを整数とし、△ｘ、△ｙを正の定数とし、前記仮想平面内に、ｘ座標がｉ△ｘであり、ｙ座標がｊ△ｙであるタイプ１の格子点(ｉ，ｊ）、ｘ座標が(ｉ＋１/２）△ｘであり、ｙ座標がｊ△ｙであるタイプ２の格子点（ｉ＋１/２，ｊ）、ｘ座標がｉ△ｘであり、ｙ座標が（ｊ＋１/２）△ｙであるタイプ３の格子点（ｉ，ｊ＋１/２）を配置し、
前記タイプ１の格子点（ｉ，ｊ）に該格子点におけるｚ方向の電界成分Ｅｚ(ｉ，ｊ）と、該電界成分Ｅｚ(ｉ，ｊ）の時間変化を計算する為の誘電率ε(ｉ，ｊ）とを割り当て、前記タイプ２の格子点（ｉ＋１/２，ｊ）に該格子点におけるｙ方向の磁界成分Ｈｙ(ｉ＋１/２，ｊ）を割り当て、前記タイプ３の格子点（ｉ，ｊ＋１/２）に該格子点におけるｘ方向の磁界成分Ｈｘ(ｉ，ｊ＋１/２）を割り当て、前記電界成分及び磁界成分の時間変化を時間領域差分法によって求める二次元電磁界シミュレーション方法であって、
前記境界に最も近い位置にある前記タイプ１の格子点を最近接格子点（ｉｍ ，ｊ）とし、前記最近接格子点（ｉｍ ，ｊ）から前記境界までの距離を｜δ｜△ｘとし、前記最近接格子点（ｉｍ ，ｊ）が前記境界上または誘電率ε１の領域内にある場合に０≦δとし、前記最近接格子点（ｉｍ ，ｊ）が誘電率ε２の領域内にある場合にδ＜０とし、
誘電率ε１の領域にあって前記最近接格子点（ｉｍ ，ｊ）に最も近い位置にある前記タイプ１の格子点を第１近接格子点（ｉ１ ，ｊ）とし、誘電率ε２の領域にあって前記最近接格子点（ｉｍ ，ｊ）に最も近い位置にある前記タイプ１の格子点を第２近接格子点（ｉ２，ｊ）とし、
式、
(εｉ１+εｉ２+εｉ３)/3 = [(1+ 2δ/3)ε１ + (1- 2δ/3)ε２]/2、（３）
及び式、
12(εｉ３-εｉ１) = (13- 6δ２)(ε２-ε１)                   （４）
を満足するεｉ１、εｉ２及びεｉ３が、それぞれ、ε(ｉ１ ，ｊ）、ε(ｉｍ ，ｊ）及びε(ｉ２，ｊ）に相等しいとし、前記最近接格子点（ｉｍ ，ｊ）と前記第１近接格子点（ｉ１ ，ｊ）と前記第２近接格子点（ｉ２，ｊ）とを除く前記タイプ１の格子点（ｉ，ｊ）に対して、該格子点（ｉ，ｊ）が誘電率ε１の領域にある場合にはε(ｉ，ｊ）＝ε１とし、該格子点（ｉ，ｊ）が誘電率ε２の領域にある場合にはε(ｉ，ｊ）＝ε２とする事を特徴とする二次元電磁界シミュレーション方法を構成する。
【００１３】
また、本発明においては、請求項３に記載のように、
前記電界成分と、該電界成分の時間変化を計算する為の誘電率と、前記ｙ方向の磁界成分と、前記ｘ方向の磁界成分とをコンピュータハードウェア資源であるメモリに数値として蓄え、該数値を引き出して計算に用いる事を特徴とする請求項１または請求項２に記載の二次元電磁界シミュレーション方法を構成する。
【００１４】
【発明の実施の形態】
本発明は、以下に説明するように、時間２次精度・空間４次精度の時間領域差分法、または、時間４次精度・空間４次精度の時間領域差分法に関する。これらの手法では、空間４次精度の差分式を用いる為、セル幅を半分にすると空間差分近似の誤差は１/１６となり、セル幅をそれ程小さくせずとも高精度の計算が可能である。従って、コンピュータのメモリ使用量が少なくて済む。
【００１５】
本発明の特徴は、境界近傍で、電界成分の時間変化を計算するために、格子点に割り当てる誘電率の値を上記式（１）及び式（２）、或いは、式（３）及び式（４）（ただし、式（３）及び式（４）は、それぞれ、式（１）及び式（２）と同じである）によって定める事にある。境界近傍で、かかる誘電率の割り当てを行う事により、境界での反射を３次精度で近似できる為、本発明の実施によって、精度良いシミュレーションが可能となる。
【００１６】
本発明の実施の形態においては、１つの平面内に透磁率μ・誘電率ε１を有する第１の媒質と透磁率μ・誘電率ε２を有する第２の媒質とが１つの直線を境界として接し、前記平面内に、前記平面に垂直な交流磁界成分を持たない電磁波が存在する場合に、前記平面と等価な仮想平面を前記平面と同一視し、仮想平面内に格子点を前記境界に平行な方向及び垂直な方向に二次元配置して二次元格子を形成し、前記電磁波の電界成分、該電界成分の時間変化を計算する為の誘電率、及び、磁界の成分を前記格子点に割り当て、このように割り当てられた物理量を用いて、時間領域差分法に基づく数値計算によって、前記電磁波の電磁界の時間変化を求める。
【００１７】
（実施例１）
次に、本発明の１つの実施例を説明する。これは、時間２次・空間４次精度の時間領域差分法として実施されるものである。
【００１８】
図３は、本実施例における、上記の仮想平面上の格子点の配置を示す図である。図において、屈折率ε１の領域（請求項１及び請求項２における第１の媒質に対応）と屈折率ε２の領域（請求項１及び請求項２における第２の媒質に対応）とが境界（図中の縦太線）において接している。境界の左側の領域の誘電率をε１とし、境界の右側の領域の誘電率をε２としている。さらに、タイプ１の格子点とタイプ２の格子点とタイプ３の格子点とを、境界に平行な方向及び垂直な方向に沿って配置する。
【００１９】
この仮想平面内に、境界と平行なｙ軸とそれに垂直なｘ軸とを設け、ｉ，ｊを整数とし、△ｘ、△ｙを正の定数とすれば、タイプ１の格子点座標は（ｉ△ｘ，ｊ△ｙ）と表され、タイプ２の格子点座標は（（ｉ＋１/２)△ｘ，ｊ△ｙ）と表され、タイプ３の格子点座標は（ｉ△ｘ，（ｊ＋１/２)△ｙ）と表される。以後、各格子点は、そのｘ座標を△ｘで割り、そのｙ座標を△ｙで割った数の組合わせで、例えば、（ｉ，ｊ）や（ｉ＋１/２，ｊ）として表す事とする。
【００２０】
本実施例においては、タイプ１の格子点（ｉ，ｊ）に該格子点におけるｚ方向の電界成分Ｅｚ(ｉ，ｊ）と該電界成分の時間変化を計算するための誘電率ε(ｉ，ｊ）を割り当て、タイプ２の格子点（ｉ＋１/２，ｊ）には該格子点におけるｙ方向の磁界成分Ｈｙ(ｉ＋１/２，ｊ）を割り当て、タイプ３の格子点（ｉ，ｊ＋１/２）には該格子点におけるｘ方向の磁界成分Ｈｘ(ｉ，ｊ＋１/２）を割り当てる。
【００２１】
境界と格子点の位置関係について、境界に最も近いタイプ１の格子点を最近接格子点５とし、この格子点においてｉ＝０とする。この最近接格子点５から境界までの距離を｜δ｜△ｘとし、最近接格子点５が境界上または誘電率ε１の領域内にある場合に０≦δとし、その他の場合にδ＜０とする。このように決められたδは｜δ｜≦１/２を満足する。
【００２２】
なお、磁界のｚ方向成分Ｈｚと、電界のｘ方向成分Ｅｘ及びｙ方向成分Ｅｙとは常に０とする。
【００２３】
表１に各格子点に割り当てられる電磁界成分を纏めて示し、表２に本発明におけるタイプ１の格子点（ｉ，ｊ）に割り当てられる誘電率ε(ｉ，ｊ）を従来の手法の場合と対比して示した。
【００２４】
【表１】
0003827973.tif 000002
【００２５】
【表２】
0003827973.tif 000003
表２において、従来の手法は時間２次・空間２次精度の場合、ＴＥモードに対して用いられる手法である（Allen Taflove 著 "Computational Electrodynamics, The Finite-Difference Time-Domain Method." Artech House (Boston) 1995. 374頁)。
【００２６】
表２において、本発明におけるタイプ１の格子点（ｉ，ｊ）への誘電率の割り当ては、次のように行われている。すなわち、
格子点（０，ｊ）（最近接格子点５）、格子点（−１，ｊ）（これを第１近接格子点６とする）、格子点（１，ｊ）（これを第２近接格子点７とする）に、式、
(εｉ１+εｉ２+εｉ３)/3 = [(1+ 2δ/3)ε１ + (1- 2δ/3)ε２]/2、（１）
及び式、
12(εｉ３-εｉ１) = (13- 6δ２)(ε２-ε１)                   （２）
を満足するεｉ２、εｉ１及びεｉ３を、それぞれ、割り当る。すなわち、ε(０，ｊ）＝εｉ２、ε(−１，ｊ）＝εｉ１、ε(１，ｊ）＝εｉ３とする。これ以外の場合には、格子点（ｉ，ｊ）が誘電率ε１の領域にある場合、すなわち、ｉ≦−２の場合にはε(ｉ，ｊ）＝ε１とし、格子点（ｉ，ｊ）が誘電率ε２の領域にある場合、すなわち、１≦ｉの場合にはε(ｉ，ｊ）＝ε２とする。
【００２７】
以上のように、電磁界成分は空間的に異なった位置に割り当てられるが、時間的にも、前記仮想平面に時間の次元を加えてなる仮想３次元空間中の異なった位置に割り当てられる。例えば、電界が、離散化された時刻ｔ０、ｔ１、ｔ２、…ｔｎ−１、ｔｎ、ｔｎ＋１、…に割り当てられるとすると、磁界はｔ−１／２、ｔ１／２、ｔ１＋１／２、ｔ２＋１／２、…ｔｎ−１／２、ｔｎ＋１／２、ｔｎ＋３／２、…と言うように、電界成分が割り当てられた時刻の中間の時刻に割り当てられる。電磁界成分のデータは、格子点と離散化された時刻の組に対応するメモリ領域に記憶される。
【００２８】
本実施例で説明する時間２次精度の時間領域差分法では、記憶された電磁界成分から次の時刻の電磁界は次の演算をする事により求められる。
【００２９】
【数１】
0003827973.tif 000004
ここで、△ｔは時間ステップ間隔、Ｈｘｎ＋１／２(ｉ，ｊ＋１/２）は格子点（ｉ，ｊ＋１/２）・時刻ｔｎ＋１／２におけるｘ方向の磁界成分、Ｈｙｎ＋１／２(ｉ＋１/２，ｊ）は格子点（ｉ＋１/２，ｊ）・時刻ｔｎ＋１におけるｙ方向の磁界成分、Ｅｚｎ(ｉ，ｊ）は格子点（ｉ，ｊ）・時刻ｔｎ＋１におけるｚ方向の電界成分、ε(ｉ，ｊ）は格子点（ｉ，ｊ）に割り当てられた誘電率である。上記の式（５)〜(７）において、偏微分式を括弧[ ]で囲んだ項は、[ ]内に書かれた偏微分値を近似する差分式である。空間４次精度の時間領域差分法では、誘電率が一様な領域において、以下の差分式を用いる。
【００３０】
【数２】
0003827973.tif 000005
誘電率が一様な場合のように、各電磁界成分が４回偏微分可能として良い場合は、これらの差分式が四次精度を与える事を各電磁界成分の空間座標ｘ、ｙについてのテイラー展開を用いて示される。一方、誘電率境界の近傍では、各電磁界成分が境界と垂直な方向には、一般には、２回偏微分可能ではない。そこで、境界近傍で、式（９）や式（１０）を用いる事が可能かどうかは不明であった。
【００３１】
本発明では、境界近傍のタイプ１の格子点に割り当てる誘電率を表２に示したように設定すれば、境界近傍でも式（８）〜式（１１）を用いる事により、高精度の計算が可能である事を見出した。
【００３２】
演算後の結果は、再び図１の記憶装置２に記憶されるが、前の時刻の電界データを次の時刻の電界データで置き換えても、或いは、前の時刻の磁界データを次の時刻の磁界データで置き換えても、以後の電磁界の計算に支障は無い。それは、式（５）〜式（７）で示されるように、電界計算と磁界の計算が交互に行われているからである。
【００３３】
図４に本実施形態の電磁界の具体的計算手順をフローチャートとして記述した。  図４に示した計算を開始する前に、
第１ステップとして、上記変数Ｈｘｎ−１／２(ｉ，ｊ＋１/２）、Ｈｙｎ−１／２(ｉ＋１/２，ｊ）及びＥｚｎ(ｉ，ｊ）を記憶する領域をコンピュータのメモリ内に確保し、
第２ステップとして、上記変数ε(ｉ，ｊ）及び定数μを該メモリに記憶させ、
第３ステップとして、△ｔ及び最大時間ステップ数ｎｍａｘを該メモリに記憶させ、時間ステップ数ｎを記憶する領域を確保する。
【００３４】
計算を開始するに当たっては、
第４ステップとして、変数Ｈｘ−１／２(ｉ，ｊ＋１/２）、Ｈｙ−１／２(ｉ＋１/２，ｊ）及びＥｚ０(ｉ，ｊ）の設定を行い、
第５ステップとして、ｎを０に設定し、
第６ステップとして、式（５）と（６）とを用いて、時間ステップｎ＋１/２における磁界成分Ｈｘｎ＋１／２(ｉ，ｊ＋１/２）とＨｙｎ＋１／２(ｉ＋１/２，ｊ）とを計算して該メモリに格納し、
第７ステップとして、式（７）を用いて、時間ステップｎ＋１における電界成分Ｅｚｎ＋１(ｉ，ｊ）を計算して該メモリに格納し、
第８ステップとして、時間ステップｎ＋１/２における磁界成分Ｈｘｎ＋１／２(ｉ，ｊ＋１/２）とＨｙｎ＋１／２(ｉ＋１/２，ｊ）と、時間ステップｎ＋１における電界成分Ｅｚｎ＋１(ｉ，ｊ）とからデータ（例えば、磁界ベクトルの方向等）を抽出して該メモリに格納し、
第９ステップとして、ｎ＋１を改めてｎとし、
第１０ステップとして、ｎ＜ｎｍａｘの正否をコンピュータに判断させ、「正」の場合には、第６〜第１０ステップを順次実行させ、「否」の場合には、
第１１ステップとして、電磁界分布を出力すること、又は、電界強度・反射率等を計算して出力すること等を行い、
第１２ステップとして、計算プログラムを終了する。
【００３５】
（実施例２）
次に、本発明のもう１つの実施例を説明する。これは、時間４次精度・空間４次精度の時間領域差分法として実施されるものである。
【００３６】
時間４次精度の時間領域差分法の時間差分式は、特開２０００−１０５２５９号公報に記述されている。本実施例の場合に、誘電率ε１の領域、誘電率ε２の領域、両領域の境界に平行に置かれたｙ軸、ｙ軸と垂直なｘ軸、二次元格子点の設定、及び、二次元格子点に割り当てられる電界成分Ｅｚ(ｉ，ｊ）、誘電率ε(ｉ，ｊ）、磁界成分Ｈｙ(ｉ＋１/２，ｊ）及び磁界成分Ｈｘ(ｉ，ｊ＋１/２）の設定は実施例１と同一である。従って、図３と表１と表２とはそのまま使用される。また、ｚ方向（ｘ方向とｙ方向とに垂直な方向）の磁界成分Ｈｚ、ｘ方向の電界成分Ｅｘとｙ方向の電界成分Ｅｙは常に０とすることも実施例１と同様である。
【００３７】
一方、時間的な割り当ては、実施例１の時間２次精度の時間領域差分法と異なり、前記仮想平面に時間の次元を加えてなる仮想３次元空間中の同じ位置に割り当てられると考える。例えば、電界が間隔△ｔで離散化された時刻ｔ０、ｔ１、ｔ２、…ｔｎ−１、ｔｎ、ｔｎ＋１、…に割り当てられるとすると、磁界も時刻ｔ０、ｔ１、ｔ２、…ｔｎ−１、ｔｎ、ｔｎ＋１、…即ち、電界成分が割り当てられた時刻に割り当てられる。このような割り当ては、基本的な計算手続きであり、これを部分的に変形し、実施例１の時間２次精度の時間領域差分法のように、電界成分と磁界成分とを異なった時刻に割り当てるような定式化も可能である。
【００３８】
電磁界成分のデータは、格子点と離散化された時刻の組に対応するメモリ領域に記憶される。さらに、時間４次精度の時間領域差分法では、記憶された電磁界成分から次の時刻の電磁界を求める上で幾つかの中間状態を経る点が実施例１の時間２次精度の時間領域差分法と異なっている。この中間状態の電磁界成分は、右上に( )で囲んだ添え字を付して表す事にする。ここでは、ｔｎからｔｎ＋１の間に、電界が３つの中間状態Ｅ（１）、Ｅ（２）、Ｅ（３）、磁界が４つの中間状態Ｈ（１）、Ｈ（２）、Ｈ（３）、Ｈ（４）を経る場合を例にとって実施例を説明する。この電界が３つの中間状態、磁界が４つの中間状態を経由する場合の具体的手続きは、特開２０００−１０５２５９号公報の実施例には示されていないが、、特開２０００−１０５２５９号公報から容易に類推できる手法であり、その有効性は既に公表されている（廣野卓夫、ルイウェイン、関俊司、「シンプレクティック・プロパゲータを用いた三次元・四次精度ＦＤＴＤスキームの安定性と数値分散」、電気情報通信学会１９９９年総合大会講演論文集、p. 11、C-1-11）。具体的には、記憶された電磁界成分から次の時刻の電磁界は次の演算をする事により求められる。
【００３９】
【数３】
0003827973.tif 000006
【００４０】
【数４】
0003827973.tif 000007
ここで、Δｔは時間ステップ間隔、ｋを整数（１≦ｋ≦４）として、Ｈｘ（ｋ）(ｉ，ｊ＋１/２）は格子点（ｉ，ｊ＋１/２）のｋ番目の中間状態のｘ方向の磁界成分、、Ｈｙ（ｋ）(ｉ＋１/２，ｊ）は格子点（ｉ＋１/２，ｊ）のｋ番目の中間状態のｙ方向の磁界成分、Ｅｚ（ｋ）(ｉ，ｊ）は格子点（ｉ，ｊ）のｋ番目の中間状態のｚ方向の電界成分である。また、ここで、ｃ１、ｃ２、ｃ３、ｃ４、ｃ５、ｄ１、ｄ２、ｄ３、ｄ４は、ｃ１＋ｃ２＋ｃ＋ｃ４＋ｃ５＝ｄ１＋ｄ２＋ｄ３＋ｄ４＝１を満足する定数であり、本実施例においては、下記の値を用いた。
ｃ１＝ｃ５＝0.17399689146541、
ｃ２＝ｃ４＝-0.12038504121430、
ｃ３＝0.89277629949778、
ｄ１＝ｄ４＝0.62337932451322、
ｄ２＝ｄ３＝-0.12337932451322
時間４次精度の差分法での１つの時間ステップ（△ｔ）内の計算を整理する為、１つの時間ステップ内の計算は何段階かのステージからなっていると考える。パラメタｃ１とｄ１とを用いる式（１２)〜式(１４）の計算を第１ステージとする。以下、パラメタｃ２とｄ２とを用いる式（１５)〜式(１７）の計算を第２ステージとし、パラメタｃ３とｄ３とを用いる式（１８)〜式(２０）の計算を第３ステージとし、パラメタｃ４とｄ４とを用いる式（２１)〜式(２３）の計算を第４ステージとし、パラメタｃ５を用いる式（２４)〜式(２５）の計算を第５ステージとする。本実施例では、全ステージ数ｍｓは５である。上記の式（１２)〜式(２５）において、偏微分式を括弧[ ]で囲んだ項は、[ ]内に書かれた偏微分値を近似する差分式である。空間４次精度の時間領域差分法では、誘電率が一様な領域において、以下の差分式を用いる。
【００４１】
【数５】
0003827973.tif 000008
誘電率が一様な場合のように、各電磁界成分が４回偏微分可能として良い場合は、これらの差分式が四次精度を与える事を各電磁界成分の平面座標ｘ、ｙについてのテイラー展開を用いて示す事ができる。一方、誘電率が不連続的に変化する境界の近傍では、各電磁界成分が境界と垂直な方向には、一般には、前記境界を挟んでの２回偏微分が可能ではない。そこで、境界近傍で、式（２６)〜式(２９）を用いる事が可能かどうかは不明であった。
【００４２】
本発明においては、境界近傍のタイプ１の格子点に割り当てる誘電率を表３に示したように設定すれば、境界近傍においても式（２６)〜式(２９）を用いる事により、高精度の計算が可能である事を見出した。
【００４３】
式（１２)〜式(２９）で示される上記の演算は、図１の記憶装置２から電磁界データを読み出し、そのデータを用いて演算装置３によって行われる。その演算の結果は、再び図１の記憶装置２に記憶されるが、前の中間状態の電界データを次の中間状態の電界データで置き換えても、或いは、前の中間状態の磁界データを次の中間状態の磁界データで置き換えても、以後の電磁界の計算に支障は無い。それは、式（１２)〜式(２５）で示されるように、電界の計算と磁界の計算とが交互に行われているからである。
【００４４】
図５に本実施例の電磁界の具体的な計算の手順を示すフローチャートを記述した。
【００４５】
図５に示した計算を開始する前に、
第１ステップとして、上記変数Ｈｘｎ(ｉ，ｊ＋１/２）、Ｈｙｎ(ｉ＋１/２，ｊ）及びＥｚｎ(ｉ，ｊ）を記憶する領域をコンピュータのメモリ内に確保し、
第２ステップとして、変数ε(ｉ，ｊ）及び定数μを該メモリに記憶させ、
第３ステップとして、定数△ｔ及び最大時間ステップ数ｎｍａｘを該メモリに記憶させ、時間ステップ数ｎを記憶する該領域をメモリ内に確保する。
【００４６】
計算を開始するに当たっては、
第４ステップとして、変数Ｈｘ０(ｉ，ｊ＋１/２）、Ｈｙ０(ｉ＋１/２，ｊ）及びＥｚ０(ｉ，ｊ）の設定を行い、
第５ステップとして、ｎを０に設定し、
第６ステップとして、式（１２）と（１３）と（１４）とを用いて、Ｈｘ（１）(ｉ，ｊ＋１/２）とＨｙ（１）(ｉ＋１/２，ｊ）とＥｚ（１）(ｉ，ｊ）とを計算して該メモリに格納し、
第７ステップとして、式（１５）と（１６）と（１７）とを用いて、Ｈｘ（２）(ｉ，ｊ＋１/２）とＨｙ（２）(ｉ＋１/２，ｊ）とＥｚ（２）(ｉ，ｊ）とを計算して該メモリに格納し、
第８ステップとして、式（１８）と（１９）と（２０）とを用いて、Ｈｘ（３）(ｉ，ｊ＋１/２）とＨｙ（３）(ｉ＋１/２，ｊ）とＥｚ（３）(ｉ，ｊ）とを計算して該メモリに格納し、
第９ステップとして、式（２１）と（２２）と（２３）とを用いて、Ｈｘ（４）(ｉ，ｊ＋１/２）とＨｙ（４）(ｉ＋１/２，ｊ）とＥｚｎ＋１(ｉ，ｊ）とを計算して該メモリに格納し、
第１０ステップとして、式（２４）と（２５）とを用いて、Ｈｘｎ＋１(ｉ，ｊ＋１/２）とＨｙｎ＋１(ｉ＋１/２，ｊ）とを計算して該メモリに格納し、
第１１ステップとして、時間ステップｎ＋１における磁界成分Ｈｘｎ＋１(ｉ，ｊ＋１/２）とＨｙｎ＋１(ｉ＋１/２，ｊ）と、電界成分Ｅｚｎ＋１(ｉ，ｊ）とからデータ（例えば、磁界ベクトルの方向等）を抽出して該メモリに格納し、
第１２ステップとして、ｎ＋１を改めてｎとし、
第１３ステップとして、ｎ＜ｎｍａｘの正否をコンピュータに判断させ、「正」の場合には、第６〜第１３ステップを順次実行させ、「否」の場合には、
第１４ステップとして、電磁界分布を出力すること、又は、電界強度・反射率等を計算して出力すること等を行い、
第１５ステップとして、計算プログラムを終了する。
【００４７】
（参考例１）
半導体スラブ導波路の端面での反射率の計算に本発明に係る二次元電磁界シミュレーション方法を適用する。
【００４８】
図６は、本参考例におけるシミュレーション状況の模式図である。図において、屈折率3.4の半導体の媒質８中に幅0.3μｍ、屈折率3.61の導波路９が形成され、導波路９は光軸に垂直な端面１０で切断されている。自由空間波長0.9μｍのＴＥ基本モード電磁振動を端面１０の手前0.25μｍに設置した入射線（incident line）で励振し、端面に向かって入射した。計算領域境界にはＰＭＬ（perfectly matched layer）吸収領域を設定し、反射率を１０−１２％以下に抑えた。セルは正方形とし、屈折率が異なる媒質の境界は総てセル（すべてのタイプの格子点を含む最小単位格子）の辺方向と平行となるよう配置し、さらに、屈折率が異なる媒質の境界上に、その境界に平行な磁界成分の計算のための格子点が存在するようにした。即ち、本発明の適用条件が満足されるように、電磁界成分の計算のための格子点を設定した。
【００４９】
シミュレーション条件Ａを通常の時間２次精度・空間２次精度の時間領域差分法、シミュレーション条件Ｂを時間２次精度・空間４次精度の時間領域差分法に本発明の差分式を用いたもの、シミュレーション条件Ｃを時間４次精度・空間４次精度の時間領域差分法に本発明の差分式を用いたものとし、前記３条件について誤差の比較を行った。
【００５０】
境界近傍の誘電率の設定に就いては、本参考例で用いた境界とセルの位置関係の場合、通常の時間２次精度・空間２次精度の時間領域差分法では、誘電率の設定法は誘電率を設定すべき格子点について、その格子点が属する領域の誘電率を使用すれば良いことは既に知られている（Allen Taflove 著 "Computational Electrodynamics, The Finite-Difference Time-Domain Method." Artech House (Boston) 1995. 374頁)。
【００５１】
本参考例で用いた境界とセルの位置関係の場合、本発明の誘電率設定法を適用するとδ＝0.5 となる。さらに、εｉ１＝ε１、εｉ２＝（23ε１＋ε２)/24、εｉ３＝（ε１＋23ε２)/24 とすれば本発明の条件式（１）、条件式（２）（或いは条件式（３）、条件式（４））を満たす誘電率設定法となる。そこでシミュレーション手法ＢとＣではεｉ１＝ε１、εｉ２＝（23ε１＋ε２)/24、εｉ３＝（ε１＋23ε２)/24として、境界近傍の誘電率を設定した。
【００５２】
総ての手法とも、セル幅を0.05μｍ、0.025μｍ、0.0125μｍと減少して、メッシュを細かくした。反射率のセル幅 0μｍの場合への外挿値は一致して38.715％となり、この外挿値を真値とみなす事ができた。
【００５３】
図７に、各手法について、反射率誤差のセル幅依存性を示した。図７では、メッシュの細かさをコア（導波路９）における光の波長（0.9μｍ/3.61）当りのセル数で表した。図７では、セル数が増加する場合の誤差の減少が、シミュレーション条件Ａではほぼ２次精度（セル数を２倍にすると誤差が１/４になる）であるが、シミュレーション条件ＢとＣでは３次精度（セル数を２倍にすると誤差が１/８になる）と４次精度（セル数を２倍にすると誤差が１/１６になる）の間の傾きを示している。表示したどのセル数においても、シミュレーション条件ＢとＣの誤差はシミュレーション条件Ａより小さく、本発明の有効性を示している。
【００５４】
（参考例２）
半導体スラブ導波路の伝搬シミュレーションに本発明に係る二次元電磁界シミュレーション方法を適用する。
【００５５】
図８はシミュレーション状況の模式図である。解析領域の辺に平行にｘ軸とｙ軸を設け、電界成分Ｅｚ(ｉ，ｊ）と誘電率ε(ｉ，ｊ）を割り当てるタイプ１の格子点（ｉ，ｊ）、磁界成分Ｈｙ(ｉ＋１/２，ｊ）を割り当てるタイプ２の格子点（ｉ＋１/２，ｊ）、電界成分Ｈｘ(ｉ，ｊ＋１/２）を割り当てるタイプ３の格子点（ｉ，ｊ＋１/２）を設定する。
【００５６】
解析領域のバックグラウンドは空気（ただし、屈折率＝１とする）とし、ｘ軸方向に幅1.5μｍ、屈折率3.2の半導体導波路１１を置く。セルの大きさを変えると空気と半導体の境界と格子点の位置関係は変わってくる。２つの境界の内の一つでは、δ＝１/２となるように導波路と格子点の位置関係を定めた。他方の境界では、δ＝０、１/４、１/２とした。
【００５７】
伝搬する光の自由空間波長は1.55μｍとした。電界成分Ｅｚ、磁界成分Ｈｘ、磁界成分Ｈｙにより表される伝搬モードは二次元ＴＥモードに相当する。本導波路１１は、自由空間波長が1.55μｍの場合、６つのＴＥモード（０次〜５次）を保持できる。本参考例では、５次モードについて、伝搬シミュレーションを行った。解析領域のｘ方向の長さを、伝搬シミュレーションを行う５次モードの整数倍とし、解析領域の周囲は周期境界条件とする。
【００５８】
シミュレーション手法Ａを、両方の境界においてδ＝１/２となるようにし、通常の時間２次精度・空間２次精度の時間領域差分法を適用した場合とする。シミュレーション手法Ｂを、片方の境界でδ＝１/２もう一方の境界ではδ＝０とし、境界近傍では本発明の誘電率設定法を用いた場合とする。シミュレーション手法Ｃを、片方の境界でδ＝１/２もう一方の境界でδ＝１/４とし、境界近傍では本発明の誘電率設定法を用いた場合とする。シミュレーション手法Ｄを、片方の境界でδ＝１/２もう一方の境界でもδ＝１/２とし、境界近傍では本発明の誘電率設定法を用いた場合とする。
【００５９】
シミュレーション手法Ａは、時間２次精度・空間２次精度の時間領域差分法を用いて二次元ＴＥモードの伝搬シュミレーションを行う事であるが、δ＝１/２の場合には境界近傍の誘電率を設定すべき格子点では、その格子点が属する領域の誘電率を用いれば良いことが知られており（Allen Taflove 著 "Computational Electrodynamics, The Finite-Difference Time-Domain Method." Artech House (Boston) 1995. 374頁）、シミュレーション手法Ａでもその設定法を用いた。シミュレーション手法Ｂ、Ｃ、Ｄにおいて境界近傍に割り当てた誘電率を表３に示した。
【００６０】
【表３】
0003827973.tif 000009
表３のεｉ１、εｉ２、εｉ３は請求項に示した式（１）、式（２）（或いは、式（３）、式（４））を満足する。種々のセル幅についてシミュレーションを行い、位相速度誤差を比較した。２次元の導波路モードは解析的に計算できるので、位相速度誤差は厳密に評価できる。
【００６１】
各手法による計算の結果を図９に示す。図９では、メッシュの細かさをコア（半導体導波路１１）内での波長当りのセル数で表した。図９は本発明の手法を用いた場合の誤差が従来法の場合より遥かに小さい事を示している。また、図９において、セル数が増加する場合の誤差の減少が、シミュレーション条件Ａではほぼ２次精度（セル数を２倍にすると誤差が１/４になる）であるが、シミュレーション条件Ｂ、Ｃ、Ｄでは三次精度（セル数を２倍にすると誤差が１/８になる）以上の傾きを示している。また、図９は、セルと境界の位置関係がδ＝１/２以外の場合でも、本発明が有効である事を示している。
【００６２】
【発明の効果】
本発明の実施によって、誘電率が異なる２つの領域が１つの直線を境界として接している場合の二次元ＴＥモード電磁界を高い精度で計算することを可能とする二次元電磁界シミュレーション方法を提供する事ができる。
【図面の簡単な説明】
【図１】電磁界シミュレーション装置の機能ブロック図である。
【図２】電磁界シミュレーション手法の概要を示すフローチャートである。
【図３】本発明の時間領域差分法で用いられる格子点の説明図である。
【図４】実施例１のフローチャートである。
【図５】実施例２のフローチャートである。
【図６】参考例１で用いたシミュレーション領域の模式図である。
【図７】参考例１の結果を示す図である。
【図８】参考例２で用いたシミュレーション領域の模式図である。
【図９】参考例２の結果を示す図である。
【符号の説明】
１…入力装置、２…記憶装置、３…演算装置、４…出力装置、５…最近接格子点、６…第１近接格子点、７…第２近接格子点、８…屈折率３.４の媒質、９…屈折率３.６１の導波路、１０…導波路の端面、１１…半導体導波路。
【図１】
0003827973.tif 000010
【図２】
0003827973.tif 000011
【図３】
0003827973.tif 000012
【図４】
0003827973.tif 000013
【図５】
0003827973.tif 000014
【図６】
0003827973.tif 000015
【図７】
0003827973.tif 000016
【図８】
0003827973.tif 000017
【図９】
0003827973.tif 000018