JP3639116B2 - Electromagnetic wave analysis apparatus and computer-readable recording medium recording electromagnetic wave analysis program - Google Patents

Description

【０００５】
【数８】

【０００６】
ここで、ｎは演算開始時点から時間軸方向に進めたステップの数であり、Δｔは時間離散間隔（１ステップで進められる時間）であり、ｉは差分格子（セル）の原点からの順番であり、Δｘは空間離散間隔（分割領域の幅）である。この微分方程式を陽解法で解く場合、時間微分に対して前進差分で近似を行い、
【０００７】
【数９】

【０００８】
という式を得る。この式（９）はｒ＝Δｔ／Δｘ² とおいて
【０００９】
【数１０】

【００１０】
と書くことができる。これにより、ｎ＋１番目のステップにおける解が、過去の時刻（ｎ番目のステップ）の解から代入計算によりただちに計算できる。
【００１１】
【発明が解決しようとする課題】
しかし、数（１０）の陽解法型差分法では、解が発散しないための安定条件として
【００１２】
【数１１】

【００１３】
を満たしていなければならない。
また、ＦＤＴＤ法においては、解が発散しないためにはＣＦＬ(Courant, Friedrich and Levy) 条件あるいはCourant 条件を満たす必要がある。
【００１４】
ＣＦＬ条件は、解が発散しないためには、
【００１５】
【数１２】

【００１６】
の安定条件を満たしていなければならないという条件である。ＦＤＴＤ法のＣＦＬ条件については、「A.Taflove,Computational Electrodynamics,MA,Artech,1995」に詳しい。なお、式（１２）は、３次元波を解析する場合のCourant 条件である。ここで、ｖは、電磁波の速度であり、Δｘ_min 、Δｙ_min 、Δｚ_min はそれぞれＸＹＺ方向の空間離散間隔の最小値である。
【００１７】
これらの制約があることにより、陽解法においては、時間離散間隔の最大値が空間離散間隔の最小値によって決定されるという性質がある。このため、微細な構造を有する解析対象の時間応答を長時間計算する場合、シミュレーションの時刻を進行させるための計算の繰り返し回数が増加し、シミュレーションに要する時間が増加するといった問題があった。
【００１８】
なお、式（７）、式（８）に示した偏微分方程式は、陰解法(implicit method) で解くこともできる。例えば、式（７）、式（８）のような１次元放物型偏微分方程式を陰解法で解く場合には、時間微分について後退差分で近似を行い、
【００１９】
【数１３】

【００２０】
を得る。この式（１３）は、陽解法の場合と同様にｒ＝Δｔ／Δｘ² とおいて
【００２１】
【数１４】

【００２２】
と書くことができる。このような式（１４）を解くのが陰解法である。この陰解法は、陽解法と異なり常に安定である。
ただし、陰解法は、以下のような連立方程式を解かなくてはｆ_i ⁿ⁺¹ の値を得ることができない。
【００２３】
【数１５】

【００２４】
ここで、ａ_i は式（１４）において「ｉ＝１，２，・・・imax」とした時の左辺第３項の定数部分（−ｒ）であり、ｂ_i は式（１４）において「ｉ＝１，２，・・・imax」とした時の左辺第２項の定数部分（１＋２ｒ）であり、ｃ_i は式（１４）において「ｉ＝１，２，・・・imax」とした時の左辺第１項の定数部分（−ｒ）である。
【００２５】
電磁波解析の過渡的な現象を解析するには、少なくとも２次元以上の偏微分方程式を解く必要があるが、２次元以上の偏微分方程式の陰解法については別途注意が必要となる。例えば２次元放物型偏微分方程式は、以下の式で表される。
【００２６】
【数１６】

【００２７】
この式（１６）を、例えばCrank-Nicolson法で解く場合には、
【００２８】
【数１７】

【００２９】
のように時間微分を近似する。この式（１７）を連立方程式によって解く場合、（ｘ方向の格子数−１）×（ｙ方向の格子数−１）元の連立１次方程式を解く必要がある。そのため、格子数が多い場合には計算量や必要な記憶容量が膨大になる。すなわち、格子を微細化すると、陽解法で解析した場合と同様に非常に時間のかかる計算となる。
【００３０】
本発明はこのような点に鑑みてなされたものであり、微細な構造を有する解析対象物の電磁波解析を高速に行うことのできる電磁波解析装置を提供することを目的とする。
【００３１】
また、本発明の他の目的は、微細な構造を有する解析対象物の電磁波解析をコンピュータに高速に行わせることのできる電磁波解析プログラムを記録したコンピュータ読み取り可能な記録媒体を提供することである。
【００３２】
【課題を解決するための手段】
本発明では上記課題を解決するために、マックスウェル方程式を時間と空間とについて差分法で解くことにより電磁波の伝搬を解析する電磁波解析装置において、時間離散間隔をΔｔとしたとき、第１の時刻における電磁波の状態に基づいて、前記第１の時刻からΔｔ経過後の第２の時刻における電磁波の状態を交互方向陰解法型の有限差分時間領域法により求めることで、電磁波の挙動を解析する解析手段、を有することを特徴とする電磁波解析装置が提供される。
【００３３】
このような電磁波解析装置によれば、解析手段により、第１の時刻における電磁波の状態に基づいて、第１の時刻からΔｔ経過後の第２の時刻における電磁波の状態が交互方向陰解法型の有限差分時間領域法により求められる。このように交互方向陰解法型の有限差分時間領域法によって得られた解は、時間離散間隔の大小によらず安定している。
【００３４】
また、上記課題を解決するために、マックスウェル方程式を時間と空間とについて差分法で解くことにより電磁波の伝搬を解析するための電磁波解析プログラムを記録したコンピュータ読み取り可能な記録媒体において、時間離散間隔をΔｔとしたとき、第１の時刻における電磁波の状態に基づいて、前記第１の時刻からΔｔ経過後の第２の時刻における電磁波の状態を交互方向陰解法型の有限差分時間領域法により求めることで、電磁波の挙動を解析する解析手段、としてコンピュータを機能させることを特徴とする電磁波解析プログラムを記録したコンピュータ読み取り可能な記録媒体が提供される。
【００３５】
この記録媒体に記録された電磁波解析プログラムをコンピュータに実行させれば、上記本発明に係る電磁波解析装置の機能がコンピュータ上で実現される。
【００３６】
【発明の実施の形態】
以下、本発明の実施の形態を図面を参照して説明する。本発明の電磁波解析装置は、ＡＤＩ(Alternating Direction Implicit)法型のＦＤＴＤ法を用いて電磁波の解析を行うことにより、解の安定性を確保しつつ高速の演算を可能としたものである。
【００３７】
図１は、本発明の原理構成図である。解析指令入力手段１は、電磁波解析に必要な初期値や境界吸収条件を指定して、解析指令を入力する。
解析手段２は、解析指令入力手段１からの解析指令を受け取ると、２次元ＴＥ(Transverse Electric) 波の伝搬を、交互方向陰解法型−差分時間領域（ＡＤＩ−ＦＤＴＤ）法によって数値解析する。ここで、ＴＥ波とは以下の式で表現される電磁波である。
【００３８】
【数１８】

【００３９】
【数１９】

【００４０】
【数２０】

【００４１】
また、ＡＤＩ−ＦＤＴＤ法は、時間離散間隔を前半と後半とに分割し、前半過程の演算と後半過程の演算とを交互に行う方法である。例えば、式（１６）の２次元放物型偏微分方程式をＡＤＩ法で解く場合には、時間微分の値を、
【００４２】
【数２１】

【００４３】
【数２２】

【００４４】
のように近似する。式（２１）が前半過程であり、左辺の時間微分は、ｘ方向のｆから見ると後退差分で近似し、ｙ方向のｆから見ると前進差分で近似している。式（２２）が後半過程であり、左辺の時間微分は、ｘ方向ｆから見ると前進差分で近似し、ｙ方向のｆから見ると後退差分で近似している。前半過程は前半過程演算手段２ａが演算し、後半過程は後半過程演算手段２ｂが演算する。なお、Δｔは時間離散間隔であり、Δｘはｘ方向の空間離散間隔であり、Δｙはｙ方向の空間離散間隔である。
【００４５】
前半過程演算手段２ａでは、後半過程演算手段２ｂから時刻「ｎΔｔ」（ｎはシミュレーション開始時からの時刻進行処理のステップ数であり、前半過程と後半過程ではそれぞれ１／２ステップずつ進められる）における電磁波の状態の情報を受け取ると、時刻「（ｎ＋１／２）Δｔ」の電磁波の状態を計算する。その際、時刻「（ｎ＋１／２）Δｔ」におけるｘ方向の電界Ｅｘ^n+1/2 を陽解法により求め、ｙ方向の電界Ｅｙ^n+1/2 とｚ方向の磁界Ｈｚ^n+1/2 とを陰解法により求める。計算結果は、後半過程演算手段２ｂに引き渡す。
【００４６】
後半過程演算手段２ｂは、前半過程演算手段２ａから時刻「（ｎ＋１／２）Δｔ」における電磁波の状態を受け取ると、時刻「（ｎ＋１）Δｔ」の電磁波の状態を計算する。その際、時刻「（ｎ＋１）Δｔ」におけるｙ方向の電界Ｅｙⁿ⁺¹ を陽解法により求め、ｘ方向の電界Ｅｘⁿ⁺¹ とｚ方向の磁界Ｈｚⁿ⁺¹ とを陰解法により求める。計算結果は、前半過程演算手段２ａに引き渡す。
【００４７】
解析結果出力手段３は、解析手段２による解析結果を出力する。
以下に、時刻「ｎΔｔ」の電磁波の状態から時刻「（ｎ＋１）Δｔ」の状態を求める際の解析手段２の計算内容について説明する。
【００４８】
まず、解析対象領域を長方形の格子に分割する。ここで、解析対象の領域のｘ方向をｉmax 個の格子に分割し、ｙ方向をｊmax 個の格子に分割したものとする。
【００４９】
図２は、格子の拡大図である。格子点はの位置は（ｉ，ｊ）のような２組の数字で表される。ｉは、原点からのｘ方向の格子の順番である。ｊは、原点からのｙ方向の格子の順番である。また、「ｉ」と「ｉ＋１」との中間を「ｉ＋１／２」と表す。ここで、格子点（ｉ，ｊ）を基準とした領域の２次元ＴＥ波の解析を行う場合、（ｉ＋１／２，ｊ）におけるｘ方向の電界、（ｉ，ｊ＋１／２）におけるｙ方向の電界、及び（ｉ＋１／２，ｊ＋１／２）におけるｚ方向の磁界の時間変化を差分法で解析する。
【００５０】
マックスウェル方程式を時間と空間について差分法で解く場合、前半過程ではＥ_x の時間微分を前進差分で近似し（以下の式（２３）に示す）、Ｅ_y の時間微分を後退差分で近似し（以下の式（２４）に示す）、Ｈ_z の時間微分はＥ_x から見れば前進差分で近似するとともに、Ｅ_y から見れば後退差分で近似する（以下の式（２５）に示す）。その結果、時刻「（ｎ＋１／２）Δｔ」の状態を求める前半過程は、以下の式で表せる。
【００５１】
【数２３】

【００５２】
【数２４】

【００５３】
【数２５】

