(57)【特許請求の範囲】
【請求項１】
マックスウェル方程式を時間と空間とについて差分法で解くことにより電磁波の伝搬を解析する電磁波解析装置において、
時間離散間隔をΔｔとしたとき、第１の時刻における電磁波の状態に基づいて、前記第１の時刻からΔｔ経過後の第２の時刻における電磁波の状態を交互方向陰解法型の有限差分時間領域法により求めることで、電磁波の挙動を解析する解析手段、
を有することを特徴とする電磁波解析装置。
【請求項２】
前記解析手段は、２次元ＴＥ波の伝搬を解析することを特徴とする請求項１記載の電磁波解析装置。
【請求項３】
前記解析手段は、電界の振動面をｘ−ｙ平面とし、前記ｘ−ｙ平面の任意の格子点（ｉ，ｊ）を基準とした差分格子の電磁波を解析する際には、格子点（ｉ，ｊ）からｘ方向への空間離散間隔をΔｘ、格子点（ｉ，ｊ）からｙ方向への空間離散間隔をΔｙとした場合に、（ｉ＋１／２，ｊ）におけるｘ方向の電界、（ｉ，ｊ＋１／２）におけるｙ方向の電界、及び（ｉ＋１／２，ｊ＋１／２）における前記ｘ−ｙ平面に垂直なｚ方向の磁界を求めることを特徴とする請求項２記載の電磁波解析装置。
【請求項４】
前記解析手段は、ｘ方向の電界に対して陽解法を用い、ｙ方向の電界、及びｚ方向の磁界に対して陰解法を用いることにより、前記第１の時刻における電磁波の状態に基づいて前記第１の時刻からΔｔ／２経過後の中間時刻における電磁波の状態を演算する前半過程演算手段と、ｙ方向の電界に対して陽解法を用い、ｘ方向の電界、及びｚ方向の磁界に対して陰解法を用いることにより、前記前半過程演算手段により算出された前記中間時刻の電磁波の状態に基づいて前記第２の時刻における電磁波の状態を演算する後半過程演算手段と、を具備することを特徴とする請求項３記載の電磁波解析装置。
【請求項５】
前記前半過程演算手段は、シミュレーションにおけるｎステップ目の時刻を前記第１の時刻、ｎ＋１ステップ目の時刻を前記第２の時とし、前記第１の時刻における  （ｉ＋１／２，ｊ）でのｘ方向の電界をＥｘi+1/2,j n 、前記第１の時刻における（ｉ，ｊ＋１／２）でのｙ方向の電界をＥｙi,j+1/2 n 、前記第１の時刻における（ｉ＋１／２，ｊ＋１／２）でのｚ方向の磁界をＨｚi+1/2,j+1/2 n 、前記中間時刻における（ｉ＋１／２，ｊ）でのｘ方向の電界をＥｘi+1/2,j n+1/2 、前記中間時刻における（ｉ，ｊ＋１／２）でのｙ方向の電界をＥｙi,j+1/2 n+1/2 、前記中間時刻における（ｉ＋１／２，ｊ＋１／２）でのｚ方向の磁界をＨｚi+1/2,j+1/2 n+1/2 、解析すべき領域の誘電率をε、解析すべき領域の透磁率をμとした場合に、格子状に分割された前記ｘ−ｙ平面の任意の格子点（ｉ，ｊ）を基準とした電磁波を
【数１】
0003639116.tif 000080
【数２】
0003639116.tif 000081
【数３】
0003639116.tif 000082
と定義することで、前記中間時刻におけるすべての格子点を基準とした電磁波の状態を演算し、
前記後半過程演算手段は、前記第２の時刻における（ｉ＋１／２，ｊ）でのｘ方向の電界をＥｘi+1/2,j n+1 、前記第２の時刻における（ｉ，ｊ＋１／２）でのｙ方向の電界をＥｙi,j+1/2 n+1 、前記第２の時刻における（ｉ＋１／２，ｊ＋１／２）でのｚ方向の磁界をＨｚi+1/2,j+1/2 n+1 とした場合に、格子状に分割された前記ｘ−ｙ平面の任意の格子点（ｉ，ｊ）を基準とした電磁波を
【数４】
0003639116.tif 000083
【数５】
0003639116.tif 000084
【数６】
0003639116.tif 000085
と定義することで、前記第２の時刻におけるすべての格子点を基準とした電磁波の状態を演算することを特徴とする請求項４記載の電磁波解析装置。
【請求項６】
マックスウェル方程式を時間と空間とについて差分法で解くことにより電磁波の伝搬を解析するための電磁波解析プログラムを記録したコンピュータ読み取り可能な記録媒体において、
時間離散間隔をΔｔとしたとき、第１の時刻における電磁波の状態に基づいて、前記第１の時刻からΔｔ経過後の第２の時刻における電磁波の状態を交互方向陰解法型の有限差分時間領域法により求めることで、電磁波の挙動を解析する解析手段、
としてコンピュータを機能させることを特徴とする電磁波解析プログラムを記録したコンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は電磁波の挙動を電子計算機を使って解析する電磁波解析装置及び電磁波解析プログラムを記録したコンピュータ読み取り可能な記録媒体に関し、特に有限差分時間領域法を用いて解析を行う電磁波解析装置及び電磁波解析プログラムを記録したコンピュータ読み取り可能な記録媒体に関する。
【０００２】
【従来の技術】
電磁波の過渡的な挙動を、電子計算機を使った数値シミュレーションによって解析する方法の１つに有限差分時間領域（ＦＤＴＤ：Finite Difference Time Domain ）法がある。この方法は、マックスウェル方程式を時間と空間とについて差分法で解く手法であり、その適用範囲の広さ等から広く用いられている手法である。
【０００３】
ところで、従来の有限差分時間領域法は陽解法(explicit method) 型差分法で解かれていた。以下に、陽解法について説明する。例えば、簡単のために以下のような１次元放物型偏微分方程式を考える。
【０００４】
【数７】
0003639116.tif 000002
【０００５】
【数８】
0003639116.tif 000003
【０００６】
ここで、ｎは演算開始時点から時間軸方向に進めたステップの数であり、Δｔは時間離散間隔（１ステップで進められる時間）であり、ｉは差分格子（セル）の原点からの順番であり、Δｘは空間離散間隔（分割領域の幅）である。この微分方程式を陽解法で解く場合、時間微分に対して前進差分で近似を行い、
【０００７】
【数９】
0003639116.tif 000004
【０００８】
という式を得る。この式（９）はｒ＝Δｔ／Δｘ2 とおいて
【０００９】
【数１０】
0003639116.tif 000005
【００１０】
と書くことができる。これにより、ｎ＋１番目のステップにおける解が、過去の時刻（ｎ番目のステップ）の解から代入計算によりただちに計算できる。
【００１１】
【発明が解決しようとする課題】
しかし、数（１０）の陽解法型差分法では、解が発散しないための安定条件として
【００１２】
【数１１】
0003639116.tif 000006
【００１３】
を満たしていなければならない。
また、ＦＤＴＤ法においては、解が発散しないためにはＣＦＬ(Courant, Friedrich and Levy) 条件あるいはCourant 条件を満たす必要がある。
【００１４】
ＣＦＬ条件は、解が発散しないためには、
【００１５】
【数１２】
0003639116.tif 000007
【００１６】
の安定条件を満たしていなければならないという条件である。ＦＤＴＤ法のＣＦＬ条件については、「A.Taflove,Computational Electrodynamics,MA,Artech,1995」に詳しい。なお、式（１２）は、３次元波を解析する場合のCourant 条件である。ここで、ｖは、電磁波の速度であり、Δｘmin 、Δｙmin 、Δｚmin はそれぞれＸＹＺ方向の空間離散間隔の最小値である。
【００１７】
これらの制約があることにより、陽解法においては、時間離散間隔の最大値が空間離散間隔の最小値によって決定されるという性質がある。このため、微細な構造を有する解析対象の時間応答を長時間計算する場合、シミュレーションの時刻を進行させるための計算の繰り返し回数が増加し、シミュレーションに要する時間が増加するといった問題があった。
【００１８】
なお、式（７）、式（８）に示した偏微分方程式は、陰解法(implicit method) で解くこともできる。例えば、式（７）、式（８）のような１次元放物型偏微分方程式を陰解法で解く場合には、時間微分について後退差分で近似を行い、
【００１９】
【数１３】
0003639116.tif 000008
【００２０】
を得る。この式（１３）は、陽解法の場合と同様にｒ＝Δｔ／Δｘ2 とおいて
【００２１】
【数１４】
0003639116.tif 000009
【００２２】
と書くことができる。このような式（１４）を解くのが陰解法である。この陰解法は、陽解法と異なり常に安定である。
ただし、陰解法は、以下のような連立方程式を解かなくてはｆi n+1 の値を得ることができない。
【００２３】
【数１５】
0003639116.tif 000010
【００２４】
ここで、ａi は式（１４）において「ｉ＝１，２，・・・imax」とした時の左辺第３項の定数部分（−ｒ）であり、ｂi は式（１４）において「ｉ＝１，２，・・・imax」とした時の左辺第２項の定数部分（１＋２ｒ）であり、ｃi は式（１４）において「ｉ＝１，２，・・・imax」とした時の左辺第１項の定数部分（−ｒ）である。
【００２５】
電磁波解析の過渡的な現象を解析するには、少なくとも２次元以上の偏微分方程式を解く必要があるが、２次元以上の偏微分方程式の陰解法については別途注意が必要となる。例えば２次元放物型偏微分方程式は、以下の式で表される。
【００２６】
【数１６】
0003639116.tif 000011
【００２７】
この式（１６）を、例えばCrank-Nicolson法で解く場合には、
【００２８】
【数１７】
0003639116.tif 000012
【００２９】
のように時間微分を近似する。この式（１７）を連立方程式によって解く場合、（ｘ方向の格子数−１）×（ｙ方向の格子数−１）元の連立１次方程式を解く必要がある。そのため、格子数が多い場合には計算量や必要な記憶容量が膨大になる。すなわち、格子を微細化すると、陽解法で解析した場合と同様に非常に時間のかかる計算となる。
【００３０】
本発明はこのような点に鑑みてなされたものであり、微細な構造を有する解析対象物の電磁波解析を高速に行うことのできる電磁波解析装置を提供することを目的とする。
【００３１】
また、本発明の他の目的は、微細な構造を有する解析対象物の電磁波解析をコンピュータに高速に行わせることのできる電磁波解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提供することである。
【００３２】
【課題を解決するための手段】
本発明では上記課題を解決するために、マックスウェル方程式を時間と空間とについて差分法で解くことにより電磁波の伝搬を解析する電磁波解析装置において、時間離散間隔をΔｔとしたとき、第１の時刻における電磁波の状態に基づいて、前記第１の時刻からΔｔ経過後の第２の時刻における電磁波の状態を交互方向陰解法型の有限差分時間領域法により求めることで、電磁波の挙動を解析する解析手段、を有することを特徴とする電磁波解析装置が提供される。
【００３３】
このような電磁波解析装置によれば、解析手段により、第１の時刻における電磁波の状態に基づいて、第１の時刻からΔｔ経過後の第２の時刻における電磁波の状態が交互方向陰解法型の有限差分時間領域法により求められる。このように交互方向陰解法型の有限差分時間領域法によって得られた解は、時間離散間隔の大小によらず安定している。
【００３４】
また、上記課題を解決するために、マックスウェル方程式を時間と空間とについて差分法で解くことにより電磁波の伝搬を解析するための電磁波解析プログラムを記録したコンピュータ読み取り可能な記録媒体において、時間離散間隔をΔｔとしたとき、第１の時刻における電磁波の状態に基づいて、前記第１の時刻からΔｔ経過後の第２の時刻における電磁波の状態を交互方向陰解法型の有限差分時間領域法により求めることで、電磁波の挙動を解析する解析手段、としてコンピュータを機能させることを特徴とする電磁波解析プログラムを記録したコンピュータ読み取り可能な記録媒体が提供される。
【００３５】
この記録媒体に記録された電磁波解析プログラムをコンピュータに実行させれば、上記本発明に係る電磁波解析装置の機能がコンピュータ上で実現される。
【００３６】
【発明の実施の形態】
以下、本発明の実施の形態を図面を参照して説明する。本発明の電磁波解析装置は、ＡＤＩ(Alternating Direction Implicit)法型のＦＤＴＤ法を用いて電磁波の解析を行うことにより、解の安定性を確保しつつ高速の演算を可能としたものである。
【００３７】
図１は、本発明の原理構成図である。解析指令入力手段１は、電磁波解析に必要な初期値や境界吸収条件を指定して、解析指令を入力する。
解析手段２は、解析指令入力手段１からの解析指令を受け取ると、２次元ＴＥ(Transverse Electric) 波の伝搬を、交互方向陰解法型−差分時間領域（ＡＤＩ−ＦＤＴＤ）法によって数値解析する。ここで、ＴＥ波とは以下の式で表現される電磁波である。
【００３８】
【数１８】
0003639116.tif 000013
【００３９】
【数１９】
0003639116.tif 000014
【００４０】
【数２０】
0003639116.tif 000015
【００４１】
また、ＡＤＩ−ＦＤＴＤ法は、時間離散間隔を前半と後半とに分割し、前半過程の演算と後半過程の演算とを交互に行う方法である。例えば、式（１６）の２次元放物型偏微分方程式をＡＤＩ法で解く場合には、時間微分の値を、
【００４２】
【数２１】
0003639116.tif 000016
【００４３】
【数２２】
0003639116.tif 000017
【００４４】
のように近似する。式（２１）が前半過程であり、左辺の時間微分は、ｘ方向のｆから見ると後退差分で近似し、ｙ方向のｆから見ると前進差分で近似している。式（２２）が後半過程であり、左辺の時間微分は、ｘ方向ｆから見ると前進差分で近似し、ｙ方向のｆから見ると後退差分で近似している。前半過程は前半過程演算手段２ａが演算し、後半過程は後半過程演算手段２ｂが演算する。なお、Δｔは時間離散間隔であり、Δｘはｘ方向の空間離散間隔であり、Δｙはｙ方向の空間離散間隔である。
【００４５】
前半過程演算手段２ａでは、後半過程演算手段２ｂから時刻「ｎΔｔ」（ｎはシミュレーション開始時からの時刻進行処理のステップ数であり、前半過程と後半過程ではそれぞれ１／２ステップずつ進められる）における電磁波の状態の情報を受け取ると、時刻「（ｎ＋１／２）Δｔ」の電磁波の状態を計算する。その際、時刻「（ｎ＋１／２）Δｔ」におけるｘ方向の電界Ｅｘn+1/2 を陽解法により求め、ｙ方向の電界Ｅｙn+1/2 とｚ方向の磁界Ｈｚn+1/2 とを陰解法により求める。計算結果は、後半過程演算手段２ｂに引き渡す。
【００４６】
後半過程演算手段２ｂは、前半過程演算手段２ａから時刻「（ｎ＋１／２）Δｔ」における電磁波の状態を受け取ると、時刻「（ｎ＋１）Δｔ」の電磁波の状態を計算する。その際、時刻「（ｎ＋１）Δｔ」におけるｙ方向の電界Ｅｙn+1 を陽解法により求め、ｘ方向の電界Ｅｘn+1 とｚ方向の磁界Ｈｚn+1 とを陰解法により求める。計算結果は、前半過程演算手段２ａに引き渡す。
【００４７】
解析結果出力手段３は、解析手段２による解析結果を出力する。
以下に、時刻「ｎΔｔ」の電磁波の状態から時刻「（ｎ＋１）Δｔ」の状態を求める際の解析手段２の計算内容について説明する。
【００４８】
まず、解析対象領域を長方形の格子に分割する。ここで、解析対象の領域のｘ方向をｉmax 個の格子に分割し、ｙ方向をｊmax 個の格子に分割したものとする。
【００４９】
図２は、格子の拡大図である。格子点はの位置は（ｉ，ｊ）のような２組の数字で表される。ｉは、原点からのｘ方向の格子の順番である。ｊは、原点からのｙ方向の格子の順番である。また、「ｉ」と「ｉ＋１」との中間を「ｉ＋１／２」と表す。ここで、格子点（ｉ，ｊ）を基準とした領域の２次元ＴＥ波の解析を行う場合、（ｉ＋１／２，ｊ）におけるｘ方向の電界、（ｉ，ｊ＋１／２）におけるｙ方向の電界、及び（ｉ＋１／２，ｊ＋１／２）におけるｚ方向の磁界の時間変化を差分法で解析する。
【００５０】
マックスウェル方程式を時間と空間について差分法で解く場合、前半過程ではＥx の時間微分を前進差分で近似し（以下の式（２３）に示す）、Ｅy の時間微分を後退差分で近似し（以下の式（２４）に示す）、Ｈz の時間微分はＥx から見れば前進差分で近似するとともに、Ｅy から見れば後退差分で近似する（以下の式（２５）に示す）。その結果、時刻「（ｎ＋１／２）Δｔ」の状態を求める前半過程は、以下の式で表せる。
【００５１】
【数２３】
0003639116.tif 000018
【００５２】
【数２４】
0003639116.tif 000019
【００５３】
【数２５】
0003639116.tif 000020
【００５４】
ここで、電界Ｅや磁界Ｈの右肩の値は、シミュレーションの時刻ステップを表している。電界Ｅや磁界Ｈの右下の値は、解析領域を格子で分割した際の格子番号である。εは、解析対象物の誘電率である。μは、解析対象物の透磁率である。
【００５５】
ところで、式（２３）は陽解法的表現となっており、式（２４）と式（２５）とは陰解法的表現となっている。つまり前半過程では、Ｅx の時間微分を前進差分で近似しているため、時刻「（ｎ＋１／２）Δｔ」におけるｘ方向の電界Ｅx n+1/2 は陽解法により求めることができる。また、前半過程では、Ｅy の時間微分を後退差分で近似しているため、時刻「（ｎ＋１／２）Δｔ」におけるｙ方向の電界Ｅｙn+1/2 は陰解法により求めなければならない。さらに、Ｈz の時間微分をＥy から見ると後退差分で近似しているので、時刻「（ｎ＋１／２）Δｔ」におけるｚ方向の磁界Ｈｚn+1/2 は陰解法により求めなければならない。
【００５６】
同様に、時刻「（ｎ＋１）Δｔ」の状態を求める後半過程は、Ｅx の時間微分を後退差分で近似し（以下の式（２６）に示す）、Ｅy の時間微分を前進差分で近似し（以下の式（２７）に示す）、Ｈz の時間微分はＥx から見れば後退差分で近似するとともに、Ｅy から見れば前進差分で近似する（以下の式（２８）に示す）。これにより、後半過程は以下の式で表せる。
【００５７】
【数２６】
0003639116.tif 000021
【００５８】
【数２７】
0003639116.tif 000022
【００５９】
【数２８】
0003639116.tif 000023
【００６０】
後半過程ではｙ方向の電界を求める式（２７）が陽解法的表現となり、その他の式（２６）、式（２８）が陰解法的表現となっている。つまり後半過程では、Ｅx の時間微分を後退差分で近似しているため、時刻「（ｎ＋１）Δｔ」におけるｘ方向の電界Ｅx n+1 は陰解法により求めなければならない。また、Ｅy の時間微分を前進差分で近似しているため、時刻「（ｎ＋１）Δｔ」におけるｙ方向の電界Ｅｙn+1 は陽解法により求めることができる。さらに、Ｈz の時間微分をＥx から見ると後退差分で近似しているので、時刻「（ｎ＋１）Δｔ」におけるｚ方向の磁界Ｈｚn+1 は陰解法により求めなければならない。
【００６１】
前半過程の式（２３）により、時刻「（ｎ＋１／２）Δｔ」のｘ方向の電界は、時刻「ｎΔｔ」の電界及び磁界の代入計算によって求めることができる。同様に、後半過程の式（２７）も簡単な代入計算で求めることができる。
【００６２】
一方、前半過程の式（２４）、式（２５）は、右辺に時刻「（ｎ＋１／２）Δｔ」の電界や磁界が含まれているため、そのままでは解けない。そこで、式（２４）と式（２５）とを連立させて磁界Ｈｚn+1/2 を消去し、下の漸化式を得る。
【００６３】
【数２９】
0003639116.tif 000024
【００６４】
後半過程の式（２６）、式（２８）も、右辺に時刻「（ｎ＋１）Δｔ」の電界と磁界が含まれるため、式（２６）と式（２８）とを連立させて磁界Ｈｚn+1 を消去し、以下の漸化式を得る。
【００６５】
【数３０】
0003639116.tif 000025
【００６６】
式（２９）の漸化式を「ｉ＝１，・・・，ｉmax 」とした場合のそれぞれについて生成すれば、以下のような連立一次方程式に帰着する。
【００６７】
【数３１】
0003639116.tif 000026
【００６８】
ここで式（３１）におけるａi は、「ｉ＝１，・・・，ｉmax 」とした場合に得られる式（２９）の左辺第１項の定数部分である。ｂi は、「ｉ＝１，・・・，ｉmax 」とした場合に得られる式（２９）の左辺第２項の定数部分である。ｃi は、「ｉ＝１，・・・，ｉmax 」とした場合に得られる式（２９）の左辺第３項の定数部分である。Ｄi n は、「ｉ＝１，・・・，ｉmax 」とした場合に得られる式（２９）の右辺の値である。また、第１行と最終行には、吸収境界条件が適用される。
【００６９】
式（３１）は、ｘ方向の格子点「ｉ」についての３項方程式(tridiagonal system of equation)であり、トーマス法等によって解くことができる。このような３項方程式を「ｉ＝１，・・・，ｉmax 」とした場合のそれぞれについて解けば、時刻「ｎ＋１／２」における全ての格子でのｙ方向の電界が求められる。ｙ方向の電界が求まれば、式（１６）によってｚ方向の磁界も求まる。
【００７０】
後半過程に関する漸化式である式（３０）も前半過程と同様に、以下のような「ｊ」について３項方程式を得ることができる。
【００７１】
【数３２】
0003639116.tif 000027
【００７２】
ここで式（３２）におけるａj は、「ｊ＝１，・・・，ｊmax 」とした場合に得られる式（３０）の左辺第１項の定数部分である。ｂj は、「ｊ＝１，・・・，ｊmax 」とした場合に得られる式（３０）の左辺第２項の定数部分である。ｃj は、「ｊ＝１，・・・，ｊmax 」とした場合に得られる式（３０）の左辺第３項の定数部分である。Ｄj n+1/2 は、「ｊ＝１，・・・，ｊmax 」とした場合に得られる式（３０）の右辺の値である。また、第１行と最終行には、吸収境界条件が適用される。
【００７３】
このように、前半過程の式（２４）と後半過程の式（２６）とのそれぞれの演算を３項方程式に帰着することができる。従って、前半過程演算手段２ａにおいては式（２３）と式（３１）の３項方程式とを解き、さらに式（２５）を解くことになる。後半過程演算手段２ｂにおいては式（２７）と式（３２）の３項方程式とを解き、それから式（２８）を解くことになる。このような演算を交互に行うことにより時刻のステップを進めていけば、２次元ＴＥ波の伝搬をＡＤＩ−ＦＤＴＤ法により解析することができる。
【００７４】
図３は、ＡＤＩ−ＦＤＴＤ法による２次元ＴＥ波解析の処理手順を示すフローチャートである。この処理をステップ番号に沿って説明する。なお、吸収境界条件、時間離散間隔、空間離散間隔等の各種初期条件は、すでに与えられているものとする。シミュレーション開始時点では、シミュレーション時刻「Ｔ＝０」である。
［Ｓ１］前半過程演算手段２ａは、まず、すでに求められているＥx n 、Ｈz n の値を代入して式（２３）を解くことで、Ｅx n+1/2 を計算する。次に、すでに求められているＥx n 、Ｅy n 、Ｈz n の値を用いて式（３１）の３項方程式を解くことで、Ｅy n+1/2 を計算する。さらに、式（３１）を解くことによって得られたＥy n+1/2 の値と既知のＥx n 、Ｈz n の値とを用いて式（２５）を解くことで、Ｈz n+1/2 を計算する。これにより、シミュレーションの時刻Ｔは、ΔＴ／２だけ進められる。
［Ｓ２］後半過程演算手段２ｂは、まず、すでに求められているＥx n+1/2 、Ｅy 1+1/2 、Ｈz n+1/2 の値を用いて式（３２）の３項方程式を解くことで、Ｅx n+1 を計算する。次に、すでに求められているＥy n+1/2 、Ｈz n+1/2 の値を用いて式（２７）を解くことで、Ｅy n+1 を計算する。さらに、式（３２）を解くことによって得られたＥx n+1 の値と既知のＥy n+1/2 、Ｈz n+1/2 の値とを用いて式（２８）を解くことで、Ｈz n+1 を計算する。これにより、シミュレーションの時刻Ｔは、ΔＴ／２だけ進められる。その結果、ステップＳ１とステップＳ２とのトータルで時間離散間隔であるΔｔだけ時刻が進められたことになる。
［Ｓ３］解析手段２は、時刻Ｔがあらかじめ設定されているシミュレーション最大時刻Ｔmax を超えているか否かを判断する。Ｔmax を超えていればシミュレーションを終了し、超えていなければステップＳ１に進みシミュレーションを継続する。
【００７５】
以上のような処理を解析手段２が実行することにより、ＡＤＩ−ＦＤＴＤ法を用いた２次元ＴＥ波の解析が行われる。
なお、以上の式は簡略化のために、電磁波を、無損失、等方性、非分散性、均質の媒質を伝搬するものと仮定し、さらに、差分格子（セル）は一様（ΔｔとΔｙがすべてのセルで等しい）であるものと仮定しているため、「ａi 」、「ｃi 」、「ａj 」、「ｃj 」の値はすべて「１」である。媒質が不均質または差分格子（セル）が一様でない場合は、式の拡張を行う必要がある。前述の式を媒質が不均質または差分格子（セル）が一様でない場合へ拡張する方法は、通常のＦＤＴＤ法で行われていた拡張方法と同じである。
【００７６】
以下に、２次元ＴＥ波の伝搬をＡＤＩ−ＦＤＴＤ法により解析した場合の解の安定性についてフーリエ級数法を用いて検討した結果を示す。
まず、時間発展に伴う誤差の成長係数ξは、時刻のステップ数ｎ’を用いて
【００７７】
【数３３】
0003639116.tif 000028
【００７８】
と表される。なお、γは複素定数である。このような誤差の成長係数を用いて、前半過程および後半過程の解をそれぞれ次のようにおく。
【００７９】
【数３４】
0003639116.tif 000029
【００８０】
【数３５】
0003639116.tif 000030
【００８１】
ここで、ξ1 n'及びξ2 n'は、それぞれ前半過程及び後半過程での時間発展にともなう誤差の成長係数である。また、虚数単位は添字とまぎらわしいため、−１の平方根で表している。すると、これらは次式で表される。
【００８２】
【数３６】
0003639116.tif 000031
【００８３】
【数３７】
0003639116.tif 000032
【００８４】
ここで、ｐ、ｑは、それぞれ
【００８５】
【数３８】
0003639116.tif 000033
【００８６】
【数３９】
0003639116.tif 000034
【００８７】
である。式（３６）、式（３７）から、全過程での成長係数ξは常に
【００８８】
【数４０】
0003639116.tif 000035
【００８９】
となる。これは、いかなる場合も安定で、しかもエネルギー散逸がないことを示している。
次に、本発明の電磁波解析装置のハードウェア構成について説明する。
【００９０】
図４は、電磁波解析装置のハードウェア構成を示す図である。この例では、コンピュータに電磁波解析プログラムを実行させることで、電磁波解析装置を実現している。
【００９１】
この電磁波解析装置は、ＣＰＵ(Central Processing Unit) １１を中心に構成されている。ＣＰＵ１１は、メモリ１２に記憶されたプログラムに基づいて電磁波解析の演算を行うとともに、バス１７を介して接続された各種機器を制御する。バス１７に接続された周辺機器は、次のようなものである。
【００９２】
ディスプレイコントローラ１３は、ＣＰＵ１１から送られてくる描画命令に従って表示画像を生成し、表示装置１３ａに送る。ディスプレイコントローラ１３に接続された表示装置１３ａは、ディスプレイコントローラ１３から送られた表示画像情報に従い、その画像を画面に表示する。
【００９３】
入力機器インタフェース１４は、キーボード１４ａやマウス１４ｂが接続されており、キーボード１４ａやマウス１４ｂからの入力信号をＣＰＵ１１へ転送する。
【００９４】
ネットワークインタフェース１５は、ＬＡＮ(Local Area Network)に接続されており、ＬＡＮを介したデータ通信を制御する。すなわち、ＣＰＵ１１から送られたデータをＬＡＮ上の他の装置へ転送するとともに、ＬＡＮを介して送られてきたデータを受け取りＣＰＵ１１に渡す。
【００９５】
ＨＤＤ(Hard Disk Drive) コントローラ１６は、接続されたハードディスク装置（ＨＤＤ）１６ａへのデータの入出力を制御する。ＨＤＤ１６ａには、ＣＰＵ１１が実行すべきシステムプログラムや、電磁波解析処理を実行する際の処理内容が記述されたプログラムが格納されている。
【００９６】
このような構成の電磁波解析装置において、電磁波解析解析を行う場合には、まずユーザがキーボード１４ａを用いて初期条件を入力する。初期条件は、ｘ方向とｙ方向との空間離散間隔、時間離散間隔、境界条件、演算領域の範囲指定等である。入力された初期条件は、入力機器インタフェース１４からＣＰＵ１１に渡され、メモリ１２に保持される。そして、ユーザが電磁波解析処理の実行指令を入力すると、その指令がＣＰＵ１１に送られる。ＣＰＵ１１は、本発明に係るＦＤＴＤ法による電磁波解析処理を行うためのプログラムをＨＤＤ１６ａからメモリ１２にロードし、そのプログラムを実行する。実行結果は、ディスプレイコントローラ１３を介して表示装置１３ａの画面に表示される。
【００９７】
以下に、具体的な数値計算例を示す。なお、以下の計算は、媒質が不均質な場合でも、または差分格子（セル）が一様でない場合でも計算できるようなプログラムにより行った。このような機能拡張の方法は、従来のＦＤＴＤ法で行われている方法と同じである。
【００９８】
図５は、第１のシミュレーションモデルを示す図である。これは、平行平板導波路２０の内部に損失性誘電体膜２１，２２がある場合の透過性について、計算したものである。損失性誘電体膜２１，２２は、２４ｍｍの間隔で２箇所に配置されている。損失性誘電体膜２１，２２は、幅０．３ｍｍのスリット２１ａ，２２ａを有し、厚さは０．５ｍｍ、比誘電率２．０、比透磁率２．０、電気伝導率１５．０Ｓ／ｍである。上下の平行平板２３，２４の表面には電気壁境界条件を設定し、導波路の終端２５，２６にはＭｕｒの吸収境界条件を設定した。なお、Ｍｕｒの吸収境界条件については、「G,Mur “absorbing Boundary Conditions for the finite-difference approximation of the time domain electromagnetic field equations,"  IEEE Trans,EMC-23,1981,pp.377-382 」に記されている。セルサイズは、損失性誘電体膜２１，２２及びスリット２１ａ，２２ａを構成する部分が０．１×０．１ｍｍ2 、それ以外は基本的に１．０×１．０ｍｍ2 としている。
【００９９】
図５の左端終端部手前の電界垂直成分をガウス波形で励振し、右端終端部手前の電界垂直成分の波形を観測して、励振波と観測波とのフーリエ変換から透過特性を算出した。なお、時間離散間隔Δｔは、本発明の手法を用いる場合は、時間波形をフーリエ変換する際のサンプリング定理による制約から３．３３ｐｓ、通常のＦＤＴＤ法を用いる場合は、ＣＦＬ条件の制約から０．２３ｐｓとした。
【０１００】
図６は、第１のシミュレーションモデルにおける計算結果を示す図である。この図は、横軸が電磁波の周波数（ＧＨｚ）であり、縦軸は導体損失（ｄＢ）である。また、図中に点線で示した「ADI FDTD」が本発明による計算結果であり、実線で示した「Normal FDTD 」が通常のＦＤＴＤ法による計算結果である。図６から分かるように、「ADI FDTD」と「Normal FDTD 」との結果は、良好に一致している。
【０１０１】
ここで、計算コストに関するデータを以下の表に示す。
【０１０２】
【表１】
0003639116.tif 000036
【０１０３】
なお、この計算は、一般的なワークステーションを用いて行った。また、表中の［  ］は、比率を表している。この比率は、「時間離散間隔( Δｔ) 」と「メモリー容量(Memory)」については「Normal FDTD 」を基準とし、「ステップ数(Step)」と「演算時間(CPU Time)」については、「ADI FDTD」を基準としている。
【０１０４】
比較の結果、「ADI FDTD」では「Normal FDTD 」に比べて、１４．５倍の時間離散間隔を設定できる。そのため、ステップ数は、１／１４．５ですむ。その結果、演算時間が１／３．８に短縮できた。メモリー容量は、若干増加し、１．３倍となったが、これは、電界成分を記憶する配列数の増加によるものである。  次に、さらにセルサイズを小さくした場合のシミュレーション結果について説明する。
【０１０５】
図７は、第２のシミュレーションモデルを示す図である。この平行平板導波路３０の損失性誘電体膜３１，３２は、１９．６４ｍｍの間隔で２箇所に配置されている。損失性誘電体膜３１，３２は幅０．３ｍｍのスリット３１ａ，３２ａを有し、厚さは０．１ｍｍ、比誘電率２．０、比透磁率２．０、電気伝導率７４．９Ｓ／ｍである。セルサイズは、損失性誘電体膜３１，３２のスリット３１ａ，３２ａを構成する部分は０．０２×０．１ｍｍ2 、それ以外は基本的に１．０×１．０ｍｍ2 である。上下の平行平板３３，３４の表面には電気壁境界条件を設定し、導波路の終端３５，３６にはＭｕｒの吸収境界条件を設定した。励起条件は第１のシミュレーションモデルと同様である。時間離散間隔Δｔは、本発明の手法を用いる場合には３．３３ｐｓとし、通常のＦＤＴＤ法を用いる場合には、ＣＦＬ条件の制約から０．０６５ｐｓとした。
【０１０６】
図８は、第２のシミュレーションの計算結果を示す図である。第２のシミュレーションでは、第１のシミュレーションの結果に比べて定量的な差違が若干増加したが、周波数特性としては良好に一致している。
【０１０７】
第２のシミュレーションの計算コストに関するデータを以下の表に示す。
【０１０８】
【表２】
0003639116.tif 000037
【０１０９】
「ADI FDTD」では、「Normal FDTD 」に比べてΔｔを５１．２倍の大きさで設定できる。そのため、時間ステップ数が１／５１．２ですむ。その結果、演算時間が１／１４に短縮できた。また、使用メモリ容量の増加は、第１のシミュレーションの場合と同様である。
【０１１０】
ここで、第２のシミュレーションの方が第１のシミュレーションよりも演算時間の短縮効果が高いことが分かる。すなわち、セルサイズが小さくなるほど、本発明の手法を適用することによる処理の高速化の効果が得られる。
【０１１１】
ところで、上記の実施の形態では２次元ＴＥ波の解析を行うための電磁波解析装置について説明したが、同様にして、２次元ＴＭ(Transverse Magnetic) 波や３次元波の伝搬を解析することもできる。なお、ＴＭ波とは、磁界ベクトルがすべての場所で波の進行方向に垂直な電磁波である。
【０１１２】
２次元ＴＭ波の伝搬を解析する場合には、格子点（ｉ，ｊ）を基準とした格子に関して、（ｉ，ｊ＋１／２）のｘ方向の磁界Ｈi,j+1/2 、（ｉ＋１／２，ｊ）のｙ方向の磁界Ｈi+1/2,j 、及び（ｉ，ｊ）のｚ方向の電界Ｅi,j を求める。そして、２次元ＴＭ波の伝搬は以下のように定式化される。
【０１１３】
時刻「（ｎ＋１／２）Δｔ」における電磁波の状態を求める前半過程は、
【０１１４】
【数４１】
0003639116.tif 000038
【０１１５】
【数４２】
0003639116.tif 000039
【０１１６】
【数４３】
0003639116.tif 000040
【０１１７】
であり、時刻「（ｎ＋１）Δｔ」における電磁波の状態を求める後半過程は、
【０１１８】
【数４４】
0003639116.tif 000041
【０１１９】
【数４５】
0003639116.tif 000042
【０１２０】
【数４６】
0003639116.tif 000043
【０１２１】
である。式（４１）、式（４５）は簡単な代入計算で解き、式（４２）と式（４４）とは、２次元ＴＥ波の場合と同様に３項方程式を解くことで、前半過程及び後半過程の電磁波の伝搬を解析する。さらに、式（４３）と式（４６）とは、それぞれ式（４２）、式（４４）の計算結果を用いて、簡単な代入式で解くことができる。
【０１２２】
また、３次元の電磁波の伝搬を解析する場合には、ｘ−ｙ−ｚの３次元空間内の格子点（ｉ，ｊ，ｋ）からのｘ方向の空間離散間隔をΔｘ、ｙ方向の空間離散間隔をΔｙ、ｚ方向の空間離散間隔をΔｚとしたとき、格子点（ｉ，ｊ，ｋ）を基準とした格子に関して、（ｉ＋１／２，ｊ，ｋ）のｘ方向の電界Ｅi+1/2,j,k 、（ｉ，ｊ＋１／２，ｋ）のｙ方向の電界Ｅi,j+1/2,k 、（ｉ，ｊ，ｋ＋１／２）のｚ方向の電界Ｅi,j,k+1/2 、（ｉ，ｊ＋１／２，ｋ＋１／２）のｘ方向の磁界Ｈi,j+1/2,k+1/2 、（ｉ＋１／２，ｊ，ｋ＋１／２）のｙ方向の磁界Ｈi+1/2,j,k+1/2 、及び（ｉ＋１／２，ｊ＋１／２，ｋ）のｚ方向の磁界Ｈi+1/2,j+1/2,k を求める。３次元の電磁波は２次元のＴＥ波とＴＭ波の重ね合わせで表現できるため、以下のように定式化される。
【０１２３】
時刻「（ｎ＋１／２）Δｔ」における電磁波の状態を求める前半過程は、
【０１２４】
【数４７】
0003639116.tif 000044
【０１２５】
【数４８】
0003639116.tif 000045
【０１２６】
【数４９】
0003639116.tif 000046
【０１２７】
【数５０】
0003639116.tif 000047
【０１２８】
【数５１】
0003639116.tif 000048
【０１２９】
【数５２】
0003639116.tif 000049
【０１３０】
であり、時刻「（ｎ＋１）Δｔ」における電磁波の状態を求める後半過程は、
【０１３１】
【数５３】
0003639116.tif 000050
【０１３２】
【数５４】
0003639116.tif 000051
【０１３３】
【数５５】
0003639116.tif 000052
【０１３４】
【数５６】
0003639116.tif 000053
【０１３５】
【数５７】
0003639116.tif 000054
【０１３６】
【数５８】
0003639116.tif 000055
【０１３７】
である。前半過程においては、式（４７）と式（５１）、式（４８）と式（５２）、及び式（４９）と式（５０）の間で漸化式を導出する。また、後半過程においては、式（５３）と式（５８）、式（５４）と式（５６）、及び式（５５）と式（５７）の間で漸化式を導出する。導出された各漸化式は、ＴＥ波と同様に３項方程式に帰着する。従って、３次元波の電磁波の伝搬も同様にシミュレーションすることができる。
【０１３８】
なお、電磁波解析装置が有すべき機能の処理内容は、コンピュータで読み取り可能な記録媒体に記録されたプログラムに記述しておくことができる。そして、このプログラムをコンピュータで実行することにより、上記処理がコンピュータで実現される。コンピュータで読み取り可能な記録媒体としては、磁気記録装置や半導体メモリ等がある。市場に流通させる場合には、ＣＤ−ＲＯＭ(Compact Disk Read Only Memory) やフロッピーディスク等の可搬型記録媒体にプログラムを格納して流通させたり、ネットワークを介して接続されたコンピュータの記憶装置に格納しておき、ネットワークを通じて他のコンピュータに転送することもできる。コンピュータで実行する際には、コンピュータ内のハードディスク装置等にプログラムを格納しておき、メインメモリにロードして実行する。
【０１３９】
以下に、本発明の特徴を項ごとに区分し、その発明を特定するために必要と認める事項を列挙する。
（第１項）  マックスウェル方程式を時間と空間とについて差分法で解くことにより電磁波の伝搬を解析する電磁波解析装置において、
時間離散間隔をΔｔとしたとき、第１の時刻における電磁波の状態に基づいて、前記第１の時刻からΔｔ経過後の第２の時刻における電磁波の状態を交互方向陰解法型の有限差分時間領域法により求めることで、電磁波の挙動を解析する解析手段、
を有することを特徴とする電磁波解析装置。
【０１４０】
（第２項）  前記解析手段は、２次元ＴＥ波の伝搬を解析することを特徴とする第１項記載の電磁波解析装置。
（第３項）  前記解析手段は、電界の振動面をｘ−ｙ平面とし、前記ｘ−ｙ平面の任意の格子点（ｉ，ｊ）を基準とした差分格子の電磁波を解析する際には、格子点（ｉ，ｊ）からｘ方向への空間離散間隔をΔｘ、格子点（ｉ，ｊ）からｙ方向への空間離散間隔をΔｙとした場合に、（ｉ＋１／２，ｊ）におけるｘ方向の電界、（ｉ，ｊ＋１／２）におけるｙ方向の電界、及び（ｉ＋１／２，ｊ＋１／２）における前記ｘ−ｙ平面に垂直なｚ方向の磁界を求めることを特徴とする第２項記載の電磁波解析装置。
【０１４１】
（第４項）  前記解析手段は、ｘ方向の電界に対して陽解法を用い、ｙ方向の電界、及びｚ方向の磁界に対して陰解法を用いることにより、前記第１の時刻における電磁波の状態に基づいて前記第１の時刻からΔｔ／２経過後の中間時刻における電磁波の状態を演算する前半過程演算手段と、ｙ方向の電界に対して陽解法を用い、ｘ方向の電界、及びｚ方向の磁界に対して陰解法を用いることにより、前記前半過程演算手段により算出された前記中間時刻の電磁波の状態に基づいて前記第２の時刻における電磁波の状態を演算する後半過程演算手段と、を具備することを特徴とする第３項記載の電磁波解析装置。
【０１４２】
（第５項）  前記前半過程演算手段は、シミュレーションにおけるｎステップ目の時刻を前記第１の時刻、ｎ＋１ステップ目の時刻を前記第２の時とし、前記第１の時刻における（ｉ＋１／２，ｊ）でのｘ方向の電界をＥｘi+1/2,j n 、前記第１の時刻における（ｉ，ｊ＋１／２）でのｙ方向の電界をＥｙi,j+1/2 n 、前記第１の時刻における（ｉ＋１／２，ｊ＋１／２）でのｚ方向の磁界をＨｚi+1/2,j+1/2 n 、前記中間時刻における（ｉ＋１／２，ｊ）でのｘ方向の電界をＥｘi+1/2,j n+1/2 、前記中間時刻における（ｉ，ｊ＋１／２）でのｙ方向の電界をＥｙi,j+1/2 n+1/2 、前記中間時刻における（ｉ＋１／２，ｊ＋１／２）でのｚ方向の磁界をＨｚi+1/2,j+1/2 n+1/2 、解析すべき領域の誘電率をε、解析すべき領域の透磁率をμとした場合に、格子状に分割された前記ｘ−ｙ平面の任意の格子点（ｉ，ｊ）を基準とした電磁波を
【０１４３】
【数５９】
0003639116.tif 000056
【０１４４】
【数６０】
0003639116.tif 000057
【０１４５】
【数６１】
0003639116.tif 000058
【０１４６】
と定義することで、前記中間時刻におけるすべての格子点を基準とした電磁波の状態を演算し、
前記後半過程演算手段は、前記第２の時刻における（ｉ＋１／２，ｊ）でのｘ方向の電界をＥｘi+1/2,j n+1 、前記第２の時刻における（ｉ，ｊ＋１／２）でのｙ方向の電界をＥｙi,j+1/2 n+1 、前記第２の時刻における（ｉ＋１／２，ｊ＋１／２）でのｚ方向の磁界をＨｚi+1/2,j+1/2 n+1 とした場合に、格子状に分割された前記ｘ−ｙ平面の任意の格子点（ｉ，ｊ）を基準とした電磁波を
【０１４７】
【数６２】
0003639116.tif 000059
【０１４８】
【数６３】
0003639116.tif 000060
【０１４９】
【数６４】
0003639116.tif 000061
【０１５０】
と定義することで、前記第２の時刻におけるすべての格子点を基準とした電磁波の状態を演算することを特徴とする第４項記載の電磁波解析装置。
（第６項）  前記解析手段は、２次元ＴＭ波の伝搬を解析することを特徴とする第１項記載の電磁波解析装置。
【０１５１】
（第７項）  前記解析手段は、磁界の振動面をｘ−ｙ平面とし、前記ｘ−ｙ平面の任意の格子点（ｉ，ｊ）を基準とした差分格子の電磁波を解析する際には、格子点（ｉ，ｊ）からｘ方向への空間離散間隔をΔｘ、格子点（ｉ，ｊ）からｙ方向への空間離散間隔をΔｙとした場合に、（ｉ，ｊ＋１／２）におけるｘ方向の磁界、（ｉ＋１／２，ｊ）におけるｙ方向の磁界、及び（ｉ，ｊ）における前記ｘ−ｙ平面に垂直なｚ方向の電界を求めることを特徴とする第６項記載の電磁波解析装置。
【０１５２】
（第８項）  前記解析手段は、ｘ方向の磁界に対して陽解法を用い、ｙ方向の磁界、及びｚ方向の電界に対して陰解法を用いることにより、前記第１の時刻における電磁波の状態に基づいて前記第１の時刻からΔｔ／２経過後の中間時刻における電磁波の状態を演算する前半過程演算手段と、ｙ方向の磁界に対して陽解法を用い、ｘ方向の磁界、及びｚ方向の電界に対して陰解法を用いることにより、前記前半過程演算手段により算出された前記中間時刻の電磁波の状態に基づいて前記第２の時刻における電磁波の状態を演算する後半過程演算手段と、を具備することを特徴とする第７項記載の電磁波解析装置。
【０１５３】
（第９項）  前記前半過程演算手段は、シミュレーションにおけるｎステップ目の時刻を前記第１の時刻、ｎ＋１ステップ目の時刻を前記第２の時とし、前記第１の時刻における（ｉ，ｊ＋１／２）でのｘ方向の磁界をＨｘi,j+1/2 n 、前記第１の時刻における（ｉ＋１／２，ｊ）でのｙ方向の磁界をＨｙi+1/2,j n 、前記第１の時刻における（ｉ，ｊ）でのｚ方向の電界をＥｚi,j n 、前記中間時刻における（ｉ，ｊ＋１／２）でのｘ方向の磁界をＨｘi,j+1/2 n+1/2 、前記中間時刻における（ｉ＋１／２，ｊ）でのｙ方向の磁界をＨｙi+1/2,j n+1/2 、前記中間時刻における（ｉ，ｊ）でのｚ方向の電界をＥｚi,j n+1/2 、解析すべき領域の誘電率をε、解析すべき領域の透磁率をμとした場合に、格子状に分割された前記ｘ−ｙ平面の任意の格子点（ｉ，ｊ）を基準とした電磁波を
【０１５４】
【数６５】
0003639116.tif 000062
【０１５５】
【数６６】
0003639116.tif 000063
【０１５６】
【数６７】
0003639116.tif 000064
【０１５７】
と定義することで、前記中間時刻におけるすべての格子点を基準とした電磁波の状態を演算し、
前記後半過程演算手段は、前記第２の時刻における（ｉ，ｊ＋１／２）でのｘ方向の磁界をＨｘi,j+1/2 n+1 、前記第２の時刻における（ｉ＋１／２，ｊ）でのｙ方向の磁界をＨｙi+1/2,j n+1 、前記第２の時刻における（ｉ，ｊ）でのｚ方向の電界をＥｚi,j n+1 とした場合に、格子状に分割された前記ｘ−ｙ平面の任意の格子点（ｉ，ｊ）を基準とした電磁波を
【０１５８】
【数６８】
0003639116.tif 000065
【０１５９】
【数６９】
0003639116.tif 000066
【０１６０】
【数７０】
0003639116.tif 000067
【０１６１】
と定義することで、前記第２の時刻におけるすべての格子点を基準とした電磁波の状態を演算することを特徴とする第８項記載の電磁波解析装置。
（第１０項）  前記解析手段は、３次元波の伝搬を解析することを特徴とする第１項記載の電磁波解析装置。
【０１６２】
（第１１項）  前記解析手段は、ｘ−ｙ−ｚの３次元空間を想定し、３次元空間の任意の格子点（ｉ，ｊ，ｋ）を基準とした差分格子の電磁波を解析する際には、格子点（ｉ，ｊ，ｋ）からｘ方向への空間離散間隔をΔｘ、格子点（ｉ，ｊ，ｋ）からｙ方向への空間離散間隔をΔｙ、格子点（ｉ，ｊ，ｋ）からｚ方向への空間離散間隔をΔｚとした場合に、（ｉ＋１／２，ｊ，ｋ）におけるｘ方向の電界、（ｉ，ｊ＋１／２，ｋ）におけるｙ方向の電界、（ｉ，ｊ，ｋ＋１／２）におけるｚ方向の電界、（ｉ，ｊ＋１／２，ｋ＋１／２）におけるｘ方向の磁界、（ｉ＋１／２，ｊ，ｋ＋１／２）におけるｙ方向の磁界、及び（ｉ＋１／２，ｊ＋１／２，ｋ）におけるｚ方向の磁界を求めることを特徴とする第１０項記載の電磁波解析装置。
【０１６３】
（第１２項）  前記解析手段は、
シミュレーションにおけるｎステップ目の時刻を前記第１の時刻、ｎ＋１ステップ目の時刻を前記第２の時とし、前記第１の時刻における（ｉ＋１／２，ｊ，ｋ）でのｘ方向の電界をＥｘi+1/2,j,k n 、前記第１の時刻における（ｉ，ｊ＋１／２，ｋ）でのｙ方向の電界をＥｙi,j+1/2,k n 、前記第１の時刻における（ｉ，ｊ，ｋ＋１／２）でのｚ方向の電界をＥｚi,j,k+1/2 n 、前記第１の時刻における（ｉ，ｊ＋１／２，ｋ＋１／２）でのｘ方向の磁界をＨｘi,j+1/2,k+1/2 n 、前記第１の時刻における（ｉ＋１／２，ｊ，ｋ＋１／２）でのｙ方向の磁界をＨｙi+1/2,j,k+1/2 n 、前記第１の時刻における（ｉ＋１／２，ｊ＋１／２，ｋ）でのｚ方向の磁界をＨｚi+1/2,j+1/2,k n 、前記中間時刻における（ｉ＋１／２，ｊ，ｋ）でのｘ方向の電界をＥｘi+1/2,j,k n+1/2 、前記中間時刻における（ｉ，ｊ＋１／２，ｋ）でのｙ方向の電界をＥｙi,j+1/2,k n+1/2 、前記中間時刻における（ｉ，ｊ，ｋ＋１／２）でのｚ方向の電界をＥｚi,j,k+1/2 n+1/2 、前記中間時刻における（ｉ，ｊ＋１／２，ｋ＋１／２）でのｘ方向の磁界をＨｘi,j+1/2,k+1/2 n+1/2 、前記中間時刻における（ｉ＋１／２，ｊ，ｋ＋１／２）でのｙ方向の磁界をＨｙi+1/2,j,k+1/2 n+1/2 、前記中間時刻における（ｉ＋１／２，ｊ＋１／２，ｋ）でのｚ方向の磁界をＨｚi+1/2,j+1/2,k n+1/2 、解析すべき領域の誘電率をε、解析すべき領域の透磁率をμとした場合に、格子状に分割された前記３次元空間の任意の格子点（ｉ，ｊ，ｋ）を基準とした電磁波を
【０１６４】
【数７１】
0003639116.tif 000068
【０１６５】
【数７２】
0003639116.tif 000069
【０１６６】
【数７３】
0003639116.tif 000070
【０１６７】
【数７４】
0003639116.tif 000071
【０１６８】
【数７５】
0003639116.tif 000072
【０１６９】
【数７６】
0003639116.tif 000073
【０１７０】
と定義することで、前記中間時刻におけるすべての格子点を基準とした電磁波の状態を演算する前半過程演算手段と、
前記中間時刻における（ｉ＋１／２，ｊ，ｋ）でのｘ方向の電界をＥｘi+1/2,j,k n+1 、前記中間時刻における（ｉ，ｊ＋１／２，ｋ）でのｙ方向の電界をＥｙi,j+1/2,k n+1 、前記中間時刻における（ｉ，ｊ，ｋ＋１／２）でのｚ方向の電界をＥｚi,j,k+1/2 n+1 、前記中間時刻における（ｉ，ｊ＋１／２，ｋ＋１／２）でのｘ方向の磁界をＨｘi,j+1/2,k+1/2 n+1 、前記中間時刻における（ｉ＋１／２，ｊ，ｋ＋１／２）でのｙ方向の磁界をＨｙi+1/2,j,k+1/2 n+1 、前記中間時刻における（ｉ＋１／２，ｊ＋１／２，ｋ）でのｚ方向の磁界をＨｚi+1/2,j+1/2,k n+1 とした場合に、格子状に分割された前記３次元空間の任意の格子点（ｉ，ｊ、ｋ）を基準とした電磁波を
【０１７１】
【数７７】
0003639116.tif 000074
【０１７２】
【数７８】
0003639116.tif 000075
【０１７３】
【数７９】
0003639116.tif 000076
【０１７４】
【数８０】
0003639116.tif 000077
【０１７５】
【数８１】
0003639116.tif 000078
【０１７６】
【数８２】
0003639116.tif 000079
【０１７７】
と定義することで、前記第２の時刻におけるすべての格子点を基準とした電磁波の状態を演算する前記後半過程演算手段と、
を具備することを特徴とする第１１項記載の電磁波解析装置。
【０１７８】
（第１３項）  マックスウェル方程式を時間と空間とについて差分法で解くことにより電磁波の伝搬を解析するための電磁波解析プログラムを記録したコンピュータ読み取り可能な記録媒体において、
時間離散間隔をΔｔとしたとき、第１の時刻における電磁波の状態に基づいて、前記第１の時刻からΔｔ経過後の第２の時刻における電磁波の状態を交互方向陰解法型の有限差分時間領域法により求めることで、電磁波の挙動を解析する解析手段、
としてコンピュータを機能させることを特徴とする電磁波解析プログラムを記録したコンピュータ読み取り可能な記録媒体。
【０１７９】
【発明の効果】
以上説明したように本発明の電磁波解析装置では、電磁波の状態を交互方向陰解法型の有限差分時間領域法により求めることで、電磁波の挙動を解析するようにしたため、空間離散間隔に依存せずに時間離散間隔を定めることが可能となる。その結果、微細な構造を有する解析対象物の電磁波解析を行う際の時間離散間隔を比較的大きくとることが可能となり、解析処理を高速に行うことができる。
【０１８０】
また、本発明の電磁波解析プログラムを記録したコンピュータ読み取り可能な記録媒体では、記録された電磁波解析プログラムをコンピュータに実行させることにより、電磁波の状態を交互方向陰解法型の有限差分時間領域法による電磁波の解析をコンピュータで高速に行うことが可能となる。
【図面の簡単な説明】
【図１】本発明の原理構成図である。
【図２】格子の拡大図である。
【図３】ＡＤＩ−ＦＤＴＤ法による２次元ＴＥ波解析の処理手順を示すフローチャートである。
【図４】電磁波解析装置のハードウェア構成を示す図である。
【図５】第１のシミュレーションモデルを示す図である。
【図６】第１のシミュレーションモデルにおける計算結果を示す図である。
【図７】第２のシミュレーションモデルを示す図である。
【図８】第２のシミュレーションの計算結果を示す図である。
【符号の説明】
１  解析指令入力手段
２  解析手段
２ａ  前半過程演算手段
２ｂ  後半過程演算手段
３  解析結果出力手段
【図１】
0003639116.tif 000086
【図２】
0003639116.tif 000087
【図３】
0003639116.tif 000088
【図４】
0003639116.tif 000089
【図５】
0003639116.tif 000090
【図６】
0003639116.tif 000091
【図７】
0003639116.tif 000092
【図８】
0003639116.tif 000093