(57)【要約】【課題】吸収境界条件としてＢｅｒｅｎｇｅｒのＰＭＬを適用した場合の解析精度を維持しつつ、計算に必要なメモリ量や計算時間を低減すること。【解決手段】ＦＤＴＤ法により解析する解析領域解析手段（Ｃ１Ｄ）と、吸収境界領域（５）で分割された領域毎に電界（Ｅ）および磁界（Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）を記憶する吸収境界電磁界記憶手段（Ｃ１E3）と、仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電界（Ｅ）または磁界（Ｈ）のいずれか一方を演算する仮想時間電磁界演算手段（Ｃ１E1a＋Ｃ１E1b）と、解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の電界（Ｅ）または磁界（Ｈ）のいずれか一方を演算する解析時間後電磁界演算手段（Ｃ１E2a＋Ｃ１E2b）と、電磁界（Ｅ，Ｈ）を表示する解析電磁界表示手段（Ｃ２）とを備えたことを特徴とする電磁界解析装置（ＰＣ）。【選択図】図２
2008234335.tif 000002
【特許請求の範囲】
【請求項１】
電界および磁界からなる電磁界をＦＤＴＤ法で解析する電磁界解析装置において、
電界および磁界を解析する対象の領域である解析領域を、電界用分割領域および前記電界用分割領域とは位置がずれた磁界用分割領域に分割し、分割された領域毎にＦＤＴＤ法により電界および磁界を解析する解析領域解析手段と、
前記解析領域の境界に設けられて電界および磁界を吸収して減衰させるための吸収境界領域を、電界用分割領域および前記電界用分割領域とは位置がずれた磁界用分割領域に分割し、前記吸収境界領域で分割された領域毎に電界および磁界の強さを多次元で記憶する吸収境界電磁界記憶手段と、
前記吸収境界領域の前記電界または磁界のいずれか一方に対応する分割領域毎に、前記多次元の特定の第１軸方向の前記分割領域の長さに対する電界または磁界のいずれか他方の空間差分と、解析時間間隔に対する電界または磁界のいずれか一方の時間差分と、に基づいて、仮想時間における電界または磁界のいずれか一方を演算する仮想時間電磁界演算手段と、
前記吸収境界領域の前記電界または磁界のいずれか一方に対応する分割領域毎に、前記第１軸方向とは独立な軸方向の前記分割領域の長さに対する磁界または電界のいずれか他方の空間差分と、解析時間間隔に対する前記仮想時間における電界または磁界のいずれか一方の時間差分とに基づいて、前記解析時間後の電界または磁界のいずれか一方を演算する解析時間後電磁界演算手段と、
前記解析された解析時間後の電磁界を表示する解析電磁界表示手段と、
を備えたことを特徴とする電磁界解析装置。
【請求項２】
前記解析領域および前記吸収境界領域の大きさ、および、解析時間間隔、誘電率、透磁率を含む解析条件パラメタを入力する解析条件入力画像を表示する解析条件入力画像表示手段と、
を備えたことを特徴とする請求項１に記載の電磁界解析装置。
【請求項３】
電界および磁界を解析する対象の領域である解析領域を、電界用分割領域および前記電界用分割領域とは位置がずれた磁界用分割領域に分割し、分割された領域毎にＦＤＴＤ法により電界および磁界を解析すると共に、
前記解析領域の境界に設けられて電界および磁界を吸収して減衰させるための吸収境界領域において、前記吸収境界領域の前記電界または磁界のいずれか一方に対応する分割領域毎に、多次元の特定の軸方向の前記分割領域の長さに対する電界または磁界のいずれか他方の空間差分と、解析時間間隔に対する電界または磁界のいずれか一方の時間差分と、に基づいて、仮想時間における電界または磁界のいずれか一方を演算し、
前記解析時間間隔に対する前記仮想時間における電界または磁界のいずれか一方の時間差分と、前記第一軸方向とは独立な軸方向の前記分割領域の長さに対する磁界または電界のいずれか他方の空間差分とに基づいて、前記解析時間後の電界または磁界のいずれか一方を演算することを特徴とする電磁界解析方法。
【請求項４】
電界および磁界からなる電磁界をＦＤＴＤ法で解析する電磁界解析装置を、
電界および磁界を解析する対象の領域である解析領域を、電界用分割領域および前記電界用分割領域とは位置がずれた磁界用分割領域に分割し、分割された領域毎にＦＤＴＤ法により電界および磁界を解析する解析領域解析手段、
前記解析領域の境界に設けられて電界および磁界を吸収して減衰させるための吸収境界領域を、電界用分割領域および前記電界用分割領域とは位置がずれた磁界用分割領域に分割し、前記吸収境界領域で分割された領域毎に電界および磁界の強さを多次元で記憶する記憶手段、
前記吸収境界領域の前記電界または磁界のいずれか一方に対応する分割領域毎に、前記多次元の特定の第１軸方向の前記分割領域の長さに対する電界または磁界のいずれか他方の空間差分と、解析時間間隔に対する電界または磁界のいずれか一方の時間差分と、に基づいて、仮想時間における電界または磁界のいずれか一方を演算する仮想時間電磁界演算手段、
前記吸収境界領域の前記電界または磁界のいずれか一方に対応する分割領域毎に、前記第１軸方向とは独立な軸方向の前記分割領域の長さに対する磁界または電界のいずれか他方の空間差分と、解析時間間隔に対する前記仮想時間における電界または磁界のいずれか一方の時間差分とに基づいて、前記解析時間後の電界または磁界のいずれか一方を演算する解析時間後電磁界演算手段、
前記解析された解析時間後の電磁界を表示する解析電磁界表示手段、
として機能させるための電磁界解析プログラム。
【請求項５】
第１要素の時間差分が第２要素の空間差分で表現され、前記第２要素の時間差分が前記第１要素の空間差分で表現される場をＦＤＴＤ法で解析する場解析装置において、
前記第１要素および前記第２要素を解析する対象の領域である解析領域を、第１要素用分割領域および前記第１要素用分割領域とは位置がずれた第２要素用分割領域に分割し、分割された領域毎にＦＤＴＤ法により前記第１要素および前記第２要素を解析する解析領域解析手段と、
前記解析領域の境界に設けられて前記第１要素および前記第２要素を吸収して減衰させるための吸収境界領域を、第１要素用分割領域および前記第１要素用分割領域とは位置がずれた第２要素用分割領域に分割し、前記吸収境界領域で分割された領域毎に前記第１要素および前記第２要素の強さを多次元で記憶する吸収境界要素記憶手段と、
前記吸収境界領域の前記第１要素または前記第２要素のいずれか一方に対応する分割領域毎に、前記多次元の特定の第１軸方向の前記分割領域の長さに対する前記第１要素または前記第２要素のいずれか他方の空間差分と、解析時間間隔に対する前記第１要素または前記第２要素のいずれか一方の時間差分と、に基づいて、仮想時間における前記第１要素または前記第２要素のいずれか一方を演算する仮想時間要素演算手段と、
前記吸収境界領域の前記第１要素または前記第２要素のいずれか一方に対応する分割領域毎に、前記第１軸方向とは独立な軸方向の前記分割領域の長さに対する前記第２要素または前記第１要素のいずれか他方の空間差分と、解析時間間隔に対する前記仮想時間における前記第１要素または前記第２要素のいずれか一方の時間差分とに基づいて、前記解析時間後の前記第１要素または前記第２要素のいずれか一方を演算する解析時間後要素演算手段と、
前記解析された解析時間後の前記場を表示する解析場表示手段と、
を備えたことを特徴とする場解析装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、予め与えられた解析領域から電磁界等の場を解析できる電磁界解析装置、電磁界解析方法、電磁界解析プログラムおよび場解析装置に関する。
【背景技術】
【０００２】
従来から、電磁波、音波、地震波等の波の伝搬の解析が行われており、特に近年、携帯電話の発達に応じて、アンテナ等から発生する電磁波が生体に与える影響の解析が盛んに行われている。このような解析は、膨大な計算量と記憶容量が必要だったため、スーパーコンピュータを使用して行われていたが、近年の計算機技術の発展に伴い、高性能なパーソナルコンピュータ等の計算機が安価で入手できるようになった。このため、前記解析実験の解析時間や金銭的な問題を低減するために、前記解析は、前記計算機上の数値シミュレーションとして行われるようになった。この結果、前記解析の数値シミュレーションは、大学等の研究機関により活発に行われることとなり、その解析手法についても、多くの研究が行われることとなった。
【０００３】
（ＦＤＴＤ法について）
図１２はＦＤＴＤ法の説明図であり、図１２ＡはＦＤＴＤ法の解析領域と解析領域セルの関係の説明図であり、図１２Ｂは図１２Ａの解析領域セルと電磁界の配置の説明図である。
アンテナ解析や電磁界解析等の数値電磁解析（computational electromagnetics）の分野における解析手法として、１９６６年にＫ．Ｓ．Ｙｅｅによって提案されたＦＤＴＤ法（Finite Difference Time Domain method）がある。
前記ＦＤＴＤ法では、図１２Ａ、図１２Ｂに示すように、電磁波等の入力波の入力波源（散乱体、アンテナ等）０１を囲むように解析領域０２をとり、前記解析領域０２全体を微小直方体である解析領域セル０３に分割する。
【０００４】
また、前記ＦＤＴＤ法では、電界Ｅの回転が磁界Ｈを作り、磁界Ｈの回転が電界Ｅを作るというマクスウェルの方程式（マクスウェル回転方程式系）が基礎方程式として利用されている。前記マクスウェル回転方程式系は、以下の式（１），（２）で示される。
ｒｏｔＥ＝∇×Ｅ＝−（∂／∂ｔ）Ｂ＝−μ（∂／∂ｔ）Ｈ                  …（１）
ｒｏｔＨ＝∇×Ｈ＝（∂／∂ｔ）Ｄ＋Ｊ＝ε（∂／∂ｔ）Ｅ＋σＥ            …（２）
ここで、ｒｏｔＨは磁界Ｈの回転（rotation）、ｒｏｔＥは電界Ｅの回転、∇はベクトル量としての微分演算子、∇×Ｈはベクトル量としての磁界Ｈと前記微分演算子∇とのベクトル積、∇×Ｅはベクトル量としての電界Ｅと前記微分演算子∇とのベクトル積、Ｄは電束密度、Ｊは電流密度、Ｂは磁束密度、μは透磁率、εは誘電率、σは導電率（電気伝導率）である。
【０００５】
まず、前記ＦＤＴＤ法のアルゴリズムであるＹｅｅアルゴリズムでは、図１２Ｂに示すように、前記解析領域セル０３を、Δｘ×Δｙ×Δｚのサイズの微小直方体とし、前記電界ＥについてのＸ軸方向、Ｙ軸方向、Ｚ軸方向のベクトル成分であるＥｘ，Ｅｙ，Ｅｚは、前記解析領域セル０３の各辺の中央部から前記各辺に沿って割り当て、前記磁界ＨについてのＸ軸方向、Ｙ軸方向、Ｚ軸方向のベクトル成分であるＨｘ，Ｈｙ，Ｈｚは、前記解析領域セル０３の各面の中心から垂直方向に割り当てる。次に、Ｙｅｅアルゴリズムでは、前記式（１），（２）の偏微分を中心差分近似で計算する。
【０００６】
ここで、前記Ｙｅｅアルゴリズムにおける中心差分近似の具体例として、前記式（１），（２）を、前記単位時間Δｔ間隔の中心差分近似の式に置換した式（３），（４）を以下に示す。
（Ｈｎ＋１／２−Ｈｎ−１／２）／Δｔ＝−（１／μ）∇×Ｅｎ             …（３）
ε（Ｅｎ−Ｅｎ−１）／Δｔ＋σ（Ｅｎ＋Ｅｎ−１）／２＝∇×Ｈｎ−１／２  …（４）
なお、図１２Ｂに示すように、前記電界Ｅおよび前記磁界Ｈが割り当てられた場合、前記電界Ｅおよび前記磁界Ｈは、時間的に交互に配置されることから、前記式（３），（４）では、（ｎ−１）・Δｔ，ｎ・Δｔの各タイミングにおける電界ＥをそれぞれＥｎ−１，Ｅｎとし、（ｎ−１／２）・Δｔ，（ｎ＋１／２）・Δｔの各タイミングにおける磁界ＨをそれぞれＨｎ−１／２，Ｈｎ＋１／２とした。また、前記式（３），（４）をＨｎ＋１／２，Ｅｎについて解いた式を式（３′），（４′）に示す。
Ｈｎ＋１／２＝Ｈｎ−１／２−（Δｔ／μ）∇×Ｅｎ                     …（３′）
Ｅｎ＝（（１−σ・Δｔ／２ε）／（１＋σ・Δｔ／２ε））Ｅｎ−１
＋（（Δｔ／ε）／（１＋σ・Δｔ／２ε））∇×Ｈｎ−１／２      …（４′）
【０００７】
また、前記式（３），（４）のベクトル積∇×Ｈｎ−１／２，∇×Ｅｎも同様に、単位長Δｘ，Δｙ，Δｚの式に置換できる。
例えば、磁界ＨのＸ軸方向のベクトル成分Ｈｘは、前記式（３′）から、以下の式（３−１′）が得られる。
Ｈｘｎ＋１／２＝Ｈｘｎ−１／２
−（Δｔ／μ）（（∂／∂ｙ）Ｅｚｎ−（∂／∂ｚ）Ｅｙｎ）…（３−１′）
ここで、図１２に示すように、前記解析領域セル０３の座標が（ｉ，ｊ＋１／２，ｋ＋１／２）の磁界Ｈのベクトル成分をＨｘ（ｉ，ｊ＋１／２，ｋ＋１／２）とし、前記座標が（ｉ，ｊ，ｋ＋１／２），（ｉ，ｊ＋１，ｋ＋１／２）の電界Ｅのベクトル成分をＥｚ（ｉ，ｊ，ｋ＋１／２），Ｅｚ（ｉ，ｊ＋１，ｋ＋１／２）とし、Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２）についてのＥｚｎの偏微分を（∂／∂ｙ）Ｅｚｎ｜（ｉ，ｊ＋１／２，ｋ＋１／２）とした場合、以下の式（５）が成立する。
【０００８】
（∂／∂ｙ）Ｅｚｎ｜（ｉ，ｊ＋１／２，ｋ＋１／２）
＝（Ｅｚｎ（ｉ，ｊ＋１，ｋ＋１／２）−Ｅｚｎ（ｉ，ｊ，ｋ＋１／２））／Δｙ
…（５）
また、前記座標が（ｉ，ｊ＋１／２，ｋ），（ｉ，ｊ＋１／２，ｋ＋１）の電界Ｅのベクトル成分をＥｙ（ｉ，ｊ＋１／２，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ＋１）とし、Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２）についてのＥｙｎの偏微分を（∂／∂ｙ）Ｅｙｎ｜（ｉ，ｊ＋１／２，ｋ＋１／２）とした場合、以下の式（６）が成立する。
（∂／∂ｚ）Ｅｙｎ｜（ｉ，ｊ＋１／２，ｋ＋１／２）
＝（Ｅｙｎ（ｉ，ｊ＋１／２，ｋ＋１）−Ｅｙｎ（ｉ，ｊ＋１／２，ｋ））／Δｚ
…（６）
【０００９】
また、その他の偏微分についても同様に変換できる。
例えば、磁界ＨのＺ軸方向のベクトル成分Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ）については、前記式（３′）から、以下の式（３−２′）が得られる。
Ｈｚｎ＋１／（ｉ＋１／２，ｊ＋１／２，ｋ）
＝Ｈｚｎ−１／２（ｉ＋１／２，ｊ＋１／２，ｋ）
−（Δｔ／μ）｛（∂／∂ｘ）Ｅｙｎ−（∂／∂ｙ）Ｅｘｎ｝        …（３−２′）
ここで、Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ）についてのＥｘ，Ｅｙに関する偏微分（∂／∂ｘ）Ｅｙｎ｜（ｉ＋１／２，ｊ＋１／２，ｋ），（∂／∂ｙ）Ｅｘｎ｜（ｉ＋１／２，ｊ＋１／２，ｋ）については、以下の式（５′），（６′）が成立する。
（∂／∂ｘ）Ｅｙｎ｜（ｉ＋１／２，ｊ＋１／２，ｋ）
＝（Ｅｙｎ（ｉ＋１，ｊ＋１／２，ｋ）
−Ｅｙｎ（ｉ，ｊ＋１／２，ｋ））／Δｘ                            …（５′）
（∂／∂ｙ）Ｅｘｎ｜（ｉ＋１／２，ｊ＋１／２，ｋ）
＝（Ｅｘｎ（ｉ＋１／２，ｊ＋１，ｋ）
−Ｅｘｎ（ｉ＋１／２，ｊ，ｋ））／Δｙ                            …（６′）
【００１０】
また、電界ＥのＸ軸方向のベクトル成分Ｅｘ（ｉ＋１／２，ｊ，ｋ）については、前記式（４′）から、以下の式（４−１′）が得られる。
Ｅｘｎ（ｉ＋１／２，ｊ，ｋ）
＝｛（１−σ・Δｔ／２ε）／（１＋σ・Δｔ／２ε）｝
×Ｅｘｎ−１（ｉ＋１／２，ｊ，ｋ）
＋｛（Δｔ／ε）／（１＋σ・Δｔ／２ε）｝
×（（∂／∂ｙ）Ｈｚｎ−１／２−（∂／∂ｚ）Ｈｙｎ−１／２）   …（４−１′）
ここで、Ｅｘ（ｉ＋１／２，ｊ，ｋ）についてのＨｙ，Ｈｚに関する偏微分（∂／∂ｘ）Ｈｙｎ−１／２｜（ｉ＋１／２，ｊ，ｋ），（∂／∂ｙ）Ｈｚｎ−１／２｜（ｉ＋１／２，ｊ，ｋ）については、以下の式（５″），（６″）が成立する。
（∂／∂ｘ）Ｈｙｎ−１／２｜（ｉ＋１／２，ｊ，ｋ）
＝（Ｈｙｎ−１／２（ｉ＋１／２，ｊ，ｋ＋１／２）
−Ｈｙｎ−１／２（ｉ＋１／２，ｊ，ｋ−１／２））／Δｘ           …（５″）
（∂／∂ｙ）Ｈｚｎ−１／２｜（ｉ＋１／２，ｊ，ｋ）
＝（Ｈｚｎ−１／２（ｉ，ｊ＋１／２，ｋ＋１／２）
−Ｈｚｎ−１／２（ｉ，ｊ＋１／２，ｋ−１／２））／Δｙ           …（６″）
【００１１】
したがって、前記解析領域０２における各座標位置、すなわち、計算対象となる解析領域セル０３の位置毎のベクトル積∇×Ｈｎ−１／２，∇×Ｅｎは、計算対象となる解析領域セル０３のベクトル量Ｅｎ−１，Ｈｎ−１／２と、計算対象となる解析領域セル０３に隣接するセル０３のベクトル量Ｅｎ−１，Ｈｎ−１／２との中心差分近似の式に基づいて計算できる。
なお、前記各単位長Δｘ，Δｙ，Δｚおよび前記単位時間Δｔについては、Ｃｏｕｒａｎｔの安定条件により、光速をｃとした場合に、以下の式（７）が成立している必要があることが知られている。
ｃΔｔ≦１／（（１／Δｘ）２＋（１／Δｙ）２＋（１／Δｚ）２）１／２    …（７）
【００１２】
したがって、前記解析領域０２における各座標位置のベクトル量Ｅ，Ｈの時間的な変化、すなわち、前記入力波源０１から発生した前記入力波の伝搬の様子は、前記ＦＤＴＤ法のＹｅｅアルゴリズムに基づいて計算できる。すなわち、前記解析領域セル０３毎に、（ｎ−１）・Δｔのタイミングで計算されて既知のベクトル量Ｅｎ−１＝（Ｅｘｎ−１，Ｅｙｎ−１，Ｅｚｎ−１），Ｈｎ−１／２＝（Ｈｘｎ−１／２，Ｈｙｎ−１／２，Ｈｚｎ−１／２）から、次のｎ・Δｔのタイミングのベクトル量Ｅｎ＝（Ｅｘｎ，Ｅｙｎ，Ｅｚｎ），Ｈｎ＋１／２＝（Ｈｘｎ＋１／２，Ｈｙｎ＋１／２，Ｈｚｎ＋１／２）を計算することができる。
【００１３】
なお、前記ＦＤＴＤ法は、例えば、Ｚ軸方向に変化のない電磁波（平面波）の場合、２次元問題として扱うことができるため、前記解析領域０２を２次元領域、すなわち、Ｘ軸とＹ軸とによる平面領域として扱うことができる。すなわち、前記平面波の解析については２次元ＦＤＴＤ法として簡略化することが可能である。
前記ＦＤＴＤ法は、アルゴリズムが簡単であること、優れた解析精度を持つこと、複雑な物質の解析や材料定数の異なる物質の解析にも適している等の理由から、前記数値電磁解析の有効な解析手法として広く用いられている。なお、前記数値電磁解析としては、前記ＦＤＴＤ法以外にも、モーメント法や有限要素法が知られている（非特許文献１参照）。
【００１４】
（吸収境界について）
図１３は吸収境界条件の説明図である。
ここで、前記ＦＤＴＤ法は、図１２に示すように、前記解析領域０２で閉じた領域（閉領域）の解析手法である。よって、記解析領域０２が完全導体（電気抵抗がゼロとなる物質）で囲まれているような前記閉領域の問題を扱う場合には、前記完全導体表面である前記解析領域０２の外面０４の電界の接線成分、すなわち、外面０４に沿って割り当てられた電界Ｅのベクトル成分（Ｅｘ，Ｅｙ，Ｅｚ）を０とすればよいが、散乱解析あるいはアンテナ解析等のように、開放された領域（開放領域）の問題を扱う場合には、前記外面０４の電界の接線成分（Ｅｘ，Ｅｙ，Ｅｚ）を０にできない。したがって、前記開放領域の問題を扱う場合には、前記外面０４に入射した電磁波は、前記外面０４の電界の接線成分（Ｅｘ，Ｅｙ，Ｅｚ）の影響を受け、前記外面０４で電磁波の反射（反射波）が発生してしまう。この結果、発生した前記反射波が前記解析領域０２に戻るため、解析における誤差の原因となるという問題がある。
【００１５】
しかし、この問題は、図１３に示すように、前記外面０４の外方に電磁波を反射させずに吸収してしまう仮想的な領域である吸収境界０５を設けて、前記解析領域０２を閉じてしまうことで解決する。具体的には、前記吸収境界０５は、前記吸収境界０５全体を前記解析領域セル０３と同様の微小直方体である吸収境界セル０６に分割された複数の層によって構成する。そして、前記吸収境界０５において、前記電磁波は、前記複数の層を通過するうちに、減衰するようにして、前記吸収境界０５は、前記電磁波を吸収できるようにすることで解決する。
【００１６】
（ＢｅｒｅｎｇｅｒのＰＭＬ（Ｂ−ＰＭＬ）について）
ここで、前記吸収境界０５では、前記Ｙｅｅアルゴリズムが適用できないため、特別な条件、いわゆる、吸収境界条件（ABC（Absorbing Boundary Condition））が必要となる。前記吸収境界０５の吸収境界条件は、これまで、数多くの提案がされている。
特に、Ｂｅｒｅｎｇｅｒによって提案されたＢｅｒｅｎｇｅｒのＰＭＬ（Berenger’s PML（Perfectly Matched Layer）、以下Ｂ−ＰＭＬと記載する）は、解析における誤差が最も少なくなる、すなわち、解析精度が最も高い吸収境界条件のうちの一つであり、広く用いられている。
【００１７】
前記吸収境界条件として前記Ｂ−ＰＭＬを適用した場合、前記吸収境界０５は、非物理的な仮想媒質により構成されることになる。具体的には、前記Ｂ−ＰＭＬでは、前記吸収境界０５において、Ｘ軸に垂直な面（ＹＺ面）についてのＸ軸，Ｙ軸，Ｚ軸方向のベクトル成分の導電率をσｘｘ，σｘｙ，σｘｚ、導磁率（導電率に対応する仮想的な値）をσｘｘ＊，σｘｙ＊，σｘｚ＊とし、Ｙ軸に垂直な面（ＺＸ面）についてのＸ軸方向、Ｙ軸方向、Ｚ軸方向のベクトル成分の導電率をσｙｘ，σｙｙ，σｙｚ、導磁率をσｙｘ＊，σｙｙ＊，σｙｚ＊とし、Ｚ軸に垂直な面（ＸＹ面）に電磁波が入射した場合のＸ軸方向、Ｙ軸方向、Ｚ軸方向のベクトル成分の導電率をσｚｘ，σｚｙ，σｚｚ、導磁率をσｚｘ＊，σｚｙ＊，σｚｚ＊とした場合に、以下の式（８−１）〜（８−３）および式（９−１），（９−２）が成立する。
【００１８】
σｘｘ／ε０＝σｘｘ＊／μ０                                        …（８−１）
σｙｙ／ε０＝σｙｙ＊／μ０                                        …（８−２）
σｚｚ／ε０＝σｚｚ＊／μ０                                        …（８−３）
σｘｙ＝σｘｚ＝σｙｚ＝σｙｘ＝σｚｙ＝σｚｘ＝０                  …（９−１）
σｘｙ＊＝σｘｚ＊＝σｙｚ＊＝σｙｘ＊＝σｚｙ＊＝σｚｘ＊＝０     …（９−２）
なお、ε０は真空中の誘電率、μ０は真空中の透磁率である。また、前記導電率σｘｘ，σｙｙ，σｚｚおよび導磁率σｘｘ＊，σｙｙ＊，σｚｚ＊は、前記外面０４からの距離に比例して大きくなるように設定されている。したがって、前記吸収境界０５の外面の前記吸収境界セル０６における前記導電率σｘｘ，σｙｙ，σｚｚおよび導磁率σｘｘ＊，σｙｙ＊，σｚｚ＊が十分に大きくなるよう設定することにより、前記吸収境界０５内で各ベクトル量Ｅ，Ｈを十分に減衰させることができる。このため、前記吸収境界０５の外面の前記吸収境界セル０６における前記導電率σｘｘ，σｙｙ，σｚｚは、最大減衰量と呼ばれることがある。
【００１９】
ここで、前記式（８−１）〜（８−３）は、インピーダンスマッチング条件と呼ばれ、前記吸収境界０５に入射した電磁波が前記吸収境界０５で反射せずに吸収するための条件である。また、前記式（９−１），（９−２）は、前記電磁波が前記吸収境界０５に斜めに入射した場合でも反射せずに吸収するための仮想的な条件である。この場合、前記吸収境界０５に入射した前記電磁波は、前記電磁波の周波数とは無関係に反射なしで前記媒質内に浸透する。また、浸透した前記電磁波は、前記導電率σｘｘ，σｙｙ，σｚｚ、前記導磁率σｘｘ＊，σｙｙ＊，σｚｚ＊が十分に大きい場合、速やかに減衰することになる。
【００２０】
また、前記Ｂ−ＰＭＬでは、前記電界Ｅおよび前記磁界ＨのＸ軸方向、Ｙ軸方向、Ｚ軸方向の各ベクトル成分Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚについて、以下の式（１０）〜（１５）が成立するものとする。
Ｅｘ＝Ｅｘｙ＋Ｅｘｚ                                                 …（１０）
Ｅｙ＝Ｅｙｚ＋Ｅｙｘ                                                 …（１１）
Ｅｚ＝Ｅｚｘ＋Ｅｚｙ                                                 …（１２）
Ｈｘ＝Ｈｘｙ＋Ｈｘｚ                                                 …（１３）
Ｈｙ＝Ｈｙｚ＋Ｈｙｘ                                                 …（１４）
Ｈｚ＝Ｈｚｘ＋Ｈｚｙ                                                 …（１５）
この場合、各ベクトル成分Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚを分解した成分Ｅｘｙ，Ｅｘｚ，Ｅｙｚ，Ｅｙｘ，Ｅｚｘ，Ｅｚｙ，Ｈｘｙ，Ｈｘｚ，Ｈｙｚ，Ｈｙｘ，Ｈｚｘ，Ｈｚｙについて、前記マクスウェル回転方程式系、すなわち、式（１），（２）が、以下の式（１０−１）〜（１５−１），（１０−２）〜（１５−２）で表現できる。
【００２１】
（∂／∂ｔ）Ｅｘｙ＋σｙｙＥｘｙ＝（∂／∂ｙ）Ｈｚ                …（１０−１）
（∂／∂ｔ）Ｅｘｚ＋σｚｚＥｘｚ＝−（∂／∂ｚ）Ｈｙ              …（１０−２）
（∂／∂ｔ）Ｅｙｚ＋σｚｚＥｙｚ＝（∂／∂ｚ）Ｈｘ                …（１１−１）
（∂／∂ｔ）Ｅｙｘ＋σｘｘＥｙｘ＝−（∂／∂ｘ）Ｈｚ              …（１１−２）
（∂／∂ｔ）Ｅｚｘ＋σｘｘＥｚｘ＝（∂／∂ｘ）Ｈｙ                …（１２−１）
（∂／∂ｔ）Ｅｚｙ＋σｙｙＥｚｙ＝−（∂／∂ｙ）Ｈｘ              …（１２−２）
（∂／∂ｔ）Ｈｘｙ＋σｙｙ＊Ｈｘｙ＝−（∂／∂ｙ）Ｅｚ            …（１３−１）
（∂／∂ｔ）Ｈｘｚ＋σｚｚ＊Ｈｘｚ＝（∂／∂ｚ）Ｅｙ              …（１３−２）
（∂／∂ｔ）Ｈｙｚ＋σｚｚ＊Ｈｙｚ＝−（∂／∂ｚ）Ｅｘ            …（１４−１）
（∂／∂ｔ）Ｈｙｘ＋σｘｘ＊Ｈｙｘ＝（∂／∂ｘ）Ｅｚ              …（１４−２）
（∂／∂ｔ）Ｈｚｘ＋σｘｘ＊Ｈｚｘ＝−（∂／∂ｘ）Ｅｙ            …（１５−１）
（∂／∂ｔ）Ｈｚｙ＋σｙｙ＊Ｈｚｙ＝（∂／∂ｙ）Ｅｘ              …（１５−２）
【００２２】
また、前記Ｂ−ＰＭＬでは、前記ＦＤＴＤ法と同様に、前記式（８−１）〜（１３−１），（８−２）〜（１３−２）の偏微分を中心差分近似で計算する。
例えば、（ｎ−１）・Δｔのタイミングで計算されて既知の電界Ｅｎ−１のＺ軸方向のベクトル成分Ｅｚｎ−１から、次のｎ・Δｔのタイミングの電界ＥｎのＺ軸方向のベクトル成分Ｅｚｎを計算する場合、まず、前記式（１２）は、以下の式（１２′）に変換される。
Ｅｚｎ＝Ｅｚｘｎ＋Ｅｚｙｎ                                         …（１２′）
そして、前記式（１２−１），（１２−２）を、前記単位時間Δｔ間隔の中心差分近似の式に変換すると、以下の式（１２−１′）〜（１０−３′）となる。
（Ｅｚｘｎ−Ｅｚｘｎ−１）／Δｔ＋σｘｘ（Ｅｚｘｎ＋Ｅｚｘｎ−１）／２
＝（∂／∂ｘ）Ｈｙｎ−１／２                                   …（１２−１′）
（Ｅｚｙｎ−Ｅｚｙｎ−１）／Δｔ＋σｙｙ（Ｅｚｙｎ＋Ｅｚｙｎ−１）／２
＝−（∂／∂ｙ）Ｈｘｎ−１／２                                 …（１２−２′）
また、前記式（１２−１′），（１２−２′）の右辺の（∂／∂ｘ）Ｈｙｎ−１／２，（∂／∂ｙ）Ｈｘｎ−１／２は、前記式（５″），（６″）により、前記単位長Δｘ，Δｙ毎の中間差分近似の式に置換できる。
【００２３】
したがって、前記式（１２−１′）に基づいて、Δｔ，σｘｘ，Ｅｚｘｎ−１，Δｘ，Ｈｙｎ−１／２からＥｚｘｎが計算でき、前記式（１２−２′）に基づいて、Δｔ，σｙｙ，Ｅｚｙｎ−１，Δｙ，Ｈｘｎ−１／２からＥｚｙｎが計算できる。さらに、前記式（１２′）に基づいて、計算されたＥｚｘｎ，ＥｚｙｎからＥｚｎが計算できる。そして、その他のベクトル成分Ｅｘｎ，Ｅｙｎ，Ｈｘｎ＋１／２，Ｈｙｎ＋１／２，Ｈｚｎ＋１／２も、前記ベクトル成分Ｅｚｎと同様に計算できる。
【００２４】
この結果、前記吸収境界０６における前記入力波の伝搬の様子を解析するために、（ｎ−１）・Δｔのタイミングで計算されて既知のベクトル量Ｅｎ−１＝（Ｅｘｎ−１，Ｅｙｎ−１，Ｅｚｎ−１），Ｈｎ−１／２＝（Ｈｘｎ−１／２，Ｈｙｎ−１／２，Ｈｚｎ−１／２）から、次のｎ・Δｔのタイミングのベクトル量Ｅｎ＝（Ｅｘｎ，Ｅｙｎ，Ｅｚｎ），Ｈｎ＋１／２＝（Ｈｘｎ＋１／２，Ｈｙｎ＋１／２，Ｈｚｎ＋１／２）を計算することができる。なお、前記Ｂ−ＰＭＬにおいて、前記式（１０）〜（１５），（１０−１）〜（１５−１），（１０−２）〜（１５−２）を利用して計算する手法（伝搬法）は、ＳＦ（split field）法と呼ばれている。
【００２５】
この結果、前記吸収境界０５における各座標位置のベクトル量Ｅ，Ｈの時間的な変化、すなわち、前記入力波源０１から発生した前記入力波の伝搬の様子は、前記ＳＦ法によって計算できる。
なお、ＦＤＴＤ法、Ｂ−ＰＭＬについては、特許文献１や非特許文献１に記載されており、公知である。
【００２６】
【特許文献１】特開２００３−２２３４２６号公報（「０００２」〜「００３１」、図１〜図５）
【非特許文献１】宇野亨著、「ＦＤＴＤ法による電磁界およびアンテナ解析」、第４版、コロナ社、２００２年１０月、ｐ．ｉ−ｉｉ，１−４５，５０−５９，８０−１０２，２３８−２４４
【発明の開示】
【発明が解決しようとする課題】
【００２７】
（従来技術の問題点）
吸収境界条件として前記Ｂ−ＰＭＬを適用し、前記ＳＦ法で前記吸収境界０５における各座標位置のベクトル量Ｅ，Ｈの時間的な変化を計算する場合には、前記吸収境界セル０６毎に、１８個の変数（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ，Ｅｙｘ，Ｅｚｘ，Ｈｙｘ，Ｈｚｘ，Ｅｘｙ，Ｅｚｙ，Ｅｘｚ，Ｅｙｚ，Ｈｘｙ，Ｈｚｙ，Ｈｘｚ，Ｈｙｚ）が必要となる。すなわち、前記ＳＦ法は、前記吸収境界セル０６毎に、前記１８個の変数を記憶するメモリ領域を確保する必要がある。
【００２８】
ここで、図１３に示すように、前記解析領域０２が一辺２０個の立方体領域、すなわち、総数が８０００個の前記解析領域セル０３の立法体を考える（２０３＝８０００）。そして、前記吸収境界０５の層が上下左右前後にそれぞれｄ層づつある場合を考える。
この場合、前記各変数に必要なメモリ量をｍとすると、前記吸収境界０５の層が一層増えたときの、計算に必要なメモリ量の増加分は、以下の式（１６）で表現できる。
［（２０＋２ｄ）２×６＋（２０＋２ｄ）×１２＋８］×１８×ｍ     …（１６）
したがって、前記ＳＦ法は、前記吸収境界０５の層が１層増える度に、増加した吸収境界セル０６の個数（（２０＋２ｄ）２×６＋（２０＋２ｄ）×１２＋８）の１８×ｍ倍ものメモリ領域を確保する必要があるため、計算に必要なメモリ量が急激に大きくなってしまうという問題があった。
【００２９】
また、前記Ｂ−ＰＭＬと数学的に同等の解析精度を有する吸収境界条件に、Ｚ．Ｓ．Ｓａｃｋｓらによって提案されたＵ−ＰＭＬ（Uniaxial PML）が知られている（非特許文献１参照）。
吸収境界条件として前記Ｕ−ＰＭＬを適用した場合には、前記Ｂ−ＰＭＬを適用した場合と比べて、計算に必要なメモリ量が約２／３倍にまで低減できること知られている。しかしながら、前記Ｕ−ＰＭＬを適用した場合には、前記Ｂ−ＰＭＬを適用した場合と比べて、計算量が多くなるため、計算時間が長くなってしまうという問題があった。
さらに、前記Ｂ−ＰＭＬおよび前記Ｕ−ＰＭＬ以外の従来公知の吸収境界条件では、計算に必要なメモリ量や計算時間が低減できるが、解析精度が低くなってしまうという問題があった。
【００３０】
本発明は、前述の事情に鑑み、次の記載内容（Ｏ01）を技術的課題とする。
（Ｏ01）吸収境界条件としてＢｅｒｅｎｇｅｒのＰＭＬを適用した場合の解析精度を維持しつつ、計算に必要なメモリ量や計算時間を低減すること。
【課題を解決するための手段】
【００３１】
前記技術的課題に鑑み、本願発明者が鋭意研究した結果、前記ＳＦ法に替わる新たな伝搬法であるＵＰＰ（uniaxial pseudo propagation、一軸擬伝搬）法を見出し、本発明をするに至った。
【００３２】
次に、前記課題を解決した本発明を説明するが、本発明の構成要素と後述の実施例の構成要素との対応を容易に理解できるようにするため、本発明の構成要素の直後に実施例の構成要素の符号をカッコ（）で囲んで付記する。
なお、本発明を後述の実施例の符号と対応させて説明する理由は、本発明の理解を容易にするためであり、本発明の範囲を実施例に限定するものではない。
【００３３】
（第１発明）
前記技術的課題を解決するために、第１発明の電磁界解析装置（ＰＣ）は、
電界（Ｅ）および磁界（Ｈ）からなる電磁界（Ｅ，Ｈ）をＦＤＴＤ法で解析する電磁界解析装置（ＰＣ）において、
電界（Ｅ）および磁界（Ｈ）を解析する対象の領域である解析領域（２）を、電界用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））および前記電界用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））とは位置がずれた磁界用分割領域（Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））に分割し、分割された領域毎にＦＤＴＤ法により電界（Ｅ）および磁界（Ｈ）を解析する解析領域解析手段（Ｃ１Ｄ）と、
前記解析領域（２）の境界に設けられて電界（Ｅ）および磁界（Ｈ）を吸収して減衰させるための吸収境界領域（５）を、電界用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））および前記電界用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））とは位置がずれた磁界用分割領域（Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））に分割し、前記吸収境界領域（５）で分割された領域毎に電界（Ｅ）および磁界（Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）を多次元で記憶する吸収境界電磁界記憶手段（Ｃ１E3）と、
前記吸収境界領域（５）の前記電界（Ｅ）または磁界（Ｈ）のいずれか一方に対応する分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２），Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））毎に、前記多次元の特定の第１軸方向の前記分割領域の長さ（Δｘ，Δｙ，Δｚ）に対する電界（Ｅ）または磁界（Ｈ）のいずれか他方の空間差分と、解析時間間隔（Δｔ）に対する電界（Ｅ）または磁界（Ｈ）のいずれか一方の時間差分と、に基づいて、仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電界（Ｅ）または磁界（Ｈ）のいずれか一方を演算する仮想時間電磁界演算手段（Ｃ１E1a＋Ｃ１E1b）と、
前記吸収境界領域（５）の前記電界（Ｅ）または磁界（Ｈ）のいずれか一方に対応する分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２），Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））毎に、前記第１軸方向とは独立な軸方向の前記分割領域の長さ（Δｘ，Δｙ，Δｚ）に対する磁界（Ｈ）または電界（Ｅ）のいずれか他方の空間差分と、解析時間間隔（Δｔ）に対する前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電界（Ｅ）または磁界（Ｈ）のいずれか一方の時間差分とに基づいて、前記解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の電界（Ｅ）または磁界（Ｈ）のいずれか一方を演算する解析時間後電磁界演算手段（Ｃ１E2a＋Ｃ１E2b）と、
前記解析された解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の電磁界（Ｅ，Ｈ）を表示する解析電磁界表示手段（Ｃ２）と、
を備えたことを特徴とする。
【００３４】
（第１発明の作用）
前記構成要件を備えた第１発明の電磁界解析装置（ＰＣ）では、解析領域解析手段（Ｃ１Ｄ）は、電界（Ｅ）および磁界（Ｈ）を解析する対象の領域である解析領域（２）を、電界用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））および前記電界用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））とは位置がずれた磁界用分割領域（Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））に分割し、分割された領域毎にＦＤＴＤ法により電界（Ｅ）および磁界（Ｈ）を解析する。
【００３５】
吸収境界電磁界記憶手段（Ｃ１E3）は、前記解析領域（２）の境界に設けられて電界（Ｅ）および磁界（Ｈ）を吸収して減衰させるための吸収境界領域（５）を、電界用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））および前記電界用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））とは位置がずれた磁界用分割領域（Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））に分割し、前記吸収境界領域（５）で分割された領域毎に電界（Ｅ）および磁界（Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）を多次元で記憶する。仮想時間電磁界演算手段（Ｃ１E1a＋Ｃ１E1b）は、前記吸収境界領域（５）の前記電界（Ｅ）または磁界（Ｈ）のいずれか一方に対応する分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２），Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））毎に、前記多次元の特定の第１軸方向の前記分割領域の長さ（Δｘ，Δｙ，Δｚ）に対する電界（Ｅ）または磁界（Ｈ）のいずれか他方の空間差分と、解析時間間隔（Δｔ）に対する電界（Ｅ）または磁界（Ｈ）のいずれか一方の時間差分と、に基づいて、仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電界（Ｅ）または磁界（Ｈ）のいずれか一方を演算する。
【００３６】
解析時間後電磁界演算手段（Ｃ１E2a＋Ｃ１E2b）は、前記吸収境界領域（５）の前記電界（Ｅ）または磁界（Ｈ）のいずれか一方に対応する分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２），Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））毎に、前記第１軸方向とは独立な軸方向の前記分割領域の長さ（Δｘ，Δｙ，Δｚ）に対する磁界（Ｈ）または電界（Ｅ）のいずれか他方の空間差分と、解析時間間隔（Δｔ）に対する前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電界（Ｅ）または磁界（Ｈ）のいずれか一方の時間差分とに基づいて、前記解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の電界（Ｅ）または磁界（Ｈ）のいずれか一方を演算する。解析電磁界表示手段（Ｃ２）は、前記解析された解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の電磁界（Ｅ，Ｈ）を表示する。
【００３７】
したがって、第１発明の電磁界解析装置（ＰＣ）は、前記仮想時間電磁界演算手段（Ｃ１E1a＋Ｃ１E1b）で仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電磁界（Ｅ，Ｈ）を演算し、演算された前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電磁界（Ｅ，Ｈ）に基づいて、解析時間後電磁界演算手段（Ｃ１E2a＋Ｃ１E2b）で前記解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の電磁界（Ｅ，Ｈ）を演算できる。この結果、第１発明の電磁界解析装置（ＰＣ）は、吸収境界条件に従来公知のＢ−ＰＭＬを適用し、前記ＳＦ法で電磁界（Ｅ，Ｈ）を演算する電磁界解析装置と比べ、解析精度を低減させずに演算ステップ数を低減することができ、前記電磁界（Ｅ，Ｈ）を演算するときの演算時間を低減することができる。
また、第１発明の電磁界解析装置（ＰＣ）は、吸収境界電磁界記憶手段（Ｃ１E3）で多次元における前記電磁界（Ｅ，Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）を記憶することにより、前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電磁界（Ｅ，Ｈ）および前記解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の電磁界（Ｅ，Ｈ）が演算できる。
【００３８】
すなわち、前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電磁界（Ｅ，Ｈ）を演算したときに、吸収境界電磁界記憶手段（Ｃ１E3）が記憶している前記電磁界（Ｅ，Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）を消去して、前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電磁界（Ｅ，Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）に更新して記憶できる。また、前記解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の電磁界（Ｅ，Ｈ）を演算したときに、吸収境界電磁界記憶手段（Ｃ１E3）が記憶している前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における前記電磁界（Ｅ，Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）を消去して、前記解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の電磁界（Ｅ，Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）に更新して記憶できる。この結果、第１発明の電磁界解析装置（ＰＣ）は、吸収境界条件に従来公知のＢ−ＰＭＬを適用し、前記ＳＦ法で電磁界（Ｅ，Ｈ）を演算する電磁界解析装置と比べ、前記電磁界（Ｅ，Ｈ）を演算するときのメモリ量を低減することができる。
【００３９】
（第１発明の形態１）
第１発明の形態１の電磁界解析装置（ＰＣ）は、前記第１発明において、
前記解析領域（２）および前記吸収境界領域（５）の大きさ（Ｉｍａｘ，Ｊｍａｘ，Ｋｍａｘ，ｄ）、および、解析時間間隔（Δｔ）、誘電率（ε）、透磁率（μ）を含む解析条件パラメタを入力する解析条件入力画像（１１）を表示する解析条件入力画像表示手段（Ｃ１A1）と、
を備えたことを特徴とする。
【００４０】
（第１発明の形態１の作用）
前記構成要件を備えた第１発明の形態１の電磁界解析装置（ＰＣ）では、解析条件入力画像表示手段（Ｃ１A1）は、前記解析領域（２）および前記吸収境界領域（５）の大きさ（Ｉｍａｘ，Ｊｍａｘ，Ｋｍａｘ，ｄ）、および、解析時間間隔（Δｔ）、誘電率（ε）、透磁率（μ）を含む解析条件パラメタを入力する解析条件入力画像（１１）を表示する。
したがって、第１発明の形態１の電磁界解析装置（ＰＣ）では、表示された前記解析条件入力画像（１１）によって、前記解析条件パラメタを入力することができる。
【００４１】
（第２発明）
前記技術的課題を解決するために、第２発明の電磁界解析方法は、
電界（Ｅ）および磁界（Ｈ）を解析する対象の領域である解析領域（２）を、電界用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））および前記電界用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））とは位置がずれた磁界用分割領域（Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））に分割し、分割された領域毎にＦＤＴＤ法により電界（Ｅ）および磁界（Ｈ）を解析すると共に、
前記解析領域（２）の境界に設けられて電界（Ｅ）および磁界（Ｈ）を吸収して減衰させるための吸収境界領域（５）において、前記吸収境界領域（５）の前記電界（Ｅ）または磁界（Ｈ）のいずれか一方に対応する分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２），Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））毎に、多次元の特定の軸方向の前記分割領域の長さ（Δｘ，Δｙ，Δｚ）に対する電界（Ｅ）または磁界（Ｈ）のいずれか他方の空間差分と、解析時間間隔（Δｔ）に対する電界（Ｅ）または磁界（Ｈ）のいずれか一方の時間差分と、に基づいて、仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電界（Ｅ）または磁界（Ｈ）のいずれか一方を演算し、
前記解析時間間隔（Δｔ）に対する前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電界（Ｅ）または磁界（Ｈ）のいずれか一方の時間差分と、前記第一軸方向とは独立な軸方向の前記分割領域の長さ（Δｘ，Δｙ，Δｚ）に対する磁界（Ｈ）または電界（Ｅ）のいずれか他方の空間差分とに基づいて、前記解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の電界（Ｅ）または磁界（Ｈ）のいずれか一方を演算することを特徴とする。
【００４２】
（第２発明の作用）
前記構成要件を備えた第２発明の電磁界解析方法は、電界（Ｅ）および磁界（Ｈ）を解析する対象の領域である解析領域（２）を、電界用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））および前記電界用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））とは位置がずれた磁界用分割領域（Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））に分割し、分割された領域毎にＦＤＴＤ法により電界（Ｅ）および磁界（Ｈ）を解析すると共に、前記解析領域（２）の境界に設けられて電界（Ｅ）および磁界（Ｈ）を吸収して減衰させるための吸収境界領域（５）において、前記吸収境界領域（５）の前記電界（Ｅ）または磁界（Ｈ）のいずれか一方に対応する分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２），Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））毎に、多次元の特定の軸方向の前記分割領域の長さ（Δｘ，Δｙ，Δｚ）に対する電界（Ｅ）または磁界（Ｈ）のいずれか他方の空間差分と、解析時間間隔（Δｔ）に対する電界（Ｅ）または磁界（Ｈ）のいずれか一方の時間差分と、に基づいて、仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電界（Ｅ）または磁界（Ｈ）のいずれか一方を演算し、前記解析時間間隔（Δｔ）に対する前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電界（Ｅ）または磁界（Ｈ）のいずれか一方の時間差分と、前記第一軸方向とは独立な軸方向の前記分割領域の長さ（Δｘ，Δｙ，Δｚ）に対する磁界（Ｈ）または電界（Ｅ）のいずれか他方の空間差分とに基づいて、前記解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の電界（Ｅ）または磁界（Ｈ）のいずれか一方を演算する。
【００４３】
したがって、第２発明の電磁界解析方法は、仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電磁界（Ｅ，Ｈ）を演算し、演算された前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電磁界（Ｅ，Ｈ）に基づいて、前記解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の電磁界（Ｅ，Ｈ）を演算できる。この結果、第２発明の電磁界解析方法は、吸収境界条件に従来公知のＢ−ＰＭＬを適用し、前記ＳＦ法で電磁界（Ｅ，Ｈ）を演算する電磁界解析方法と比べ、解析精度を低減させずに演算ステップ数を低減することができ、前記電磁界（Ｅ，Ｈ）を演算するときの演算時間を低減することができる。
【００４４】
また、第２発明の電磁界解析方法は、前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電磁界（Ｅ，Ｈ）を演算したときに、演算に必要であった前記電磁界（Ｅ，Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）を、前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電磁界（Ｅ，Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）に更新できる。また、前記解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の電磁界（Ｅ，Ｈ）を演算したときに、演算に必要であった前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における前記電磁界（Ｅ，Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）を、前記解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の電磁界（Ｅ，Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）に更新できる。この結果、第２発明の電磁界解析方法は、吸収境界条件に従来公知のＢ−ＰＭＬを適用し、前記ＳＦ法で電磁界（Ｅ，Ｈ）を演算する電磁界解析方法と比べ、前記電磁界（Ｅ，Ｈ）を演算するときのメモリ量を低減することができる。
【００４５】
（第３発明）
前記技術的課題を解決するために、第３発明の電磁界解析プログラム（ＡＰ１）は、
電界（Ｅ）および磁界（Ｈ）からなる電磁界（Ｅ，Ｈ）をＦＤＴＤ法で解析する電磁界解析装置（ＰＣ）を、
電界（Ｅ）および磁界（Ｈ）を解析する対象の領域である解析領域（２）を、電界用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））および前記電界用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））とは位置がずれた磁界用分割領域（Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））に分割し、分割された領域毎にＦＤＴＤ法により電界（Ｅ）および磁界（Ｈ）を解析する解析領域解析手段（Ｃ１Ｄ）、
前記解析領域（２）の境界に設けられて電界（Ｅ）および磁界（Ｈ）を吸収して減衰させるための吸収境界領域（５）を、電界用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））および前記電界用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））とは位置がずれた磁界用分割領域（Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））に分割し、前記吸収境界領域（５）で分割された領域毎に電界（Ｅ）および磁界（Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）を多次元で記憶する吸収境界電磁界記憶手段（Ｃ１E3）、
前記吸収境界領域（５）の前記電界（Ｅ）または磁界（Ｈ）のいずれか一方に対応する分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２），Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））毎に、前記多次元の特定の第１軸方向の前記分割領域の長さ（Δｘ，Δｙ，Δｚ）に対する電界（Ｅ）または磁界（Ｈ）のいずれか他方の空間差分と、解析時間間隔（Δｔ）に対する電界（Ｅ）または磁界（Ｈ）のいずれか一方の時間差分と、に基づいて、仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電界（Ｅ）または磁界（Ｈ）のいずれか一方を演算する仮想時間電磁界演算手段（Ｃ１E1a＋Ｃ１E1b）、
前記吸収境界領域（５）の前記電界（Ｅ）または磁界（Ｈ）のいずれか一方に対応する分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２），Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））毎に、前記第１軸方向とは独立な軸方向の前記分割領域の長さ（Δｘ，Δｙ，Δｚ）に対する磁界（Ｈ）または電界（Ｅ）のいずれか他方の空間差分と、解析時間間隔（Δｔ）に対する前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電界（Ｅ）または磁界（Ｈ）のいずれか一方の時間差分とに基づいて、前記解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の電界（Ｅ）または磁界（Ｈ）のいずれか一方を演算する解析時間後電磁界演算手段（Ｃ１E2a＋Ｃ１E2b）、
前記解析された解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の電磁界（Ｅ，Ｈ）を表示する解析電磁界表示手段（Ｃ２）、
として機能させる。
【００４６】
（第３発明の作用）
前記構成要件を備えた第３発明の電磁界解析プログラム（ＡＰ１）では、解析領域解析手段（Ｃ１Ｄ）は、電界（Ｅ）および磁界（Ｈ）を解析する対象の領域である解析領域（２）を、電界用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））および前記電界用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））とは位置がずれた磁界用分割領域（Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））に分割し、分割された領域毎にＦＤＴＤ法により電界（Ｅ）および磁界（Ｈ）を解析する。
【００４７】
吸収境界電磁界記憶手段（Ｃ１E3）は、前記解析領域（２）の境界に設けられて電界（Ｅ）および磁界（Ｈ）を吸収して減衰させるための吸収境界領域（５）を、電界用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））および前記電界用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））とは位置がずれた磁界用分割領域（Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））に分割し、前記吸収境界領域（５）で分割された領域毎に電界（Ｅ）および磁界（Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）を多次元で記憶する。仮想時間電磁界演算手段（Ｃ１E1a＋Ｃ１E1b）は、前記吸収境界領域（５）の前記電界（Ｅ）または磁界（Ｈ）のいずれか一方に対応する分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２），Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））毎に、前記多次元の特定の第１軸方向の前記分割領域の長さ（Δｘ，Δｙ，Δｚ）に対する電界（Ｅ）または磁界（Ｈ）のいずれか他方の空間差分と、解析時間間隔（Δｔ）に対する電界（Ｅ）または磁界（Ｈ）のいずれか一方の時間差分と、に基づいて、仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電界（Ｅ）または磁界（Ｈ）のいずれか一方を演算する。
【００４８】
解析時間後電磁界演算手段（Ｃ１E2a＋Ｃ１E2b）は、前記吸収境界領域（５）の前記電界（Ｅ）または磁界（Ｈ）のいずれか一方に対応する分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２），Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））毎に、前記第１軸方向とは独立な軸方向の前記分割領域の長さ（Δｘ，Δｙ，Δｚ）に対する磁界（Ｈ）または電界（Ｅ）のいずれか他方の空間差分と、解析時間間隔（Δｔ）に対する前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電界（Ｅ）または磁界（Ｈ）のいずれか一方の時間差分とに基づいて、前記解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の電界（Ｅ）または磁界（Ｈ）のいずれか一方を演算する。解析電磁界表示手段（Ｃ２）は、前記解析された解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の電磁界（Ｅ，Ｈ）を表示する。
【００４９】
したがって、第３発明の電磁界解析プログラム（ＡＰ１）は、前記仮想時間電磁界演算手段（Ｃ１E1a＋Ｃ１E1b）で仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電磁界（Ｅ，Ｈ）を演算し、演算された前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電磁界（Ｅ，Ｈ）に基づいて、解析時間後電磁界演算手段（Ｃ１E2a＋Ｃ１E2b）で前記解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の電磁界（Ｅ，Ｈ）を演算できる。この結果、第３発明の電磁界解析プログラム（ＡＰ１）は、吸収境界条件に従来公知のＢ−ＰＭＬを適用し、前記ＳＦ法で電磁界（Ｅ，Ｈ）を演算する電磁界解析プログラムと比べ、解析精度を低減させずに演算ステップ数を低減することができ、前記電磁界（Ｅ，Ｈ）を演算するときの前記電磁界解析装置（ＰＣ）における演算時間を低減することができる。
また、第３発明の電磁界解析プログラム（ＡＰ１）は、吸収境界電磁界記憶手段（Ｃ１E3）で多次元における前記電磁界（Ｅ，Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）を記憶することにより、前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電磁界（Ｅ，Ｈ）および前記解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の電磁界（Ｅ，Ｈ）が演算できる。
【００５０】
すなわち、前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電磁界（Ｅ，Ｈ）を演算したときに、吸収境界電磁界記憶手段（Ｃ１E3）が記憶している前記電磁界（Ｅ，Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）を消去して、前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における電磁界（Ｅ，Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）に更新して記憶できる。また、前記解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の電磁界（Ｅ，Ｈ）を演算したときに、吸収境界電磁界記憶手段（Ｃ１E3）が記憶している前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における前記電磁界（Ｅ，Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）を消去して、前記解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の電磁界（Ｅ，Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）に更新して記憶できる。この結果、第３発明の電磁界解析プログラム（ＡＰ１）は、吸収境界条件に従来公知のＢ−ＰＭＬを適用し、前記ＳＦ法で電磁界（Ｅ，Ｈ）を演算する電磁界解析装置と比べ、前記電磁界（Ｅ，Ｈ）を演算するときの前記電磁界解析装置（ＰＣ）におけるメモリ量を低減することができる。
【００５１】
（第４発明）
前記技術的課題を解決するために、第４発明の場解析装置（ＰＣ）は、
第１要素（Ｅ）の時間差分が第２要素（Ｈ）の空間差分で表現され、前記第２要素（Ｈ）の時間差分が前記第１要素（Ｅ）の空間差分で表現される場（Ｅ，Ｈ）をＦＤＴＤ法で解析する場解析装置（ＰＣ）において、
前記第１要素（Ｅ）および前記第２要素（Ｈ）を解析する対象の領域である解析領域（２）を、第１要素用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））および前記第１要素用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））とは位置がずれた第２要素用分割領域（Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））に分割し、分割された領域毎にＦＤＴＤ法により前記第１要素（Ｅ）および前記第２要素（Ｈ）を解析する解析領域解析手段（Ｃ１Ｄ）と、
前記解析領域（２）の境界に設けられて前記第１要素（Ｅ）および前記第２要素（Ｈ）を吸収して減衰させるための吸収境界領域（５）を、第１要素用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））および前記第１要素用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））とは位置がずれた第２要素用分割領域（Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））に分割し、前記吸収境界領域（５）で分割された領域毎に前記第１要素（Ｅ）および前記第２要素（Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）を多次元で記憶する吸収境界要素記憶手段（Ｃ１E3）と、
前記吸収境界領域（５）の前記第１要素（Ｅ）または前記第２要素（Ｈ）のいずれか一方に対応する分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２），Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））毎に、前記多次元の特定の第１軸方向の前記分割領域の長さ（Δｘ，Δｙ，Δｚ）に対する前記第１要素（Ｅ）または前記第２要素（Ｈ）のいずれか他方の空間差分と、解析時間間隔（Δｔ）に対する前記第１要素（Ｅ）または前記第２要素（Ｈ）のいずれか一方の時間差分と、に基づいて、仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における前記第１要素（Ｅ）または前記第２要素（Ｈ）のいずれか一方を演算する仮想時間要素演算手段（Ｃ１E1a＋Ｃ１E1b）と、
前記吸収境界領域（５）の前記第１要素（Ｅ）または前記第２要素（Ｈ）のいずれか一方に対応する分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２），Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））毎に、前記第１軸方向とは独立な軸方向の前記分割領域の長さ（Δｘ，Δｙ，Δｚ）に対する前記第２要素（Ｈ）または前記第１要素（Ｅ）のいずれか他方の空間差分と、解析時間間隔（Δｔ）に対する前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における前記第１要素（Ｅ）または前記第２要素（Ｈ）のいずれか一方の時間差分とに基づいて、前記解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の前記第１要素（Ｅ）または前記第２要素（Ｈ）のいずれか一方を演算する解析時間後要素演算手段（Ｃ１E2a＋Ｃ１E2b）と、
前記解析された解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の前記場（Ｅ，Ｈ）を表示する解析場表示手段（Ｃ２）と、
を備えたことを特徴とする。
【００５２】
（第４発明の作用）
前記構成要件を備えた第４発明の場解析装置（ＰＣ）では、解析領域解析手段（Ｃ１Ｄ）は、前記第１要素（Ｅ）および前記第２要素（Ｈ）を解析する対象の領域である解析領域（２）を、第１要素用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））および前記第１要素用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））とは位置がずれた第２要素用分割領域（Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））に分割し、分割された領域毎にＦＤＴＤ法により前記第１要素（Ｅ）および前記第２要素（Ｈ）を解析する。
【００５３】
吸収境界要素記憶手段（Ｃ１E3）は、前記解析領域（２）の境界に設けられて前記第１要素（Ｅ）および前記第２要素（Ｈ）を吸収して減衰させるための吸収境界領域（５）を、第１要素用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））および前記第１要素用分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２））とは位置がずれた第２要素用分割領域（Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））に分割し、前記吸収境界領域（５）で分割された領域毎に前記第１要素（Ｅ）および前記第２要素（Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）を多次元で記憶する。仮想時間要素演算手段（Ｃ１E1a＋Ｃ１E1b）は、前記吸収境界領域（５）の前記第１要素（Ｅ）または前記第２要素（Ｈ）のいずれか一方に対応する分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２），Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））毎に、前記多次元の特定の第１軸方向の前記分割領域の長さ（Δｘ，Δｙ，Δｚ）に対する前記第１要素（Ｅ）または前記第２要素（Ｈ）のいずれか他方の空間差分と、解析時間間隔（Δｔ）に対する前記第１要素（Ｅ）または前記第２要素（Ｈ）のいずれか一方の時間差分と、に基づいて、仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における前記第１要素（Ｅ）または前記第２要素（Ｈ）のいずれか一方を演算する。
【００５４】
解析時間後要素演算手段（Ｃ１E2a＋Ｃ１E2b）は、前記吸収境界領域（５）の前記第１要素（Ｅ）または前記第２要素（Ｈ）のいずれか一方に対応する分割領域（Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２），Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ））毎に、前記第１軸方向とは独立な軸方向の前記分割領域の長さ（Δｘ，Δｙ，Δｚ）に対する前記第２要素（Ｈ）または前記第１要素（Ｅ）のいずれか他方の空間差分と、解析時間間隔（Δｔ）に対する前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における前記第１要素（Ｅ）または前記第２要素（Ｈ）のいずれか一方の時間差分とに基づいて、前記解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の前記第１要素（Ｅ）または前記第２要素（Ｈ）のいずれか一方を演算する。
解析場表示手段（Ｃ２）は、前記解析された解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の前記場（Ｅ，Ｈ）を表示する。
【００５５】
したがって、第４発明の場解析装置（ＰＣ）は、前記仮想時間要素演算手段（Ｃ１E1a＋Ｃ１E1b）で仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における場（Ｅ，Ｈ）を演算し、演算された前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における場（Ｅ，Ｈ）に基づいて、解析時間後要素演算手段（Ｃ１E2a＋Ｃ１E2b）で前記解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の場（Ｅ，Ｈ）を演算できる。この結果、第４発明の場解析装置（ＰＣ）は、吸収境界条件に従来公知のＢ−ＰＭＬを適用し、前記ＳＦ法で場（Ｅ，Ｈ）を演算する場解析装置と比べ、解析精度を低減させずに演算ステップ数を低減することができ、前記場（Ｅ，Ｈ）を演算するときの演算時間を低減することができる。
また、第４発明の場解析装置（ＰＣ）は、吸収境界要素記憶手段（Ｃ１E3）で多次元における前記場（Ｅ，Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）を記憶することにより、前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における場（Ｅ，Ｈ）および前記解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の場（Ｅ，Ｈ）が演算できる。
【００５６】
すなわち、前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における場（Ｅ，Ｈ）を演算したときに、吸収境界要素記憶手段（Ｃ１E3）が記憶している前記場（Ｅ，Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）を消去して、前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における場（Ｅ，Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）に更新して記憶できる。また、前記解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の場（Ｅ，Ｈ）を演算したときに、吸収境界電磁界記憶手段（Ｃ１E3）が記憶している前記仮想時間（ｎ′Δｔ，（ｎ＋１／２）′Δｔ）における前記場（Ｅ，Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）を消去して、前記解析時間後（ｎ・Δｔ，（ｎ＋１／２）Δｔ）の場（Ｅ，Ｈ）の強さ（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ）に更新して記憶できる。この結果、第４発明の場解析装置（ＰＣ）は、吸収境界条件に従来公知のＢ−ＰＭＬを適用し、前記ＳＦ法で場（Ｅ，Ｈ）を演算する場解析装置と比べ、前記場（Ｅ，Ｈ）を演算するときのメモリ量を低減することができる。
【発明の効果】
【００５７】
前述の本発明の電磁界解析装置、電磁界解析方法、電磁界解析プログラムおよび場解析装置は、下記の効果（Ｅ01）を奏する。
（Ｅ01）吸収境界条件としてＢｅｒｅｎｇｅｒのＰＭＬを適用した場合の解析精度を維持しつつ、計算に必要なメモリ量や計算時間を低減することができる。
【発明を実施するための最良の形態】
【００５８】
次に図面を参照しながら、本発明の実施の形態の具体例（実施例）を説明するが、本発明は以下の実施例に限定されるものではない。また、以下の図面を使用した説明において、理解の容易のために説明に必要な部材以外の図示は適宜省略されている。
なお、以後の説明の理解を容易にするために、図面において、前後方向をＸ軸方向、左右方向をＹ軸方向、上下方向をＺ軸方向とし、矢印Ｘ，−Ｘ，Ｙ，−Ｙ，Ｚ，−Ｚで示す方向または示す側をそれぞれ、前方、後方、右方、左方、上方、下方、または、前側、後側、右側、左側、上側、下側とする。
また、図中、「○」の中に「・」が記載されたものは紙面の裏から表に向かう矢印を意味し、「○」の中に「×」が記載されたものは紙面の表から裏に向かう矢印を意味するものとする。
【実施例１】
【００５９】
図１は本発明の実施例１の電磁界解析装置の説明図である。
図１において、本発明の実施例１の電磁界解析装置（場演算装置、電磁場演算装置）としてのクライアントパソコン（パーソナルコンピュータ）ＰＣは、コンピュータ装置により構成されており、コンピュータ本体Ｈ１と、ディスプレイ（情報表示画面）Ｈ２と、キーボードＨ３やマウスＨ４等の入力装置、図示しないＨＤドライブ（ハードディスクドライブ）、ＣＤドライブ（コンパクトディスクドライブ）等により構成されている。
【００６０】
（実施例１の制御部の説明）
（クライアントパソコンＰＣの制御部の説明）
図２は実施例１の電磁界解析装置としてのクライアントパソコンの制御部分が備えている各機能をブロック図（機能ブロック図）で示した図である。
図２において、クライアントパソコンＰＣのコンピュータ本体Ｈ１は、外部との信号の入出力および入出力信号レベルの調節等を行うＩ／Ｏ（入出力インターフェース）、必要な処理を行うためのプログラムおよびデータ等が記憶されたＲＯＭ（リードオンリーメモリ）、必要なデータを一時的に記憶するためのＲＡＭ（ランダムアクセスメモリ）、ハードディスクやＲＯＭ等に記憶されたプログラムに応じた処理を行うＣＰＵ（中央演算処理装置）、ならびにクロック発振器等を有している。
前記構成のクライアントパソコンＰＣは、前記ハードディスクやＲＯＭ等に記憶されたプログラムを実行することにより種々の機能を実現することができる。
【００６１】
前記クライアントパソコンＰＣの前記ハードディスクドライブには、クライアントパソコンの基本動作を制御する基本ソフト（オペレーティングシステム）ＯＳや、アプリケーションプログラムとしての電磁界解析表示プログラムＡＰ１、その他の図示しないソフトウェア（文書作成ソフトウェアや製図ソフトウェア等）が記憶されている。
【００６２】
（電磁界解析表示プログラムＡＰ１）
電磁界解析表示プログラム（電磁界解析プログラム）ＡＰ１は、下記の機能手段（プログラムモジュール）を有する。
図３は実施例１の電磁界解析表示プログラムＡＰ１により解析される解析領域および吸収境界の説明図であり、図３Ａは解析領域および吸収境界全体の説明図、図３Ｂは解析領域および吸収境界の各セルの拡大説明図である。
Ｃ１：電磁界解析手段
電磁界解析手段Ｃ１は、解析開始判別手段Ｃ１Ａと，解析条件パラメタ記憶手段Ｃ１Ｂと、入力波源計算手段Ｃ１Ｃと、解析領域解析手段Ｃ１Ｄと、吸収境界解析手段Ｃ１Ｅと、解析結果出力手段Ｃ１Ｆと、解析終了判別手段Ｃ１Ｇとを有し、入力波源１（図３Ａ参照）から発生する電磁波、すなわち、電磁界（場、電磁場）を解析する電磁界解析処理を実行する。
【００６３】
図３Ａ、図３Ｂにおいて、実施例１の前記入力波源１は、Ｚ軸方向の成分のみの電界Ｅ＝（０，０，Ｅｚ）に振幅を与える点波源により構成されている。すなわち、実施例１では、前記入力波源１の電界ＥのＺ軸方向の成分Ｅｚについて、以下の式（１７）が成立する。
Ｅｚ＝Ｅ０×ｓｉｎ（ω０ｔ）                                        …（１７）
なお、Ｅ０はＥｚの振幅、ω０はＥｚの角周波数である。
また、前記入力波源１は、前記入力波源１から発生する電磁波（入力波）を解析する解析領域２で囲まれており、前記解析領域２は、微小直方体である解析領域セル３により分割されている。
また、前記解析領域２の外方は、前記入力波源１の前記電磁波を吸収する吸収境界（吸収境界領域）５で囲まれており、前記吸収境界５は、前記解析領域セル３と同様の吸収境界セル６により分割されている。
【００６４】
また、実施例１では、Ｘ軸方向、Ｙ軸方向、Ｚ軸方向の単位長をそれぞれΔｘ，Δｙ，Δｚとし、前記セル３，６は、図３Ｂで示すように、Δｘ×Δｙ×Δｚのサイズの微小直方体とする。また、前記セル３，６には、前記電界（第１要素）ＥについてのＸ軸方向、Ｙ軸方向、Ｚ軸方向のベクトル成分であるＥｘ，Ｅｙ，Ｅｚが前記セル３，６の各辺の中央部から前記各辺に沿って割り当てられ、前記磁界（第２要素）ＨについてのＸ軸方向、Ｙ軸方向、Ｚ軸方向のベクトル成分であるＨｘ，Ｈｙ，Ｈｚが前記セル３，６の各面の中心から垂直方向に割り当てられている。すなわち、整数ｉ，ｊ，ｋを用いて、前記セル３，６の各頂点の座標を（ｉ，ｊ，ｋ）で表現した場合、前記各頂点（ｉ，ｊ，ｋ）におけるベクトル量Ｅ（ｉ，ｊ，ｋ），Ｈ（ｉ，ｊ，ｋ）の各ベクトル成分Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚは、前記各頂点（ｉ，ｊ，ｋ）に最も近い点で割り当てられる。すなわち、前記各ベクトル成分Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚは、それぞれＥｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２），Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ）で表現できる。
【００６５】
なお、実施例１の前記入力波源（点波源）１は、前記解析領域２の重心位置に配置されているため、重心位置の座標を（ｉＧ，ｊＧ，ｋＧ）とした場合、前記式（１７）は以下の式（１７′）で表現できる。
Ｅｚ（ｉＧ，ｊＧ，ｋＧ＋１／２）＝Ｅ０×ｓｉｎ（ω０ｔ）            …（１７′）
また、図３Ｂに示すように、前記電界Ｅおよび前記磁界Ｈが割り当てられた場合、前記電界Ｅおよび前記磁界Ｈは、時間的に交互に配置されることから、実施例１では、（ｎ−１）・Δｔ，ｎ・Δｔの各タイミングにおける電界ＥをそれぞれＥｎ−１，Ｅｎとし、（ｎ−１／２）・Δｔ，（ｎ＋１／２）・Δｔの各タイミングにおける磁界ＨをそれぞれＨｎ−１／２，Ｈｎ＋１／２と表現する。
【００６６】
なお、実施例１の前記解析領域２は、Ｘ軸方向、Ｙ軸方向、Ｚ軸方向の長さがそれぞれＩｍａｘ，Ｊｍａｘ，Ｋｍａｘの直方体領域とする。すなわち、前記解析領域２は、総数が（Ｉｍａｘ／Δｘ）×（Ｊｍａｘ／Δｙ）×（Ｋｍａｘ／Δｚ）個の前記解析領域セル３の直方体とする。また、前記解析領域２は、図３Ａで示すように、上下左右前後にそれぞれｄ層の吸収境界セル６で囲まれているものとする。
また、実施例１では、前記解析領域２の誘電率をε、透磁率をμとする。また、前記吸収境界５の導電率をσｘｘ，σｘｙ，σｘｚ，σｙｘ，σｙｙ，σｙｚ，σｚｘ，σｚｙ，σｚｚ、導磁率をσｘｘ＊，σｘｙ＊，σｘｚ＊，σｙｘ＊，σｙｙ＊，σｙｚ＊，σｚｘ＊，σｚｙ＊，σｚｚ＊とし、前記式（８−１′）〜（８−３′）および式（９−１′），（９−２′）が成立するものとする。
【００６７】
また、実施例１では、前記導電率σｘｘ，σｙｙ，σｚｚについては、各層ごとにσｘｘ＝σｙｙ＝σｚｚが成立するものとし、前記導電率σｘｘ，σｙｙ，σｚｚの最大減衰量をσｍａｘとする。なお、前記最大減衰量σｍａｘは、前記吸収境界５の外面となる吸収境界セル６、すなわち、ｄ層目の吸収境界セル６の導電率σｘｘ，σｙｙ，σｚｚの値である。よって、実施例１では、前記式（８−１）〜（８−３）および前記最大減衰量σｍａｘに基づいて、全ての吸収境界セル０６における導電率をσｘｘ，σｙｙ，σｚｚおよび導磁率をσｘｘ＊，σｙｙ＊，σｚｚ＊が設定される。
さらに、実施例１では、単位時間をΔｔ、前記単位時間Δｔ毎に解析する解析回数の最大値をＮｍａｘとする。
【００６８】
前記入力波源１のＥｚ（ｉＧ，ｊＧ，ｋＧ＋１／２）の振幅Ｅ０および角周波数ω０、前記セル３，６の各辺の長さΔｘ，Δｙ，Δｚ、前記解析領域２の各辺の長さＩｍａｘ，Ｊｍａｘ，Ｋｍａｘ、前記吸収境界５の層の値ｄ、前記解析領域２の誘電率εおよび透磁率μ、前記吸収境界５の導電率σｘｘ，σｙｙ，σｚｚおよび導磁率σｘｘ＊，σｙｙ＊，σｚｚ＊を設定するための前記最大減衰量σｍａｘ、前記単位時間Δｔ、前記解析回数の最大値Ｎｍａｘ等とにより実施例１の解析条件パラメタが構成されている。
【００６９】
Ｃ１Ａ：解析開始判別手段
解析開始判別手段Ｃ１Ａは、解析条件入力画像表示手段Ｃ１A1を有し、前記電磁界解析処理を開始するか否かを判別する。実施例１の前記解析開始判別手段Ｃ１Ａは、ユーザにより前記解析条件パラメタが入力されたか否かを判別することにより、前記電磁界解析処理を開始するか否かを判別する。
図４は実施例１の電磁界解析表示プログラムが表示する解析条件パラメタ入力画像の説明図である。
Ｃ１A1：解析条件入力画像表示手段
解析条件入力画像表示手段Ｃ１A1は、ユーザが前記解析条件パラメタを入力するための解析条件入力画像１１を前記ディスプレイＨ２に表示する。
図４において、実施例１の前記解析条件入力画像１１は、前記入力波源１のＥｚ（ｉＧ，ｊＧ，ｋＧ＋１／２）の振幅Ｅ０を入力するための入力波源振幅入力欄１２と、前記入力波源１のＥｚ（ｉＧ，ｊＧ，ｋＧ＋１／２）の角周波数ω０を入力するための入力波源角周波数入力欄１３と、前記解析領域２の各辺の長さＩｍａｘ，Ｊｍａｘ，Ｋｍａｘを入力するための解析領域サイズ入力欄１４，１６，１７と、前記セル３，６の各辺の長さΔｘ，Δｙ，Δｚを入力するためのセルサイズ入力欄１８，１９，２１とを有している。
【００７０】
また、前記解析条件入力画像１１は、ユーザが前記解析領域２の誘電率εを入力するための誘電率入力欄２２と、前記解析領域２の透磁率μを入力するための透磁率入力欄２３と、前記吸収境界５の層の値ｄを入力するための吸収境界層入力欄２４と、前記吸収境界５の導電率σｘｘ，σｙｙ，σｚｚおよび導磁率σｘｘ＊，σｙｙ＊，σｚｚ＊を設定するために必要となる導電率σｘｘ，σｙｙ，σｚｚの最大減衰量σｍａｘを入力するための最大源推量入力欄２６と、前記単位時間Δｔを入力するための単位時間入力欄２７と、前記解析回数ｎの最大値Ｎｍａｘを入力するための解析回数最大値入力欄２８と、入力した内容で前記電磁界解析処理を開始するための解析開始画像２９（「解析開始」アイコン参照）とを有する。
なお、実施例１の前記解析開始判別手段Ｃ１Ａは、具体的には、前記解析条件入力画像１１の前記各入力欄１２〜２８に適切な値が入力され、且つ、前記解析開始画像２９が選択されたか否かを判別することにより、前記電磁界解析処理を開始するか否かを判別する。
【００７１】
Ｃ１Ｂ：解析条件パラメタ記憶手段
解析条件パラメタ記憶手段Ｃ１Ｂは、前記解析開始判別手段Ｃ１Ａで入力された前記解析条件パラメタを記憶する。
Ｃ１Ｃ：入力波源計算手段
入力波源計算手段Ｃ１Ｃは、各タイミングの前記入力波源１の電界を計算する。実施例１の前記入力波源計算手段Ｃ１Ｃは、前記式（１７′）に基づいて、各タイミングのＥｚ（ｉＧ，ｊＧ，ｋＧ＋１／２）を計算する。具体的には、前記入力波源計算手段Ｃ１Ｃは、（ｎ−１）・ΔｔのタイミングのＥｚｎ−１（ｉＧ，ｊＧ，ｋＧ＋１／２）を、以下の式（１７″）に基づいて計算する。
Ｅｚｎ−１（ｉＧ，ｊＧ，ｋＧ＋１／２）
＝Ｅ０×ｓｉｎ（ω０（ｎ−１）・Δｔ）                           …（１７″）
Ｃ１Ｄ：解析領域解析手段
解析領域解析手段Ｃ１Ｄは、解析領域電界演算手段Ｃ１D1と、解析領域磁界演算手段Ｃ１D2と、解析領域電磁界記憶手段Ｃ１D3とを有し、前記入力波源１から発生した前記電磁波を前記解析領域２で解析する解析領域解析処理を実行する。実施例１の前記解析領域解析手段Ｃ１Ｄは、前記ＦＤＴＤ法のＹｅｅアルゴリズムに基づいて、（ｎ−１）・Δｔのタイミングで計算されて既知のベクトル量Ｅｎ−１＝（Ｅｘｎ−１，Ｅｙｎ−１，Ｅｚｎ−１），Ｈｎ−１／２＝（Ｈｘｎ−１／２，Ｈｙｎ−１／２，Ｈｚｎ−１／２）から、次のｎ・Δｔのタイミングのベクトル量Ｅｎ＝（Ｅｘｎ，Ｅｙｎ，Ｅｚｎ），Ｈｎ＋１／２＝（Ｈｘｎ＋１／２，Ｈｙｎ＋１／２，Ｈｚｎ＋１／２）を計算する前記解析領域解析処理を実行する。
【００７２】
なお、前記入力波源１が前記点波源により構成されている場合、前記解析領域２の各ベクトル成分Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚは、前記式（３），（４）に基づいて計算できることが知られている。したがって、前記解析領域解析手段Ｃ１Ｄは、前記式（３−１′），（５），（６）に基づいて、Ｈｘｎ−１／２（ｉ，ｊ＋１／２，ｋ＋１／２），Δｔ，μ，Δｙ，Ｅｚｎ（ｉ，ｊ，ｋ＋１／２），Ｅｚｎ（ｉ，ｊ＋１，ｋ＋１／２），Δｚ，Ｅｙｎ（ｉ，ｊ＋１／２，ｋ），Ｅｙｎ（ｉ，ｊ＋１／２，ｋ＋１）からＨｘｎ＋１／２（ｉ，ｊ＋１／２，ｋ＋１／２）を計算する。また、前記解析領域解析手段Ｃ１Ｄは、その他のベクトル成分Ｅｘｎ（ｉ＋１／２，ｊ，ｋ），Ｅｙｎ（ｉ，ｊ＋１／２，ｋ），Ｅｚｎ（ｉ，ｊ，ｋ＋１／２），Ｈｙｎ＋１／２（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚｎ＋１／２（ｉ＋１／２，ｊ＋１／２，ｋ）も、Ｈｘｎ＋１／２（ｉ，ｊ＋１／２，ｋ＋１／２）と同様に計算する。
【００７３】
具体的には、実施例１の前記解析領域解析手段Ｃ１Ｄでは、前記解析領域２の全ての座標（ｉ，ｊ，ｋ）において、Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２），Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ）は、ｎ＝０の場合に０である（Ｅｘ０＝Ｅｙ０＝Ｅｚ０＝Ｈｘ０＋１／２＝Ｈｙ０＋１／２＝Ｈｚ０＋１／２＝０）。また、前記解析領域解析手段Ｃ１Ｄは、ｎ＝１，２，…，Ｎｍａｘの場合には、前記式（３），（４）に基づいて、前記入力波源計算手段Ｃ１Ｃで計算されたＥｚｎ−１（ｉＧ，ｊＧ，ｋＧ＋１／２）の値に応じた全てのＥｘｎ，Ｅｙｎ，Ｅｚｎ，Ｈｘｎ＋１／２，Ｈｙｎ＋１／２，Ｈｚｎ＋１／２の値が計算される。
【００７４】
Ｃ１D1：解析領域電界計算手段
解析領域電界演算手段Ｃ１D1は、前記解析領域２において、（ｎ−１）・Δｔのタイミングで計算された電界Ｅｎ−１＝（Ｅｘｎ−１，Ｅｙｎ−１，Ｅｚｎ−１）等から、次のｎ・Δｔのタイミングの電界Ｅｎ＝（Ｅｘｎ，Ｅｙｎ，Ｅｚｎ）を計算する解析領域電界計算処理を実行する。実施例１の前記解析領域電界演算手段Ｃ１D1は、前記式（４）に基づいて、前記解析領域２の全てのＥｘｎ（ｉ＋１／２，ｊ，ｋ），Ｅｙｎ（ｉ，ｊ＋１／２，ｋ），Ｅｚｎ（ｉ，ｊ，ｋ＋１／２）を計算する。
Ｃ１D2：解析領域磁界計算手段
解析領域磁界解析手段Ｃ１D2は、前記解析領域２において、（ｎ−１／２）・Δｔのタイミングで計算された磁界Ｈｎ−１／２＝（Ｈｘｎ−１／２，Ｈｙｎ−１／２，Ｈｚｎ−１／２）等から、次の（ｎ＋１／２）・Δｔのタイミングの磁界Ｈｎ＋１／２＝（Ｈｘｎ＋１／２，Ｈｙｎ＋１／２，Ｈｚｎ＋１／２）を計算する解析領域磁界計算処理を実行する。実施例１の前記解析領域磁界演算手段Ｃ１D2は、前記式（３）に基づいて、前記解析領域２の全てのＨｘｎ＋１／２（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙｎ＋１／２（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚｎ＋１／２（ｉ＋１／２，ｊ＋１／２，ｋ）を計算する。
【００７５】
Ｃ１D3：解析領域電磁界記憶手段（解析領域要素記憶手段）
解析領域電磁界記憶手段Ｃ１D3は、解析領域電界一時記憶手段Ｃ１D3aと、解析領域磁界一時記憶手段Ｃ１D3bとを有し、前記解析領域２の電磁界、すなわち、ベクトル成分Ｅ＝（Ｅｘ，Ｅｙ，Ｅｚ），Ｈ＝（Ｈｘ，Ｈｙ，Ｈｚ）を記憶する。
Ｃ１D3a：解析領域電界一時記憶手段
解析領域電界一時記憶手段Ｃ１D3aは、前記解析領域２の電界Ｅ＝（Ｅｘ，Ｅｙ，Ｅｚ）を一時記憶する。実施例１の前記解析領域電界一時記憶手段Ｃ１D3aは、（ｎ−１）・Δｔのタイミングでは、全てのＥｘｎ−１（ｉ＋１／２，ｊ，ｋ），Ｅｙｎ−１（ｉ，ｊ＋１／２，ｋ），Ｅｚｎ−１（ｉ，ｊ，ｋ＋１／２）を一時記憶し、次のｎ・Δｔのタイミングでは、前記解析領域電界演算手段Ｃ１D1で計算された全てのＥｘｎ（ｉ＋１／２，ｊ，ｋ），Ｅｙｎ（ｉ，ｊ＋１／２，ｋ），Ｅｚｎ（ｉ，ｊ，ｋ＋１／２）に更新して一時記憶する。
Ｃ１D3b：解析領域磁界一時記憶手段
解析領域磁界一時記憶手段Ｃ１D3bは、前記解析領域２の磁界Ｈ＝（Ｈｘ，Ｈｙ，Ｈｚ）を一時記憶する。実施例１の前記解析領域磁界一時記憶手段Ｃ１D3bは、（ｎ−１／２）・Δｔのタイミングでは、全てのＨｘｎ−１／２（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙｎ−１／２（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚｎ−１／２（ｉ＋１／２，ｊ＋１／２，ｋ）を一時記憶し、次の（ｎ＋１／２）・Δｔのタイミングでは、前記解析領域磁界演算手段Ｃ１D2で計算された全てのＨｘｎ＋１／２（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙｎ＋１／２（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚｎ＋１／２（ｉ＋１／２，ｊ＋１／２，ｋ）に更新して一時記憶する。
【００７６】
Ｃ１Ｅ：吸収境界解析手段
吸収境界解析手段Ｃ１Ｅは、吸収境界電界演算手段Ｃ１E1と、吸収境界磁界演算手段Ｃ１E2と、吸収境界電磁界記憶手段Ｃ１E3とを有し、前記入力波源１から発生した前記電磁波を前記吸収境界５で解析する吸収境界解析処理を実行する。実施例１の前記吸収境界解析手段Ｃ１Ｅは、後述するＵＰＰ法に基づいて、（ｎ−１）・Δｔのタイミングで計算されて既知のベクトル量Ｅｎ−１＝（Ｅｘｎ−１，Ｅｙｎ−１，Ｅｚｎ−１），Ｈｎ−１／２＝（Ｈｘｎ−１／２，Ｈｙｎ−１／２，Ｈｚｎ−１／２）から、次のｎ・Δｔのタイミングのベクトル量Ｅｎ＝（Ｅｘｎ，Ｅｙｎ，Ｅｚｎ），Ｈｎ＋１／２＝（Ｈｘｎ＋１／２，Ｈｙｎ＋１／２，Ｈｚｎ＋１／２）を計算する前記吸収境界解析処理を実行する。
【００７７】
（ＵＰＰ（uniaxial pseudo propagation、一軸擬伝搬）法の概念）
ここで、本願発明者が提案する、前記ＳＦ法に替わる新たな伝搬法を、以降、ＵＰＰ（uniaxial pseudo propagation、一軸擬伝搬）法と呼ぶこととする。
前記ＵＰＰ法は、従来公知のＣＩＰ（Cubic-Interpolated Pseudo particle）法を応用した新たな伝搬法である。ＣＩＰ法とは、電磁波等の解析とは異なる分野で、主に成層流体の移流を解析する数値シミュレーションの分野において研究されてきた数値解析の手法であり、近年では、乱流、気泡流、多孔流、希薄流体、粘弾性流体の分野から、航空、船舶、自動車、化学反応、結晶化、天文、原子力、半導体、プラズマ、電力等の幅広い分野で利用されている。なお、前記ＣＩＰ法については、特開２００２−３５１８５７公報等に記載されており、公知である。
【００７８】
図５はＣＩＰ法におけるdirectional splittingの説明図である。
次に、前記ＣＩＰ法の多次元化処理に利用されるdirectional splittingについて簡単に説明する。例えば、図５に示すように、ある入力波ｆについて、現在のタイミング（ｎ−１）・Δｔの前記入力波ｆの位置をｆｎ−１とし、次のタイミングｎ・Δｔの前記入力波ｆの位置をｆｎとした場合に、ｆｎ−１からｆｎをＸＹ平面上で計算する場合を考える。
このとき、前記directional splittingでは、前記入力波ｆがｆｎ−１からｆｎに進行したときの変化量をＸ軸方向の成分とＹ軸方向の成分とに分離する。すなわち、入力波ｆのベクトルのＸ軸方向の成分と、Ｙ軸方向の成分とが独立で計算できるように分離する。
【００７９】
この場合、仮想的なタイミングをｎ′・Δｔとし、前記仮想的なタイミングｎ′・Δｔの前記入力波ｆの位置をｆｎ′とした場合、例えば、ｆｎ′をｆｎ−１がｆｎに進行したときの変化量のＸ軸方向の成分だけ変化した位置とすれば、ｆｎ−１からｆｎ′を計算し、計算されたｆｎ′から前記変化量のＹ軸方向の成分だけ変化した位置を計算して、ｆｎを計算することができる。よって、前記directional splittingでは、入力波ｆが現在のタイミング（ｎ−１）・Δｔの位置ｆｎ−１から、次のタイミングｎ・Δｔの位置ｆｎを計算する場合、以下の式（１７），（１８）を計算する。
【００８０】
（∂／∂ｔ）ｆ＝−ｕ（∂／∂ｘ）ｆ（ｆｎ−１からｆｎ′に関して）      …（１８）
（∂／∂ｔ）ｆ＝−ｖ（∂／∂ｙ）ｆ（ｆｎ′からｆｎに関して）         …（１９）
ここで、ｕは入力波ｆのＸ軸方向の速度、ｖは入力波ｆのＹ軸方向の速度である。
また、前記directional splittingでは、ｆｎ′をｆｎ−１がｆｎに進行したときの変化量のＹ軸方向の成分だけ変化した位置とすることも可能である。すなわち、以下の式（１８′），（１９′）を計算した結果は、前記式（１８），（１９）を計算した結果と同等である。
（∂／∂ｔ）ｆ＝−ｖ（∂／∂ｙ）ｆ（ｆｎ−１からｆｎ′に関して）    …（１８′）
（∂／∂ｔ）ｆ＝−ｕ（∂／∂ｘ）ｆ（ｆｎ′からｆｎに関して）       …（１９′）
すなわち、前記ＣＩＰ法のdirectional splittingでは、ｆｎ−１からｆｎを計算する場合、前記式（１８），（１９）または式（１８′），（１９′）に基づいて、ｆｎ−１のＸ軸方向（またはＹ軸方向）の偏微分から仮想的な位置ｆｎ′を計算し、ｆｎ′のＹ軸方向（Ｘ軸方向）の偏微分から仮想的な位置ｆｎを計算する。
【００８１】
（ＵＰＰ法の説明）
ＵＰＰ法は、前記ＣＩＰ法のdirectional splittingの概念を利用し、前記マクスウェル回転方程式系（式（１），（２）参照）に適用した新たな伝搬法である。なお、前記吸収境界５における前記入力波の伝搬の様子を解析するために前記ＵＰＰ法を使用する吸収境界条件を、以降、ＵＰＰ−ＰＭＬと呼ぶこととする。
なお、前記ＵＰＰ−ＰＭＬでは、前記吸収境界５において、前記Ｂ−ＰＭＬと同等の導電率および導磁率を有するものとする。すなわち、前記ＵＰＰ−ＰＭＬでは、前記導電率σｘｘ，σｘｙ，σｘｚ，σｙｘ，σｙｙ，σｙｚ，σｚｘ，σｚｙ，σｚｚおよび前記導磁率σｘｘ＊，σｘｙ＊，σｘｚ＊，σｙｘ＊，σｙｙ＊，σｙｚ＊，σｚｘ＊，σｚｙ＊，σｚｚ＊について、前記Ｂ−ＰＭＬと同様に、前記インピーダンスマッチング条件（式（８−１）〜（８−３）参照）および前記電磁波が前記吸収境界０５に斜めに入射した場合でも反射せずに吸収するための仮想的な条件（式（９−１），（９−２）参照）を満たす。
【００８２】
（マクスウェル回転方程式系のオペレータ分割）
この場合、前記ＵＰＰ−ＰＭＬでは、前記マクスウェル回転方程式系（式（１），（２）参照）を、以下の式（２０）〜（２５）で表現できる。
（∂／∂ｔ）Ｅｘ
＝（∂／∂ｙ）Ｈｚ−σｙｙＥｘ−（∂／∂ｚ）Ｈｙ−σｚｚＥｘ        …（２０）
（∂／∂ｔ）Ｅｙ
＝（∂／∂ｚ）Ｈｘ−σｚｚＥｙ−（∂／∂ｘ）Ｈｚ−σｘｘＥｙ        …（２１）
（∂／∂ｔ）Ｅｚ
＝（∂／∂ｘ）Ｈｙ−σｘｘＥｚ−（∂／∂ｙ）Ｈｘ−σｙｙＥｚ        …（２２）
（∂／∂ｔ）Ｈｘ
＝−（∂／∂ｙ）Ｅｚ−σｙｙ＊Ｈｘ＋（∂／∂ｚ）Ｅｙ−σｚｚ＊Ｈｘ  …（２３）
（∂／∂ｔ）Ｈｙ
＝−（∂／∂ｚ）Ｅｘ−σｚｚ＊Ｈｙ＋（∂／∂ｘ）Ｅｚ−σｘｘ＊Ｈｙ  …（２４）
（∂／∂ｔ）Ｈｚ
＝−（∂／∂ｘ）Ｅｙ−σｘｘ＊Ｈｚ＋（∂／∂ｙ）Ｅｘ−σｙｙ＊Ｈｚ  …（２５）
【００８３】
ここで、ＵＰＰ法で、（ｎ−１）・Δｔのタイミングで計算されて既知の電界Ｅｎ−１のＺ軸方向のベクトル成分Ｅｚｎ−１から、次のｎ・Δｔのタイミングの電界ＥｎのＺ軸方向のベクトル成分Ｅｚｎを計算する場合について説明する。
この場合、ＵＰＰ法は、前記directional splittingの概念を利用して、前記式（２２）の右辺について、Ｘ軸方向の成分と、Ｙ軸方向の成分とを分離し、以下の式（２２−１），（２２−２）で表現する。
（∂／∂ｔ）Ｅｚ（ｘ成分）
＝（∂／∂ｘ）Ｈｙ−σｘｘＥｚ                                    …（２２−１）
（∂／∂ｔ）Ｅｚ（ｙ成分）
＝−（∂／∂ｙ）Ｈｘ−σｙｙＥｚ                                  …（２２−２）
【００８４】
（ＵＰＰ法の適用）
ここで、前記仮想的なタイミングｎ′・Δｔの電界をＥｎ′、前記仮想的なタイミングｎ′・Δｔの電界Ｅｎ′のＺ軸方向のベクトル成分をＥｚｎ′とする。前記ＵＰＰ法では、例えば、前記仮想的なタイミングｎ′・Δｔの電界Ｅｎ′は、前記ベクトル成分Ｅｚｎ′のみ変化するものと仮定し、前記次のタイミングｎ・Δｔの電界Ｅｎは、前記ベクトル成分Ｅｚｎのみ変化するものと仮定する。この場合、前記式（２２−１），（２２−２）が、以下の式（２２−１′）（２２−２′）で表現できる。
（∂／∂ｔ）Ｅｚ｜ｎｘ
＝（（∂／∂ｘ）Ｈｙｎ−１／２−σｘｘＥｚｎｘ）               …（２２−１′）
（∂／∂ｔ）Ｅｚ｜ｎｙ
＝（（∂／∂ｙ）Ｈｘｎ−１／２−σｘｘＥｚｎｙ）               …（２２−２′）
【００８５】
ここで、ｎｘは（ｎ−１）・Δｔからｎ′・Δｔまでの仮想的な時間、ｎｙはｎ′・Δｔからｎ・Δｔまでの仮想的な時間である。
また、ＵＰＰ法では、前記式（２２−１′），（２２−２′）について、左辺の（∂／∂ｔ）Ｅｚを時間差分、右辺のＥｚｎｘ，Ｅｚｎｙを中心差分した、中心差分近似の式（２２−１″），（２２−２″）に変換する。
（Ｅｚｎ′−Ｅｚｎ−１）／Δｔ
＝（∂／∂ｘ）Ｈｙｎ−１／２−σｘｘ（Ｅｚｎ′＋Ｅｚｎ−１）／２…（２２−１″）
（Ｅｚｎ−Ｅｚｎ′）／Δｔ
＝−（∂／∂ｙ）Ｈｘｎ−１／２−σｙｙ（Ｅｚｎ＋Ｅｚｎ′）／２…（２２−２″）
【００８６】
また、前記式（２２−１″），（２２−２″）の右辺の（∂／∂ｘ）Ｈｙｎ−１／２，（∂／∂ｙ）Ｈｘｎ−１／２は、前記式（５″），（６″）によって、前記単位長Δｘ，Δｙの式に置換できる。
したがって、前記式（２２−１″）に基づいて、Δｔ，σｘｘ，Ｅｚｎ−１，Ｈｙｎ−１／２からＥｚｎ′が計算でき、前記式（２２−２″）に基づいて、Δｔ，σｙｙ，Ｅｚｎ′，Ｈｙｎ−１／２からＥｚｎが計算できる。すなわち、Ｅｚｎ−１から前記式（２２−１″）よりＥｚｎ′を計算し、Ｅｚｎ′から前記式（２２−２″）よりＥｚｎを計算できる。
【００８７】
また、Ｈｚｎ＋１／２についても同様に、仮想のタイミングを（ｎ＋１／２）′・Δｔとした場合に、前記式（２５−１），（２５−２）が、以下の式（２５−１″）（２５−２″）で表現できる。
（Ｈｚ（ｎ＋１／２）′−Ｈｚｎ−１／２）／Δｔ
＝−（∂／∂ｘ）Ｅｙｎ−σｘｘ＊（Ｈｚ（ｎ＋１／２）′＋Ｈｚｎ−１／２）／２
…（２５−１″）
（Ｈｚｎ＋１／２−Ｈｚ（ｎ＋１／２）′）／Δｔ
＝（∂／∂ｙ）Ｅｘｎ−σｙｙ＊（Ｈｚｎ＋１／２＋Ｈｚ（ｎ＋１／２）′）／２
…（２５−２″）
また、前記式（２５−１″），（２５−２″）の右辺の（∂／∂ｘ）Ｅｙｎ，（∂／∂ｙ）Ｅｘｎは、前記式（５′），（６′）によって、前記単位長Δｘ，Δｙの式に置換できる。
【００８８】
そして、その他のベクトル成分Ｅｘｎ，Ｅｙｎ，Ｈｘｎ＋１／２，Ｈｙｎ＋１／２も、前記ベクトル成分Ｅｚｎ，Ｈｚｎ＋１／２と同様に計算できる。
したがって、前記ＵＰＰ法は、前記式（２０）〜（２５）に基づいて、前記吸収境界５における前記入力波の伝搬の様子を解析することができる。すなわち、（ｎ−１）・Δｔのタイミングで計算されて既知のベクトル量Ｅｎ−１＝（Ｅｘｎ−１，Ｅｙｎ−１，Ｅｚｎ−１），Ｈｎ−１／２＝（Ｈｘｎ−１／２，Ｈｙｎ−１／２，Ｈｚｎ−１／２）から、ｎ′・Δｔのベクトル量を経て、次のｎ・Δｔのタイミングのベクトル量Ｅｎ＝（Ｅｘｎ，Ｅｙｎ，Ｅｚｎ），Ｈｎ＋１／２＝（Ｈｘｎ＋１／２，Ｈｙｎ＋１／２，Ｈｚｎ＋１／２）を計算することができる。
この結果、前記吸収境界５における各座標位置のベクトル量Ｅ，Ｈの時間的な変化、すなわち、前記入力波源１から発生した前記入力波の伝搬の様子は、前記ＵＰＰ法によって計算できる。
【００８９】
Ｃ１E1：吸収境界電界演算手段
吸収境界電界演算手段Ｃ１E1は、仮想時間電界演算手段Ｃ１E1aと、解析時間後電界演算手段Ｃ１E1bとを有し、前記吸収境界５において、（ｎ−１）・Δｔのタイミングで計算された電界Ｅｎ−１＝（Ｅｘｎ−１，Ｅｙｎ−１，Ｅｚｎ−１）等から、ｎ′・Δｔのタイミングの電界Ｅｎ′＝（Ｅｘｎ′，Ｅｙｎ′，Ｅｚｎ′）を計算し、ｎ′・Δｔのタイミングの電界Ｅｎ′＝（Ｅｘｎ′，Ｅｙｎ′，Ｅｚｎ′）等から、次のｎ・Δｔのタイミングの電界Ｅｎ＝（Ｅｘｎ，Ｅｙｎ，Ｅｚｎ）を計算する吸収境界電界計算処理を実行する。実施例１の前記吸収境界電界演算手段Ｃ１E1は、前記式（２０）〜（２２）に基づいて、Ｅｘｎ（ｉ＋１／２，ｊ，ｋ），Ｅｙｎ（ｉ，ｊ＋１／２，ｋ），Ｅｚｎ（ｉ，ｊ，ｋ＋１／２）を計算する（式（２２−１″），（２２−２″），（５″），（６″）参照）。
Ｃ１E1a：仮想時間電界演算手段（仮想時間要素演算手段）
仮想時間電界演算手段Ｃ１E1aは、前記吸収境界５において、（ｎ−１）・Δｔのタイミングで計算された電界Ｅｎ−１＝（Ｅｘｎ−１，Ｅｙｎ−１，Ｅｚｎ−１）等から、ｎ′・Δｔのタイミングの電界Ｅｎ′＝（Ｅｘｎ′，Ｅｙｎ′，Ｅｚｎ′）を計算する（式（２２−１″），（５″）参照）。
Ｃ１E1b：解析時間後電界演算手段（解析時間後演算記憶手段）
解析時間後電界演算手段Ｃ１E1bは、前記吸収境界５において、ｎ′・Δｔのタイミングの電界Ｅｎ′＝（Ｅｘｎ′，Ｅｙｎ′，Ｅｚｎ′）等から、次のｎ・Δｔのタイミングの電界Ｅｎ＝（Ｅｘｎ，Ｅｙｎ，Ｅｚｎ）を計算する（式（２２−２″），（６″）参照）。
【００９０】
Ｃ１E2：吸収境界磁界演算手段
吸収境界磁界演算手段Ｃ１E2は、仮想時間磁界演算手段Ｃ１E2aと、解析時間後磁界演算手段Ｃ１E2bとを有し、前記吸収境界５において、（ｎ−１／２）・Δｔのタイミングで計算された磁界Ｈｎ−１／２＝（Ｈｘｎ−１／２，Ｈｙｎ−１／２，Ｈｚｎ−１／２）等から、（ｎ＋１／２）′・Δｔのタイミングの磁界Ｈ（ｎ＋１／２）′＝（Ｈｘ（ｎ＋１／２）′，Ｈｙ（ｎ＋１／２）′，Ｈｚ（ｎ＋１／２）′）を計算し、（ｎ＋１／２）′・Δｔのタイミングの磁界Ｈ（ｎ＋１／２）′＝（Ｈｘ（ｎ＋１／２）′，Ｈｙ（ｎ＋１／２）′，Ｈｚ（ｎ＋１／２）′）等から、次の（ｎ＋１／２）・Δｔのタイミングの磁界Ｈｎ＋１／２＝（Ｈｘｎ＋１／２，Ｈｙｎ＋１／２，Ｈｚｎ＋１／２）を計算する吸収境界磁界計算処理を実行する。実施例１の前記吸収境界磁界演算手段Ｃ１E2は、前記式（２３）〜（２５）に基づいて、Ｈｘｎ＋１／２（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙｎ＋１／２（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚｎ＋１／２（ｉ＋１／２，ｊ＋１／２，ｋ）を計算する（式（２５−１″），（２５−２″），（５′），（６′）参照）。
Ｃ１E2a：仮想時間磁界演算手段（仮想時間要素演算手段）
仮想時間磁界演算手段Ｃ１E2aは、前記吸収境界５において、（ｎ−１／２）・Δｔのタイミングで計算された磁界Ｈｎ−１／２＝（Ｈｘｎ−１／２，Ｈｙｎ−１／２，Ｈｚｎ−１／２）等から、（ｎ＋１／２）′・Δｔのタイミングの磁界Ｈ（ｎ＋１／２）′＝（Ｈｘ（ｎ＋１／２）′，Ｈｙ（ｎ＋１／２）′，Ｈｚ（ｎ＋１／２）′）を計算する（式（２５−１″），（５′）参照）。
Ｃ１E2b：解析時間後磁界演算手段（解析時間後要素演算手段）
解析時間後磁界演算手段Ｃ１E2bは、（ｎ＋１／２）′・Δｔのタイミングの磁界Ｈ（ｎ＋１／２）′＝（Ｈｘ（ｎ＋１／２）′，Ｈｙ（ｎ＋１／２）′，Ｈｚ（ｎ＋１／２）′）等から、次の（ｎ＋１／２）・Δｔのタイミングの磁界Ｈｎ＋１／２＝（Ｈｘｎ＋１／２，Ｈｙｎ＋１／２，Ｈｚｎ＋１／２）を計算する（式（２５−２″），（６′）参照）。
【００９１】
Ｃ１E3：吸収境界電磁界記憶手段（吸収境界要素記憶手段）
吸収境界電磁界記憶手段Ｃ１E3は、吸収境界電界一時記憶手段Ｃ１E3aと、吸収境界磁界一時記憶手段Ｃ１E3bとを有し、前記吸収境界５の電磁界、すなわち、ベクトル成分Ｅ＝（Ｅｘ，Ｅｙ，Ｅｚ），Ｈ＝（Ｈｘ，Ｈｙ，Ｈｚ）を記憶する。
Ｃ１E3a：吸収境界電界一時記憶手段
吸収境界電界一時記憶手段Ｃ１E3aは、前記吸収境界５の電界Ｅ＝（Ｅｘ，Ｅｙ，Ｅｚ）を一時記憶する。実施例１の前記吸収境界電界一時記憶手段Ｃ１E3aは、（ｎ−１）・Δｔのタイミングでは、全てのＥｘｎ−１（ｉ＋１／２，ｊ，ｋ），Ｅｙｎ−１（ｉ，ｊ＋１／２，ｋ），Ｅｚｎ−１（ｉ，ｊ，ｋ＋１／２）を一時記憶し、ｎ′・Δｔのタイミングでは、全てのＥｘｎ′（ｉ＋１／２，ｊ，ｋ），Ｅｙｎ′（ｉ，ｊ＋１／２，ｋ），Ｅｚｎ′（ｉ，ｊ，ｋ＋１／２）に更新して一時記憶し、次のｎ・Δｔのタイミングでは、前記吸収境界電界演算手段Ｃ１E1で計算された全てのＥｘｎ（ｉ＋１／２，ｊ，ｋ），Ｅｙｎ（ｉ，ｊ＋１／２，ｋ），Ｅｚｎ（ｉ，ｊ，ｋ＋１／２）に更新して一時記憶する。
Ｃ１E3b：吸収境界磁界一時記憶手段
吸収境界磁界一時記憶手段Ｃ１E3bは、前記吸収境界５の磁界Ｈ＝（Ｈｘ，Ｈｙ，Ｈｚ）を一時記憶する。実施例１の前記吸収境界磁界一時記憶手段Ｃ１E3bは、（ｎ−１／２）・Δｔのタイミングでは、全てのＨｘｎ−１／２（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙｎ−１／２（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚｎ−１／２（ｉ＋１／２，ｊ＋１／２，ｋ）を一時記憶し、（ｎ＋１／２）′・Δｔのタイミングでは、全てのＨｘ（ｎ＋１／２）′（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｎ＋１／２）′（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｎ＋１／２）′（ｉ＋１／２，ｊ＋１／２，ｋ）に更新して一時記憶し、次の（ｎ＋１／２）・Δｔのタイミングでは、前記吸収境界磁界演算手段Ｃ１E2で計算された全てのＨｘｎ＋１／２（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙｎ＋１／２（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚｎ＋１／２（ｉ＋１／２，ｊ＋１／２，ｋ）に更新して一時記憶する。
【００９２】
Ｃ１Ｆ：解析結果出力手段
解析結果出力手段Ｃ１Ｆは、前記解析領域解析手段Ｃ１Ｄおよび前記吸収境界解析手段Ｃ１Ｅによる電磁界解析処理の解析結果を出力する。実施例１の前記解析結果出力手段Ｃ１Ｆは、前記解析結果を解析結果ファイルに出力する。具体的には、前記解析結果出力手段Ｃ１Ｆは、前記単位時間Δｔ毎に、前記解析領域電界一時記憶手段Ｃ１D3a、前記解析領域磁界一時記憶手段Ｃ１D3b、前記吸収境界電界一時記憶手段Ｃ１E3a、前記吸収境界磁界一時記憶手段Ｃ１E3bに一時記憶された全てのＥｘｎ（ｉ＋１／２，ｊ，ｋ），Ｅｙｎ（ｉ，ｊ＋１／２，ｋ），Ｅｚｎ（ｉ，ｊ，ｋ＋１／２），Ｈｘｎ＋１／２（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙｎ＋１／２（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚｎ＋１／２（ｉ＋１／２，ｊ＋１／２，ｋ）を、前記解析結果ファイルに一括出力する。すなわち、実施例１の前記解析結果ファイルには、ｎ・Δｔ（ｎ＝１，２，…，Ｎｍａｘ）毎に更新された全てＥｘｎ（ｉ＋１／２，ｊ，ｋ），Ｅｙｎ（ｉ，ｊ＋１／２，ｋ），Ｅｚｎ（ｉ，ｊ，ｋ＋１／２），Ｈｘｎ＋１／２（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙｎ＋１／２（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚｎ＋１／２（ｉ＋１／２，ｊ＋１／２，ｋ）が一括出力されている。
【００９３】
Ｃ１Ｇ：解析終了判別手段
解析終了判別手段Ｃ１Ｇは、前記電磁界解析処理を終了するか否かを判別する。実施例１の前記解析終了判別手段Ｃ１Ｇは、演算回数ｎ（ｎ＝１，２，…，Ｎｍａｘ）が、ユーザにより入力された前記解析回数の最大値Ｎｍａｘまで計数されたか否かを判別することにより、前記電磁界解析処理を終了するか否かを判別する。
図６は実施例１の解析結果表示画像の説明図であり、解析結果を動画で表示したときの一例を示す図である。
Ｃ２：解析電磁界表示手段（解析場表示手段）
解析電磁界表示手段Ｃ２は、解析結果動画表示手段Ｃ２Ａを有し、前記電磁界解析処理により出力された解析結果、すなわち、前記解析された電磁界を表示する。実施例１の解析電磁界表示手段Ｃ２は、図６に示すように、解析結果表示画像３１をディスプレイＨ２に表示し、前記解析結果出力手段Ｃ１Ｅにより出力された解析結果ファイルに基づいて、前記解析領域解析手段Ｃ１Ｄおよび前記吸収境界解析手段Ｃ１Ｅによる電磁界解析処理の解析結果を前記解析結果表示画像３１に表示する。
【００９４】
Ｃ２Ａ：解析結果動画表示手段
解析結果動画表示手段Ｃ２Ａは、前記解析結果出力手段Ｃ１Ｅにより出力された解析結果ファイルに基づいて、前記解析領域解析手段Ｃ１Ｄおよび前記吸収境界解析手段Ｃ１Ｅによる電磁界解析処理の解析結果を動画を作成し、作成された前記解析結果の動画を前記解析結果表示画像３１に表示する（図６参照）。
実施例１では、前記解析結果の動画は、ｎ・Δｔ（ｎ＝１，２，…，Ｎｍａｘ）の各タイミングで計算された前記解析領域２の全て座標（ｉ，ｊ，ｋ）の電界Ｅｎ（ｉ，ｊ，ｋ）＝（Ｅｘｎ（ｉ＋１／２，ｊ，ｋ），Ｅｙｎ（ｉ，ｊ＋１／２，ｋ），Ｅｚｎ（ｉ，ｊ，ｋ＋１／２））に基づいて、等電位線を描いた画像を作成し、各タイミングの等電位線の画像を連続表示することにより作成される。なお、実施例１では、前記解析結果の動画を前記解析領域２の電界Ｅに基づいて作成するが、これに限定されず、例えば、前記解析領域２の全て座標（ｉ，ｊ，ｋ）の磁界Ｈｎ＋１／２（ｉ，ｊ，ｋ）＝（Ｈｘｎ＋１／２（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙｎ＋１／２（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚｎ＋１／２（ｉ＋１／２，ｊ＋１／２，ｋ））に基づいて、等磁位線の画像を連続表示した動画を作成し、前記解析結果の動画とすることも可能である。
【００９５】
（実施例１のフローチャートの説明）
次に、実施例１のクライアントパソコンＰＣの処理の流れをフローチャートを使用して説明する。なお、前記電磁界解析表示プログラムＡＰ１の電磁気解析結果表示手段Ｃ２の処理については、解析結果ファイルに基づいて、解析結果を動画に変換して表示するだけであるため、フローチャートによる図示及び詳細な説明を省略する。
【００９６】
（電磁界解析処理のフローチャートの説明）
図７は本発明の実施例１の電磁界解析装置の電磁界解析処理のフローチャートである。
図７のフローチャートの各ＳＴ（ステップ）の処理は、クライアントパソコンＰＣの電磁界解析表示プログラムＡＰ１に従って行われる。また、この処理は、クライアントパソコンＰＣの他の各種処理と並行して実行される。
図７に示すフローチャートは電磁界解析表示プログラムＡＰ１の起動により開始される。
【００９７】
図７のＳＴ１において、ユーザにより前記解析条件入力画像１１の前記各入力欄１２〜２８に適切な値が入力され、且つ、前記解析開始画像２９が選択されたか否かを判別することにより、前記電磁界解析処理を開始するか否かを判別する（図４参照）。イエス（Ｙ）の場合はＳＴ２に移り、ノー（Ｎ）の場合はＳＴ１を繰り返す。
ＳＴ２において、以下の（１）〜（３）の処理を実行し、ＳＴ３に移る。
（１）入力された解析条件パラメタに基づいて、入力波源１の電界ＥのＺ軸方向のベクトル成分Ｅｚ（ｉＧ，ｊＧ，ｋＧ＋１／２）の振幅Ｅ０および角周波数ω０、セル３，６の各辺の長さΔｘ，Δｙ，Δｚ、解析領域２の各辺の長さＩｍａｘ，Ｊｍａｘ，Ｋｍａｘ、吸収境界５の層の値ｄ、解析領域２の誘電率εおよび透磁率μ、吸収境界５の導電率σｘｘ，σｙｙ，σｚｚおよび導磁率σｘｘ＊，σｙｙ＊，σｚｚ＊、単位時間Δｔ、解析回数ｎの最大値Ｎｍａｘをセットする。
（２）解析領域２および吸収境界５の全ての座標（ｉ，ｊ，ｋ）において、Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２），Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ）をリセットする。すなわち、
Ｅｘ（ｉ＋１／２，ｊ，ｋ）＝Ｅｘ０（ｉ＋１／２，ｊ，ｋ）＝０，
Ｅｙ（ｉ，ｊ＋１／２，ｋ）＝Ｅｙ０（ｉ，ｊ＋１／２，ｋ）＝０，
Ｅｚ（ｉ，ｊ，ｋ＋１／２）＝Ｅｚ０（ｉ，ｊ，ｋ＋１／２）＝０，
Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２）＝Ｈｘ１／２（ｉ，ｊ＋１／２，ｋ＋１／２）＝０，
Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２）＝Ｈｙ１／２（ｉ＋１／２，ｊ，ｋ＋１／２）＝０，
Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ）＝Ｈｚ１／２（ｉ＋１／２，ｊ＋１／２，ｋ）＝０
とする。
（３）解析回数ｎに１をセットする（ｎ＝１）。
【００９８】
ＳＴ３において、以下の（１），（２）の処理を実行し、ＳＴ４に移る。
（１）入力波源１の座標（ｉＧ，ｊＧ，ｋＧ）の電界Ｅ（ｉＧ，ｊＧ，ｋＧ）について（Ｅ＝（０，０，Ｅｚ））、（ｎ−１）・ΔｔのタイミングのＥｚｎ−１（ｉＧ，ｊＧ，ｋＧ＋１／２）を、式（１７″）に基づいて計算する。すなわち、
Ｅｚｎ−１（ｉＧ，ｊＧ，ｋＧ＋１／２）＝Ｅ０×ｓｉｎ（ω０（ｎ−１）・Δｔ）
を計算する。
（２）Ｅｚ（ｉＧ，ｊＧ，ｋＧ＋１／２）を、計算されたＥｚｎ−１（ｉＧ，ｊＧ，ｋＧ＋１／２）に更新して一時記憶する。すなわち、
Ｅｚ（ｉＧ，ｊＧ，ｋＧ＋１／２）＝Ｅｚｎ−１（ｉＧ，ｊＧ，ｋＧ＋１／２）
とする。
【００９９】
ＳＴ４において、解析領域２の全ての座標（ｉ，ｊ，ｋ）の電界Ｅ（ｉ，ｊ，ｋ）について、（ｎ−１）・Δｔのタイミングで計算された電界Ｅｎ−１＝（Ｅｘｎ−１，Ｅｙｎ−１，Ｅｚｎ−１）等から、次のｎ・Δｔのタイミングの電界Ｅｎ＝（Ｅｘｎ，Ｅｙｎ，Ｅｚｎ）を計算する解析領域電界計算処理を実行する。なお、解析領域電界計算処理については、従来公知のＦＤＴＤ法における電界Ｅｎ＝（Ｅｘｎ，Ｅｙｎ，Ｅｚｎ）の計算をする処理と同様であるため、詳細な説明を省略する（式（４）参照）。そして、ＳＴ５に移る。
ＳＴ５において、吸収境界５の全ての座標（ｉ，ｊ，ｋ）の電界Ｅ（ｉ，ｊ，ｋ）について、（ｎ−１）・Δｔのタイミングで計算された電界Ｅｎ−１＝（Ｅｘｎ−１，Ｅｙｎ−１，Ｅｚｎ−１）等から、ｎ′・Δｔのタイミングの電界Ｅｎ′＝（Ｅｘｎ′，Ｅｙｎ′，Ｅｚｎ′）を計算し、ｎ′・Δｔのタイミングの電界Ｅｎ′＝（Ｅｘｎ′，Ｅｙｎ′，Ｅｚｎ′）等から、次のｎ・Δｔのタイミングの電界Ｅｎ＝（Ｅｘｎ，Ｅｙｎ，Ｅｚｎ）を計算する吸収境界電界計算処理（後述する図８のサブルーチン参照）を実行する。そして、ＳＴ６に移る。
【０１００】
ＳＴ６において、解析領域２の全ての座標（ｉ，ｊ，ｋ）の磁界Ｈ（ｉ，ｊ，ｋ）について、（ｎ−１／２）・Δｔのタイミングで計算された磁界Ｈｎ−１／２＝（Ｈｘｎ−１／２，Ｈｙｎ−１／２，Ｈｚｎ−１／２）等から、次の（ｎ＋１／２）・Δｔのタイミングの磁界Ｈｎ＋１／２＝（Ｈｘｎ＋１／２，Ｈｙｎ＋１／２，Ｈｚｎ＋１／２）を計算する解析領域磁界計算処理を実行する。なお、解析領域磁界計算処理については、従来公知のＦＤＴＤ法における磁界Ｈｎ＋１／２＝（Ｈｘｎ＋１／２，Ｈｙｎ＋１／２，Ｈｚｎ＋１／２）の計算をする処理と同様であるため、詳細な説明を省略する（式（３）参照）。そして、ＳＴ７に移る。
ＳＴ７において、吸収境界５の全ての座標（ｉ，ｊ，ｋ）の磁界Ｈ（ｉ，ｊ，ｋ）について、（ｎ−１／２）・Δｔのタイミングで計算された磁界Ｈｎ−１／２＝（Ｈｘｎ−１／２，Ｈｙｎ−１／２，Ｈｚｎ−１／２）等から、（ｎ＋１／２）′・Δｔのタイミングの磁界Ｈ（ｎ＋１／２）′＝（Ｈｘ（ｎ＋１／２）′，Ｈｙ（ｎ＋１／２）′，Ｈｚ（ｎ＋１／２）′）を計算し、（ｎ＋１／２）′・Δｔのタイミングの磁界Ｈ（ｎ＋１／２）′＝（Ｈｘ（ｎ＋１／２）′，Ｈｙ（ｎ＋１／２）′，Ｈｚ（ｎ＋１／２）′）等から、次の（ｎ＋１／２）・Δｔのタイミングの磁界Ｈｎ＋１／２＝（Ｈｘｎ＋１／２，Ｈｙｎ＋１／２，Ｈｚｎ＋１／２）を計算する吸収境界磁界計算処理（後述する図９のサブルーチン参照）を実行する。そして、ＳＴ８に移る。
【０１０１】
ＳＴ８において、解析領域２および吸収境界５の全ての座標（ｉ，ｊ，ｋ）において、各一時記憶手段Ｃ１D3，Ｃ１D4，Ｃ１E3，Ｃ１E4でＥｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２），Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ）として一時記憶されたＥｘｎ（ｉ＋１／２，ｊ，ｋ），Ｅｙｎ（ｉ，ｊ＋１／２，ｋ），Ｅｚｎ（ｉ，ｊ，ｋ＋１／２），Ｈｘｎ＋１／２（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙｎ＋１／２（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚｎ＋１／２（ｉ＋１／２，ｊ＋１／２，ｋ）を解析結果ファイルに一括出力する。そして、ＳＴ９に移る。
ＳＴ９において、解析回数ｎに＋１を加算する（ｎ＝ｎ＋１）。そして、ＳＴ１０に移る。
ＳＴ１０において、解析回数ｎが最大値Ｎｍａｘより大きくなったか否かを判別する。イエス（Ｙ）の場合はＳＴ１に戻り、ノー（Ｎ）の場合はＳＴ３に戻る。
【０１０２】
（吸収境界電界計算処理のフローチャートの説明）
図８は吸収境界電界計算処理のフローチャートであり、図７のＳＴ５のサブルーチンの説明図である。
図８のＳＴ１０１において、以下の（１），（２）の処理を実行し、ＳＴ１０２に移る。
（１）吸収境界５の各座標（ｉ，ｊ，ｋ）において、吸収境界電界一時記憶手段Ｃ１E3aでＥｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２）として一時記憶されたＥｘｎ−１（ｉ＋１／２，ｊ，ｋ），Ｅｙｎ−１（ｉ，ｊ＋１／２，ｋ），Ｅｚｎ−１（ｉ，ｊ，ｋ＋１／２）と、Δｘ〜Δｚ，Δｔ，σｘｘ〜σｚｚ，Ｈｘｎ−１／２〜Ｈｚｎ−１／２とから、Ｅｘｎ′（ｉ＋１／２，ｊ，ｋ），Ｅｙｎ′（ｉ，ｊ＋１／２，ｋ），Ｅｚｎ′（ｉ，ｊ，ｋ＋１／２）を計算する。
例えば、Ｅｚｎ′（ｉ，ｊ，ｋ＋１／２）は、式（２２−１″），（５″）から得られる以下の式（２６）をＥｚｎ′（ｉ，ｊ，ｋ＋１／２）について解くことにより計算できる。
｛Ｅｚｎ′（ｉ，ｊ，ｋ＋１／２）−Ｅｚｎ−１（ｉ，ｊ，ｋ＋１／２）｝／Δｔ
＝｛Ｈｙｎ−１／２（ｉ＋１／２，ｊ，ｋ＋１／２）
−Ｈｙｎ−１／２（ｉ−１／２，ｊ，ｋ＋１／２）｝／Δｘ
−σｘｘ｛Ｅｚｎ′（ｉ，ｊ，ｋ＋１／２）＋Ｅｚｎ−１（ｉ，ｊ，ｋ＋１／２）｝／２
…（２６）
また、Ｅｘｎ′（ｉ＋１／２，ｊ，ｋ），Ｅｙｎ′（ｉ，ｊ＋１／２，ｋ）もＥｚｎ′（ｉ，ｊ，ｋ＋１／２）と同様に計算できる。
（２）吸収境界電界一時記憶手段Ｃ１E3aでＥｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２）として一時記憶されているＥｘｎ−１（ｉ＋１／２，ｊ，ｋ），Ｅｙｎ−１（ｉ，ｊ＋１／２，ｋ），Ｅｚｎ−１（ｉ，ｊ，ｋ＋１／２）を消去して、計算されたＥｘｎ′（ｉ＋１／２，ｊ，ｋ），Ｅｙｎ′（ｉ，ｊ＋１／２，ｋ），Ｅｚｎ′（ｉ，ｊ，ｋ＋１／２）に更新して一時記憶する。すなわち、
Ｅｘ（ｉ＋１／２，ｊ，ｋ）＝Ｅｘｎ′（ｉ＋１／２，ｊ，ｋ），
Ｅｙ（ｉ，ｊ＋１／２，ｋ）＝Ｅｙｎ′（ｉ，ｊ＋１／２，ｋ），
Ｅｚ（ｉ，ｊ，ｋ＋１／２）＝Ｅｚｎ′（ｉ，ｊ，ｋ＋１／２）
とする。
【０１０３】
ＳＴ１０２において、以下の（１），（２）の処理を実行し、ＳＴ１０３に移る。
（１）吸収境界電界一時記憶手段Ｃ１E3aでＥｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２）として一時記憶されたＥｚｎ′（ｉ＋１／２，ｊ，ｋ），Ｅｙｎ′（ｉ，ｊ＋１／２，ｋ），Ｅｚｎ′（ｉ，ｊ，ｋ＋１／２）と、Δｘ〜Δｚ，Δｔ，σｘｘ〜σｚｚ，Ｈｘｎ−１／２〜Ｈｚｎ−１／２とから、Ｅｚｎ（ｉ＋１／２，ｊ，ｋ），Ｅｙｎ（ｉ，ｊ＋１／２，ｋ），Ｅｚｎ（ｉ，ｊ，ｋ＋１／２）を計算する。
例えば、Ｅｚｎ（ｉ，ｊ，ｋ＋１／２）は、式（２２−２″），（６″）から得られる以下の式（２７）をＥｚｎ（ｉ，ｊ，ｋ＋１／２）について解くことにより計算できる。
｛Ｅｚｎ（ｉ，ｊ，ｋ＋１／２）−Ｅｚｎ′（ｉ，ｊ，ｋ＋１／２）｝／Δｔ
＝｛Ｈｘｎ−１／２（ｉ，ｊ＋１／２，ｋ＋１／２）
−Ｈｘｎ−１／２（ｉ，ｊ−１／２，ｋ＋１／２）｝／Δｙ
−σｙｙ｛Ｅｚｎ（ｉ，ｊ，ｋ＋１／２）＋Ｅｚｎ′（ｉ，ｊ，ｋ＋１／２）｝／２
…（２７）
また、Ｅｘｎ（ｉ＋１／２，ｊ，ｋ），Ｅｙｎ（ｉ，ｊ＋１／２，ｋ）もＥｚｎ（ｉ，ｊ，ｋ＋１／２）と同様に計算できる。
（２）吸収境界電界一時記憶手段Ｃ１E3aでＥｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２）として一時記憶されているＥｚｎ′（ｉ＋１／２，ｊ，ｋ），Ｅｙｎ′（ｉ，ｊ＋１／２，ｋ），Ｅｚｎ′（ｉ，ｊ，ｋ＋１／２）を消去して、計算されたＥｘｎ（ｉ＋１／２，ｊ，ｋ），Ｅｙｎ（ｉ，ｊ＋１／２，ｋ），Ｅｚｎ（ｉ，ｊ，ｋ＋１／２）に更新して一時記憶する。すなわち、
Ｅｘ（ｉ＋１／２，ｊ，ｋ）＝Ｅｘｎ（ｉ＋１／２，ｊ，ｋ），
Ｅｙ（ｉ，ｊ＋１／２，ｋ）＝Ｅｙｎ（ｉ，ｊ＋１／２，ｋ），
Ｅｚ（ｉ，ｊ，ｋ＋１／２）＝Ｅｚｎ（ｉ，ｊ，ｋ＋１／２）
とする。
【０１０４】
ＳＴ１０３において、吸収境界５の全ての座標（ｉ，ｊ，ｋ）において、Ｅｘｎ（ｉ＋１／２，ｊ，ｋ），Ｅｙｎ（ｉ，ｊ＋１／２，ｋ），Ｅｚｎ（ｉ，ｊ，ｋ＋１／２）が計算されたか否かを判別する。ノー（Ｎ）の場合はＳＴ１０１に戻り、イエス（Ｙ）の場合は吸収境界電界計算処理を終了して図７の電磁界解析処理に戻る。
【０１０５】
（吸収境界磁界計算処理のフローチャートの説明）
図９は吸収境界磁界計算処理のフローチャートであり、図７のＳＴ７のサブルーチンの説明図である。
図９のＳＴ２０１において、以下の（１），（２）の処理を実行し、ＳＴ２０２に移る。
（１）吸収境界５の各座標（ｉ，ｊ，ｋ）において、吸収境界磁界一時記憶手段Ｃ１E3bでＨｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ）として一時記憶されたＨｘｎ−１／２（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙｎ−１／２（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚｎ−１／２（ｉ＋１／２，ｊ＋１／２，ｋ）と、Δｘ〜Δｚ，Δｔ，σｘｘ*〜σｚｚ*，Ｅｘｎ〜Ｅｚｎとから、Ｈｘ（ｎ＋１／２）′（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｎ＋１／２）′（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｎ＋１／２）′（ｉ＋１／２，ｊ＋１／２，ｋ）を計算する。
例えば、Ｈｚ（ｎ＋１／２）′（ｉ＋１／２，ｊ＋１／２，ｋ）は、式（２５−１″），（５′）から得られる以下の式（２６′）をＨｚ（ｎ＋１／２）′（ｉ＋１／２，ｊ＋１／２，ｋ）について解くことにより計算できる。
｛Ｈｚ（ｎ＋１／２）′（ｉ＋１／２，ｊ＋１／２，ｋ）
−Ｈｚｎ−１／２（ｉ＋１／２，ｊ＋１／２，ｋ）｝／Δｔ
＝−｛Ｅｙｎ（ｉ＋１，ｊ＋１／２，ｋ）−Ｅｙｎ（ｉ，ｊ＋１／２，ｋ）｝／Δｘ
−σｘｘ＊｛Ｈｚ（ｎ＋１／２）′（ｉ＋１／２，ｊ＋１／２，ｋ）
＋Ｈｚｎ−１／２（ｉ＋１／２，ｊ＋１／２，ｋ）｝／２   …（２６′）
また、Ｈｘ（ｎ＋１／２）′（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｎ＋１／２）′（ｉ＋１／２，ｊ，ｋ＋１／２）もＨｚ（ｎ＋１／２）′（ｉ＋１／２，ｊ＋１／２，ｋ）と同様に計算できる。
（２）吸収境界磁界一時記憶手段Ｃ１E3bでＨｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ）として一時記憶されているＨｘｎ−１／２（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙｎ−１／２（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ）を消去して、計算されたＨｘ（ｎ＋１／２）′（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｎ＋１／２）′（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｎ＋１／２）′（ｉ＋１／２，ｊ＋１／２，ｋ）に更新して一時記憶する。すなわち、
Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２）
＝Ｈｘ（ｎ＋１／２）′（ｉ，ｊ＋１／２，ｋ＋１／２），
Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２）
＝Ｈｙ（ｎ＋１／２）′（ｉ＋１／２，ｊ，ｋ＋１／２），
Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ）
＝Ｈｚ（ｎ＋１／２）′（ｉ＋１／２，ｊ＋１／２，ｋ）
とする。
【０１０６】
ＳＴ２０２において、以下の（１），（２）の処理を実行し、ＳＴ２０３に移る。
（１）吸収境界磁界一時記憶手段Ｃ１E3bでＨｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ）として一時記憶されたＨｘ（ｎ＋１／２）′（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｎ＋１／２）′（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｎ＋１／２）′（ｉ＋１／２，ｊ＋１／２，ｋ）と、Δｘ〜Δｚ，Δｔ，σｘｘ*〜σｚｚ*，Ｅｘｎ〜Ｅｚｎとから、Ｈｘｎ＋１／２（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙｎ＋１／２（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚｎ＋１／２（ｉ＋１／２，ｊ＋１／２，ｋ）を計算する。
例えば、Ｈｘｎ＋１／２（ｉ，ｊ＋１／２，ｋ＋１／２）は、式（２５−２″），（６′）から得られる以下の式（２７′）をＨｘｎ＋１／２（ｉ，ｊ＋１／２，ｋ＋１／２）について解くことにより計算できる。
｛Ｈｚｎ＋１／２（ｉ＋１／２，ｊ＋１／２，ｋ）
−Ｈｚ（ｎ＋１／２）′（ｉ＋１／２，ｊ＋１／２，ｋ）｝／Δｔ
＝｛Ｅｘｎ（ｉ＋１／２，ｊ＋１，ｋ）−Ｅｘｎ（ｉ＋１／２，ｊ，ｋ）｝／Δｙ
−σｙｙ＊｛Ｈｚｎ＋１／２（ｉ＋１／２，ｊ＋１／２，ｋ）
＋Ｈｚ（ｎ＋１／２）′（ｉ＋１／２，ｊ＋１／２，ｋ）｝／２ 
…（２７′）
また、Ｈｙｎ＋１／２（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚｎ＋１／２（ｉ＋１／２，ｊ＋１／２，ｋ）もＨｘｎ＋１／２（ｉ，ｊ＋１／２，ｋ＋１／２）と同様に計算できる。
（２）吸収境界磁界一時記憶手段Ｃ１E3bでＨｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ）として一時記憶されているＨｘ（ｎ＋１／２）′（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｎ＋１／２）′（ｉ＋１／２，ｊ，ｋ＋１／２）を消去して、計算されたＨｘｎ＋１／２（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙｎ＋１／２（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚｎ＋１／２（ｉ＋１／２，ｊ＋１／２，ｋ）に更新して一時記憶する。すなわち、
Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２）＝Ｈｘｎ＋１／２（ｉ，ｊ＋１／２，ｋ＋１／２），
Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２）＝Ｈｙｎ＋１／２（ｉ＋１／２，ｊ，ｋ＋１／２），
Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ）＝Ｈｚｎ＋１／２（ｉ＋１／２，ｊ＋１／２，ｋ）
とする。
【０１０７】
ＳＴ２０３において、吸収境界５の全ての座標（ｉ，ｊ，ｋ）において、Ｈｘｎ＋１／２（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙｎ＋１／２（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚｎ＋１／２（ｉ＋１／２，ｊ＋１／２，ｋ）が計算されたか否かを判別する。ノー（Ｎ）の場合はＳＴ２０１に戻り、イエス（Ｙ）の場合は吸収境界磁界計算処理を終了して図７の電磁界解析処理に戻る。
【０１０８】
（実施例１の作用）
前記構成を備えた実施例１の電磁界解析装置としてのクライアントパソコンＰＣでは、アプリケーションプログラムとしての電磁界解析表示プログラムＡＰ１を起動した状態で、前記解析条件パラメタ入力画面１１（図３参照）が表示され、ユーザにより前記各入力欄１２〜２８に前記解析条件パラメタが入力され、前記解析開始画像２９が選択された場合に（図５のＳＴ１参照）、前記解析領域２および前記吸収境界５の全ての座標（ｉ，ｊ，ｋ）のベクトル量Ｅ（ｉ，ｊ，ｋ），Ｈ（ｉ，ｊ，ｋ）の計算、すなわち、電磁波の解析が、前記単位時間Δｔ毎に実行される（図５のＳＴ３〜ＳＴ７、図６のＳＴ１０１〜ＳＴ１０３、図７のＳＴ２０１〜ＳＴ２０３参照）。また、解析回数ｎが最大値Ｎｍａｘに達した場合には、前記電磁波の解析が終了し、解析結果が前記解析結果表示画像３１に表示される（図６参照）。
【０１０９】
このとき、実施例１では、前記吸収境界５における解析で、吸収境界条件に前記ＵＰＰ−ＰＭＬが適用されており、前記ベクトル量Ｅ（ｉ，ｊ，ｋ），Ｈ（ｉ，ｊ，ｋ）のベクトル成分Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２），Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ）が前記ＵＰＰ法によって計算される（図６のＳＴ１０１〜ＳＴ１０３，図７のＳＴ２０１〜ＳＴ２０３参照）。
なお、以降の説明より、各ベクトル成分Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２），Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ）における座標表記を、簡単のため省略し、Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚと記載する。
【０１１０】
ここで、従来公知の前記ＳＦ法では、（ｎ−１）・ΔｔのタイミングのＥｘｎ−１，Ｅｙｎ−１，Ｅｚｎ−１，Ｈｘｎ−１／２，Ｈｙｎ−１／２，Ｈｚｎ−１／２から、次のｎ・ΔｔのタイミングのＥｘｎ，Ｅｙｎ，Ｅｚｎ，Ｈｘｎ＋１／２，Ｈｙｎ＋１／２，Ｈｚｎ＋１／２を計算する場合、前記式（１０）〜（１５），（１０−１）〜（１５−１），（１０−２）〜（１５−２）を利用して計算する。例えば、Ｅｘｎ−１からＥｚｎを計算する場合、
（１）Ｅｚｘｎ−１，Ｈｙｎ−１／２，σｘｘ，ΔｔからＥｚｘｎを式（１２−１′）で計算し、
（２）Ｅｚｙｎ−１，Ｈｘｎ−１／２，σｙｙ，ΔｔからＥｚｙｎを式（１２−２′）で計算し、
（３）Ｅｚｘｎ，ＥｚｙｎからＥｚｎを式（１２′）で計算する必要があった。
すなわち、前記各ベクトル成分Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚにつき、（１）〜（３）の３つのステップで計算する必要があった。
【０１１１】
これに対し、前記ＵＰＰ法では、前記directional splittingの概念を利用して、前記式（２０）〜（２５）の右辺について、２つの軸方向の成分の式に分離し（式（２２−１），（２２−２）参照）、分離された２つの式を利用して、（ｎ−１）・ΔｔのタイミングのＥｘｎ−１，Ｅｙｎ−１，Ｅｚｎ−１，Ｈｘｎ−１／２，Ｈｙｎ−１／２，Ｈｚｎ−１／２から仮想的なタイミングであるｎ′・ΔｔのタイミングのＥｘｎ′，Ｅｙｎ′，Ｅｚｎ′，Ｈｘ（ｎ＋１／２）′，Ｈｙ（ｎ＋１／２）′，Ｈｚ（ｎ＋１／２）′を計算し、Ｅｘｎ′，Ｅｙｎ′，Ｅｚｎ′，Ｈｘ（ｎ＋１／２）′，Ｈｙ（ｎ＋１／２）′，Ｈｚ（ｎ＋１／２）′から次のｎ・ΔｔのタイミングのＥｘｎ，Ｅｙｎ，Ｅｚｎ，Ｈｘｎ＋１／２，Ｈｙｎ＋１／２，Ｈｚｎ＋１／２を計算することができる（式（２６），（２７）および式（２６′），（２７′）参照）。例えば、Ｅｘｎ−１からＥｚｎを計算する場合、
（１）Ｅｚｎ−１，Ｈｙｎ−１／２，σｘｘ，ΔｔからＥｚｎ′を式（２２−１″）で計算し、
（２）Ｅｚｎ′，Ｈｘｎ−１／２，σｙｙ，ΔｔからＥｚｎを式（２２−２″）で計算する。
すなわち、前記各ベクトル成分Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚにつき、（１），（２）の２つのステップで計算できる。したがって、前記ＵＰＰ法は、前記ＳＦ法と比べ、各ベクトル成分Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚを計算するために必要なステップ数を、前記単位時間Δｔ毎にそれぞれ１ステップずつ減らすことができる。
この結果、実施例１の電磁界解析装置としてのクライアントパソコンＰＣは、吸収境界条件として前記Ｂ−ＰＭＬを適用した従来の電磁界解析装置と比べ、計算時間を低減することができる。
【０１１２】
また、前記ＵＰＰ法では、前記単位時間Δｔ毎のベクトル量Ｅ（ｉ，ｊ，ｋ），Ｈ（ｉ，ｊ，ｋ）を計算するために、Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚの全６個の変数用の一時記憶領域（メモリ）があれば、（ｎ−１）・ΔｔのタイミングのＥｘｎ−１，Ｅｙｎ−１，Ｅｚｎ−１，Ｈｘｎ−１／２，Ｈｙｎ−１／２，Ｈｚｎ−１／２、仮想的なタイミングであるｎ′・ΔｔのタイミングのＥｘｎ′，Ｅｙｎ′，Ｅｚｎ′，Ｈｘ（ｎ＋１／２）′，Ｈｙ（ｎ＋１／２）′，Ｈｚ（ｎ＋１／２）′、次のｎ・ΔｔのタイミングのＥｘｎ，Ｅｙｎ，Ｅｚｎ，Ｈｘｎ＋１／２，Ｈｙｎ＋１／２，Ｈｚｎ＋１／２という具合に、計算結果を順次更新してゆくことができる（図６のＳＴ１０１の（２），ＳＴ１０２の（２）および図６のＳＴ２０１の（２），ＳＴ２０２の（２）参照）。
これに対し、前記ＳＦ法で同様の計算する場合には、１８個の変数（Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ，Ｅｙｘ，Ｅｚｘ，Ｈｙｘ，Ｈｚｘ，Ｅｘｙ，Ｅｚｙ，Ｅｘｚ，Ｅｙｚ，Ｈｘｙ，Ｈｚｙ，Ｈｘｚ，Ｈｙｚ）が必要であった。
したがって、吸収境界条件に前記ＵＰＰ−ＰＭＬを適用した場合には、前記Ｂ−ＰＭＬを適用した場合と比べ、吸収境界条件に前記変数の数を１／３（＝６／１８）にすることができる。
【０１１３】
なお、実施例１において、各変数に必要なメモリ量をｍとし、前記吸収境界５が１層増えた場合に必要なメモリ量は、以下の式（１６′）で表現できる。
［（（Ｉｍａｘ／Δｘ）×（Ｊｍａｘ／Δｙ）＋２ｄ）×２
＋（（Ｊｍａｘ／Δｙ）×（Ｋｍａｘ／Δｚ）＋２ｄ）×２
＋（（Ｋｍａｘ／Δｚ）×（Ｉｍａｘ／Δｘ）＋２ｄ）×２
＋（（Ｉｍａｘ／Δｘ）＋２ｄ）×４
＋（（Ｊｍａｘ／Δｙ）＋２ｄ）×４
＋（（Ｋｍａｘ／Δｚ）＋２ｄ）×４
＋８］×６×ｍ                                                      …（１６′）
【０１１４】
したがって、実施例１の電磁界解析装置としてのクライアントパソコンＰＣでは、吸収境界条件として前記ＵＰＰ−ＰＭＬを適用しているため、前記Ｂ−ＰＭＬを適用した場合と比べ、吸収境界５におけるＥｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚの計算に必要なメモリ量を１／３まで低減できる。また、実施例１の電磁界解析装置としてのクライアントパソコンＰＣでは、前記吸収境界５における電磁波の解析をする場合に必要なメモリ総量は、Ｂ−ＰＭＬの場合に比べて、１／３にすることができ、前記吸収境界５の層が多くなると、その効果も大きくなる。また、従来公知のＵ−ＰＭＬを適用した場合には、前記Ｂ−ＰＭＬを適用した場合と比べ、メモリ量を２／３にできるため、実施例１の前記クライアントパソコンＰＣは、前記Ｕ−ＰＭＬを適用した場合と比べても、吸収境界５におけるＥｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚの計算に必要なメモリ量を約１／２（＝（１／３）／（２／３））まで低減できる。
【０１１５】
（実験例）
ここで、実施例１の前記クライアントパソコンＰＣが、吸収境界条件として前記ＵＰＰ−ＰＭＬを適用したことにより、計算コスト（計算に必要なメモリ量および計算時間）や、解析精度がどのように影響するかを調べるために、次の実験例１および比較例１を準備した。
（実験例１）
実験例１の電磁界解析装置は、実施例１のクライアントパソコンＰＣと同様の構成である。すなわち、実験例１の電磁界解析装置は、実施例１における前記吸収境界解析手段Ｃ１Ｅが前記ＵＰＰ−ＰＭＬのＵＰＰ法に基づいて構成されている。
（比較例１）
比較例１の電磁界解析装置は、従来公知のＢ−ＰＭＬを吸収境界条件として適用し、前記ＳＦ法により吸収境界５における前記ベクトル成分Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚの計算を実行する。すなわち、比較例１の電磁界解析装置は、実施例１における前記吸収境界解析手段Ｃ１Ｅが前記ＵＰＰ−ＰＭＬのＵＰＰ法に基づいて構成されておらず、前記Ｂ−ＰＭＬのＳＦ法に基づいて構成されている。なお、比較例１の電磁界解析装置は、前記吸収境界解析手段Ｃ１Ｅ以外については前記実験例１と同様の構成である。
【０１１６】
（実験条件）
また、前記実験例１および前記比較例１についての実験条件は以下の通りである。
（解析条件パラメタについて）
まず、実験条件の解析条件パラメタを、以下の（１）〜（８）に示す。
（１）入力波源１の電界Ｅ＝（０，０，Ｅｚ）のＺ軸方向のベクトル成分Ｅｚ（ｉＧ，ｊＧ，ｋＧ＋１／２）については、
Ｅ０＝１［Ｖ／ｍ］，
ω０＝２π×３×１０９［ｒａｄ／ｓ］
が設定される。すなわち、Ｅｚ（ｉＧ，ｊＧ，ｋＧ＋１／２）は、
Ｅｚ（ｉＧ，ｊＧ，ｋＧ＋１／２）＝１×ｓｉｎ（（２π×３×１０９）ｔ）
で表現できる振幅（Ｅ０）が１［Ｖ／ｍ］、周波数（ω０／２π）が３［ＧＨｚ］の点波源、いわゆる、強制励振源である。
（２）前記セル３，６の単位長Δｘ，Δｙ，Δｚについては、
Δｘ＝Δｙ＝Δｚ＝５×１０−３［ｍ］
が設定される。すなわち、前記セル３，６は、各辺の長さが全て５［ｍｍ］の微小立方体である。
【０１１７】
（３）解析領域２の各辺の長さＩｍａｘ，Ｊｍａｘ，Ｋｍａｘについては、
Ｉｍａｘ＝Ｊｍａｘ＝Ｋｍａｘ＝０．５［ｍ］
が設定される。すなわち、解析領域２は、各辺の長さが全て０．５［ｍ］の立方体領域であり、総数が１００×１００×１００個（＝（０．５／（５×１０−３））３）の解析領域セル３の立方体領域である。
（４）解析領域２の誘電率εおよび透磁率μについては、
ε＝ε０＝８．８５４×１０−１２［Ｆ／ｍ］，
μ＝μ０＝４π×１０−７［Ｈ／ｍ］
が設定される。すなわち、解析領域２内は真空である。
（５）吸収境界５の層の値ｄについては、
ｄ＝２０
が設定される。すなわち、解析領域２の外方は、２０層の吸収境界５に覆われており、解析領域２および吸収境界５は、総数が１４０×１４０×１４０個（＝（１００＋２×２０）３）のセル３，６の立方体領域である。
【０１１８】
（６）吸収境界５の最大減衰量σｍａｘについては、
σｍａｘ＝０．３［Ｓ／ｍ］
が設定される。すなわち、各層ごとの吸収境界セル６の導電率σｘｘ，σｙｙ，σｚｚ（σｘｘ＝σｙｙ＝σｚｚ）および導磁率σｘｘ＊，σｙｙ＊，σｚｚ＊（σｘｘ＊＝σｙｙ＊＝σｚｚ＊）は、σｍａｘ＝０．３［Ｓ／ｍ］および式（８−１）〜（８−３）に基づいて設定される。
（７）単位時間Δｔについては、
Δｔ＝２［ｐｓ］（２×１０−１２［ｓ］）
が設定される。すなわち、Ｅｚ（ｉＧ，ｊＧ，ｋＧ＋１／２）の１周期である、
１／（３×１０９）［ｓ］
に対して、
｛１／（３×１０９）｝／（２×１０−１２）＝１／６×１０３［回］
の前記電磁界解析処理が行われる。
（８）解析回数の最大値Ｎｍａｘは実験結果を得るために十分に大きな値とする必要がある。したがって、例えば、最大値Ｎｍａｘについて、
Ｎｍａｘ＝１／６×１０３×５［回］
が設定される。すなわち、Ｅｚ（ｉＧ，ｊＧ，ｋＧ＋１／２）の５周期にあたる約１．７［ｎｓ］の解析を実行する（５×１／（３×１０９）＝１．６６６…×１０−９）。
【０１１９】
（実験１（コスト評価））
実験１では、前記解析条件パラメタ（１）〜（８）に基づいて、前記実験例１および前記比較例１の電磁界解析装置で、それぞれ前記電磁界解析処理（図５のＳＴ２〜ＳＴ１０、図６のＳＴ１０１〜ＳＴ１０３、図７のＳＴ２０１〜ＳＴ２０３参照）を実行したときの、計算時に必要なメモリ量と計算時間とを測定した。
（実験結果１（コスト評価））
図１０は実験例１の実験１の実験結果の説明図であり、比較例１の計算に必要なメモリ量および計算時間を１としたときの実験例１の値を示した棒グラフである。
図１０に示すように、実験例１の電磁界解析装置は、実験条件の解析条件パラメタ（１）〜（８）で前記電磁界解析処理を実行した場合、比較例１の電磁界解析装置と比べ、計算に必要なメモリ量が４５［％］、計算時間が９０［％］で計算できたことがわかる。すなわち、計算に必要なメモリ量が５５［％］、計算時間が１０［％］低減できたことがわかる。これは、ＦＤＴＤ法で計算する前記解析領域２では、実験例１と比較例１とは、同じメモリ量、計算ステップであるが、前記吸収境界５では、メモリ量が１／３、計算ステップ数が１ステップ減ったことに依っている。
【０１２０】
（実験２（ＰＭＬ境界面反射量測定））
実験２では、前記実験１の後に得られた実験例１および比較例１の解析結果ファイルから、以下の式（２８）で示されるＰＭＬ境界面からの反射量Ｒ（ｔ）［ｄＢ］を測定した。
Ｒ（ｔ）
＝（１／Ｖ）
×［∫Ｖ｛（Ｅｚ（ｔ）−Ｅｚｉｎｃ（ｔ））２／（Ｅｚｉｎｃ（ｔ））２｝ｄＶ］１／２
…（２８）
ここで、Ｖは解析領域２内の全ての座標（ｉ，ｊ，ｋ）を対象とする積分経路、Ｅｚ（ｔ）は実験開始からの時間ｔ［ｎｓ］における解析領域２の任意の座標（ｉ，ｊ，ｋ）の電界Ｅ（ｉ，ｊ，ｋ）のベクトル成分Ｅｚ（ｉ，ｊ，ｋ＋１／２）の値、Ｅｚｉｎｃ（ｔ）は実験開始からの時間ｔ［ｎｓ］における入力波源１の座標（ｉＧ，ｊＧ，ｋＧ）の電界Ｅ（ｉＧ，ｊＧ，ｋＧ）のベクトル成分Ｅｚ（ｉＧ，ｊＧ，ｋＧ＋１／２）の値である。
【０１２１】
（実験結果２（ＰＭＬ境界面反射量測定））
図１１は実験例１の実験２の実験結果の説明図であり、横軸に反射量を−２５０［ｄＢ］から−１５０［ｄＢ］の範囲で取り、縦軸に実験開始からの時間ｔ［ｎｓ］を０．６［ｎｓ］から１．４［ｎｓ］の範囲で取った実験例１と比較例１との実験結果のグラフである。
図１１に示すように、実験例１および比較例１の反射量は、０．６［ｎｓ］，０．８［ｎｓ］の時には、およそ−２５０［ｄＢ］，−１７５［ｄＢ］でほぼ等しいことがわかる。また、１［ｎｓ］，１．２［ｎｓ］，１．４［ｎｓ］の時には、実験例１と比較例１との反射量に誤差が生じていることがわかる。しかしながら、１［ｎｓ］の時には、実施例１がおよそ−１５２［ｄＢ］で比較例１がおよそ−１５４［ｄＢ］であるため、約２［ｄＢ］の微小な誤差（２％弱）である。また、１．２［ｎｓ］の時には、実施例１がおよそ−１５８［ｄＢ］で比較例１がおよそ−１６５［ｄＢ］であり、約７［ｄＢ］の誤差（５％弱）、１．４［ｎｓ］の時には、実施例１がおよそ−１６０［ｄＢ］で比較例１がおよそ−１６８［ｄＢ］であり、約８［ｄＢ］の誤差（約５％弱）である。したがって、前記反射量Ｒ（ｔ）は、実験例１と比較例１とを比較しても、ほとんど誤差が発生していないことがわかる。すなわち、実験例１と比較例１とは、ほぼ同等の解析精度を有することがわかる。
【０１２２】
したがって、前記実施例１の結果、実験例１の電磁界解析装置は、比較例１の電磁界解析装置と比べて、前記反射量Ｒ（ｔ）に基づいて測定できる解析精度がほぼ同等であり、且つ、計算コスト、すなわち、計算に必要なメモリ量や計算時間を低減できることがわかる。この結果、実施例１の電磁界解析装置としてのクライアントパソコンＰＣは、吸収境界条件として前記ＵＰＰ−ＰＭＬを適用することにより、吸収境界条件として従来公知のＢ−ＰＭＬを適用した場合の解析精度を維持しつつ、計算に必要なメモリ量や計算時間を低減することができる。
【０１２３】
（変更例）
以上、本発明の実施例を詳述したが、本発明は、前記実施例に限定されるものではなく、特許請求の範囲に記載された本発明の要旨の範囲内で、種々の変更を行うことが可能である。本発明の変更例（Ｈ01）〜（Ｈ04）を下記に例示する。
（Ｈ01）本発明の実施例１では、前記解析条件入力画像１１（図４参照）の誘電率入力欄２２と透磁率入力欄２３とによって、解析領域２の誘電率εと透磁率μとをそれぞれ一つ選択できるように設計したが、これに限定されず、例えば、解析領域２内の解析領域セル３毎に誘電率εや透磁率μを設定できるように変更することも可能である。例えば、障害物や遮蔽物、人体等が存在する場合に、対応する解析領域セル３に適切な誘電率εと透磁率μとを設定して、前記電磁界解析処理が実行できる。
【０１２４】
（Ｈ02）本発明の実施例１では、前記解析領域２および前記吸収境界５は、直交座標系の一つである、いわゆる、カーテシアン座標系により構成されていたが、これに限定されず、例えば、極座標系や円筒座標系で構成した場合でも、前記ＵＰＰ法により、前記ベクトル量Ｅ，Ｈの計算を行うことが可能である。また、本発明の実施例１では、前記解析領域２および前記吸収境界５を３次元空間として前記電磁界解析処理を実行したが、３次元に限定されず、４次元以上の多次元での解析も可能である。また、例えば、Ｚ軸方向に変化のない電磁波（平面波）の場合、２次元平面（ＸＹ平面）で前記電磁界解析処理を実行することも可能である。
（Ｈ03）本発明の実施例１では、前記入力波源１は、Ｚ軸方向の成分のみの電界Ｅ＝（０，０，Ｅｚ）に振幅を与える点波源、いわゆる、励振波源により構成されているが、これに限定されず、ＦＤＴＤ法で解析可能な全ての波源を適用可能である。例えば、前記入力波源１を散乱界、２次波源、電流源、磁流源等に置き換えることも可能である。なお、電流源（ＪＤ）、磁流源（ＪＢ）等に置き換えて前記電磁界解析処理を実行する場合、前記マクスウェル回転方程式系（式（１），（２）参照）に、これらの値ＪＤ，ＪＢを加えた方程式系に基づいてベクトル成分Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚをそれぞれ計算する必要がある。
【０１２５】
（Ｈ04）本発明の実施例１では、前記入力波源１から発生する前記入力波は電磁波であり、電界と磁界を解析したが、これに限定されず、ＦＤＴＤ法で解析可能であれば、例えば、音波、弾性波等を解析することも可能である。
例えば、弾性波の一例としての３次元弾性体中の超音波の伝搬について、以下の式（２９）〜（３２）で表すことができる場合がある。
ρ（∂／∂ｔ）ｖｘ
＝（∂／∂ｘ）τｘｘ＋（∂／∂ｙ）τｘｙ＋（∂／∂ｚ）τｘｚ          …（２９）
（∂／∂ｔ）τｘｘ＋γｘｘτｘｘ
＝（λ＋２μ）（∂／∂ｘ）ｖｘ＋（∂／∂ｙ）ｖｙ＋（∂／∂ｚ）ｖｚ   …（３０）
（∂／∂ｔ）τｘｙ＋ａｘｙτｘｙ
＝μ（（∂／∂ｙ）ｖｘ＋（∂／∂ｘ）ｖｙ）                           …（３１）
（∂／∂ｔ）τｘｚ＋ａｘｚτｘｚ
＝μ（（∂／∂ｚ）ｖｘ＋（∂／∂ｘ）ｖｚ）                           …（３２）
ここで、ｖｘ，ｖｙ，ｖｚはＸ軸方向、Ｙ軸方向、Ｚ軸方向における粒子速度、τｘｘは垂直応力、τｘｙ，τｘｚは剪断応力であり、ａｘｘ，ａｘｙ，ａｘｚは応力τｘｘ，τｘｙ，τｘｚに対する損失を表す抵抗係数であり、ρは媒質密度、λ，μはラメの定数である。
すなわち、粒子速度ｖｘの時間微分は、式（２９）で示すように、応力τｘｘ，τｘｙ，τｘｚの空間微分で表現でき、且つ、応力τｘｘ，τｘｙ，τｘｚの時間微分は、式（３０）〜（３２）で示すように、粒子速度ｖｘ，ｖｙ，ｖｚの空間微分で表現できる。なお、ｖｙ，ｖｚについても、応力τｙｘ，τｙｙ，τｙｚ，τｚｘ，τｚｙ，τｚｚと抵抗係数ａｙｘ，ａｙｙ，ａｙｚ，ａｚｘ，ａｚｙ，ａｚｚとを用いて、式（２９）〜（３２）と同様の式を示すことができる。この場合、電界と磁界に代えて、応力と速度について前記ＦＤＴＤ法が適用でき、前記弾性波を前記解析領域２および前記吸収境界５で解析できる。
【図面の簡単な説明】
【０１２６】
【図１】図１は本発明の実施例１の電磁界解析装置の説明図である。
【図２】図２は実施例１の電磁界解析装置としてのクライアントパソコンの制御部分が備えている各機能をブロック図（機能ブロック図）で示した図である。
【図３】図３は実施例１の電磁界解析表示プログラムＡＰ１により解析される解析領域および吸収境界の説明図であり、図３Ａは解析領域および吸収境界全体の説明図、図３Ｂは解析領域および吸収境界の各セルの拡大説明図である。
【図４】図４は実施例１の電磁界解析表示プログラムが表示する解析条件パラメタ入力画像の説明図である。
【図５】図５はＣＩＰ法におけるdirectional splittingの説明図である。
【図６】図６は実施例１の解析結果表示画像の説明図であり、解析結果を動画で表示したときの一例を示す図である。
【図７】図７は本発明の実施例１の電磁界解析装置の電磁界解析処理のフローチャートである。
【図８】図８は吸収境界電界計算処理のフローチャートであり、図７のＳＴ５のサブルーチンの説明図である。
【図９】図９は吸収境界磁界計算処理のフローチャートであり、図７のＳＴ７のサブルーチンの説明図である。
【図１０】図１０は実験例１の実験１の実験結果の説明図であり、比較例１の計算に必要なメモリ量および計算時間を１としたときの実験例１の値を示した棒グラフである。
【図１１】図１１は実験例１の実験２の実験結果の説明図であり、横軸に反射量を−２５０［ｄＢ］から−１５０［ｄＢ］の範囲で取り、縦軸に実験開始からの時間ｔ［ｎｓ］を０．６［ｎｓ］から１．４［ｎｓ］の範囲で取った実験例１と比較例１との実験結果のグラフである。
【図１２】図１２はＦＤＴＤ法の説明図であり、図１２ＡはＦＤＴＤ法の解析領域と解析領域セルの関係の説明図であり、図１２Ｂは図１２Ａの解析領域セルと電磁界の配置の説明図である。
【図１３】図１３は吸収境界条件の説明図である。
【符号の説明】
【０１２７】
Δｔ…解析時間間隔、Δｘ，Δｙ，Δｚ…分割領域の長さ、ε…誘電率、μ…透磁率、ＡＰ１…電磁界解析プログラム、Ｃ１A1…解析条件入力画像表示手段、Ｃ１Ｄ…解析領域解析手段、Ｃ１E1a＋Ｃ１E1b…仮想時間電磁界演算手段、Ｃ１E2a＋Ｃ１E2b…解析時間後電磁界演算手段、Ｃ１E3…吸収境界電磁界記憶手段、吸収境界要素記憶手段、Ｃ２…解析電磁界表示手段、解析場表示手段、Ｅ，Ｈ…電磁界、場、Ｅ…電界、第１要素、Ｅｘ，Ｅｙ，Ｅｚ，Ｈｘ，Ｈｙ，Ｈｚ…電界および磁界の強さ、場の強さ、Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２）…電界用分割領域、第１要素分割領域、Ｅｘ（ｉ＋１／２，ｊ，ｋ），Ｅｙ（ｉ，ｊ＋１／２，ｋ），Ｅｚ（ｉ，ｊ，ｋ＋１／２），Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ）…分割領域、ｎ・Δｔ，（ｎ＋１／２）Δｔ…解析時間後、ｎ′Δｔ，（ｎ＋１／２）′Δｔ…仮想時間、Ｈ…磁界、第２要素、Ｈｘ（ｉ，ｊ＋１／２，ｋ＋１／２），Ｈｙ（ｉ＋１／２，ｊ，ｋ＋１／２），Ｈｚ（ｉ＋１／２，ｊ＋１／２，ｋ）…磁界用分割領域、第２要素分割領域、Ｉｍａｘ，Ｊｍａｘ，Ｋｍａｘ，ｄ…解析領域および吸収境界領域の大きさ、ＰＣ…電磁界解析装置、場解析装置、１１…解析条件入力画像、２…解析領域、５…吸収境界領域。
【図１】
2008234335.tif 000003
【図２】
2008234335.tif 000004
【図３】
2008234335.tif 000005
【図４】
2008234335.tif 000006
【図５】
2008234335.tif 000007
【図７】
2008234335.tif 000008
【図８】
2008234335.tif 000009
【図９】
2008234335.tif 000010
【図１０】
2008234335.tif 000011
【図１１】
2008234335.tif 000012
【図１２】
2008234335.tif 000013
【図１３】
2008234335.tif 000014
【図６】
2008234335.tif 000015