【００５４】
ここで、電界Ｅや磁界Ｈの右肩の値は、シミュレーションの時刻ステップを表している。電界Ｅや磁界Ｈの右下の値は、解析領域を格子で分割した際の格子番号である。εは、解析対象物の誘電率である。μは、解析対象物の透磁率である。
【００５５】
ところで、式（２３）は陽解法的表現となっており、式（２４）と式（２５）とは陰解法的表現となっている。つまり前半過程では、Ｅ_x の時間微分を前進差分で近似しているため、時刻「（ｎ＋１／２）Δｔ」におけるｘ方向の電界Ｅx ^n+1/2 は陽解法により求めることができる。また、前半過程では、Ｅ_y の時間微分を後退差分で近似しているため、時刻「（ｎ＋１／２）Δｔ」におけるｙ方向の電界Ｅｙ^n+1/2 は陰解法により求めなければならない。さらに、Ｈ_z の時間微分をＥ_y から見ると後退差分で近似しているので、時刻「（ｎ＋１／２）Δｔ」におけるｚ方向の磁界Ｈｚ^n+1/2 は陰解法により求めなければならない。
【００５６】
同様に、時刻「（ｎ＋１）Δｔ」の状態を求める後半過程は、Ｅ_x の時間微分を後退差分で近似し（以下の式（２６）に示す）、Ｅ_y の時間微分を前進差分で近似し（以下の式（２７）に示す）、Ｈ_z の時間微分はＥ_x から見れば後退差分で近似するとともに、Ｅ_y から見れば前進差分で近似する（以下の式（２８）に示す）。これにより、後半過程は以下の式で表せる。
【００５７】
【数２６】

【００５８】
【数２７】

【００５９】
【数２８】

【００６０】
後半過程ではｙ方向の電界を求める式（２７）が陽解法的表現となり、その他の式（２６）、式（２８）が陰解法的表現となっている。つまり後半過程では、Ｅ_x の時間微分を後退差分で近似しているため、時刻「（ｎ＋１）Δｔ」におけるｘ方向の電界Ｅx ⁿ⁺¹ は陰解法により求めなければならない。また、Ｅ_y の時間微分を前進差分で近似しているため、時刻「（ｎ＋１）Δｔ」におけるｙ方向の電界Ｅｙⁿ⁺¹ は陽解法により求めることができる。さらに、Ｈ_z の時間微分をＥ_x から見ると後退差分で近似しているので、時刻「（ｎ＋１）Δｔ」におけるｚ方向の磁界Ｈｚⁿ⁺¹ は陰解法により求めなければならない。
【００６１】
前半過程の式（２３）により、時刻「（ｎ＋１／２）Δｔ」のｘ方向の電界は、時刻「ｎΔｔ」の電界及び磁界の代入計算によって求めることができる。同様に、後半過程の式（２７）も簡単な代入計算で求めることができる。
【００６２】
一方、前半過程の式（２４）、式（２５）は、右辺に時刻「（ｎ＋１／２）Δｔ」の電界や磁界が含まれているため、そのままでは解けない。そこで、式（２４）と式（２５）とを連立させて磁界Ｈｚ^n+1/2 を消去し、下の漸化式を得る。
【００６３】
【数２９】

【００６４】
後半過程の式（２６）、式（２８）も、右辺に時刻「（ｎ＋１）Δｔ」の電界と磁界が含まれるため、式（２６）と式（２８）とを連立させて磁界Ｈｚⁿ⁺¹ を消去し、以下の漸化式を得る。
【００６５】
【数３０】

【００６６】
式（２９）の漸化式を「ｉ＝１，・・・，ｉmax 」とした場合のそれぞれについて生成すれば、以下のような連立一次方程式に帰着する。
【００６７】
【数３１】

【００６８】
ここで式（３１）におけるａ_i は、「ｉ＝１，・・・，ｉmax 」とした場合に得られる式（２９）の左辺第１項の定数部分である。ｂ_i は、「ｉ＝１，・・・，ｉmax 」とした場合に得られる式（２９）の左辺第２項の定数部分である。ｃ_i は、「ｉ＝１，・・・，ｉmax 」とした場合に得られる式（２９）の左辺第３項の定数部分である。Ｄ_i ⁿ は、「ｉ＝１，・・・，ｉmax 」とした場合に得られる式（２９）の右辺の値である。また、第１行と最終行には、吸収境界条件が適用される。
【００６９】
式（３１）は、ｘ方向の格子点「ｉ」についての３項方程式(tridiagonal system of equation)であり、トーマス法等によって解くことができる。このような３項方程式を「ｉ＝１，・・・，ｉmax 」とした場合のそれぞれについて解けば、時刻「ｎ＋１／２」における全ての格子でのｙ方向の電界が求められる。ｙ方向の電界が求まれば、式（１６）によってｚ方向の磁界も求まる。
【００７０】
後半過程に関する漸化式である式（３０）も前半過程と同様に、以下のような「ｊ」について３項方程式を得ることができる。
【００７１】
【数３２】

【００７２】
ここで式（３２）におけるａ_j は、「ｊ＝１，・・・，ｊmax 」とした場合に得られる式（３０）の左辺第１項の定数部分である。ｂ_j は、「ｊ＝１，・・・，ｊmax 」とした場合に得られる式（３０）の左辺第２項の定数部分である。ｃ_j は、「ｊ＝１，・・・，ｊmax 」とした場合に得られる式（３０）の左辺第３項の定数部分である。Ｄ_j ^n+1/2 は、「ｊ＝１，・・・，ｊmax 」とした場合に得られる式（３０）の右辺の値である。また、第１行と最終行には、吸収境界条件が適用される。
【００７３】
このように、前半過程の式（２４）と後半過程の式（２６）とのそれぞれの演算を３項方程式に帰着することができる。従って、前半過程演算手段２ａにおいては式（２３）と式（３１）の３項方程式とを解き、さらに式（２５）を解くことになる。後半過程演算手段２ｂにおいては式（２７）と式（３２）の３項方程式とを解き、それから式（２８）を解くことになる。このような演算を交互に行うことにより時刻のステップを進めていけば、２次元ＴＥ波の伝搬をＡＤＩ−ＦＤＴＤ法により解析することができる。
【００７４】
図３は、ＡＤＩ−ＦＤＴＤ法による２次元ＴＥ波解析の処理手順を示すフローチャートである。この処理をステップ番号に沿って説明する。なお、吸収境界条件、時間離散間隔、空間離散間隔等の各種初期条件は、すでに与えられているものとする。シミュレーション開始時点では、シミュレーション時刻「Ｔ＝０」である。
［Ｓ１］前半過程演算手段２ａは、まず、すでに求められているＥ_x ⁿ 、Ｈ_z ⁿ の値を代入して式（２３）を解くことで、Ｅ_x ^n+1/2 を計算する。次に、すでに求められているＥ_x ⁿ 、Ｅ_y ⁿ 、Ｈ_z ⁿ の値を用いて式（３１）の３項方程式を解くことで、Ｅ_y ^n+1/2 を計算する。さらに、式（３１）を解くことによって得られたＥ_y ^n+1/2 の値と既知のＥ_x ⁿ 、Ｈ_z ⁿ の値とを用いて式（２５）を解くことで、Ｈ_z ^n+1/2 を計算する。これにより、シミュレーションの時刻Ｔは、ΔＴ／２だけ進められる。
［Ｓ２］後半過程演算手段２ｂは、まず、すでに求められているＥ_x ^n+1/2 、Ｅ_y ^1+1/2 、Ｈ_z ^n+1/2 の値を用いて式（３２）の３項方程式を解くことで、Ｅ_x ⁿ⁺¹ を計算する。次に、すでに求められているＥ_y ^n+1/2 、Ｈ_z ^n+1/2 の値を用いて式（２７）を解くことで、Ｅ_y ⁿ⁺¹ を計算する。さらに、式（３２）を解くことによって得られたＥ_x ⁿ⁺¹ の値と既知のＥ_y ^n+1/2 、Ｈ_z ^n+1/2 の値とを用いて式（２８）を解くことで、Ｈ_z ⁿ⁺¹ を計算する。これにより、シミュレーションの時刻Ｔは、ΔＴ／２だけ進められる。その結果、ステップＳ１とステップＳ２とのトータルで時間離散間隔であるΔｔだけ時刻が進められたことになる。
［Ｓ３］解析手段２は、時刻Ｔがあらかじめ設定されているシミュレーション最大時刻Ｔ_max を超えているか否かを判断する。Ｔ_max を超えていればシミュレーションを終了し、超えていなければステップＳ１に進みシミュレーションを継続する。
【００７５】
以上のような処理を解析手段２が実行することにより、ＡＤＩ−ＦＤＴＤ法を用いた２次元ＴＥ波の解析が行われる。
なお、以上の式は簡略化のために、電磁波を、無損失、等方性、非分散性、均質の媒質を伝搬するものと仮定し、さらに、差分格子（セル）は一様（ΔｔとΔｙがすべてのセルで等しい）であるものと仮定しているため、「ａ_i 」、「ｃ_i 」、「ａ_j 」、「ｃ_j 」の値はすべて「１」である。媒質が不均質または差分格子（セル）が一様でない場合は、式の拡張を行う必要がある。前述の式を媒質が不均質または差分格子（セル）が一様でない場合へ拡張する方法は、通常のＦＤＴＤ法で行われていた拡張方法と同じである。
【００７６】
以下に、２次元ＴＥ波の伝搬をＡＤＩ−ＦＤＴＤ法により解析した場合の解の安定性についてフーリエ級数法を用いて検討した結果を示す。
まず、時間発展に伴う誤差の成長係数ξは、時刻のステップ数ｎ’を用いて
【００７７】
【数３３】

【００７８】
と表される。なお、γは複素定数である。このような誤差の成長係数を用いて、前半過程および後半過程の解をそれぞれ次のようにおく。
【００７９】
【数３４】

【００８０】
【数３５】

【００８１】
ここで、ξ₁ ^n'及びξ₂ ^n'は、それぞれ前半過程及び後半過程での時間発展にともなう誤差の成長係数である。また、虚数単位は添字とまぎらわしいため、−１の平方根で表している。すると、これらは次式で表される。
【００８２】
【数３６】

【００８３】
【数３７】

【００８４】
ここで、ｐ、ｑは、それぞれ
【００８５】
【数３８】

【００８６】
【数３９】

【００８７】
である。式（３６）、式（３７）から、全過程での成長係数ξは常に
【００８８】
【数４０】

【００８９】
となる。これは、いかなる場合も安定で、しかもエネルギー散逸がないことを示している。
次に、本発明の電磁波解析装置のハードウェア構成について説明する。
【００９０】
図４は、電磁波解析装置のハードウェア構成を示す図である。この例では、コンピュータに電磁波解析プログラムを実行させることで、電磁波解析装置を実現している。
【００９１】
この電磁波解析装置は、ＣＰＵ(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 」に記されている。セルサイズは、損失性誘電体膜２１，２２及びスリット２１ａ，２２ａを構成する部分が０．１×０．１ｍｍ² 、それ以外は基本的に１．０×１．０ｍｍ² としている。
【００９９】
図５の左端終端部手前の電界垂直成分をガウス波形で励振し、右端終端部手前の電界垂直成分の波形を観測して、励振波と観測波とのフーリエ変換から透過特性を算出した。なお、時間離散間隔Δｔは、本発明の手法を用いる場合は、時間波形をフーリエ変換する際のサンプリング定理による制約から３．３３ｐｓ、通常のＦＤＴＤ法を用いる場合は、ＣＦＬ条件の制約から０．２３ｐｓとした。
【０１００】
図６は、第１のシミュレーションモデルにおける計算結果を示す図である。この図は、横軸が電磁波の周波数（ＧＨｚ）であり、縦軸は導体損失（ｄＢ）である。また、図中に点線で示した「ADI FDTD」が本発明による計算結果であり、実線で示した「Normal FDTD 」が通常のＦＤＴＤ法による計算結果である。図６から分かるように、「ADI FDTD」と「Normal FDTD 」との結果は、良好に一致している。
【０１０１】
ここで、計算コストに関するデータを以下の表に示す。
【０１０２】
【表１】

【０１０３】
なお、この計算は、一般的なワークステーションを用いて行った。また、表中の［ ］は、比率を表している。この比率は、「時間離散間隔( Δｔ) 」と「メモリー容量(Memory)」については「Normal FDTD 」を基準とし、「ステップ数(Step)」と「演算時間(CPU Time)」については、「ADI FDTD」を基準としている。
【０１０４】
比較の結果、「ADI FDTD」では「Normal FDTD 」に比べて、１４．５倍の時間離散間隔を設定できる。そのため、ステップ数は、１／１４．５ですむ。その結果、演算時間が１／３．８に短縮できた。メモリー容量は、若干増加し、１．３倍となったが、これは、電界成分を記憶する配列数の増加によるものである。 次に、さらにセルサイズを小さくした場合のシミュレーション結果について説明する。
【０１０５】
図７は、第２のシミュレーションモデルを示す図である。この平行平板導波路３０の損失性誘電体膜３１，３２は、１９．６４ｍｍの間隔で２箇所に配置されている。損失性誘電体膜３１，３２は幅０．３ｍｍのスリット３１ａ，３２ａを有し、厚さは０．１ｍｍ、比誘電率２．０、比透磁率２．０、電気伝導率７４．９Ｓ／ｍである。セルサイズは、損失性誘電体膜３１，３２のスリット３１ａ，３２ａを構成する部分は０．０２×０．１ｍｍ² 、それ以外は基本的に１．０×１．０ｍｍ² である。上下の平行平板３３，３４の表面には電気壁境界条件を設定し、導波路の終端３５，３６にはＭｕｒの吸収境界条件を設定した。励起条件は第１のシミュレーションモデルと同様である。時間離散間隔Δｔは、本発明の手法を用いる場合には３．３３ｐｓとし、通常のＦＤＴＤ法を用いる場合には、ＣＦＬ条件の制約から０．０６５ｐｓとした。
【０１０６】
図８は、第２のシミュレーションの計算結果を示す図である。第２のシミュレーションでは、第１のシミュレーションの結果に比べて定量的な差違が若干増加したが、周波数特性としては良好に一致している。
【０１０７】
第２のシミュレーションの計算コストに関するデータを以下の表に示す。
【０１０８】
【表２】

【０１０９】
「ADI FDTD」では、「Normal FDTD 」に比べてΔｔを５１．２倍の大きさで設定できる。そのため、時間ステップ数が１／５１．２ですむ。その結果、演算時間が１／１４に短縮できた。また、使用メモリ容量の増加は、第１のシミュレーションの場合と同様である。
【０１１０】
ここで、第２のシミュレーションの方が第１のシミュレーションよりも演算時間の短縮効果が高いことが分かる。すなわち、セルサイズが小さくなるほど、本発明の手法を適用することによる処理の高速化の効果が得られる。
【０１１１】
ところで、上記の実施の形態では２次元ＴＥ波の解析を行うための電磁波解析装置について説明したが、同様にして、２次元ＴＭ(Transverse Magnetic) 波や３次元波の伝搬を解析することもできる。なお、ＴＭ波とは、磁界ベクトルがすべての場所で波の進行方向に垂直な電磁波である。
【０１１２】
２次元ＴＭ波の伝搬を解析する場合には、格子点（ｉ，ｊ）を基準とした格子に関して、（ｉ，ｊ＋１／２）のｘ方向の磁界Ｈ_i,j+1/2 、（ｉ＋１／２，ｊ）のｙ方向の磁界Ｈ_i+1/2,j 、及び（ｉ，ｊ）のｚ方向の電界Ｅ_i,j を求める。そして、２次元ＴＭ波の伝搬は以下のように定式化される。
【０１１３】
時刻「（ｎ＋１／２）Δｔ」における電磁波の状態を求める前半過程は、
【０１１４】
【数４１】

【０１１５】
【数４２】

【０１１６】
【数４３】

【０１１７】
であり、時刻「（ｎ＋１）Δｔ」における電磁波の状態を求める後半過程は、
【０１１８】
【数４４】

【０１１９】
【数４５】

【０１２０】
【数４６】

【０１２１】
である。式（４１）、式（４５）は簡単な代入計算で解き、式（４２）と式（４４）とは、２次元ＴＥ波の場合と同様に３項方程式を解くことで、前半過程及び後半過程の電磁波の伝搬を解析する。さらに、式（４３）と式（４６）とは、それぞれ式（４２）、式（４４）の計算結果を用いて、簡単な代入式で解くことができる。
【０１２２】
また、３次元の電磁波の伝搬を解析する場合には、ｘ−ｙ−ｚの３次元空間内の格子点（ｉ，ｊ，ｋ）からのｘ方向の空間離散間隔をΔｘ、ｙ方向の空間離散間隔をΔｙ、ｚ方向の空間離散間隔をΔｚとしたとき、格子点（ｉ，ｊ，ｋ）を基準とした格子に関して、（ｉ＋１／２，ｊ，ｋ）のｘ方向の電界Ｅ_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} を求める。３次元の電磁波は２次元のＴＥ波とＴＭ波の重ね合わせで表現できるため、以下のように定式化される。
【０１２３】
時刻「（ｎ＋１／２）Δｔ」における電磁波の状態を求める前半過程は、
【０１２４】
【数４７】

【０１２５】
【数４８】

【０１２６】
【数４９】

【０１２７】
【数５０】

【０１２８】
【数５１】

【０１２９】
【数５２】

【０１３０】
であり、時刻「（ｎ＋１）Δｔ」における電磁波の状態を求める後半過程は、
【０１３１】
【数５３】

【０１３２】
【数５４】

【０１３３】
【数５５】

【０１３４】
【数５６】

【０１３５】
【数５７】

【０１３６】
【数５８】

【０１３７】
である。前半過程においては、式（４７）と式（５１）、式（４８）と式（５２）、及び式（４９）と式（５０）の間で漸化式を導出する。また、後半過程においては、式（５３）と式（５８）、式（５４）と式（５６）、及び式（５５）と式（５７）の間で漸化式を導出する。導出された各漸化式は、ＴＥ波と同様に３項方程式に帰着する。従って、３次元波の電磁波の伝搬も同様にシミュレーションすることができる。
【０１３８】
なお、電磁波解析装置が有すべき機能の処理内容は、コンピュータで読み取り可能な記録媒体に記録されたプログラムに記述しておくことができる。そして、このプログラムをコンピュータで実行することにより、上記処理がコンピュータで実現される。コンピュータで読み取り可能な記録媒体としては、磁気記録装置や半導体メモリ等がある。市場に流通させる場合には、ＣＤ−ＲＯＭ(Compact Disk Read Only Memory) やフロッピーディスク等の可搬型記録媒体にプログラムを格納して流通させたり、ネットワークを介して接続されたコンピュータの記憶装置に格納しておき、ネットワークを通じて他のコンピュータに転送することもできる。コンピュータで実行する際には、コンピュータ内のハードディスク装置等にプログラムを格納しておき、メインメモリにロードして実行する。
【０１３９】
以下に、本発明の特徴を項ごとに区分し、その発明を特定するために必要と認める事項を列挙する。
（第１項） マックスウェル方程式を時間と空間とについて差分法で解くことにより電磁波の伝搬を解析する電磁波解析装置において、
時間離散間隔をΔｔとしたとき、第１の時刻における電磁波の状態に基づいて、前記第１の時刻からΔｔ経過後の第２の時刻における電磁波の状態を交互方向陰解法型の有限差分時間領域法により求めることで、電磁波の挙動を解析する解析手段、
を有することを特徴とする電磁波解析装置。
【０１４０】
（第２項） 前記解析手段は、２次元ＴＥ波の伝搬を解析することを特徴とする第１項記載の電磁波解析装置。
（第３項） 前記解析手段は、電界の振動面をｘ−ｙ平面とし、前記ｘ−ｙ平面の任意の格子点（ｉ，ｊ）を基準とした差分格子の電磁波を解析する際には、格子点（ｉ，ｊ）からｘ方向への空間離散間隔をΔｘ、格子点（ｉ，ｊ）からｙ方向への空間離散間隔をΔｙとした場合に、（ｉ＋１／２，ｊ）におけるｘ方向の電界、（ｉ，ｊ＋１／２）におけるｙ方向の電界、及び（ｉ＋１／２，ｊ＋１／２）における前記ｘ−ｙ平面に垂直なｚ方向の磁界を求めることを特徴とする第２項記載の電磁波解析装置。
【０１４１】
（第４項） 前記解析手段は、ｘ方向の電界に対して陽解法を用い、ｙ方向の電界、及びｚ方向の磁界に対して陰解法を用いることにより、前記第１の時刻における電磁波の状態に基づいて前記第１の時刻からΔｔ／２経過後の中間時刻における電磁波の状態を演算する前半過程演算手段と、ｙ方向の電界に対して陽解法を用い、ｘ方向の電界、及びｚ方向の磁界に対して陰解法を用いることにより、前記前半過程演算手段により算出された前記中間時刻の電磁波の状態に基づいて前記第２の時刻における電磁波の状態を演算する後半過程演算手段と、を具備することを特徴とする第３項記載の電磁波解析装置。
【０１４２】
（第５項） 前記前半過程演算手段は、シミュレーションにおけるｎステップ目の時刻を前記第１の時刻、ｎ＋１ステップ目の時刻を前記第２の時とし、前記第１の時刻における（ｉ＋１／２，ｊ）でのｘ方向の電界をＥｘ_i+1/2,j ⁿ 、前記第１の時刻における（ｉ，ｊ＋１／２）でのｙ方向の電界をＥｙ_i,j+1/2 ⁿ 、前記第１の時刻における（ｉ＋１／２，ｊ＋１／２）でのｚ方向の磁界をＨｚ_i+1/2,j+1/2 ⁿ 、前記中間時刻における（ｉ＋１／２，ｊ）でのｘ方向の電界をＥｘ_i+1/2,j ^n+1/2 、前記中間時刻における（ｉ，ｊ＋１／２）でのｙ方向の電界をＥｙ_i,j+1/2 ^n+1/2 、前記中間時刻における（ｉ＋１／２，ｊ＋１／２）でのｚ方向の磁界をＨｚ_i+1/2,j+1/2 ^n+1/2 、解析すべき領域の誘電率をε、解析すべき領域の透磁率をμとした場合に、格子状に分割された前記ｘ−ｙ平面の任意の格子点（ｉ，ｊ）を基準とした電磁波を
【０１４３】
【数５９】

【０１４４】
【数６０】

【０１４５】
【数６１】

【０１４６】
と定義することで、前記中間時刻におけるすべての格子点を基準とした電磁波の状態を演算し、
前記後半過程演算手段は、前記第２の時刻における（ｉ＋１／２，ｊ）でのｘ方向の電界をＥｘ_i+1/2,j ⁿ⁺¹ 、前記第２の時刻における（ｉ，ｊ＋１／２）でのｙ方向の電界をＥｙ_i,j+1/2 ⁿ⁺¹ 、前記第２の時刻における（ｉ＋１／２，ｊ＋１／２）でのｚ方向の磁界をＨｚ_i+1/2,j+1/2 ⁿ⁺¹ とした場合に、格子状に分割された前記ｘ−ｙ平面の任意の格子点（ｉ，ｊ）を基準とした電磁波を
【０１４７】
【数６２】

【０１４８】
【数６３】

【０１４９】
【数６４】

【０１５０】
と定義することで、前記第２の時刻におけるすべての格子点を基準とした電磁波の状態を演算することを特徴とする第４項記載の電磁波解析装置。
（第６項） 前記解析手段は、２次元ＴＭ波の伝搬を解析することを特徴とする第１項記載の電磁波解析装置。
【０１５１】
（第７項） 前記解析手段は、磁界の振動面をｘ−ｙ平面とし、前記ｘ−ｙ平面の任意の格子点（ｉ，ｊ）を基準とした差分格子の電磁波を解析する際には、格子点（ｉ，ｊ）からｘ方向への空間離散間隔をΔｘ、格子点（ｉ，ｊ）からｙ方向への空間離散間隔をΔｙとした場合に、（ｉ，ｊ＋１／２）におけるｘ方向の磁界、（ｉ＋１／２，ｊ）におけるｙ方向の磁界、及び（ｉ，ｊ）における前記ｘ−ｙ平面に垂直なｚ方向の電界を求めることを特徴とする第６項記載の電磁波解析装置。
【０１５２】
（第８項） 前記解析手段は、ｘ方向の磁界に対して陽解法を用い、ｙ方向の磁界、及びｚ方向の電界に対して陰解法を用いることにより、前記第１の時刻における電磁波の状態に基づいて前記第１の時刻からΔｔ／２経過後の中間時刻における電磁波の状態を演算する前半過程演算手段と、ｙ方向の磁界に対して陽解法を用い、ｘ方向の磁界、及びｚ方向の電界に対して陰解法を用いることにより、前記前半過程演算手段により算出された前記中間時刻の電磁波の状態に基づいて前記第２の時刻における電磁波の状態を演算する後半過程演算手段と、を具備することを特徴とする第７項記載の電磁波解析装置。
【０１５３】
（第９項） 前記前半過程演算手段は、シミュレーションにおけるｎステップ目の時刻を前記第１の時刻、ｎ＋１ステップ目の時刻を前記第２の時とし、前記第１の時刻における（ｉ，ｊ＋１／２）でのｘ方向の磁界をＨｘ_i,j+1/2 ⁿ 、前記第１の時刻における（ｉ＋１／２，ｊ）でのｙ方向の磁界をＨｙ_i+1/2,j ⁿ 、前記第１の時刻における（ｉ，ｊ）でのｚ方向の電界をＥｚ_i,j ⁿ 、前記中間時刻における（ｉ，ｊ＋１／２）でのｘ方向の磁界をＨｘ_i,j+1/2 ^n+1/2 、前記中間時刻における（ｉ＋１／２，ｊ）でのｙ方向の磁界をＨｙ_i+1/2,j ^n+1/2 、前記中間時刻における（ｉ，ｊ）でのｚ方向の電界をＥｚ_i,j ^n+1/2 、解析すべき領域の誘電率をε、解析すべき領域の透磁率をμとした場合に、格子状に分割された前記ｘ−ｙ平面の任意の格子点（ｉ，ｊ）を基準とした電磁波を
【０１５４】
【数６５】

【０１５５】
【数６６】

【０１５６】
【数６７】

【０１５７】
と定義することで、前記中間時刻におけるすべての格子点を基準とした電磁波の状態を演算し、
前記後半過程演算手段は、前記第２の時刻における（ｉ，ｊ＋１／２）でのｘ方向の磁界をＨｘ_i,j+1/2 ⁿ⁺¹ 、前記第２の時刻における（ｉ＋１／２，ｊ）でのｙ方向の磁界をＨｙ_i+1/2,j ⁿ⁺¹ 、前記第２の時刻における（ｉ，ｊ）でのｚ方向の電界をＥｚ_i,j ⁿ⁺¹ とした場合に、格子状に分割された前記ｘ−ｙ平面の任意の格子点（ｉ，ｊ）を基準とした電磁波を
【０１５８】
【数６８】

【０１５９】
【数６９】

【０１６０】
【数７０】

【０１６１】
と定義することで、前記第２の時刻におけるすべての格子点を基準とした電磁波の状態を演算することを特徴とする第８項記載の電磁波解析装置。
（第１０項） 前記解析手段は、３次元波の伝搬を解析することを特徴とする第１項記載の電磁波解析装置。
【０１６２】
（第１１項） 前記解析手段は、ｘ−ｙ−ｚの３次元空間を想定し、３次元空間の任意の格子点（ｉ，ｊ，ｋ）を基準とした差分格子の電磁波を解析する際には、格子点（ｉ，ｊ，ｋ）からｘ方向への空間離散間隔をΔｘ、格子点（ｉ，ｊ，ｋ）からｙ方向への空間離散間隔をΔｙ、格子点（ｉ，ｊ，ｋ）からｚ方向への空間離散間隔をΔｚとした場合に、（ｉ＋１／２，ｊ，ｋ）におけるｘ方向の電界、（ｉ，ｊ＋１／２，ｋ）におけるｙ方向の電界、（ｉ，ｊ，ｋ＋１／２）におけるｚ方向の電界、（ｉ，ｊ＋１／２，ｋ＋１／２）におけるｘ方向の磁界、（ｉ＋１／２，ｊ，ｋ＋１／２）におけるｙ方向の磁界、及び（ｉ＋１／２，ｊ＋１／２，ｋ）におけるｚ方向の磁界を求めることを特徴とする第１０項記載の電磁波解析装置。
【０１６３】
（第１２項） 前記解析手段は、
シミュレーションにおけるｎステップ目の時刻を前記第１の時刻、ｎ＋１ステップ目の時刻を前記第２の時とし、前記第１の時刻における（ｉ＋１／２，ｊ，ｋ）でのｘ方向の電界をＥｘ_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} ⁿ 、前記中間時刻における（ｉ＋１／２，ｊ，ｋ）でのｘ方向の電界をＥｘ_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 、解析すべき領域の誘電率をε、解析すべき領域の透磁率をμとした場合に、格子状に分割された前記３次元空間の任意の格子点（ｉ，ｊ，ｋ）を基準とした電磁波を
【０１６４】
【数７１】

【０１６５】
【数７２】

【０１６６】
【数７３】

【０１６７】
【数７４】

【０１６８】
【数７５】

【０１６９】
【数７６】

【０１７０】
と定義することで、前記中間時刻におけるすべての格子点を基準とした電磁波の状態を演算する前半過程演算手段と、
前記中間時刻における（ｉ＋１／２，ｊ，ｋ）でのｘ方向の電界をＥｘ_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} ⁿ⁺¹ とした場合に、格子状に分割された前記３次元空間の任意の格子点（ｉ，ｊ、ｋ）を基準とした電磁波を
【０１７１】
【数７７】

【０１７２】
【数７８】

【０１７３】
【数７９】

【０１７４】
【数８０】

【０１７５】
【数８１】

【０１７６】
【数８２】

【０１７７】
と定義することで、前記第２の時刻におけるすべての格子点を基準とした電磁波の状態を演算する前記後半過程演算手段と、
を具備することを特徴とする第１１項記載の電磁波解析装置。
【０１７８】
（第１３項） マックスウェル方程式を時間と空間とについて差分法で解くことにより電磁波の伝搬を解析するための電磁波解析プログラムを記録したコンピュータ読み取り可能な記録媒体において、
時間離散間隔をΔｔとしたとき、第１の時刻における電磁波の状態に基づいて、前記第１の時刻からΔｔ経過後の第２の時刻における電磁波の状態を交互方向陰解法型の有限差分時間領域法により求めることで、電磁波の挙動を解析する解析手段、
としてコンピュータを機能させることを特徴とする電磁波解析プログラムを記録したコンピュータ読み取り可能な記録媒体。
【０１７９】
【発明の効果】
以上説明したように本発明の電磁波解析装置では、電磁波の状態を交互方向陰解法型の有限差分時間領域法により求めることで、電磁波の挙動を解析するようにしたため、空間離散間隔に依存せずに時間離散間隔を定めることが可能となる。その結果、微細な構造を有する解析対象物の電磁波解析を行う際の時間離散間隔を比較的大きくとることが可能となり、解析処理を高速に行うことができる。
【０１８０】
また、本発明の電磁波解析プログラムを記録したコンピュータ読み取り可能な記録媒体では、記録された電磁波解析プログラムをコンピュータに実行させることにより、電磁波の状態を交互方向陰解法型の有限差分時間領域法による電磁波の解析をコンピュータで高速に行うことが可能となる。
【図面の簡単な説明】
【図１】本発明の原理構成図である。
【図２】格子の拡大図である。
【図３】ＡＤＩ−ＦＤＴＤ法による２次元ＴＥ波解析の処理手順を示すフローチャートである。
【図４】電磁波解析装置のハードウェア構成を示す図である。
【図５】第１のシミュレーションモデルを示す図である。
【図６】第１のシミュレーションモデルにおける計算結果を示す図である。
【図７】第２のシミュレーションモデルを示す図である。
【図８】第２のシミュレーションの計算結果を示す図である。
【符号の説明】
１ 解析指令入力手段
２ 解析手段
２ａ 前半過程演算手段
２ｂ 後半過程演算手段
３ 解析結果出力手段[0001]
BACKGROUND OF THE INVENTION
The present invention relates to an electromagnetic wave analysis apparatus that analyzes the behavior of electromagnetic waves using an electronic computer and a computer-readable recording medium that records an electromagnetic wave analysis program, and in particular, an electromagnetic wave analysis apparatus and an electromagnetic wave analysis that perform analysis using a finite difference time domain method The present invention relates to a computer-readable recording medium on which a program is recorded.
[0002]
[Prior art]
One of the methods for analyzing the transient behavior of electromagnetic waves by numerical simulation using an electronic computer is a finite difference time domain (FDTD) method. This method is a method of solving the Maxwell equation with respect to time and space by a difference method, and is a method widely used because of its wide range of application.
[0003]
By the way, the conventional finite difference time domain method has been solved by an explicit method type difference method. Below, the explicit method will be described. For example, consider the following one-dimensional parabolic partial differential equation for simplicity.
[0004]
[Expression 7]

[0005]
[Equation 8]

[0006]
Here, n is the number of steps advanced in the time axis direction from the calculation start time, Δt is a time discrete interval (time advanced in one step), and i is an order from the origin of the difference grid (cell). Yes, Δx is the spatial discrete interval (width of the divided region). When solving this differential equation with the explicit method, approximation with a forward difference to the time derivative,
[0007]
[Equation 9]

[0008]
The following formula is obtained. This equation (9) is r = Δt / Δx²Anyway
[0009]
[Expression 10]

[0010]
Can be written. As a result, the solution at the (n + 1) th step can be immediately calculated from the solution at the past time (nth step) by substitution calculation.
[0011]
[Problems to be solved by the invention]
However, in the explicit difference type method of number (10), as a stable condition for the solution not to diverge
[0012]
## EQU11 ##

[0013]
Must be met.
In the FDTD method, the CFL (Courant, Friedrich and Levy) condition or the Courant condition must be satisfied so that the solution does not diverge.
[0014]
The CFL condition is that the solution does not diverge
[0015]
[Expression 12]

[0016]
It is a condition that the stability condition of must be satisfied. Details of the CFL conditions of the FDTD method are described in “A. Taflove, Computational Electrodynamics, MA, Artech, 1995”. Equation (12) is a Courant condition for analyzing a three-dimensional wave. Where v is the velocity of the electromagnetic wave and Δx_min, Δy_min, Δz_minAre the minimum values of the spatial discrete intervals in the XYZ directions.
[0017]
Due to these restrictions, the explicit method has the property that the maximum value of the time discrete interval is determined by the minimum value of the spatial discrete interval. For this reason, when the time response of the analysis target having a fine structure is calculated for a long time, there is a problem that the number of repetitions of the calculation for advancing the simulation time increases and the time required for the simulation increases.
[0018]
Note that the partial differential equations shown in Equations (7) and (8) can also be solved by an implicit method. For example, when solving a one-dimensional parabolic partial differential equation such as Equation (7) and Equation (8) by implicit method, approximation with a backward difference for time differentiation,
[0019]
[Formula 13]

[0020]
Get. This equation (13) is expressed as r = Δt / Δx as in the explicit method.²Anyway
[0021]
[Expression 14]

[0022]
Can be written. Solving such equation (14) is an implicit method. Unlike the explicit method, this implicit method is always stable.
However, the implicit method must solve the following simultaneous equations f_i ^{n + 1}Can't get the value of
[0023]
[Expression 15]

[0024]
Where a_iIs the constant part (−r) of the third term on the left side when “i = 1, 2,..., Imax” in equation (14), b_iIs the constant part (1 + 2r) of the second term on the left side when “i = 1, 2,..., Imax” in equation (14), and c_iIs the constant part (−r) of the first term on the left side when “i = 1, 2,..., Imax” in equation (14).
[0025]
In order to analyze a transient phenomenon of electromagnetic wave analysis, it is necessary to solve at least a two-dimensional partial differential equation, but separate attention is required for an implicit method of a two-dimensional or higher partial differential equation. For example, a two-dimensional parabolic partial differential equation is expressed by the following equation.
[0026]
[Expression 16]

[0027]
When solving this equation (16) by the Crank-Nicolson method, for example,
[0028]
[Expression 17]

[0029]
The time derivative is approximated as follows. When solving the equation (17) by simultaneous equations, it is necessary to solve the simultaneous linear equations of (number of lattices in the x direction−1) × (number of lattices in the y direction−1). Therefore, when the number of grids is large, the calculation amount and the necessary storage capacity become enormous. That is, if the lattice is refined, the calculation takes a very long time as in the case of analysis by the explicit method.
[0030]
The present invention has been made in view of such a point, and an object thereof is to provide an electromagnetic wave analysis apparatus capable of performing an electromagnetic wave analysis of an analysis object having a fine structure at high speed.
[0031]
Another object of the present invention is to provide a computer-readable recording medium on which an electromagnetic wave analysis program capable of causing a computer to perform electromagnetic wave analysis of an analysis object having a fine structure at high speed is recorded.
[0032]
[Means for Solving the Problems]
In the present invention, in order to solve the above-mentioned problem, in the electromagnetic wave analysis apparatus for analyzing the propagation of electromagnetic waves by solving the Maxwell equation with respect to time and space by the difference method, the first time is set when the time discrete interval is Δt. Analysis of electromagnetic wave behavior by obtaining the state of the electromagnetic wave at the second time after the lapse of Δt from the first time by the alternating direction implicit finite difference time domain method based on the state of the electromagnetic wave at There is provided an electromagnetic wave analyzing apparatus characterized by comprising means.
[0033]
According to such an electromagnetic wave analysis apparatus, the state of the electromagnetic wave at the second time after the lapse of Δt from the first time is based on the alternating direction implicit method based on the state of the electromagnetic wave at the first time. It is obtained by the finite difference time domain method. Thus, the solution obtained by the alternating direction implicit method type finite difference time domain method is stable regardless of the size of the time discrete interval.
[0034]
In order to solve the above problem, a computer-readable recording medium recording an electromagnetic wave analysis program for analyzing propagation of electromagnetic waves by solving the Maxwell equation with respect to time and space by a difference method, a time discrete interval Is set to Δt, the state of the electromagnetic wave at the second time after the lapse of Δt from the first time is obtained by an alternating direction implicit finite difference time domain method based on the state of the electromagnetic wave at the first time. Thus, a computer-readable recording medium recording an electromagnetic wave analysis program characterized by causing a computer to function as an analysis means for analyzing electromagnetic wave behavior is provided.
[0035]
When the computer executes the electromagnetic wave analysis program recorded on the recording medium, the function of the electromagnetic wave analysis apparatus according to the present invention is realized on the computer.
[0036]
DETAILED DESCRIPTION OF THE INVENTION
Hereinafter, embodiments of the present invention will be described with reference to the drawings. The electromagnetic wave analysis apparatus of the present invention enables high-speed computation while ensuring solution stability by analyzing electromagnetic waves using an ADI (Alternating Direction Implicit) type FDTD method.
[0037]
FIG. 1 is a principle configuration diagram of the present invention. The analysis command input means 1 designates initial values and boundary absorption conditions necessary for electromagnetic wave analysis and inputs an analysis command.
When the analysis unit 2 receives the analysis command from the analysis command input unit 1, the analysis unit 2 numerically analyzes the propagation of the two-dimensional TE (Transverse Electric) wave by the alternating direction implicit method-differential time domain (ADI-FDTD) method. Here, the TE wave is an electromagnetic wave expressed by the following equation.
[0038]
[Expression 18]

[0039]
[Equation 19]

[0040]
[Expression 20]

[0041]
The ADI-FDTD method is a method in which the time discrete interval is divided into the first half and the second half, and the first half process and the second half process are alternately performed. For example, when solving the two-dimensional parabolic partial differential equation of Equation (16) by the ADI method,
[0042]
[Expression 21]

[0043]
[Expression 22]

[0044]
It approximates as follows. Equation (21) is the first half process, and the time derivative on the left side is approximated by a backward difference when viewed from f in the x direction, and approximated by a forward difference when viewed from f in the y direction. Equation (22) is the latter half process, and the time derivative on the left side is approximated by a forward difference when viewed from the x direction f, and approximated by a backward difference when viewed from f in the y direction. The first half process is calculated by the first half process calculation means 2a, and the second half process is calculated by the second half process calculation means 2b. Note that Δt is a time discrete interval, Δx is a spatial discrete interval in the x direction, and Δy is a spatial discrete interval in the y direction.
[0045]
In the first half process calculating means 2a, the time “nΔt” from the second half process calculating means 2b (where n is the number of steps of the time progression process from the start of the simulation, and the first half process and the second half process are each advanced by 1/2 step). When information on the state of the electromagnetic wave is received, the state of the electromagnetic wave at time “(n + 1/2) Δt” is calculated. At that time, the electric field Ex in the x direction at the time “(n + 1/2) Δt”.^{n + 1/2}Is obtained by the explicit method, and the electric field Ey in the y direction^{n + 1/2}And z-direction magnetic field Hz^{n + 1/2}Is obtained by the implicit method. The calculation result is delivered to the second half process calculation means 2b.
[0046]
When the second half process calculating means 2b receives the state of the electromagnetic wave at the time “(n + 1/2) Δt” from the first half process calculating means 2a, the second half process calculating means 2b calculates the state of the electromagnetic wave at the time “(n + 1) Δt”. At that time, the electric field Ey in the y direction at time “(n + 1) Δt”^{n + 1}Is obtained by an explicit solution, and the electric field Ex in the x direction^{n + 1}And z-direction magnetic field Hz^{n + 1}Is obtained by the implicit method. The calculation result is delivered to the first half process calculation means 2a.
[0047]
The analysis result output unit 3 outputs the analysis result obtained by the analysis unit 2.
The calculation contents of the analysis means 2 when determining the state of the time “(n + 1) Δt” from the state of the electromagnetic wave at the time “nΔt” will be described below.
[0048]
First, the analysis target area is divided into rectangular grids. Here, it is assumed that the x direction of the region to be analyzed is divided into imax lattices and the y direction is divided into jmax lattices.
[0049]
FIG. 2 is an enlarged view of the lattice. The position of the grid point is represented by two sets of numbers such as (i, j). i is the order of the grid in the x direction from the origin. j is the order of the grid in the y direction from the origin. Further, an intermediate between “i” and “i + 1” is represented as “i + 1/2”. Here, when analyzing the two-dimensional TE wave in the region with the lattice point (i, j) as a reference, the electric field in the x direction at (i + 1/2, j) and the y direction in (i, j + 1/2) The time change of the electric field and the magnetic field in the z direction at (i + 1/2, j + 1/2) is analyzed by the difference method.
[0050]
When solving Maxwell's equations for time and space using the finite difference method,_xIs approximated by a forward difference (shown in equation (23) below), and E_yIs approximated by a backward difference (shown in equation (24) below), and H_zIs the time derivative of E_xFrom the above, it is approximated by forward difference and E_yFrom the above, it is approximated by the backward difference (shown in the following equation (25)). As a result, the first half of the process for obtaining the state of time “(n + 1/2) Δt” can be expressed by the following equation.
[0051]
[Expression 23]

[0052]
[Expression 24]

[0053]
[Expression 25]

[0054]
Here, the right shoulder values of the electric field E and the magnetic field H represent the time step of the simulation. The lower right values of the electric field E and the magnetic field H are lattice numbers when the analysis region is divided by the lattice. ε is the dielectric constant of the object to be analyzed. μ is the magnetic permeability of the object to be analyzed.
[0055]
By the way, Expression (23) is explicit expression, and Expression (24) and Expression (25) are implicit expression. In other words, in the first half process, E_xIs approximated by a forward difference, the electric field Ex in the x direction at the time “(n + 1/2) Δt”.^{n + 1/2}Can be obtained by explicit method. In the first half of the process, E_yIs approximated by a backward difference, the electric field Ey in the y direction at time “(n + 1/2) Δt”.^{n + 1/2}Must be obtained by implicit method. In addition, H_zThe time derivative of E_ySince it is approximated by a backward difference when viewed from, the magnetic field Hz in the z direction at time “(n + 1/2) Δt”^{n + 1/2}Must be obtained by implicit method.
[0056]
Similarly, the latter half of the process for obtaining the state of time “(n + 1) Δt” is E_xIs approximated by a backward difference (shown in equation (26) below), and E_yIs approximated by a forward difference (shown in equation (27) below) and H_zIs the time derivative of E_xFrom the above, it is approximated by the backward difference and E_yFrom the above, it is approximated by a forward difference (shown in the following equation (28)). As a result, the latter half of the process can be expressed by the following equation.
[0057]
[Equation 26]

[0058]
[Expression 27]

[0059]
[Expression 28]

[0060]
In the latter half of the process, the expression (27) for obtaining the electric field in the y direction is explicit expression, and the other expressions (26) and (28) are implicit expression. In other words, E_xIs approximated by a backward difference, the electric field Ex in the x direction at time “(n + 1) Δt”^{n + 1}Must be obtained by implicit method. E_y, The electric field Ey in the y direction at time “(n + 1) Δt”^{n + 1}Can be obtained by explicit method. In addition, H_zThe time derivative of E_xSince it is approximated by a backward difference when viewed from, the magnetic field Hz in the z direction at time “(n + 1) Δt”^{n + 1}Must be obtained by implicit method.
[0061]
The electric field in the x direction at the time “(n + 1/2) Δt” can be obtained by substituting the electric field and the magnetic field at the time “nΔt” according to the expression (23) in the first half process. Similarly, equation (27) in the latter half of the process can be obtained by simple substitution calculation.
[0062]
On the other hand, the expressions (24) and (25) in the first half process cannot be solved as they are because the right side includes the electric field and magnetic field at time “(n + 1/2) Δt”. Therefore, the equation (24) and the equation (25) are combined and the magnetic field Hz^{n + 1/2}To get the following recurrence formula.
[0063]
[Expression 29]

[0064]
Since the expressions (26) and (28) in the latter half of the process also include an electric field and a magnetic field at time “(n + 1) Δt” on the right side, the expressions (26) and (28) are combined to generate a magnetic field Hz.^{n + 1}Is eliminated, and the following recurrence formula is obtained.
[0065]
[30]

[0066]
If the recurrence formula of the equation (29) is generated for each of the cases where “i = 1,..., Imax”, it results in the following simultaneous linear equations.
[0067]
[31]

[0068]
Where a in equation (31)_iIs the constant part of the first term on the left side of equation (29) obtained when “i = 1,..., Imax”. b_iIs the constant part of the second term on the left side of equation (29) obtained when “i = 1,..., Imax”. c_iIs the constant part of the third term on the left side of equation (29) obtained when “i = 1,..., Imax”. D_i ⁿIs the value on the right side of Equation (29) obtained when “i = 1,..., Imax”. In addition, the absorption boundary condition is applied to the first row and the last row.
[0069]
Equation (31) is a tridiagonal system of equation for the lattice point “i” in the x direction, and can be solved by the Thomas method or the like. If such a ternary equation is solved for each of “i = 1,..., Imax”, the electric field in the y direction in all the lattices at time “n + 1/2” is obtained. If the electric field in the y direction is obtained, the magnetic field in the z direction can also be obtained from Equation (16).
[0070]
Equation (30), which is a recurrence formula for the second half process, can also obtain a ternary equation for “j” as follows, similarly to the first half process.
[0071]
[Expression 32]

[0072]
Where a in equation (32)_jIs the constant part of the first term on the left side of equation (30) obtained when “j = 1,..., Jmax”. b_jIs the constant part of the second term on the left side of Equation (30) obtained when “j = 1,..., Jmax”. c_jIs the constant part of the third term on the left side of Equation (30) obtained when “j = 1,..., Jmax”. D_j ^{n + 1/2}Is the value on the right side of Equation (30) obtained when “j = 1,..., Jmax”. In addition, the absorption boundary condition is applied to the first row and the last row.
[0073]
In this way, the respective calculations of the expression (24) in the first half process and the expression (26) in the second half process can be reduced to a ternary equation. Therefore, in the first half process calculation means 2a, the ternary equation of Expression (23) and Expression (31) is solved, and further, Expression (25) is solved. In the second half process calculation means 2b, the ternary equations of the equations (27) and (32) are solved, and then the equation (28) is solved. If the time step is advanced by alternately performing such calculation, the propagation of the two-dimensional TE wave can be analyzed by the ADI-FDTD method.
[0074]
FIG. 3 is a flowchart showing a processing procedure of two-dimensional TE wave analysis by the ADI-FDTD method. This process will be described along with step numbers. It is assumed that various initial conditions such as absorption boundary conditions, time discrete intervals, and spatial discrete intervals have already been given. At the start of the simulation, the simulation time is “T = 0”.
[S1] The first half process calculation means 2a first calculates E which has already been obtained._x ⁿ, H_z ⁿBy substituting the value of and solving equation (23), E_x ^{n + 1/2}Calculate Next, E that has already been sought_x ⁿ, E_y ⁿ, H_z ⁿBy solving the ternary equation of equation (31) using the value of_y ^{n + 1/2}Calculate Furthermore, E obtained by solving equation (31)_y ^{n + 1/2}Value and known E_x ⁿ, H_z ⁿBy solving the equation (25) using the value of_z ^{n + 1/2}Calculate As a result, the simulation time T is advanced by ΔT / 2.
[S2] The second half process calculation means 2b first calculates E which has already been obtained._x ^{n + 1/2}, E_y ^{1 + 1/2}, H_z ^{n + 1/2}By solving the ternary equation of equation (32) using the value of_x ^{n + 1}Calculate Next, E that has already been sought_y ^{n + 1/2}, H_z ^{n + 1/2}By solving equation (27) using the value of_y ^{n + 1}Calculate Furthermore, E obtained by solving equation (32)_x ^{n + 1}Value and known E_y ^{n + 1/2}, H_z ^{n + 1/2}By solving the equation (28) using the value of_z ^{n + 1}Calculate As a result, the simulation time T is advanced by ΔT / 2. As a result, the time is advanced by Δt which is a total time discrete interval between step S1 and step S2.
[S3] The analyzing means 2 uses the simulation maximum time T at which the time T is preset._maxIt is determined whether or not it exceeds. T_maxIf it exceeds, the simulation is terminated. If not, the process proceeds to step S1 and the simulation is continued.
[0075]
When the analysis unit 2 executes the above processing, the analysis of the two-dimensional TE wave using the ADI-FDTD method is performed.
For the sake of simplification, the above equation assumes that the electromagnetic wave propagates through a lossless, isotropic, non-dispersive, and homogeneous medium, and the difference grating (cell) is uniform (Δt and Since Δy is equal in all cells)_i"," C_i"," A_j"," C_jThe values of “” are all “1”. If the medium is inhomogeneous or the difference grid (cell) is not uniform, the equation needs to be expanded. The method of extending the above equation to the case where the medium is inhomogeneous or the difference lattice (cell) is not uniform is the same as the expansion method performed in the normal FDTD method.
[0076]
Below, the result of having examined using the Fourier series method about the stability of the solution at the time of analyzing the propagation of a two-dimensional TE wave by the ADI-FDTD method is shown.
First, the growth factor ξ of the error with time development is calculated by using the step number n ′ of time.
[0077]
[Expression 33]

[0078]
It is expressed. Note that γ is a complex constant. Using the error growth coefficient, the solutions of the first half process and the second half process are set as follows.
[0079]
[Expression 34]

[0080]
[Expression 35]

[0081]
Where ξ₁ ^{n '}And ξ₂ ^{n '}Is the growth factor of the error with time development in the first half process and the second half process, respectively. Moreover, since the imaginary unit is misleading as a subscript, it is represented by the square root of -1. Then, these are expressed by the following equation.
[0082]
[Expression 36]

[0083]
[Expression 37]

[0084]
Here, p and q are respectively
[0085]
[Formula 38]

[0086]
[39]

[0087]
It is. From equation (36) and equation (37), the growth coefficient ξ in the whole process is always
[0088]
[Formula 40]

[0089]
It becomes. This indicates that in any case it is stable and there is no energy dissipation.
Next, the hardware configuration of the electromagnetic wave analysis apparatus of the present invention will be described.
[0090]
FIG. 4 is a diagram illustrating a hardware configuration of the electromagnetic wave analysis device. In this example, an electromagnetic wave analysis apparatus is realized by causing a computer to execute an electromagnetic wave analysis program.
[0091]
This electromagnetic wave analyzing apparatus is configured around a CPU (Central Processing Unit) 11. The CPU 11 performs an electromagnetic wave analysis calculation based on the program stored in the memory 12 and controls various devices connected via the bus 17. The peripheral devices connected to the bus 17 are as follows.
[0092]
The display controller 13 generates a display image according to a drawing command sent from the CPU 11 and sends it to the display device 13a. The display device 13a connected to the display controller 13 displays the image on the screen according to the display image information sent from the display controller 13.
[0093]
The input device interface 14 is connected to a keyboard 14a and a mouse 14b, and transfers input signals from the keyboard 14a and the mouse 14b to the CPU 11.
[0094]
The network interface 15 is connected to a LAN (Local Area Network) and controls data communication via the LAN. That is, the data sent from the CPU 11 is transferred to another device on the LAN, and the data sent via the LAN is received and passed to the CPU 11.
[0095]
An HDD (Hard Disk Drive) controller 16 controls input / output of data to / from a connected hard disk device (HDD) 16a. The HDD 16a stores a system program to be executed by the CPU 11 and a program in which processing contents when executing electromagnetic wave analysis processing are described.
[0096]
In the electromagnetic wave analysis apparatus having such a configuration, when performing electromagnetic wave analysis analysis, the user first inputs initial conditions using the keyboard 14a. The initial conditions are a spatial discrete interval between the x direction and the y direction, a time discrete interval, a boundary condition, a range specification of a calculation region, and the like. The input initial condition is transferred from the input device interface 14 to the CPU 11 and held in the memory 12. When the user inputs an execution command for electromagnetic wave analysis processing, the command is sent to the CPU 11. The CPU 11 loads a program for performing electromagnetic wave analysis processing by the FDTD method according to the present invention from the HDD 16a to the memory 12, and executes the program. The execution result is displayed on the screen of the display device 13 a via the display controller 13.
[0097]
A specific numerical calculation example is shown below. In addition, the following calculation was performed by the program which can be calculated even when a medium is inhomogeneous or a difference grid (cell) is not uniform. Such a function expansion method is the same as the method performed in the conventional FDTD method.
[0098]
FIG. 5 is a diagram illustrating a first simulation model. This is a calculation of the transparency when the lossy dielectric films 21 and 22 are inside the parallel plate waveguide 20. The lossy dielectric films 21 and 22 are arranged at two locations at intervals of 24 mm. The lossy dielectric films 21 and 22 have slits 21a and 22a having a width of 0.3 mm, the thickness is 0.5 mm, the relative permittivity is 2.0, the relative permeability is 2.0, and the electrical conductivity is 15.0S. / M. Electric wall boundary conditions were set on the surfaces of the upper and lower parallel plates 23 and 24, and Mur absorption boundary conditions were set on the terminal ends 25 and 26 of the waveguide. Mur's absorption boundary conditions are described in “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. Has been. The cell size is 0.1 × 0.1 mm at the portions constituting the lossy dielectric films 21 and 22 and the slits 21a and 22a.²Other than that, it is basically 1.0 × 1.0mm²It is said.
[0099]
The electric field vertical component in front of the left end terminal portion in FIG. 5 was excited with a Gaussian waveform, the waveform of the electric field vertical component in front of the right end terminal portion was observed, and the transmission characteristics were calculated from the Fourier transform of the excitation wave and the observed wave. Note that the time discrete interval Δt is 3.33 ps due to the sampling theorem when Fourier transforming the time waveform when the method of the present invention is used, and is 0.3 due to the limitation of the CFL condition when the normal FDTD method is used. It was 23 ps.
[0100]
FIG. 6 is a diagram illustrating a calculation result in the first simulation model. In this figure, the horizontal axis represents the frequency (GHz) of the electromagnetic wave, and the vertical axis represents the conductor loss (dB). Further, “ADI FDTD” indicated by a dotted line in the figure is a calculation result according to the present invention, and “Normal FDTD” indicated by a solid line is a calculation result by a normal FDTD method. As can be seen from FIG. 6, the results of “ADI FDTD” and “Normal FDTD” agree well.
[0101]
Here, the data relating to the calculation cost is shown in the following table.
[0102]
[Table 1]

[0103]
This calculation was performed using a general workstation. [] In the table represents the ratio. This ratio is based on “Normal FDTD” for “time discrete interval (Δt)” and “memory capacity”, and for “step number” and “CPU time”, Based on ADI FDTD.
[0104]
As a result of the comparison, “ADI FDTD” can set a time discrete interval 14.5 times that of “Normal FDTD”. Therefore, the number of steps can be 1 / 14.5. As a result, the calculation time could be shortened to 1 / 3.8. The memory capacity slightly increased to 1.3 times, which is due to the increase in the number of arrays for storing electric field components. Next, simulation results when the cell size is further reduced will be described.
[0105]
FIG. 7 is a diagram illustrating a second simulation model. The lossy dielectric films 31 and 32 of the parallel plate waveguide 30 are arranged at two places with an interval of 19.64 mm. The lossy dielectric films 31 and 32 have slits 31a and 32a having a width of 0.3 mm, the thickness is 0.1 mm, the relative permittivity is 2.0, the relative permeability is 2.0, and the electric conductivity is 74.9 S /. m. The cell size is 0.02 × 0.1 mm for the portions constituting the slits 31a and 32a of the lossy dielectric films 31 and 32.²Other than that, basically 1.0 x 1.0mm²It is. Electrical wall boundary conditions were set on the surfaces of the upper and lower parallel plates 33 and 34, and Mur absorption boundary conditions were set on the waveguide ends 35 and 36. Excitation conditions are the same as in the first simulation model. The time discrete interval Δt is set to 3.33 ps when the method of the present invention is used, and is set to 0.065 ps when a normal FDTD method is used due to the restriction of the CFL condition.
[0106]
FIG. 8 is a diagram illustrating a calculation result of the second simulation. In the second simulation, the quantitative difference slightly increased as compared with the result of the first simulation, but the frequency characteristics agree well.
[0107]
Data relating to the calculation cost of the second simulation is shown in the following table.
[0108]
[Table 2]

[0109]
In “ADI FDTD”, Δt can be set 51.2 times larger than “Normal FDTD”. Therefore, the number of time steps is 1 / 51.2. As a result, the calculation time was reduced to 1/14. The increase in the used memory capacity is the same as in the first simulation.
[0110]
Here, it can be seen that the second simulation has a higher calculation time reduction effect than the first simulation. That is, as the cell size decreases, the effect of increasing the processing speed by applying the method of the present invention can be obtained.
[0111]
In the above embodiment, the electromagnetic wave analysis apparatus for analyzing the two-dimensional TE wave has been described. Similarly, the propagation of a two-dimensional TM (Transverse Magnetic) wave or a three-dimensional wave can be analyzed. . The TM wave is an electromagnetic wave whose magnetic field vector is perpendicular to the traveling direction of the wave at all locations.
[0112]
When analyzing the propagation of a two-dimensional TM wave, a magnetic field H in the x direction of (i, j + 1/2) with respect to a lattice with reference to the lattice point (i, j)._{i, j + 1/2}, (I + 1/2, j) y-direction magnetic field H_{i + 1/2, j}, And (i, j) z-direction electric field E_{i, j}Ask for. The propagation of the two-dimensional TM wave is formulated as follows.
[0113]
The first half of the process of obtaining the state of the electromagnetic wave at time “(n + 1/2) Δt”
[0114]
[Expression 41]

[0115]
[Expression 42]

[0116]
[Equation 43]

[0117]
And the latter half of the process for obtaining the state of the electromagnetic wave at time “(n + 1) Δt” is
[0118]
(44)

[0119]
[Equation 45]

[0120]
[Equation 46]

[0121]
It is. Equations (41) and (45) are solved by a simple substitution calculation. Equations (42) and (44) are solved by solving a ternary equation as in the case of a two-dimensional TE wave. Analyze the propagation of electromagnetic waves in the process. Furthermore, the equations (43) and (46) can be solved by simple substitution equations using the calculation results of the equations (42) and (44), respectively.
[0122]
When analyzing the propagation of a three-dimensional electromagnetic wave, the spatial discrete interval in the x direction from the lattice point (i, j, k) in the xyz three-dimensional space is defined as a space in the Δx, y direction. When the discrete interval is Δy and the spatial discrete interval in the z direction is Δz, the electric field E in the x direction of (i + 1/2, j, k) with respect to the lattice with the lattice point (i, j, k) as a reference._{i + 1/2, j, k}, (I, j + 1/2, k) electric field E in the y direction_{i, j + 1/2, k}, (I, j, k + 1/2) z-direction electric field E_{i, j, k + 1/2}, (I, j + 1/2, k + 1/2) x-direction magnetic field H_{i, j + 1/2, k + 1/2}, (I + 1/2, j, k + 1/2) y-direction magnetic field H_{i + 1/2, j, k + 1/2}, And (i + 1/2, j + 1/2, k) z-direction magnetic field H_{i + 1/2, j + 1/2, k}Ask for. Since a three-dimensional electromagnetic wave can be expressed by a superposition of a two-dimensional TE wave and a TM wave, it is formulated as follows.
[0123]
The first half of the process of obtaining the state of the electromagnetic wave at time “(n + 1/2) Δt”
[0124]
[Equation 47]

[0125]
[Formula 48]

[0126]
[Formula 49]

[0127]
[Equation 50]

[0128]
[Formula 51]

[0129]
[Formula 52]

[0130]
And the latter half of the process for obtaining the state of the electromagnetic wave at time “(n + 1) Δt” is
[0131]
[53]

[0132]
[Formula 54]

[0133]
[Expression 55]

[0134]
[Expression 56]

[0135]
[Equation 57]

[0136]
[Formula 58]

[0137]
It is. In the first half process, a recurrence formula is derived between Formula (47) and Formula (51), Formula (48) and Formula (52), and Formula (49) and Formula (50). Further, in the latter half of the process, a recurrence formula is derived between Formula (53) and Formula (58), Formula (54) and Formula (56), and Formula (55) and Formula (57). Each derived recurrence formula results in a ternary equation similar to the TE wave. Therefore, the propagation of the three-dimensional electromagnetic wave can be similarly simulated.
[0138]
Note that the processing contents of the functions that the electromagnetic wave analysis apparatus should have can be described in a program recorded on a computer-readable recording medium. Then, by executing this program on a computer, the above processing is realized by the computer. Examples of the computer-readable recording medium include a magnetic recording device and a semiconductor memory. When distributing to the market, store the program in a portable recording medium such as a CD-ROM (Compact Disk Read Only Memory) or floppy disk, or store it in a computer storage device connected via a network. In addition, it can be transferred to another computer through the network. When executed by a computer, the program is stored in a hard disk device or the like in the computer, loaded into the main memory and executed.
[0139]
Below, the features of the present invention are classified for each item, and the matters recognized as necessary for specifying the invention are listed.
(1st term) In the electromagnetic wave analysis device for analyzing the propagation of electromagnetic waves by solving the Maxwell equation for time and space by the difference method,
When the time discrete interval is Δt, based on the state of the electromagnetic wave at the first time, the state of the electromagnetic wave at the second time after the lapse of Δt from the first time, the alternating direction implicit finite difference time region Analytical means for analyzing the behavior of electromagnetic waves by
An electromagnetic wave analysis device characterized by comprising:
[0140]
(Second Item) The electromagnetic wave analyzing apparatus according to the first item, wherein the analyzing means analyzes propagation of a two-dimensional TE wave.
(3rd term) When the analysis means analyzes the electromagnetic wave of the differential grid with the vibration plane of the electric field as the xy plane and the arbitrary grid point (i, j) on the xy plane as a reference. When the spatial discrete interval in the x direction from the lattice point (i, j) is Δx, and the spatial discrete interval in the y direction from the lattice point (i, j) is Δy, x in (i + 1/2, j) A second term characterized in that an electric field in a direction, an electric field in a y direction at (i, j + 1/2), and a magnetic field in a z direction perpendicular to the xy plane at (i + 1/2, j + 1/2) are obtained. The electromagnetic wave analysis apparatus of description.
[0141]
(Item 4) The analysis means uses an explicit method for the electric field in the x direction, and uses an implicit method for the electric field in the y direction and the magnetic field in the z direction, so that the state of the electromagnetic wave at the first time Based on the first half process calculating means for calculating the state of the electromagnetic wave at the intermediate time after the lapse of Δt / 2 from the first time, and using the explicit method for the y-direction electric field, the x-direction electric field, and the z-direction electric field A second half process computing means for computing the state of the electromagnetic wave at the second time based on the state of the electromagnetic wave at the intermediate time calculated by the first half process computing means by using an implicit method for the magnetic field. 4. The electromagnetic wave analysis apparatus according to claim 3, wherein
[0142]
(Section 5) The first-half process calculation means sets the time at the n-th step in the simulation as the first time, the time at the (n + 1) -th step as the second time, and (i + 1/2, the electric field in the x direction in j)_{i + 1/2, j} ⁿ, The electric field in the y direction at (i, j + 1/2) at the first time is expressed as Ey_{i, j + 1/2} ⁿ, The magnetic field in the z direction at (i + 1/2, j + 1/2) at the first time is Hz_{i + 1/2, j + 1/2} ⁿ, The electric field in the x direction at (i + 1/2, j) at the intermediate time is expressed as Ex_{i + 1/2, j} ^{n + 1/2}, The electric field in the y direction at (i, j + 1/2) at the intermediate time is expressed as Ey_{i, j + 1/2} ^{n + 1/2}, The magnetic field in the z direction at (i + 1/2, j + 1/2) at the intermediate time is Hz_{i + 1/2, j + 1/2} ^{n + 1/2}When the dielectric constant of the region to be analyzed is ε and the magnetic permeability of the region to be analyzed is μ, an arbitrary lattice point (i, j) on the xy plane divided into a lattice shape is used as a reference. Electromagnetic wave
[0143]
[Formula 59]

[0144]
[Expression 60]

[0145]
[Equation 61]

[0146]
By defining the electromagnetic wave state with reference to all the lattice points at the intermediate time,
The second half process calculating means calculates the electric field in the x direction at (i + 1/2, j) at the second time as Ex._{i + 1/2, j} ^{n + 1}, The electric field in the y direction at (i, j + 1/2) at the second time is expressed as Ey_{i, j + 1/2} ^{n + 1}, The magnetic field in the z direction at (i + 1/2, j + 1/2) at the second time in Hz_{i + 1/2, j + 1/2} ^{n + 1}The electromagnetic wave with reference to an arbitrary lattice point (i, j) on the xy plane divided into a lattice shape.
[0147]
[62]

[0148]
[Equation 63]

[0149]
[Expression 64]

[0150]
The electromagnetic wave analysis apparatus according to claim 4, wherein the state of the electromagnetic wave is calculated based on all the lattice points at the second time.
(Sixth Item) The electromagnetic wave analyzing apparatus according to the first item, wherein the analyzing unit analyzes propagation of a two-dimensional TM wave.
[0151]
(Seventh Item) When the analysis means analyzes the electromagnetic wave of the difference grating with the vibration plane of the magnetic field as the xy plane and the arbitrary lattice point (i, j) on the xy plane as a reference. When the spatial discrete interval in the x direction from the lattice point (i, j) is Δx, and the spatial discrete interval in the y direction from the lattice point (i, j) is Δy, x in (i, j + 1/2) 7. An electromagnetic wave analysis according to claim 6, wherein a magnetic field in the direction, a magnetic field in the y direction at (i + 1/2, j), and an electric field in the z direction perpendicular to the xy plane at (i, j) are obtained. apparatus.
[0152]
(Item 8) The analysis means uses an explicit method for the magnetic field in the x direction and uses an implicit method for the magnetic field in the y direction and the electric field in the z direction, so that the state of the electromagnetic wave at the first time Based on the first half process calculating means for calculating the state of the electromagnetic wave at the intermediate time after the lapse of Δt / 2 from the first time, and using the explicit method for the y-direction magnetic field, the x-direction magnetic field, and the z-direction magnetic field A second half process computing means for computing the state of the electromagnetic wave at the second time based on the state of the electromagnetic wave at the intermediate time calculated by the first half process computing means by using an implicit method for the electric field; The electromagnetic wave analysis apparatus according to claim 7, wherein
[0153]
(Section 9) The first half process calculating means sets the time at the nth step in the simulation as the first time, the time at the (n + 1) th step as the second time, and (i, j + 1 //) at the first time. The magnetic field in the x direction in 2) is expressed as Hx_{i, j + 1/2} ⁿ, The magnetic field in the y direction at (i + 1/2, j) at the first time is Hy_{i + 1/2, j} ⁿ, The electric field in the z direction at (i, j) at the first time is expressed as Ez_{i, j} ⁿ, The magnetic field in the x direction at (i, j + 1/2) at the intermediate time is represented by Hx_{i, j + 1/2} ^{n + 1/2}, The magnetic field in the y direction at (i + 1/2, j) at the intermediate time is Hy_{i + 1/2, j} ^{n + 1/2}, The electric field in the z direction at (i, j) at the intermediate time is expressed as Ez_{i, j} ^{n + 1/2}When the dielectric constant of the region to be analyzed is ε and the magnetic permeability of the region to be analyzed is μ, an arbitrary lattice point (i, j) on the xy plane divided into a lattice shape is used as a reference. Electromagnetic wave
[0154]
[Equation 65]

[0155]
[66]

[0156]
[Expression 67]

[0157]
By calculating the state of the electromagnetic wave with reference to all the lattice points at the intermediate time,
The second half process calculating means calculates the magnetic field in the x direction at (i, j + 1/2) at the second time as Hx._{i, j + 1/2} ^{n + 1}, The magnetic field in the y direction at (i + 1/2, j) at the second time is Hy_{i + 1/2, j} ^{n + 1}, The electric field in the z direction at (i, j) at the second time is expressed as Ez_{i, j} ^{n + 1}The electromagnetic wave with reference to an arbitrary lattice point (i, j) on the xy plane divided into a lattice shape.
[0158]
[Equation 68]

[0159]
[Equation 69]

[0160]
[Equation 70]

[0161]
The electromagnetic wave analysis apparatus according to claim 8, wherein the electromagnetic wave state is calculated based on all lattice points at the second time by defining as follows.
(Item 10) The electromagnetic wave analysis device according to item 1, wherein the analysis means analyzes propagation of a three-dimensional wave.
[0162]
(Section 11) The analysis means assumes a three-dimensional space of xyz, and analyzes an electromagnetic wave of a difference grid based on an arbitrary grid point (i, j, k) in the three-dimensional space. Includes a spatial discrete interval from the lattice point (i, j, k) in the x direction by Δx, a spatial discrete interval from the lattice point (i, j, k) in the y direction by Δy, and a lattice point (i, j, k). k) to z-direction spatial discrete interval Δz, x-direction electric field at (i + 1/2, j, k), y-direction electric field at (i, j + 1/2, k), (i, j, k + 1/2) in the z direction, (i, j + 1/2, k + 1/2) in the x direction, (i + 1/2, j, k + 1/2) in the y direction, and (i + 1 / The electromagnetic wave analysis apparatus according to claim 10, wherein a magnetic field in the z direction at 2, j + 1/2, k) is obtained.
[0163]
(Section 12) The analysis means includes:
The time at the nth step in the simulation is the first time, the time at the (n + 1) th step is the second time, and the electric field in the x direction at (i + 1/2, j, k) at the first time is Ex._{i + 1/2, j, k} ⁿ, The electric field in the y direction at (i, j + 1/2, k) at the first time is expressed as Ey_{i, j + 1/2, k} ⁿ, The electric field in the z direction at (i, j, k + 1/2) at the first time is expressed as Ez_{i, j, k + 1/2} ⁿ, The magnetic field in the x direction at (i, j + 1/2, k + 1/2) at the first time is represented by Hx_{i, j + 1/2, k + 1/2} ⁿ, The magnetic field in the y direction at (i + 1/2, j, k + 1/2) at the first time is Hy_{i + 1/2, j, k + 1/2} ⁿ, The magnetic field in the z direction at (i + 1/2, j + 1/2, k) at the first time is Hz_{i + 1/2, j + 1/2, k} ⁿ, The electric field in the x direction at (i + 1/2, j, k) at the intermediate time is expressed as Ex_{i + 1/2, j, k} ^{n + 1/2}, The electric field in the y direction at (i, j + 1/2, k) at the intermediate time is expressed as Ey_{i, j + 1/2, k} ^{n + 1/2}, The electric field in the z direction at (i, j, k + 1/2) at the intermediate time is expressed as Ez_{i, j, k + 1/2} ^{n + 1/2}, The magnetic field in the x direction at (i, j + 1/2, k + 1/2) at the intermediate time is represented by Hx_{i, j + 1/2, k + 1/2} ^{n + 1/2}, The magnetic field in the y direction at (i + 1/2, j, k + 1/2) at the intermediate time is Hy_{i + 1/2, j, k + 1/2} ^{n + 1/2}, The magnetic field in the z direction at (i + 1/2, j + 1/2, k) at the intermediate time is Hz_{i + 1/2, j + 1/2, k} ^{n + 1/2}When the dielectric constant of the region to be analyzed is ε and the magnetic permeability of the region to be analyzed is μ, an arbitrary lattice point (i, j, k) in the three-dimensional space divided into a lattice is used as a reference. Electromagnetic waves
[0164]
[Equation 71]

[0165]
[Equation 72]

[0166]
[Equation 73]

[0167]
[Equation 74]

[0168]
[Expression 75]

[0169]
[76]

[0170]
By defining the first half process calculation means for calculating the state of the electromagnetic wave with reference to all the lattice points at the intermediate time,
The electric field in the x direction at (i + 1/2, j, k) at the intermediate time is expressed as Ex_{i + 1/2, j, k} ^{n + 1}, The electric field in the y direction at (i, j + 1/2, k) at the intermediate time is expressed as Ey_{i, j + 1/2, k} ^{n + 1}, The electric field in the z direction at (i, j, k + 1/2) at the intermediate time is expressed as Ez_{i, j, k + 1/2} ^{n + 1}, The magnetic field in the x direction at (i, j + 1/2, k + 1/2) at the intermediate time is represented by Hx_{i, j + 1/2, k + 1/2} ^{n + 1}, The magnetic field in the y direction at (i + 1/2, j, k + 1/2) at the intermediate time is Hy_{i + 1/2, j, k + 1/2} ^{n + 1}, The magnetic field in the z direction at (i + 1/2, j + 1/2, k) at the intermediate time is Hz_{i + 1/2, j + 1/2, k} ^{n + 1}The electromagnetic wave based on an arbitrary lattice point (i, j, k) in the three-dimensional space divided into a lattice shape.
[0171]
[77]

[0172]
[Formula 78]

[0173]
[79]

[0174]
[80]

[0175]
[Formula 81]

[0176]
[Formula 82]

[0177]
By defining the latter half process calculating means for calculating the state of the electromagnetic wave with reference to all lattice points at the second time,
The electromagnetic wave analysis apparatus according to claim 11, comprising:
[0178]
(Section 13) In a computer-readable recording medium recording an electromagnetic wave analysis program for analyzing propagation of electromagnetic waves by solving Maxwell's equations for time and space by a difference method,
When the time discrete interval is Δt, based on the state of the electromagnetic wave at the first time, the state of the electromagnetic wave at the second time after the lapse of Δt from the first time, the alternating direction implicit finite difference time region Analytical means for analyzing the behavior of electromagnetic waves by
A computer-readable recording medium recording an electromagnetic wave analysis program characterized by causing the computer to function as
[0179]
【The invention's effect】
As described above, in the electromagnetic wave analysis apparatus of the present invention, the behavior of the electromagnetic wave is analyzed by obtaining the state of the electromagnetic wave by the alternating direction implicit finite difference time domain method, so that it does not depend on the spatial discrete interval. It is possible to determine time discrete intervals. As a result, it is possible to take a relatively large time discrete interval when performing electromagnetic wave analysis of an analysis object having a fine structure, and analysis processing can be performed at high speed.
[0180]
Further, in the computer-readable recording medium on which the electromagnetic wave analysis program of the present invention is recorded, the electromagnetic wave state is changed to the electromagnetic wave by the alternating direction implicit finite difference time domain method by causing the computer to execute the recorded electromagnetic wave analysis program. Can be analyzed at high speed by a computer.
[Brief description of the drawings]
FIG. 1 is a principle configuration diagram of the present invention.
FIG. 2 is an enlarged view of a lattice.
FIG. 3 is a flowchart showing a processing procedure for two-dimensional TE wave analysis by an ADI-FDTD method.
FIG. 4 is a diagram illustrating a hardware configuration of the electromagnetic wave analysis device.
FIG. 5 is a diagram showing a first simulation model.
FIG. 6 is a diagram showing a calculation result in the first simulation model.
FIG. 7 is a diagram illustrating a second simulation model.
FIG. 8 is a diagram illustrating a calculation result of a second simulation.
[Explanation of symbols]
1 Analysis command input means
2 Analysis means
2a First half process calculation means
2b Second half process calculation means
3 Analysis result output means

Claims (6)

In an electromagnetic wave analysis device that analyzes the propagation of electromagnetic waves by solving the Maxwell equation for time and space by the difference method,
When the time discrete interval is Δt, based on the state of the electromagnetic wave at the first time, the state of the electromagnetic wave at the second time after the lapse of Δt from the first time, the alternating direction implicit finite difference time domain Analytical means for analyzing the behavior of electromagnetic waves by
An electromagnetic wave analysis device characterized by comprising: The electromagnetic wave analysis apparatus according to claim 1, wherein the analysis unit analyzes propagation of a two-dimensional TE wave. The analysis means uses the lattice point (i) when analyzing the electromagnetic wave of the differential grating with the vibration plane of the electric field as the xy plane and the arbitrary lattice point (i, j) on the xy plane as a reference. , J) to x direction spatial discrete interval Δx, and assuming that the spatial discrete interval from grid point (i, j) to y direction is Δy, the electric field in x direction at (i + 1/2, j), ( 3. The electromagnetic wave analysis apparatus according to claim 2, wherein an electric field in the y direction at i, j + 1/2) and a magnetic field in the z direction perpendicular to the xy plane at (i + 1/2, j + 1/2) are obtained. . The analysis means uses an explicit method for the electric field in the x direction and uses an implicit method for the electric field in the y direction and the magnetic field in the z direction, thereby determining the first time based on the state of the electromagnetic wave at the first time. The first half process calculating means for calculating the state of the electromagnetic wave at the intermediate time after the lapse of Δt / 2 from the time of 1 and the explicit method for the electric field in the y direction, and the A second-half process calculating means for calculating the state of the electromagnetic wave at the second time based on the state of the electromagnetic wave at the intermediate time calculated by the first-half process calculating means by using a solution. The electromagnetic wave analysis apparatus according to claim 3. The first-half process calculation means sets the time at the nth step in the simulation as the first time, the time at the (n + 1) th step as the second time, and x at (i + 1/2, j) at the first time. The electric field in the direction is Ex _{i + 1/2, j} ⁿ , and the electric field in the y direction at (i, j + 1/2) at the first time is Ey _{i, j + 1/2} ⁿ , at the first time. The magnetic field in the z direction at (i + 1/2, j + 1/2) is Hz _{i + 1/2, j + 1/2} ⁿ , and the electric field in the x direction at (i + 1/2, j) at the intermediate time is Ex _{i. + 1/2, j} ^{n + 1/2} , the electric field in the y direction at (i, j + 1/2) at the intermediate time is represented by Ey _{i, j + 1/2} ^{n + 1/2} , (i + 1 at the intermediate time) / 2, j + 1/2), the magnetic field in the z direction is Hz _{i + 1/2, j + 1/2} ^{n + 1/2} , the dielectric constant of the region to be analyzed is ε, and the magnetic permeability of the region to be analyzed is In the case of μ, the above divided into a lattice shape An electromagnetic wave based on an arbitrary lattice point (i, j) on the xy plane

By calculating the state of the electromagnetic wave with reference to all the lattice points at the intermediate time,
The latter half process calculating means sets the electric field in the x direction at (i + 1/2, j) at the second time to Ex _{i + 1/2, j} ^{n + 1} , and (i, j + 1 // at the second time). The electric field in the y direction at 2) is Ey _{i, j + 1/2} ^{n + 1} , and the magnetic field in the z direction at (i + 1/2, j + 1/2) at the second time is Hz _{i + 1/2, When j + 1/2} ^{n + 1} , an electromagnetic wave based on an arbitrary lattice point (i, j) on the xy plane divided in a lattice shape is used.

The electromagnetic wave analysis apparatus according to claim 4, wherein the state of the electromagnetic wave is calculated based on all the lattice points at the second time. In a computer-readable recording medium recording an electromagnetic wave analysis program for analyzing propagation of electromagnetic waves by solving the Maxwell equation by time and space by a difference method,
When the time discrete interval is Δt, based on the state of the electromagnetic wave at the first time, the state of the electromagnetic wave at the second time after the lapse of Δt from the first time, the alternating direction implicit finite difference time domain Analytical means for analyzing the behavior of electromagnetic waves by
A computer-readable recording medium recording an electromagnetic wave analysis program characterized by causing the computer to function as

Images