(57)【要約】【課題】  対象とする金属の導体表面の解析に際し、ＦＤＴＤ法における定式化を行い、広帯域に渡って近似精度の高い導体表面上の電界（Ｅt）を算出することができる導体表面解析プログラムを提供することを課題とする。【解決手段】  単位等価回路として、各段ごとにレジスタンス（Ｒi）およびインダクタンス（Ｌi）を含むＲＬ単位等価回路を構成した上で、ＦＤＴＤ法における定式化を行うべく、前記インダクタンス（Ｌi）を、所定長さを有し先端が短絡されて一対となった無損失伝送線路に置換して構成される単位等価回路を用い、前記単位等価回路中の前記伝送回路の単位等価回路入力インピーダンスＺin(ω，ｉ)において回路定数（Ｒ，Ｌ）を入力することにより、前記無損失伝送線路の特性インピーダンス（Ｚ（Ｌ，ｉ））を選定し、導体表面上の電界（Ｅt）を算出すべくコンピュータを機能させる。【選択図】  図２
2005233841.tif 000002
【特許請求の範囲】
【請求項１】
対象とする金属の導体表面インピーダンス（ＺS（ω）  ω：角振動数）を近似すべく、単位等価回路を複数段（ｐ段）接続した等価回路であって、前記単位等価回路として、各段（ｉ＝１，２，…，ｐ）ごとにレジスタンス（Ｒi）およびインダクタンス（Ｌi）を含むＲＬ単位等価回路を構成した上で、ＦＤＴＤ法における定式化を行うべく、前記インダクタンス（Ｌi）を、所定長さを有し先端が短絡されて一対となった無損失伝送線路に置換して構成される単位等価回路を用い、前記単位等価回路中の前記伝送回路の入力インピーダンス（Ｚin(ω，ｉ)）において回路定数（Ｒ，Ｌ）を入力することにより、前記無損失伝送線路の特性インピーダンス（Ｚ（Ｌ，ｉ））を選定し、導体表面上の電界（Ｅt）を算出すべくコンピュータを機能させるプログラムであって、
コンピュータを、
前記インダクタンス（Ｌi）と前記特性インピーダンス（Ｚ（Ｌ，ｉ））との関係を示す第１式を予め記憶する第１式記憶手段、
所定時刻（（ｎ＋１／２）Δｔ  ｎ：演算ステップ数（整数））における前記単位等価回路に印加される電圧（Ｖn+1/2）を前記特性インピーダンス（Ｚ（Ｌ，ｉ））および前記時刻（（ｎ＋１／２）Δｔ）以前の時刻における前記単位等価回路に印加される電圧および前記単位等価回路に流れる電流を用いて示す第２式を予め記憶する第２式記憶手段、
前記回路定数である前記単位等価回路ごとのレジスタンス（Ｒi）およびインダクタンス（Ｌi）をそれぞれ入力する入力手段、
前記入力手段から入力された前記インダクタンス（Ｌi）を読み出し、前記第１式記憶手段に記憶された前記第１式に代入し、前記特性インピーダンス（Ｚ（Ｌ，ｉ））を演算する特性インピーダンス演算手段、
前記特性インピーダンス演算手段において演算された前記特性インピーダンス（Ｚ（Ｌ，ｉ））および前記第２式記憶手段に記憶された前記第２式を読み出し、前記導体表面上の電界（Ｅt）を、前記所定時刻（（ｎ＋１／２）Δｔ）の前記演算ステップ数（ｎ：整数）の変化に伴い繰り返し演算する導体表面解析手段、および
前記導体表面解析手段で算出された前記電界（Ｅt）の結果を出力する出力手段
として機能させ、
前記導体表面解析手段は、
コンピュータを、
前記入力手段から入力された前記レジスタンス（Ｒi）を読み出し、前記特性インピーダンス演算手段において演算された前記特性インピーダンス（Ｚ（Ｌ，ｉ））とともに、前記第２式記憶手段に記憶された前記第２式に代入し、前記所定時刻（（ｎ＋１／２）Δｔ）における前記単位等価回路に印加される電圧（Ｖn+1/2）を演算する単位等価回路電圧演算手段、および
前記単位等価回路電圧演算手段において、演算された前記電圧（Ｖn+1/2）を前記複数段（ｐ段）足し合わせて前記導体表面上の電界（Ｅt）を演算する電界演算手段
として機能させることを特徴とする導体表面解析プログラム。
【請求項２】
前記第１式は、以下の関数
2005233841.tif 000088
で定義され、
前記第２式は、以下の関数
2005233841.tif 000089
で定義されることを特徴とする請求項１記載の導体表面解析プログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、対象とする金属の導体表面インピーダンス（ＺS（ω））を近似すべく、単位等価回路を複数段（ｐ段）接続した等価回路であって、前記単位等価回路として、各段（ｉ＝１，２，…，ｐ）ごとにレジスタンス（Ｒi）およびインダクタンス（Ｌi）を含むＲＬ単位等価回路を構成した上で、ＦＤＴＤ法における定式化を行うべく、前記インダクタンス（Ｌi）を、所定長さを有し先端が短絡されて一対となった無損失伝送線路に置換して構成される単位等価回路を用い、前記単位等価回路中の伝送線路の単位等価回路入力インピーダンス（Ｚin(ω，ｉ)）において回路定数（Ｒ，Ｌ）を入力することにより、前記無損失伝送線路の特性インピーダンス（Ｚ（Ｌ，ｉ））を選定し、広帯域に渡って近似精度の高い導体表面上の電界（Ｅt）を算出することができる導体表面解析プログラムに関する。
【背景技術】
【０００２】
時間領域の電磁界解析法である差分時間領域法（ＦＤＴＤ法：例えば、非特許文献１参照）、伝送線路行列法（ＴＬＭ法：例えば、非特許文献２参照）、あるいは空間回路網法（ＳＮＭ：例えば、非特許文献３参照）等は、電磁波回路の設計・開発における重要なツールとして注目されている。近年の無線通信における使用周波数帯域の高周波化・広帯域化に伴い、回路に用いる金属の有限な導電率に起因する損失およびその周波数分散は無視することができなくなっており、電磁界解析でそれらを考慮できる導体処理がいくつか提案されてきた。
【０００３】
ＦＤＴＤ法には、Frequency Dependent ＦＤＴＤ法（例えば、非特許文献４，５参照）がある。この方法は、周波数分散性媒質の内部の電磁界を計算するものであるが、金属内部では電磁波が急激に減衰するため、内部の格子間隔を導体外部に比べ非常に細かく設定しなければならず、膨大な計算機メモリおよび計算時間が伴うという問題があった。
【０００４】
そこで、表面インピーダンス境界条件（Surface Impedance Boundary Condition：ＳＩＢＣ）を用いた処理法が提案されている。これは、電磁波が金属表面のごく浅いところまでしか浸透せず、金属表面のごく近い電磁界だけが外部の電磁界に影響することを考慮して、金属表面の電磁界のみ計算するものである。本来、表面インピーダンスは周波数領域の概念である。導体の表面インピーダンスは周波数の関数であり、時間領域では導体表面の電磁界と導体表面インピーダンスとの関係は畳み込み積分となる。これは、導体表面の過去の電磁界を全て蓄積しておかなければならず、さらに計算時間も著しく増加する（例えば、非特許文献６参照）。
【０００５】
時間領域においてＰｒｏｎｙ近似を用いる方法（例えば、非特許文献７，８参照）があるが、時間領域での近似誤差がどの周波数成分にどれくらいの影響を与えるのかが容易に把握できない。そのため、解析の対象として考えている周波数帯域内で近似誤差の影響を少なくするのが難しい。また、非特許文献８で述べられているＳＩＢＣを実用的観点から検討した結果、ＦＤＴＤ法の演算ステップ数が増加するにつれ、精度の良い近似式を得ることが急速に困難となることが確認されている（例えば、非特許文献９参照）。
【０００６】
さらに、周波数領域で近似を用いる方法もある（例えば、非特許文献１０〜１５参照）。これらは、近似に用いる独自の係数（例えば、近似の次数等）を変えることによって特定の誤差範囲で周波数分散性表面インピーダンスを広帯域近似可能であるが、その中で高精度近似できる周波数帯域は方法独自の係数によって自動的に決まる。つまり、ユーザは解析対象とする周波数帯域の導体表面インピーダンスを高精度近似したくてもできない場合が生じるという問題があった。
【非特許文献１】K.S.Yee, "Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media," IEEE Trans. Antennas & Propag., vol.AP-14, no.3, pp.302-307, May 1966.
【非特許文献２】S.Akhtarzad and P.B.Johns, "Solution of Maxwall's equations in three space dimensions and time by the t.l.m.method of numerical analysis," IEE Proc., vod.122, no.12, pp.1344-13348, Dec.1975.
【非特許文献３】吉田則信，深井一郎，福岡醇一，“電磁界の節点方程式による過渡解析”，信学論（Ｂ），vol.J63-B, no.9, pp.876-883, Sept. 1980.
【非特許文献４】R. J. Luebbers, F. Hunsberger, K. S. Kunz, R. B. Standler, M. Schneider, "A frequency-dependent finite-difference time-domain formulation for dispersive materials," IEEE Trans. Electromagnetic Compatibility, vol.32, pp.222-227, Aug. 1990.
【非特許文献５】R. J. Luebbers, F. Hunsberger, K. S. Kunz, "A frequency-dependent finite-difference time-domain formulation for transient propagation in plasma," IEEE Trans. Antennas & Propag., vol.39, no.1, pp.29-34, Jan. 1991.
【非特許文献６】F. M. Tesche, "On the inclusion of loss in time-domain solutions of electromagnetic interaction problem," IEEE Trans. Electromagnetic. Compat., vol.32, pp.1-4, Feb. 1990.
【非特許文献７】J. G. Maloney, G. S. Smith, "The use of surface impedance concepts in the finite-difference time-domain method," IEEE Trans. Antennas & Propag., vol.40, no.1, pp.38-48, Jan. 1992.
【非特許文献８】J. H. Beggs, R. J. Luebbers, K. S. Yee, K. S. Kunz, "Finite-difference time-domain implementation of surface impedance boundary conditions," IEEE Trans. Antennas & Propag., vol.40, no.1, pp.49-56, Jan. 1992. See also "Corrections," IEEE Trans. Antennas & Propag., vol.41, no.1, pp.118, Jan. 1993.
【非特許文献９】高田政明，飯田幸雄，“FDTD法における導体表面インピーダンス境界条件についての検討”，信学技報，MW2002-13,April 2002.
【非特許文献１０】C. F. Lee, R. T. Shin, J. A. Kong, "Time domain modeling of impedance boundary conditions," IEEE Trans. Microwave Theory & Tech., vol.40, no.9, pp.1847-1850, Sept. 1992.
【非特許文献１１】K. S. Oh, J. E. Schutt-Aine, "An efficient implementation of surface impedance boundary conditions for the finite-difference time-domain method," IEEE Trans. Antennas & Propag., vol.43, no.7, pp.660-666, July 1995.
【非特許文献１２】Y. Nishioka, O. Maeshima, T. Uno, S. Adachi, "Fdtd implementation of surface impedance boundary condition for dispersive layer backed by perfect conductor," IEICE Trans. Electron., vol.E81-C, no.12, pp.1902-1904, Dec. 1998.
【非特許文献１３】C. W. Penney, R. J. Luebbers, J. W. Schuster, "Scattering from coated targets using a frequency-dependent, surface impedance boundary condition in fdtd," IEEE Trans. Antennas & Propag., vol.44, no.4, pp.434-443, April 1996.
【非特許文献１４】M. C. Marcysiak, W. Gwarek, M. Sypniewski, "A simple and effective approach to fd-td modeling of structures including lossy metals," APMC'98, pp.991-993, 1998.
【非特許文献１５】本谷智宏，並木武文，伊藤公一，“導体損失を考慮したFDTD法による伝送線路の解析”，信学技報，MW99-215, Feb. 2000.
【発明の開示】
【発明が解決しようとする課題】
【０００７】
そこで、本発明は、上記従来技術の問題点を解決するために、対象とする金属の導体表面の解析に際し、ＦＤＴＤ法における定式化を行い、広帯域に渡って近似精度の高い導体表面上の電界（Ｅt）を算出することができる導体表面解析プログラムを提供することを課題とする。
【課題を解決するための手段】
【０００８】
本発明に係る導体表面解析プログラムは、上記課題を解決するためになされたものであり、対象とする金属の導体表面インピーダンス（ＺS（ω）  ω：角振動数）を近似すべく、単位等価回路を複数段（ｐ段）接続した等価回路であって、前記単位等価回路として、各段（ｉ＝１，２，…，ｐ）ごとにレジスタンス（Ｒi）およびインダクタンス（Ｌi）を含むＲＬ単位等価回路を構成した上で、ＦＤＴＤ法における定式化を行うべく、前記インダクタンス（Ｌi）を、所定長さを有し先端が短絡されて一対となった無損失伝送線路に置換して構成される単位等価回路を用い、前記単位等価回路中の前記伝送回路の入力インピーダンス（Ｚin(ω，ｉ)）において回路定数（Ｒ，Ｌ）を入力することにより、前記無損失伝送線路の特性インピーダンス（Ｚ（Ｌ，ｉ））を選定し、導体表面上の電界（Ｅt）を算出すべくコンピュータを機能させるプログラムであって、コンピュータを、前記インダクタンス（Ｌi）と前記特性インピーダンス（Ｚ（Ｌ，ｉ））との関係を示す第１式を予め記憶する第１式記憶手段、所定時刻（（ｎ＋１／２）Δｔ  ｎ：演算ステップ数（整数））における前記単位等価回路に印加される電圧（Ｖn+1/2）を前記特性インピーダンス（Ｚ（Ｌ，ｉ））および前記時刻（（ｎ＋１／２）Δｔ）以前の時刻における前記単位等価回路に印加される電圧および前記単位等価回路に流れる電流を用いて示す第２式を予め記憶する第２式記憶手段、前記回路定数である前記単位等価回路ごとのレジスタンス（Ｒi）およびインダクタンス（Ｌi）をそれぞれ入力する入力手段、前記入力手段から入力された前記インダクタンス（Ｌi）を読み出し、前記第１式記憶手段に記憶された前記第１式に代入し、前記特性インピーダンス（Ｚ（Ｌ，ｉ））を演算する特性インピーダンス演算手段、前記特性インピーダンス演算手段において演算された前記特性インピーダンス（Ｚ（Ｌ，ｉ））および前記第２式記憶手段に記憶された前記第２式を読み出し、前記導体表面上の電界（Ｅt）を、前記所定時刻（（ｎ＋１／２）Δｔ）の前記演算ステップ数（ｎ：整数）の変化に伴い繰り返し演算する導体表面解析手段、および前記導体表面解析手段で算出された前記電界（Ｅt）の結果を出力する出力手段として機能させ、前記導体表面解析手段は、コンピュータを、前記入力手段から入力された前記レジスタンス（Ｒi）を読み出し、前記特性インピーダンス演算手段において演算された前記特性インピーダンス（Ｚ（Ｌ，ｉ））とともに、前記第２式記憶手段に記憶された前記第２式に代入し、前記所定時刻（（ｎ＋１／２）Δｔ）における前記単位等価回路に印加される電圧（Ｖn+1/2）を演算する単位等価回路電圧演算手段、および前記単位等価回路電圧演算手段において、演算された前記電圧（Ｖn+1/2）を前記複数段（ｐ段）足し合わせて前記導体表面上の電界（Ｅt）を演算する電界演算手段として機能させることを特徴とするものである。
【０００９】
好ましくは、前記第１式は、以下の関数
2005233841.tif 000003
で定義され、前記第２式は、以下の関数
2005233841.tif 000004
で定義されるように構成される。
【発明の効果】
【００１０】
本発明の導体表面解析プログラムによれば、対象とする金属の導体表面インピーダンス（ＺS（ω））を近似すべく、単位等価回路を複数段（ｐ段）接続した等価回路であって、前記単位等価回路として、各段（ｉ＝１，２，…，ｐ）ごとにレジスタンス（Ｒi）およびインダクタンス（Ｌi）を含むＲＬ単位等価回路を構成した上で、前記インダクタンス（Ｌi）を、所定長さを有し先端が短絡されて一対となった無損失伝送線路に置換して構成される単位等価回路を用い、前記単位等価回路中の伝送線路の入力インピーダンス（Ｚin(ω，ｉ) ）において回路定数（Ｒ，Ｌ）を入力することにより、前記無損失伝送線路の特性インピーダンス（Ｚ（Ｌ，ｉ））を選定し、広帯域に渡って近似精度の高い導体表面上の電界（Ｅt）を算出することができる。
【発明を実施するための最良の形態】
【００１１】
以下、図面等を参照して、本発明の一実施形態について説明する。
【００１２】
なお、以下においては、具体的な導体表面解析プログラムの説明（フローチャート等を用いた説明）の前に、解析原理を成す各種関数の定義式、仮定の定め方等について説明する。
【００１３】
周波数領域において、導体表面に平行な電界Etおよび磁界Ｈtは、一般に以下の数１で表される。
【数１】
2005233841.tif 000005
【００１４】
数１において、ｎは導体表面外向きの単位法線ベクトル、ωは角周波数、εは導体の誘電率、μは導体の透磁率、σは導体の導電率である。また、ＺS(ω)は、磁気損失の無い導体の表面インピーダンスであり、導体は、実用上の周波数において、σ／ωε≫１の条件を満たすため、以下の数２のように置き換えることができる。
【数２】
2005233841.tif 000006
【００１５】
数１における時間領域表示は、導体表面インピーダンスＺSと磁界Ｈtとの畳み込みである。この畳み込みの計算を行うためには、導体表面上の過去の磁界をすべて記憶しておく必要があり、計算機で演算する場合において多大な容量のメモリを必要とする。この問題を回避するために、数２のＺS(ω）を以下の手順で近似する。
【００１６】
＜ＲＬ定式化＞
図１は本実施形態において導体表面インピーダンスＺs(ω)を近似するＲＬ等価回路の概略図である。本実施形態においては、図１に示すように、１つの抵抗１１と１つのコイル１２とを並列接続した回路（以下、ＲＬ単位等価回路１０という）を複数段（ｐ段）接続して、ＲＬ等価回路１が構成されている。
【００１７】
そして、このＲＬ等価回路１の等価回路インピーダンスＺ’in(ω)は、以下の数３にて定義される。
【数３】
2005233841.tif 000007
【００１８】
ここで、未知数Ri，Ｌi（ｉ＝１，２，…，ｐ）は、Ri＞０，Ｌi＞０であり、広帯域周波数において、以下の数４の関係を満たすものとする。
【数４】
2005233841.tif 000008
【００１９】
本実施形態における上記未知数Ri，Ｌi値の導出手順は、ｐ≧２の場合、全て同様の手順である。そこで、ここでは、ｐ＝５の場合を例にして導出手順を説明する。
【００２０】
まず、任意の角周波数ω0＞０において、以下の数５が成り立つための関係式（以下「ＲＬ関係式」という。）を求める。
【数５】
2005233841.tif 000009
【００２１】
上記数５は、すなわち、以下の数６に置き換えられる。
【数６】
2005233841.tif 000010
【００２２】
さて、上記数６を成立させるためのＲＬ関係式を導出するためには、この数６の左辺は実部と虚部が等しいため、まず、右辺の実部と虚部を等しくするための条件を求める。数６の右辺を実部と虚部とに整理すると、以下の数７となる。
【数７】
2005233841.tif 000011
【００２３】
ここで、上記数７の実部と虚部とを等しくするための一条件は、以下の数８〜数１２である。
【数８】
2005233841.tif 000012
【数９】
2005233841.tif 000013
【数１０】
2005233841.tif 000014
【数１１】
2005233841.tif 000015
【数１２】
2005233841.tif 000016
【００２４】
そして、上記数８を展開して整理すると、以下の数１３となる。
【数１３】
2005233841.tif 000017
【００２５】
本実施形態においては、ω0＞０，Ri＞０，Ｌi＞０（ｉ＝１，２，…，ｐ）と定義していることもあり、上式を満たす以下の数１４の条件を用いることとする。
【数１４】
2005233841.tif 000018
【００２６】
よって、以下の数１５が得られる。
【数１５】
2005233841.tif 000019
【００２７】
さらに、上記数１５は数１２を満足する。同様に、上記数９からは、以下の数１６が得られる。そして、この数１６は、数１１を満足する。
【数１６】
2005233841.tif 000020
【００２８】
同様に、上記数１０からは、以下の数１７が得られる。
【数１７】
2005233841.tif 000021
【００２９】
以上のことより、上記数７の実部と虚部とを等しくする条件である数１５〜数１７が得られる。ここでは、これらの条件（数１５〜数１７）を、以下の数１８にまとめることとする。
【数１８】
2005233841.tif 000022
【００３０】
次に、上記数６の左辺と右辺を等しくする。上記数７に数１８を代入すると、以下の数１９となり、実部と虚部とが等しい。
【数１９】
2005233841.tif 000023
【００３１】
したがって、上記数６を成り立たせるためには、以下の数２０の関係を維持すればよい。
【数２０】
2005233841.tif 000024
【００３２】
そして、本実施形態においては、未知数α，β，γ，δを用いて、Ｒ1〜Ｒ5を以下の数２１のように定義した。
【数２１】
2005233841.tif 000025
【００３３】
ここで、未知数α，β，γ，δは、Ri＞０（ｉ＝１，２，…，５）となる値、すなわち、α＞０，β＞０，γ＞０，δ＞０の値にする必要がある。上記数２０に数２１を適用すると、数２０の左辺はＲ1と未知数α，β，γ，δで表現でき、以下の数２２となる。
【数２２】
2005233841.tif 000026
【００３４】
以上の手順により、ｐ＝５（ＲＬ単位等価回路を５段接続した場合）において、上記数６を成り立たせるＲＬ関係式、数１８，数２１，数２２が得られる。なお、必要に応じて、上記数８〜数１２，数１４，数２１の条件を変えることで、上記数１８，数２１，数２２とは異なるＲＬ関係式を得ることも可能である。つまり、上述した導出手順によれば、ＲＬ単位等価回路を５段接続して構成された等価回路におけるＲＬ関係式が以下の数２３および数２４（数１８と同様）のように求められる。
【数２３】
2005233841.tif 000027
【数２４】
2005233841.tif 000028
【００３５】
同様の導出手順によれば、ＲＬ単位等価回路を２段接続して構成された等価回路（ｐ＝２）におけるＲＬ関係式は、以下の数２５および数２６のように求められる。
【数２５】
2005233841.tif 000029
【数２６】
2005233841.tif 000030
【００３６】
また、ＲＬ単位等価回路を３段接続して構成された等価回路（ｐ＝３）におけるＲＬ関係式は、以下の数２７および数２８のように求められる。
【数２７】
2005233841.tif 000031
【数２８】
2005233841.tif 000032
【００３７】
さらに、ＲＬ単位等価回路を４段接続して構成された等価回路（ｐ＝４）におけるＲＬ関係式は、以下の数２９および数３０のように求められる。
【数２９】
2005233841.tif 000033
【数３０】
2005233841.tif 000034
【００３８】
以上のようにして導出されて定義されたＲＬ関係式を利用して、本実施形態に係る回路定数解析プログラムは構成されており、これによって、図１に示すように、対象とする金属の導体表面インピーダンスＺs(ω)を近似すべく等価回路を構成して、この等価回路のインピーダンスＺ’in(ω)と導体表面インピーダンスＺs(ω)とを広帯域にわたって高精度に近似させることができる。
【００３９】
なお、上記実施形態においては、等価回路が２段接続された状態から、５段接続された状態までの計算を順番に行い、適切な値が得られた段階で演算が終了する場合について説明したが、本発明はこの構成に限定されず、必要に応じて、この段数については様々な設定が可能である。例えば、データ入力時に、予め何段の等価回路について演算を行うのかを入力すべく構成してもよい。また、例えば、５段のみの設定を行い、全てについて５段接続された等価回路についての演算を行うように構成してもよい。
【００４０】
＜ＦＤＴＤ定式化＞
以上のように求められたレジスタンスＲおよびインダクタンスＬを用いて導体表面上の電界Ｅtを微分方程式を用いて直接的に求めることも可能であるが、複雑かつ安定度があまりよくない演算をしなければならない。これを回避するために、本実施形態においては、単位等価回路中のインダクタンスＬを所定長さを有し先端が短絡されて一対となった無損失伝送線路に置換した等価回路として再構成し、演算する。
【００４１】
図２は本実施形態における導体表面上の電界Ｅtを演算するための伝送線路等価回路の概略図である。本実施形態においては、図２に示すように、図１におけるコイル１２の代わりに先端を短絡した所定長さｘの無損失伝送線路１３を備えた単位等価回路２０が複数段（ｐ段）接続され、伝送線路等価回路２を構成している。
【００４２】
そして、この単位等価回路２０中の伝送線路の入力インピーダンス（Ｚin(ω，ｉ) ただし、ｉ＝１，２，…，ｐ）は、以下の数３１にて定義される。ここで、βは伝搬定数である。
【数３１】
2005233841.tif 000035
【００４３】
一方、電磁波がこの線路を一往復する時間をＦＤＴＤ法における時間ステップΔｔに等しく選ぶことによって、以下の数３２の関係が成立する。ここで、ｖpは、電磁波の位相速度である。
【数３２】
2005233841.tif 000036
【００４４】
数３１および数３２より、ｘを消去して以下の数３３を得る。
【数３３】
2005233841.tif 000037
【００４５】
ωｔがπより十分小さいとき、上記数３３は、下記数３４のように表せる。
【数３４】
2005233841.tif 000038
ただし、
【数３５】
2005233841.tif 000039
【００４６】
ここで、数３５は、図１のコイル１２を図２の短絡伝送線路１３に置換することにより生じる誤差である。
【００４７】
そして、この単位等価回路２０の入力インピーダンスＺin(ω，ｉ)を図１のコイル１２のインダクタンスＬiで置き換えると下記数３６のように表せる。
【数３６】
2005233841.tif 000040
【００４８】
上記数３４、数３５および数３６より、誤差α(ω)が理想的に０であるとして、下記数３７で表される前記インダクタンス（Ｌi）と前記特性インピーダンスＺ（Ｌ，ｉ）との関係を示す第１式が得られ、この数３７（第１式）より特性インピーダンスＺ（Ｌ，ｉ）が選定される。
【数３７】
2005233841.tif 000041
【００４９】
なお、数３２より、数３８が得られる。
【数３８】
2005233841.tif 000042
【００５０】
以上より、上記数３８で表される長さｘを有し、数３７の特性インピーダンスＺ（Ｌ，ｉ）を有する短絡伝送線路１３は、近似的にインダクタンスＬiの等価物として動作する。この場合の誤差は、数３６で与えられるが、ΔｔがＦＤＴＤ法における時間ステップで与えられているため、この誤差α（ω）は、実際十分小さいものと言える。
【００５１】
続いて、上記特性インピーダンスＺ（Ｌ，ｉ）を用いて導体表面上の電界を求める。まず、上記と同様、電磁波の一往復時間がＦＤＴＤ法における時間ステップΔｔに等しい図２の短絡伝送線路１３において、時刻ｔ＝ｎΔｔ（ｎ：演算ステップ数（整数））における短絡伝送線路１３の入力端子電流をＩinとし、端子電圧をＶinとする。また、入射波（図２において右向きに進む）の電圧をｖinとする。
【００５２】
ここで、短絡伝送線路１３の入力端から入り込み短絡伝送線路１３を伝搬する電磁波は、当該短絡伝送線路１３の先端にて反射し、Δｔ秒後に入力端子まで到達する。前記先端は、短絡されているため、先端における反射により、電磁波電圧が逆相となる。
【００５３】
このとき、時刻ｎΔｔにおいて、短絡伝送線路１３の入力端子には、入射波（図２において右向きに進む）の電圧ｖinおよび反射波（図２において左向きに進む）の電圧ｖin-1が存在するため、入力端子電圧Ｖinは、以下の数３９となる。
【数３９】
2005233841.tif 000043
【００５４】
一方、短絡伝送線路１３を流れる電流は、先端において同相反射するため、以下の数４０となる。
【数４０】
2005233841.tif 000044
【００５５】
また、一般に、伝送線路においては、以下の数４１の関係が成立する。
【数４１】
2005233841.tif 000045
【００５６】
数３９〜数４１より、以下の数４２を得る。
【数４２】
2005233841.tif 000046
【００５７】
さらに、数４０の両辺にＺ（Ｌ，ｉ）を掛け、数４１および数４２を代入して以下の数４３を得る。
【数４３】
2005233841.tif 000047
【００５８】
また、図２において、レジスタンスＲiについての端子間電圧Ｖinおよび入力電流Ｉinの関係は、以下の数４４で示される。
【数４４】
2005233841.tif 000048
【００５９】
ここで、数４３を数４４に代入して整理すると、以下の数４５を得る。
【数４５】
2005233841.tif 000049
【００６０】
一方、中心差分近似において成立する関係式は、以下の数４６で得られる。
【数４６】
2005233841.tif 000050
【００６１】
上記数４６を用いて数４５を整理すると、以下の数４７が得られる。
【数４７】
2005233841.tif 000051
【００６２】
また、上記数４４を数４６を用いて整理することにより、以下の数４８が得られる。
【数４８】
2005233841.tif 000052
【００６３】
この数４７および数４８は、本発明における第２式であり、所定時刻（ｎ＋１／２）Δｔにおける前記単位等価回路に印加される電圧Ｖn+1/2を前記特性インピーダンスＺ（Ｌ，ｉ）および前記時刻（ｎ＋１／２）Δｔ以前の時刻における前記単位等価回路に印加される電圧および前記単位等価回路に流れる電流を用いて示すものである。特に、数４７は、時刻（ｎ＋１／２）Δｔにおける短絡伝送線路１３の端子間電圧がそれ以前の時刻における各パラメータによって演算されることを意味しており、当該時刻（ｎ＋１／２）Δｔにおいては、数４７の各パラメータは、既知の値を有していることとなるため、当該時刻（ｎ＋１／２）Δｔにおける短絡伝送線路１３の端子間電圧も数４７によって求められることとなる。
【００６４】
したがって、導体表面上の電界Ｅtは、以下の数４９で得られる式により定式化することができる。
【数４９】
2005233841.tif 000053
【００６５】
上記数４７〜数４９の各係数は、ＦＤＴＤのセルサイズおよび時間ステップΔｔは含まれていない。したがって、本実施形態の導体表面解析を用いることによる解析精度は、上記数３６で示される近似誤差α（ω）のみに影響されることとなり、上述の通り、この近似誤差は、十分小さいと言えるため、本実施形態の導体表面解析は、非常に高い解析精度を達成することができる。
【００６６】
次に、上述した解析原理に基づいた本実施形態に係る導体表面解析プログラムについて、フローチャートを用いて具体的に説明する。
【００６７】
＜回路定数解析プログラム＞
まず、図１におけるＲＬ定数化のための回路定数解析プログラムについて説明する。
【００６８】
図３は本実施形態に係る回路定数解析プログラムのメインルーチンのフローチャートである。
【００６９】
本実施形態に係る回路定数解析プログラムにおいては、まず、対象とする金属の導体表面インピーダンスＺs(ω)を近似するためにＲＬ単位等価回路１０（図１参照）を２段接続した等価回路を用いての解析処理が行われる。そして、予め与えられた許容範囲内にて導体表面インピーダンスＺs(ω)が近似されているか否かが判断される（ステップＳ２０１〜Ｓ２０８）。２段接続した等価回路にて、許容範囲内での近似が行われている場合には、その際の回路定数であるレジスタンスＲとリアクタンスＬとが出力される。
【００７０】
一方、２段接続された等価回路を用いた解析の結果、許容範囲内での近似がなされていない場合には、３段以上に接続された等価回路を用いての解析処理（ステップＳ２１２〜Ｓ２１８（３段），ステップＳ２２２〜Ｓ２２８（４段），ステップＳ２３２〜Ｓ２３８（５段））が行われることとなる。
【００７１】
本実施形態においては、２段接続にて近似できなかった場合には３段接続にて解析を行い、３段接続にて近似できなかった場合には４段接続にて解析を行い、４段接続にて近似できなかった場合には５段接続にて解析を行う。つまり、本実施形態においては、等価回路の段数ｐを増やすことで、より高精度に広帯域近似を行うことができることとなる。ただし、段数が増えることで処理時間も増えるため、本実施形態においては、処理効率を図るため、より少ない段数から解析を開始すべく構成されている。
【００７２】
図３に示すように、本実施形態においては、まず、コンピュータのキーボードあるいはマウス等の入力手段を用いて、データの入力処理が行われる（ステップＳ２０１）。このステップＳ２０１では、具体的に、導体の導電率σと、導体表面インピーダンスを等価近似したい周波数帯域の中心周波数ｆ0＝ω0／２πと、等価近似したい周波数帯域の下限周波数ｆLと、等価近似したい周波数帯域の上限周波数ｆHと、許容範囲を指定する変数ΔＲ％が入力される。この「ΔＲ」は、ｐ段等価回路の入力インピーダンスの実部Ｚin(ω,p)と、導体表面インピーダンスの実部Ｚs(ω)との差の近似したい周波数帯域内（ｆL〜ｆH）における許容最大値を示したものである。なお、この中心周波数ｆ0＝ω0／２πにおいて、導体表面インピーダンスＺs(ω)とｐ段の等価回路の入力インピーダンスＺin(ω,p)とが等しくなる。また、それぞれの実部の角周波数に対する微係数は、以下の数５０および数５１にて定義される。
【数５０】
2005233841.tif 000054
【数５１】
2005233841.tif 000055
【００７３】
また、本実施形態においては、ＲＬ関係式を構成する係数（α等）をある範囲内にわたってあるステップで変化させて解析を行う必要があり、さらに、近似したい周波数帯域内（ｆL〜ｆH）も、いくつかに等分して解析を行う必要がある。そこで、本実施形態においては、上記入力データとは別に、以下の数式（数５２，数５３）の値がコンピュータ内に設定されている。なお、これらの値は、予め設定される場合に限定されず、必要に応じて、先のステップＳ２０１で入力されるように構成してもよい。
【数５２】
2005233841.tif 000056
（ここで、「Ｎａ」は、係数（α等）をあるステップで変化させる際の範囲を示す。また、「ａ1」および「ａ2」は、係数（α等）を変化させる際のステップを定めるために用いられるものである。）
【数５３】
2005233841.tif 000057
（ここで、「Ｎｆ」は、近似したい周波数帯域内（ｆL〜ｆH）の分割数を示す。また、「ｆS」は、分割された各角周波数を示す。）
【００７４】
データ入力（Ｓ２０１）が行われると、次いで、サブルーチン２１の処理が行われる（ステップＳ２０２）。図４はサブルーチン２１のフローチャートを示したものである。
【００７５】
図３によれば、２段接続の等価回路についての解析処理を行うべく、先に述べたＲＬ関係式である数２５，数２６、および数２，数３に基づく微係数（数５０，数５１）の演算処理等が行われる（ステップＳ２０２１）。ここで係数αは、図４のステップＳ２０２１内に示すように、上記数５２と「ｉα」とを用いた指数関数「α＝ａ1（ａ2／ａ1）＾（ｉα／Ｎａ）」として与えられる。「ｉα」は、本実施形態においては、「１〜５００」までの変数として与えられる。
【００７６】
より具体的に説明すると、このステップＳ２０２１においては、ある範囲内にてあるステップで与えられる「α」ごとに、数２５，数２６を演算すると共に、各角周波数（Ｎｆ個の周波数）において以下の数５４を演算する。つまり、等価回路インピーダンスの実部の角周波数に対する微係数「Ｓin(ω,p)」と、導体表面インピーダンスの実部の角周波数に対する微係数「Ｓs(ω)」との差の絶対値を演算する。
【数５４】
2005233841.tif 000058
【００７７】
次いで、本実施形態においては、Ｎｆ個の周波数に対する「ΔＳ（ω，ｐ）」の総和が最小となる「α」の値を検出する（ステップＳ２０２２）。つまり、「ΔＳ（ω，ｐ）」の総和が最小となる「α」の値を得ることができる「ｉα」の値を検出する。このときの「ｉα」の値を「ｉαmin」と定義する。
【００７８】
次いで、検出した最小の箇所（「ｉαmin」にて得られる「α」の箇所）の付近について、より小さいステップにて係数αを定義する（ステップＳ２０２３）。具体的には、以下の数５５にて定義する。本実施形態においては、以下のようにより小さいステップにて定義された係数αによって、次工程以降の処理が行われる。
【数５５】
2005233841.tif 000059
【００７９】
本実施形態においては、上記数５４の演算を終えて、このサブルーチン２１の処理を終了して、図３のメインルーチンにリターンする。
【００８０】
サブルーチン２１の処理（Ｓ２０２）が終了した後、メインルーチンにおいては、「ｉＴ＝１」と定義されて（ステップＳ２０３）、サブルーチン２２の処理（ステップＳ２０４）が開始される。図５はサブルーチン２２のフローチャートを示したものである。このサブルーチン２２においては、入力データおよびサブルーチン２１にて得られた結果に基づいて、回路定数Ｒ，Ｌの解析を行う。ここでは、サブルーチン２１にて検出した極小となるところの付近について、より小さいステップで係数α等を変化させて、極小点をさらに特定していく。
【００８１】
図５によれば、まず、あらためてより狭い範囲に設定されたα（α5〜α6）についてのステップ幅（ｄａ）を、以下の数５６のように定める（ステップＳ２０４１）。
【数５６】
2005233841.tif 000060
【００８２】
次いで、ステップＳ２０４２においては、先に述べたＲＬ関係式である数２５，数２６、および数２，数３に基づく微係数（数５０，数５１）、加えてΔＲ（ω，ｐ）の演算処理等が行われる。「ΔＲ（ω，ｐ）」は、図４中に記載の通り、導体表面インピーダンス（Ｚs(ω)）の実部と、等価回路インピーダンス（Ｚin(ω,p)）の実部との差の絶対値を示している。ここで係数αは、図５のステップＳ２０４２内に示すように、上記数５６と「ｉα」と係数α5とを用いた関数「α＝α5＋（ｉα×ｄａ）」として与えられる。「ｉα」は、ここでは、「１〜５００」までの変数として与えられる。
【００８３】
より具体的に説明すると、このステップＳ２０４２においては、ある範囲内にてあるステップで与えられる「α」ごとに、数２５，数２６を演算すると共に、各角周波数（Ｎｆ個の周波数）についての数３４等を演算する。つまり、等価回路インピーダンスの実部の角周波数に対する微係数「Ｓin(ω,p)」と、導体表面インピーダンスの実部の角周波数に対する微係数「Ｓs(ω)」との差の絶対値、およびΔＲ（ω，ｐ）を演算する。
【００８４】
次いで、Ｎｆ個の周波数に対する「ΔＳ（ω，ｐ）」の総和が最小となる「α」の値を検出する（ステップＳ２０４３）。つまり、「ΔＳ（ω，ｐ）」の総和が最小となる「α」の値を得ることができる「ｉα」の値を検出する。このときの「ｉα」の値を「ｉαmin」と定義する。
【００８５】
次いで、検出した最小の箇所（「ｉαmin」にて得られる「α」の箇所）の付近について、より小さいステップにて係数αを再度定義する（ステップＳ２０４４）。具体的には、以下の数５７にて定義する。ここで再定義されたより小さいステップ範囲の係数αは、後述の処理にて、必要とする精度の近似が得られなかった場合（後述するステップＳ２０４７にて「Ｎｏ」と判断された場合）に使用されることとなる。
【数５７】
2005233841.tif 000061
【００８６】
次いで、Ｎｆ個の周波数の中で「ΔＲ（ω，ｐ）」が最大となる周波数「ｆmax」の値を検出する（ステップＳ２０４５）。そして、この周波数「ｆmax」のときに得られるインピーダンスの差「ΔＲ（ω，ｐ）max」と、導体表面インピーダンスの実部「Ｚs(ω)max」との比「ｒ」を演算する（ステップＳ２０４６）。この「ｒ」の値が小さければ小さい程、等価回路インピーダンス（Ｚin(ω,p)）が導体表面インピーダンス（Ｚs(ω)）に近似していることとなる。ここで、「ΔＲ（ω，ｐ）が最大となる周波数」の箇所がピックアップされるのは、この箇所が最も近似状態が好ましくない部分であって、ここが後述する許容範囲内に存在することが確認されれば、他の箇所は当然の如く許容範囲内にあることは明らかだからである。
【００８７】
次いで、演算された「ｒ」と、先に入力された許容範囲「ΔＲ」との対比判断処理が行われる（ステップＳ２０４７）。演算された「ｒ」が許容範囲「ΔＲ」よりも小さい場合（Ｓ２０４７にて「Ｙｅｓ」）には、「ｉｈ＝１」と定義されて（ステップＳ２０４８）、このサブルーチン２２の処理が終了して、図３のメインルーチンにリターンする。一方、演算された「ｒ」が許容範囲「ΔＲ」よりも大きい場合（Ｓ２０４７にて「Ｎｏ」）には、「α5＝α5m」、「α6＝α6m」、「ｉｈ＝０」と定義されて（ステップＳ２０４９）、このサブルーチン２２の処理が終了して、図３のメインルーチンにリターンする。
【００８８】
サブルーチン２２の処理（Ｓ２０４）が終了した後、メインルーチンにおいては、「ｉｈ＝１」であるか否かが判断される（ステップＳ２０５）。つまり、サブルーチン２２の処理において、許容範囲（ΔＲ）を満たす結果が得られたか否かが判断される。
【００８９】
許容範囲を満たす結果が得られた場合（Ｓ２０５にて「Ｙｅｓ」）には、得られた結果としての回路定数Ｒ1，Ｒ2，Ｌ1，Ｌ2の値が、コンピュータのディスプレイ、あるいはコンピュータに接続された出力装置等に出力され（ステップＳ２０６）、本実施形態に係る回路定数解析プログラムを用いた回路定数の解析処理が終了する。
【００９０】
一方、許容範囲を満たす結果が得られなかった場合（Ｓ２０５にて「Ｎｏ」）には、図５のステップＳ２０４９にて新たに設定された係数α5，α6を用いて再びサブルーチン２２の処理が行われることとなる。なお、ここで、サブルーチン２２が再び行われる場合には、その前に「ｉＴ≧３」であるか否かが判断される（ステップＳ２０７）。「ｉＴ」は、ステップＳ２０３において「ｉＴ＝１」と定義されているから、サブルーチン２２が一度行われた後には、「１」のままである。よって、ここではステップＳ２０７で「Ｎｏ」と判断されて、「ｉＴ」の値が「１」加算されてから（ステップＳ２０８）、２回目のサブルーチン２２の処理が行われることとなる。
【００９１】
すなわち、本実施形態においては、この２段接続の等価回路を用いた解析処理の結果が許容範囲内とならなかった場合であっても、サブルーチン２２の処理が３回まで繰り返し行われることとなる。そして、サブルーチン２２の処理を３回繰り返しても許容範囲内の結果を得られなければ（Ｓ２０７にて「Ｙｅｓ」）、次いで、３段接続の等価回路を用いた解析処理（Ｓ２１２〜Ｓ２１８）が行われる。
【００９２】
このように構成されているのは、サブルーチン２２の処理を４回以上繰り返して（すなわち、係数αの範囲をさらに細かく設定して）２段接続の等価回路を用いた処理を行うよりも、３段接続の等価回路を用いる方が、効率的で且つ高精度の解析結果を得ることができるからである。
【００９３】
上述したように、サブルーチン２２の処理を３回繰り返しても許容範囲（ΔＲ）内の結果が得られない場合（Ｓ２０７にて「Ｙｅｓ」）には、次いで、サブルーチン３１の処理が行われる（ステップＳ２１２）。図６はサブルーチン３１のフローチャートを示したものである。
【００９４】
図６によれば、３段接続の等価回路についての解析処理を行うべく、先に述べたＲＬ関係式である数２７，数２８、および数２，数３に基づく微係数（数５０，数５１）の演算処理等が行われる（ステップＳ２１２１）。ここで係数α，βは、図６のステップＳ２１２１内に示すように、上記数５２と「ｉα」，「ｉβ」とを用いた指数関数「α＝ａ1（ａ2／ａ1）＾（ｉα／Ｎａ）」，「β＝ａ1（ａ2／ａ1）＾（ｉβ／Ｎａ）」として与えられる。「ｉα」，「ｉβ」は、本実施形態においては、「１〜５００」までの変数として与えられる。
【００９５】
より具体的に説明すると、２段接続の場合と略同様に、このステップＳ２１２１においては、ある範囲内にてあるステップで与えられる「α」，「β」ごとに、数２７，数２８を演算すると共に、各角周波数（Ｎｆ個の周波数）において上記数５４を演算する。つまり、等価回路インピーダンスの実部の角周波数に対する微係数「Ｓin(ω,p)」と、導体表面インピーダンスの実部の角周波数に対する微係数「Ｓs(ω)」との差の絶対値を演算する。
【００９６】
次いで、Ｎｆ個の周波数に対する「ΔＳ（ω，ｐ）」の総和が最小となる「α」，「β」の値を検出する（ステップＳ２１２２）。つまり、「ΔＳ（ω，ｐ）」の総和が最小となる「α」，「β」の値を得ることができる「ｉα」，「ｉβ」の値を検出する。このときのそれぞれの値を「ｉαmin」，「ｉβmin」と定義する。
【００９７】
次いで、検出した最小の箇所（「ｉαmin」，「ｉβmin」にて得られる「α」，「β」の箇所）の付近について、より小さいステップにてそれぞれの係数を定義する（ステップＳ２１２３）。具体的には、以下の数５８，数５９にて定義する。本実施形態においては、以下のようにより小さいステップにて定義された係数α，βによって、次工程以降の処理が行われる。
【数５８】
2005233841.tif 000062
【数５９】
2005233841.tif 000063
【００９８】
本実施形態においては、上記数５８，数５９の演算を終えて、このサブルーチン３１の処理を終了して、図３のメインルーチンにリターンする。
【００９９】
サブルーチン３１の処理（Ｓ２１２）が終了した後、メインルーチンにおいては、「ｉＴ＝１」と定義されて（ステップＳ２１３）、サブルーチン３２の処理（ステップＳ２１４）が開始される。図７はサブルーチン３２のフローチャートを示したものである。
【０１００】
このサブルーチン３２においては、入力データおよびサブルーチン３１にて得られた結果に基づいて、回路定数Ｒ，Ｌの解析を行う。ここでは、サブルーチン３１にて検出した極小となるところの付近について、より小さいステップで係数α，β等を変化させて、極小点をさらに特定していく。
【０１０１】
図７によれば、まず、あらためてより狭い範囲に設定されたα（α5〜α6），β（β5〜β6）についてのステップ幅ｄａ，ｄｂを、以下の数６０のように定める（ステップＳ２１４１）。
【数６０】
2005233841.tif 000064
【０１０２】
次いで、ステップＳ２１４２においては、先に述べたＲＬ関係式である数２７，数２８、および数２，数３に基づく微係数（数５０，数５１）、加えてΔＲ（ω，ｐ）の演算処理等が行われる。「ΔＲ（ω，ｐ）」は、先にも述べた通り、導体表面インピーダンス（Ｚs(ω)）の実部と、等価回路インピーダンス（Ｚin(ω,p)）の実部との差の絶対値を示している。ここで係数α，βは、図７のステップＳ２１４２内に示すように、上記数４１と「ｉα」，「ｉβ」と係数α5，β5とを用いた関数「α＝α5＋（ｉα×ｄａ）」，「β＝β5＋（ｉβ×ｄｂ）」として与えられる。「ｉα」，「ｉβ」は、ここでは、「１〜５００」までの変数として与えられる。
【０１０３】
より具体的に説明すると、このステップＳ２１４２においては、ある範囲内にてあるステップで与えられる「α」，「β」ごとに、数２７，数２８を演算すると共に、各角周波数（Ｎｆ個の周波数）についての数５４等を演算する。つまり、等価回路インピーダンスの実部の角周波数に対する微係数「Ｓin(ω,p)」と、導体表面インピーダンスの実部の角周波数に対する微係数「Ｓs(ω)」との差の絶対値、およびΔＲ（ω，ｐ）を演算する。
【０１０４】
次いで、Ｎｆ個の周波数に対する「ΔＳ（ω，ｐ）」の総和が最小となる「α」，「β」の値を検出する（ステップＳ２１４３）。つまり、「ΔＳ（ω，ｐ）」の総和が最小となる「α」，「β」の値を得ることができる「ｉα」，「ｉβ」の値を検出する。このときの「ｉα」，「ｉβ」の値を「ｉαmin」，「ｉβmin」と定義する。
【０１０５】
次いで、検出した最小の箇所（例えば、「ｉαmin」にて得られる「α」の箇所）の付近について、より小さいステップにて係数α，βを再度定義する（ステップＳ２１４４）。具体的には、以下の数６１，数６２にて定義する。ここで再定義されたより小さいステップ範囲の係数α，βは、後述の処理にて、必要とする精度の近似が得られなかった場合（後述するステップＳ２１４７にて「Ｎｏ」と判断された場合）に使用されることとなる。
【数６１】
2005233841.tif 000065
【数６２】
2005233841.tif 000066
【０１０６】
次いで、Ｎｆ個の周波数の中で「ΔＲ（ω，ｐ）」が最大となる周波数「ｆmax」の値を検出する（ステップＳ２１４５）。そして、この周波数「ｆmax」のときに得られるインピーダンスの差「ΔＲ（ω，ｐ）max」と、導体表面インピーダンスの実部「Ｚs(ω)max」との比「ｒ」を演算する（ステップＳ２１４６）。この「ｒ」の値が小さければ小さい程、等価回路インピーダンスＺin(ω,p)が導体表面インピーダンスＺs(ω)に近似していることとなる。ここで、「ΔＲ（ω，ｐ）が最大となる周波数」の箇所がピックアップされるのは、この箇所が最も近似状態が好ましくない部分であって、ここが後述する許容範囲内に存在することが確認されれば、他の箇所は当然の如く許容範囲内にあることは明らかだからである。
【０１０７】
次いで、演算された「ｒ」と、先に入力された許容範囲「ΔＲ」との対比判断処理が行われる（ステップＳ２１４７）。演算された「ｒ」が許容範囲「ΔＲ」よりも小さい場合（Ｓ２１４７にて「Ｙｅｓ」）には、「ｉｈ＝１」と定義されて（ステップＳ２１４８）、このサブルーチン３２の処理が終了して、図３のメインルーチンにリターンする。一方、演算された「ｒ」が許容範囲「ΔＲ」よりも大きい場合（Ｓ２１４７にて「Ｎｏ」）には、「α5＝α5m，α6＝α6m」、「β5＝β5m，β6＝β6m」、「ｉｈ＝０」と定義されて（ステップＳ２１４９）、このサブルーチン３２の処理が終了して、図３のメインルーチンにリターンする。
【０１０８】
サブルーチン３２の処理（Ｓ２１４）が終了した後、メインルーチンにおいては、「ｉｈ＝１」であるか否かが判断される（ステップＳ２１５）。つまり、サブルーチン３２の処理において、許容範囲（ΔＲ）を満たす結果が得られたか否かが判断される。
【０１０９】
許容範囲を満たす結果が得られた場合（Ｓ２１５にて「Ｙｅｓ」）には、得られた結果としての回路定数Ｒ1，Ｒ2，Ｒ3，Ｌ1，Ｌ2，Ｌ3の値が、コンピュータのディスプレイ、あるいはコンピュータに接続された出力装置等に出力され（ステップＳ２０６）、本実施形態に係る回路定数解析プログラムを用いた回路定数の解析処理が終了する。
【０１１０】
一方、許容範囲を満たす結果が得られなかった場合（Ｓ２１５にて「Ｎｏ」）には、図６のステップＳ２１４９にて新たに設定された係数α5，α6、β5，β6を用いて再びサブルーチン３２の処理が行われることとなる。なお、ここで、サブルーチン３２が再び行われる場合には、２段接続の場合と同様に、その前に「ｉＴ≧３」であるか否かが判断される（ステップＳ２１７）。「ｉＴ」は、ステップＳ２１３において「ｉＴ＝１」と定義されているから、サブルーチン３２が一度行われた後には、「１」のままである。よって、ここではステップＳ２１７で「Ｎｏ」と判断されて、「ｉＴ」の値が「１」加算されてから（ステップＳ２１８）、２回目のサブルーチン２２の処理が行われることとなる。
【０１１１】
つまり、本実施形態においては、２段接続の場合と同様に、この３段接続の等価回路を用いた解析処理の結果が許容範囲内とならなかった場合であっても、サブルーチン３２の処理が３回まで繰り返し行われることとなる。そして、サブルーチン３２の処理を３回繰り返しても許容範囲内の結果を得られなければ（Ｓ２１７にて「Ｙｅｓ」）、次いで、４段接続の等価回路を用いた解析処理（Ｓ２２２〜Ｓ２２８）が行われる。このように構成されているのは、２段接続のところで述べたのと同様の理由からである。
【０１１２】
上述したように、サブルーチン３２の処理を３回繰り返しても許容範囲（ΔＲ）内の結果が得られない場合（Ｓ２１７にて「Ｙｅｓ」）には、次いで、サブルーチン４１の処理が行われる（ステップＳ２２２）。図８はサブルーチン４１のフローチャートを示したものである。
【０１１３】
図８によれば、４段接続の等価回路についての解析処理を行うべく、先に述べたＲＬ関係式である数２９，数３０、および数２，数３に基づく微係数（数５０，数５１）の演算処理等が行われる（ステップＳ２２２１）。ここで係数α，β，γは、図８のステップＳ２２２１内に示すように、上記数５２と「ｉα」，「ｉβ」，「ｉγ」とを用いた指数関数「α＝ａ1（ａ2／ａ1）＾（ｉα／Ｎａ）」，「β＝ａ1（ａ2／ａ1）＾（ｉβ／Ｎａ）」，「γ＝ａ1（ａ2／ａ1）＾（ｉγ／Ｎａ）」として与えられる。「ｉα」，「ｉβ」，「ｉγ」は、本実施形態においては、「１〜５００」までの変数として与えられる。
【０１１４】
より具体的に説明すると、２段接続等の場合と略同様に、このステップＳ２２２１においては、ある範囲内にてあるステップで与えられる「α」，「β」，「γ」ごとに、数２９，数３０を演算すると共に、各角周波数（Ｎｆ個の周波数）において上記数５４を演算する。つまり、等価回路インピーダンスの実部の角周波数に対する微係数「Ｓin(ω,p)」と、導体表面インピーダンスの実部の角周波数に対する微係数「Ｓs(ω)」との差の絶対値を演算する。
【０１１５】
次いで、Ｎｆ個の周波数に対する「ΔＳ（ω，ｐ）」の総和が最小となる「α」，「β」，「γ」の値を検出する（ステップＳ２２２２）。つまり、「ΔＳ（ω，ｐ）」の総和が最小となる「α」，「β」，「γ」の値を得ることができる「ｉα」，「ｉβ」，「ｉγ」の値を検出する。このときのそれぞれの値を「ｉαmin」，「ｉβmin」，「ｉγmin」と定義する。
【０１１６】
次いで、検出した最小の箇所（「ｉαmin」，「ｉβmin」，「ｉγmin」にて得られる「α」，「β」，「γ」の箇所）の付近について、より小さいステップにてそれぞれの係数を定義する（ステップＳ２２２３）。具体的には、以下の数６３，数６４，数６５にて定義する。本実施形態においては、以下のようにより小さいステップにて定義された係数α，β，γによって、次工程以降の処理が行われる。
【数６３】
2005233841.tif 000067
【数６４】
2005233841.tif 000068
【数６５】
2005233841.tif 000069
【０１１７】
本実施形態においては、上記数６３〜数６５の演算を終えて、このサブルーチン４１の処理を終了して、図３のメインルーチンにリターンする。
【０１１８】
サブルーチン４１の処理（Ｓ２２２）が終了した後、メインルーチンにおいては、「ｉＴ＝１」と定義されて（ステップＳ２２３）、サブルーチン４２の処理（ステップＳ２２４）が開始される。図９はサブルーチン４２のフローチャートを示したものである。
【０１１９】
このサブルーチン４２においては、入力データおよびサブルーチン４１にて得られた結果に基づいて、回路定数Ｒ，Ｌの解析を行う。ここでは、サブルーチン４１にて検出した極小となるところの付近について、より小さいステップで係数α，β，γ等を変化させて、極小点をさらに特定していく。
【０１２０】
図９によれば、まず、あらためてより狭い範囲に設定されたα（α5〜α6），β（β5〜β6），γ（γ5〜γ6）についてのステップ幅ｄａ，ｄｂ，ｄｇを、以下の数６６のように定める（ステップＳ２２４１）。
【数６６】
2005233841.tif 000070
【０１２１】
次いで、ステップＳ２２４２においては、先に述べたＲＬ関係式である数２９，数３０、および数２，数３に基づく微係数（数５０，数５１）、加えてΔＲ（ω，ｐ）の演算処理等が行われる。「ΔＲ（ω，ｐ）」は、先にも述べた通り、導体表面インピーダンスＺs(ω)の実部と、等価回路インピーダンスＺin(ω,p)の実部との差の絶対値を示している。ここで係数α，β，γは、図９のステップＳ２２４２内に示すように、上記数６６と「ｉα」，「ｉβ」，「ｉγ」と係数α5，β5，γ5とを用いた関数「α＝α5＋（ｉα×ｄａ）」，「β＝β5＋（ｉβ×ｄｂ）」，「γ＝γ5＋（ｉγ×ｄｇ）」として与えられる。「ｉα」，「ｉβ」，「ｉγ」は、ここでは、「１〜５００」までの変数として与えられる。
【０１２２】
より具体的に説明すると、このステップＳ２２４２においては、ある範囲内にてあるステップで与えられる「α」，「β」，「γ」ごとに、数２９，数３０を演算すると共に、各角周波数（Ｎｆ個の周波数）についての数５４等を演算する。つまり、等価回路インピーダンスの実部の角周波数に対する微係数「Ｓin(ω,p)」と、導体表面インピーダンスの実部の角周波数に対する微係数「Ｓs(ω)」との差の絶対値、およびΔＲ（ω，ｐ）を演算する。
【０１２３】
次いで、Ｎｆ個の周波数に対する「ΔＳ（ω，ｐ）」の総和が最小となる「α」，「β」，「γ」の値を検出する（ステップＳ２２４３）。つまり、「ΔＳ（ω，ｐ）」の総和が最小となる「α」，「β」，「γ」の値を得ることができる「ｉα」，「ｉβ」，「ｉγ」の値を検出する。このときの「ｉα」，「ｉβ」，「ｉγ」の値を「ｉαmin」，「ｉβmin」，「ｉγmin」と定義する。
【０１２４】
次いで、検出した最小の箇所（例えば、「ｉαmin」にて得られる「α」の箇所）の付近について、より小さいステップにて係数α，β，γを再度定義する（ステップＳ２２４４）。具体的には、以下の数６７〜数６９にて定義する。ここで再定義されたより小さいステップ範囲の係数α，β，γは、後述の処理にて、必要とする精度の近似が得られなかった場合（後述するステップＳ２２４７にて「Ｎｏ」と判断された場合）に使用されることとなる。
【数６７】
2005233841.tif 000071
【数６８】
2005233841.tif 000072
【数６９】
2005233841.tif 000073
【０１２５】
次いで、Ｎｆ個の周波数の中で「ΔＲ（ω，ｐ）」が最大となる周波数「ｆmax」の値を検出する（ステップＳ２２４５）。そして、この周波数「ｆmax」のときに得られるインピーダンスの差「ΔＲ（ω，ｐ）max」と、導体表面インピーダンスの実部「Ｚs(ω)max」との比「ｒ」を演算する（ステップＳ２２４６）。この「ｒ」の値が小さければ小さい程、等価回路インピーダンスＺin(ω,p)が導体表面インピーダンスＺs(ω)に近似していることとなる。ここで、「ΔＲ（ω，ｐ）が最大となる周波数」の箇所がピックアップされるのは、この箇所が最も近似状態が好ましくない部分であって、ここが後述する許容範囲内に存在することが確認されれば、他の箇所は当然の如く許容範囲内にあることは明らかだからである。
【０１２６】
次いで、演算された「ｒ」と、先に入力された許容範囲「ΔＲ」との対比判断処理が行われる（ステップＳ２２４７）。演算された「ｒ」が許容範囲「ΔＲ」よりも小さい場合（Ｓ２２４７にて「Ｙｅｓ」）には、「ｉｈ＝１」と定義されて（ステップＳ２２４８）、このサブルーチン４２の処理が終了して、図３のメインルーチンにリターンする。一方、演算された「ｒ」が許容範囲「ΔＲ」よりも大きい場合（Ｓ２２４７にて「Ｎｏ」）には、「α5＝α5m，α6＝α6m」、「β5＝β5m，β6＝β6m」、「γ5＝γ5m，γ6＝γ6m」、「ｉｈ＝０」と定義されて（ステップＳ２２４９）、このサブルーチン４２の処理が終了して、図３のメインルーチンにリターンする。
【０１２７】
サブルーチン４２の処理（Ｓ２２４）が終了した後、メインルーチンにおいては、「ｉｈ＝１」であるか否かが判断される（ステップＳ２２５）。つまり、サブルーチン４２の処理において、許容範囲（ΔＲ）を満たす結果が得られたか否かが判断される。
【０１２８】
許容範囲を満たす結果が得られた場合（Ｓ２２５にて「Ｙｅｓ」）には、得られた結果としての回路定数Ｒ1〜Ｒ4，Ｌ1〜Ｌ4の値が、コンピュータのディスプレイ、あるいはコンピュータに接続された出力装置等に出力され（ステップＳ２０６）、本実施形態に係る回路定数解析プログラムを用いた回路定数の解析処理が終了する。
【０１２９】
一方、許容範囲を満たす結果が得られなかった場合（Ｓ２２５にて「Ｎｏ」）には、図９のステップＳ２２４９にて新たに設定された係数α5，α6、β5，β6、γ5，γ6を用いて再びサブルーチン４２の処理が行われることとなる。なお、ここで、サブルーチン４２が再び行われる場合には、２段接続等の場合と同様に、その前に「ｉＴ≧３」であるか否かが判断される（ステップＳ２２７）。「ｉＴ」は、ステップＳ２２３において「ｉＴ＝１」と定義されているから、サブルーチン４２が一度行われた後には、「１」のままである。よって、ここではステップＳ２２７で「Ｎｏ」と判断されて、「ｉＴ」の値が「１」加算されてから（ステップＳ２２８）、２回目のサブルーチン２２の処理が行われることとなる。
【０１３０】
つまり、本実施形態においては、２段接続等の場合と同様に、この４段接続の等価回路を用いた解析処理の結果が許容範囲内とならなかった場合であっても、サブルーチン４２の処理が３回まで繰り返し行われることとなる。そして、サブルーチン４２の処理を３回繰り返しても許容範囲内の結果を得られなければ（Ｓ２２７にて「Ｙｅｓ」）、次いで、５段接続の等価回路を用いた解析処理（Ｓ２３２〜Ｓ２３９）が行われる。このように構成されているのは、２段接続のところで述べたのと同様の理由からである
【０１３１】
上述したように、サブルーチン４２の処理を３回繰り返しても許容範囲（ΔＲ）内の結果が得られない場合（Ｓ２２７にて「Ｙｅｓ」）には、次いで、サブルーチン５１の処理が行われる（ステップＳ２３２）。図１０はサブルーチン５１のフローチャートを示したものである。
【０１３２】
図１０によれば、５段接続の等価回路についての解析処理を行うべく、先に述べたＲＬ関係式である数２３，数２４、および数２，数３に基づく微係数（数５０，数５１）の演算処理等が行われる（ステップＳ２３２１）。ここで係数α，β，γ，δは、図１０のステップＳ２３２１内に示すように、上記数５２と「ｉα」，「ｉβ」，「ｉγ」，「ｉδ」とを用いた指数関数「α＝ａ1（ａ2／ａ1）＾（ｉα／Ｎａ）」，「β＝ａ1（ａ2／ａ1）＾（ｉβ／Ｎａ）」，「γ＝ａ1（ａ2／ａ1）＾（ｉγ／Ｎａ）」，「δ＝ａ1（ａ2／ａ1）＾（ｉδ／Ｎａ）」として与えられる。「ｉα」，「ｉβ」，「ｉγ」，「ｉδ」は、本実施形態においては、「１〜５００」までの変数として与えられる。
【０１３３】
より具体的に説明すると、２段接続等の場合と略同様に、このステップＳ２３２１においては、ある範囲内にてあるステップで与えられる「α」，「β」，「γ」，「δ」ごとに、数２３，数２４を演算すると共に、各角周波数（Ｎｆ個の周波数）において上記数５４を演算する。つまり、等価回路インピーダンスの実部の角周波数に対する微係数「Ｓin(ω,p)」と、導体表面インピーダンスの実部の角周波数に対する微係数「Ｓs(ω)」との差の絶対値を演算する。
【０１３４】
次いで、Ｎｆ個の周波数に対する「ΔＳ（ω，ｐ）」の総和が最小となる「α」，「β」，「γ」，「δ」の値を検出する（ステップＳ２３２２）。つまり、「ΔＳ（ω，ｐ）」の総和が最小となる「α」，「β」，「γ」，「δ」の値を得ることができる「ｉα」，「ｉβ」，「ｉγ」，「ｉδ」の値を検出する。このときのそれぞれの値を「ｉαmin」，「ｉβmin」，「ｉγmin」，「ｉδmin」と定義する。
【０１３５】
次いで、検出した最小の箇所（「ｉαmin」，「ｉβmin」，「ｉγmin」，「ｉδmin」にて得られる「α」，「β」，「γ」，「δ」の箇所）の付近について、より小さいステップにてそれぞれの係数を定義する（ステップＳ２３２３）。具体的には、以下の数７０〜数７３にて定義する。本実施形態においては、以下のようにより小さいステップにて定義された係数α，β，γ，δによって、次工程以降の処理が行われる。
【数７０】
2005233841.tif 000074
【数７１】
2005233841.tif 000075
【数７２】
2005233841.tif 000076
【数７３】
2005233841.tif 000077
【０１３６】
本実施形態においては、上記数７０〜数７３の演算を終えて、このサブルーチン５１の処理を終了して、図３のメインルーチンにリターンする。
【０１３７】
サブルーチン５１の処理（Ｓ２３２）が終了した後、メインルーチンにおいては、「ｉＴ＝１」と定義されて（ステップＳ２３３）、サブルーチン５２の処理（ステップＳ２３４）が開始される。図１１はサブルーチン５２のフローチャートを示したものである。このサブルーチン５２においては、入力データおよびサブルーチン５１にて得られた結果に基づいて、回路定数Ｒ，Ｌの解析を行う。ここでは、サブルーチン５１にて検出した極小となるところの付近について、より小さいステップで係数α，β，γ，δ等を変化させて、極小点をさらに特定していく。
【０１３８】
図１１によれば、まず、あらためてより狭い範囲に設定されたα（α5〜α6），β（β5〜β6），γ（γ5〜γ6），δ（δ5〜δ6）についてのステップ幅ｄａ，ｄｂ，ｄｇ，ｄｄを、以下の数７４のように定める（ステップＳ２３４１）。
【数７４】
2005233841.tif 000078
【０１３９】
次いで、ステップＳ２３４２においては、先に述べたＲＬ関係式である数２３，数２４、および数２，数３に基づく微係数（数５０，数５１）、加えてΔＲ（ω，ｐ）の演算処理等が行われる。「ΔＲ（ω，ｐ）」は、先にも述べた通り、導体表面インピーダンスＺs(ω)の実部と、等価回路インピーダンスＺin(ω,p)の実部との差の絶対値を示している。ここで係数α，β，γ，δは、図１１のステップＳ２３４２内に示すように、上記数７４と「ｉα」，「ｉβ」，「ｉγ」，「ｉδ」と係数α5，β5，γ5，δ5とを用いた関数「α＝α5＋（ｉα×ｄａ）」，「β＝β5＋（ｉβ×ｄｂ）」，「γ＝γ5＋（ｉγ×ｄｇ）」，「δ＝δ5＋（ｉδ×ｄｄ）」として与えられる。「ｉα」，「ｉβ」，「ｉγ」，「ｉδ」は、ここでは、「１〜５００」までの変数として与えられる。
【０１４０】
より具体的に説明すると、このステップＳ２３４２においては、ある範囲内にてあるステップで与えられる「α」，「β」，「γ」，「δ」ごとに、数２３，数２４を演算すると共に、各角周波数（Ｎｆ個の周波数）についての数５４等を演算する。つまり、等価回路インピーダンスの実部の角周波数に対する微係数「Ｓin(ω,p)」と、導体表面インピーダンスの実部の角周波数に対する微係数「Ｓs(ω)」との差の絶対値、およびΔＲ（ω，ｐ）を演算する。
【０１４１】
次いで、Ｎｆ個の周波数に対する「ΔＳ（ω，ｐ）」の総和が最小となる「α」，「β」，「γ」，「δ」の値を検出する（ステップＳ２３４３）。つまり、「ΔＳ（ω，ｐ）」の総和が最小となる「α」，「β」，「γ」，「δ」の値を得ることができる「ｉα」，「ｉβ」，「ｉγ」，「ｉδ」の値を検出する。このときの「ｉα」，「ｉβ」，「ｉγ」，「ｉδ」の値を「ｉαmin」，「ｉβmin」，「ｉγmin」，「ｉδmin」と定義する。
【０１４２】
次いで、検出した最小の箇所（例えば、「ｉαmin」にて得られる「α」の箇所）の付近について、より小さいステップにて係数α，β，γ，δを再度定義する（ステップＳ２３４４）。具体的には、以下の数７５〜数７８にて定義する。ここで再定義されたより小さいステップ範囲の係数α，β，γ，δは、後述の処理にて、必要とする精度の近似が得られなかった場合（後述するステップＳ２３４７にて「Ｎｏ」と判断された場合）に使用されることとなる。
【数７５】
2005233841.tif 000079
【数７６】
2005233841.tif 000080
【数７７】
2005233841.tif 000081
【数７８】
2005233841.tif 000082
【０１４３】
次いで、Ｎｆ個の周波数の中で「ΔＲ（ω，ｐ）」が最大となる周波数「ｆmax」の値を検出する（ステップＳ２３４５）。そして、この周波数「ｆmax」のときに得られるインピーダンスの差「ΔＲ（ω，ｐ）max」と、導体表面インピーダンスの実部「Ｚs(ω)max」との比「ｒ」を演算する（ステップＳ２３４６）。この「ｒ」の値が小さければ小さい程、等価回路インピーダンスＺin(ω,p)が導体表面インピーダンスＺs(ω)に近似していることとなる。ここで、「ΔＲ（ω，ｐ）が最大となる周波数」の箇所がピックアップされるのは、この箇所が最も近似状態が好ましくない部分であって、ここが後述する許容範囲内に存在することが確認されれば、他の箇所は当然の如く許容範囲内にあることは明らかだからである。
【０１４４】
次いで、演算された「ｒ」と、先に入力された許容範囲「ΔＲ」との対比判断処理が行われる（ステップＳ２３４７）。演算された「ｒ」が許容範囲「ΔＲ」よりも小さい場合（Ｓ２３４７にて「Ｙｅｓ」）には、「ｉｈ＝１」と定義されて（ステップＳ２３４８）、このサブルーチン５２の処理が終了して、図３のメインルーチンにリターンする。一方、演算された「ｒ」が許容範囲「ΔＲ」よりも大きい場合（Ｓ２３４７にて「Ｎｏ」）には、「α5＝α5m，α6＝α6m」、「β5＝β5m，β6＝β6m」、「γ5＝γ5m，γ6＝γ6m」、「δ5＝δ5m，δ6＝δ6m」、「ｉｈ＝０」と定義されて（ステップＳ２３４９）、このサブルーチン５２の処理が終了して、図３のメインルーチンにリターンする。
【０１４５】
サブルーチン５２の処理（Ｓ２３４）が終了した後、メインルーチンにおいては、「ｉｈ＝１」であるか否かが判断される（ステップＳ２３５）。つまり、サブルーチン５２の処理において、許容範囲（ΔＲ）を満たす結果が得られたか否かが判断される。
【０１４６】
許容範囲を満たす結果が得られた場合（Ｓ２３５にて「Ｙｅｓ」）には、得られた結果としての回路定数Ｒ1〜Ｒ5，Ｌ1〜Ｌ5の値が、コンピュータのディスプレイ、あるいはコンピュータに接続された出力装置等に出力され（ステップＳ２０６）、本実施形態に係る回路定数解析プログラムを用いた回路定数の解析処理が終了する。
【０１４７】
一方、許容範囲を満たす結果が得られなかった場合（Ｓ２３５にて「Ｎｏ」）には、図１１のステップＳ２３４９にて新たに設定された係数α5，α6、β5，β6、γ5，γ6、δ5，δ6を用いて再びサブルーチン５２の処理が行われることとなる。なお、ここで、サブルーチン５２が再び行われる場合には、２段接続等の場合と同様に、その前に「ｉＴ≧３」であるか否かが判断される（ステップＳ２３７）。「ｉＴ」は、ステップＳ２３３において「ｉＴ＝１」と定義されているから、サブルーチン５２が一度行われた後には、「１」のままである。よって、ここではステップＳ２３７で「Ｎｏ」と判断されて、「ｉＴ」の値が「１」加算されてから（ステップＳ２３８）、２回目のサブルーチン２２の処理が行われることとなる。
【０１４８】
つまり、本実施形態においては、２段接続等の場合と同様に、この５段接続の等価回路を用いた解析処理の結果が許容範囲内とならなかった場合であっても、サブルーチン５２の処理が３回まで繰り返し行われることとなる。そして、サブルーチン５２の処理を３回繰り返しても許容範囲内の結果を得られなければ（Ｓ２３７にて「Ｙｅｓ」）、「最大５段では許容範囲を満たす解析結果を得ることができませんでした。条件を変更して下さい。」等のメッセージを表示して（ステップＳ２３９）、一連の処理を終了すべく構成されている。
【０１４９】
本実施形態に係る回路定数解析プログラムは、５段接続の等価回路を用いた解析を３回繰り返しても、許容範囲（ΔＲ）内の結果が得られない場合には、上述したように、「メッセージ」を表示して、入力条件の変更を促すが、実用的な範囲においては、ほとんどの場合は、この回路定数解析プログラムの範囲内（すなわち５段接続までの範囲内）で、高精度近似を実現することが可能である。
【０１５０】
本実施形態に係る回路定数解析プログラムは、以上のように構成され機能するため、次のような効果を得ることができる。すなわち、１）本実施形態に係る回路定数解析プログラムによれば、単位等価回路の接続段数に応じて、それぞれ、未知数Ｒ，Ｌを数２３〜数３０のように定義しているため、それぞれの係数α，β，γ，δを調整することで、解析の対象として考えている周波数帯域内の導体表面インピーダンスを高精度に近似可能である。また、２）単位等価回路の接続段数ｐを増加させることで、さらに精度よく広帯域近似が可能である。さらに、３）本実施形態に係る回路定数解析プログラムによれば、単位等価回路の接続段数が少ないものから順番に演算を行い、要求された許容範囲を満足した段階で解析を終了すべく構成されているため、効率よく短時間で解析を行うことができる。
【０１５１】
＜導体表面解析プログラム＞
以下、本実施形態に係る導体表面解析プログラムについて、フローチャートを用いて具体的に説明する。
【０１５２】
図１２は本実施形態の導体表面解析プログラムを実行処理するコンピュータの概念図である。本実施形態の導体表面解析プログラムを実行処理するコンピュータは、キーボード、マウス等の入力手段５０、前記第１式を記憶する第１式記憶手段６１と前記第２式を記憶する第２式記憶手段６２とを含み、本実施形態の導体表面解析プログラム自身をも記憶するハードディスク等の記憶装置６０、特性インピーダンスＺ（Ｌ，ｉ）を演算する特性インピーダンス演算手段７１と導体表面解析手段７２とを含み、導体表面解析プログラムに基づいて種々の演算処理を実行するＣＰＵ７０、および演算結果を表示するディスプレイ、プリンタ等の出力手段８０を具備し、それぞれがバス線９０に接続されている。さらに、導体表面解析手段７２は、単位等価回路２０の端子間電圧Ｖin+1/2を演算する単位等価回路電圧演算手段７２１と導体表面の電界Ｅtを演算する電界演算手段７２２とを含んでいる（つまり、単位等価回路電圧演算手段７２１および電界演算手段７２２もＣＰＵ７０により実現される）。
【０１５３】
また、図示しないが、入力手段５０から入力された回路定数を入力データファイルとして保存または一時記憶したり、各演算手段において演算された演算結果をファイル保存するためのＲＡＭやハードディスク等、一般的なコンピュータの具備する構成は具備しているものとする（すなわち、本実施形態の導体表面解析プログラムは、一般的なパーソナルコンピュータ等のコンピュータで実行可能である）。さらに、入力手段５０は、ＦＤＴＤ法によって解析対象の電磁界の時間変化を時々刻々計算するために、導体表面解析に必要な位置に関する情報を設定する手段を兼ねるものであり、ＣＰＵ７０はＦＤＴＤ計算に必要な電界計算・電磁計算を行う手段を兼ねるものである。
【０１５４】
図１３は本実施形態の導体表面解析プログラムのメインフローチャートである。
【０１５５】
まず、入力手段５０により上記回路定数解析プログラム等で演算された回路定数Ｒ，Ｌを設定・入力する（ステップＳ１０１）。この回路定数の入力に際しては、伝送線路等価回路２の段数ｐを別途入力することとしてもよいし、回路定数Ｒ，Ｌをそれぞれ入力した数（ｐ個ずつ）を検出し、当該数をもって伝送線路等価回路２の段数ｐを判断（一時記憶）することとしてもよい。さらに、導体表面解析に必要な位置に関する情報（図１５において後述するような格子点の位置や導体の性質（例えば、銅、鉄等）等）を設定し（ステップＳ１０２）、プログラムをスタートさせる。
【０１５６】
プログラムのスタートを受けたＣＰＵ７０は、特性インピーダンス演算手段７１として機能し、記憶装置６０の第1式記憶手段６１に記憶されている上述の第1式（数３７）を読み出し、この第1式にステップＳ１１０において入力されたインダクタンスＬiを代入し、Ｚ（Ｌ，ｉ）を演算する（ステップＳ１０３）。この際、ｐ段分すべてのＬiについてＺ（Ｌ，ｉ）を演算しておき、一時記憶しておく。
【０１５７】
続いて、ＣＰＵ７０は、ＦＤＴＤ計算における演算ステップ数ｎの初期値としてｎ＝０を設定し（ステップＳ１０４）、時刻をＴ＝（ｎ＋１／２）Δｔだけ進める（ステップＳ１０５）。
【０１５８】
さらに、ＣＰＵ７０は、ステップＳ１０２において設定された導体表面解析の必要な位置以外の位置におけるＦＤＴＤ法に必要な電界を演算する（ステップＳ１０６）。具体的には、後述する図１５における導体壁面を除く空間部分の電界を演算する。このような演算はＦＤＴＤ法における既存の方法に従うものである。その上で、ステップＳ１０２において設定された導体表面解析の必要な位置における電界を演算する（ステップＳ１０７）。これについては、図１４を用いて後述する。その後、ステップＳ１０６，Ｓ１０７において演算された電界を用いて、ＦＤＴＤ法に必要な磁界を演算する（ステップＳ１０８）。この演算についてもＦＤＴＤ法における既存の方法に従うものである。なお、演算された磁界は、次回（ｎ＝ｎ＋１）において電界の演算を繰り返し行うためにも必要となる。
【０１５９】
ここで、ｎ＝ｎ＋１とし（ステップＳ１０９）、ＦＤＴＤ計算における演算ステップ数ｎが解析ステップ数ｎmaxより大きいか否かを判断する（ステップＳ１１０）。大きければ（ステップＳ１１０でＹｅｓ）、電磁界の演算を終了し、結果を出力する（ステップＳ１１１）。一方、演算ステップ数ｎが解析ステップ数ｎmaxより小さい場合は、演算ステップ数ｎが解析ステップ数ｎmaxを上回るまで、時刻ＴにΔｔを加算して、電磁界計算を順次繰り返す（ステップＳ１０５〜ステップＳ１１０：ＣＰＵ７０が導体表面解析手段７２として機能する）。
【０１６０】
ここで、ステップＳ１０７の電界計算について図１４を用いて、より詳しく説明する。図１４は導体表面解析プログラムのサブフローチャートである。サブプログラムのスタートを受けたＣＰＵ７０は、初期値として以下の数７９を設定する（ステップＳ１０７１）。
【数７９】
2005233841.tif 000083
【０１６１】
次に、等価回路２を流れる電流Ｉnの値として、時刻ｎにおける磁界Ｈnの値が設定される。この磁界については、図１５を用いて補足説明する。図１５は導体表面解析プログラムを適用する場合の一例としての導体表面付近の様子を示す図である。
【０１６２】
本例は、ステップＳ１０２において、図１５に示されるような導体表面を設定し、ｘｙ面の導体表面に平行なｘ方向電界Ｅxのうち、導体表面にある電界Ｅx2を導体表面プログラムによって求めるものである。この場合、Ｈnとして図中のｙ方向磁界Ｈyの時刻ｎにおける磁界Ｈynを電流Ｉnとして設定する。
【０１６３】
つまり、導体表面に平行な電界Ｅx2を導体表面解析プログラムで解析する場合に用いる磁界は、導体表面に平行で、求めたい電界Ｅxと対をなす電磁界の電力が導体の中へ入り込む成分を有するものであり、原理的には求めたい電界Ｅx2と同じ位置のものである。すなわち、用いる磁界成分は、一般に、求めたい電界成分と対をなし、ポインティング電力の向きが導体内方向の磁界成分である。ここで、一般的なＦＤＴＤ法の場合、磁界成分は、電界成分と同じ位置にはないので、図１５に示すように、求めたい電界Ｅx2の位置に最も近い磁界Ｈyを合理的に妥当な磁界値として用いるべきである。図１５に例示するように、用いる磁界の位置と電界の位置とは、ＦＤＴＤ法におけるセルサイズの半分だけ離れている。これは、他の磁界の位置に比べて単に近いというだけでなく、導体平面近くでは、導体面に平行な磁界成分における導体面に垂直な方向磁界変化がほぼ０となる性質を有するという、合理的理由があるからである（例えば、飯田幸雄、森田正信：“ＦＤＴＤ法による空洞解析でのＱ値の計算法と曲面導体壁の取扱い”，電子情報通信学会論文誌(C-I)，vol.J77-c-I,no.9,pp.496-503,Sept.1994参照）。
【０１６４】
以上のような電流Inの設定の後、ＣＰＵ７０が単位等価回路電圧演算手段７２１として機能し、記憶装置６０の第２式記憶手段６２に記憶されている上述の第２式（数４７および数４８）を読み出し、これらの第２式にステップＳ１０１において入力されたレジスタンスＲ1、ステップＳ１０３において演算された特性インピーダンスＺ（Ｌ，ｉ）をそれぞれ代入し、単位等価回路２０の端子間電圧Ｖ1n+1/2および単位等価回路２０に流れる電流Ｉ1nを演算する（ステップＳ１０７３）。演算結果となるＶ1n+1/2は、数４７に示されるように、それ以前の時刻（ｎ−１／２）Δｔ，（ｎ−３／２）Δｔ，（ｎ−１）Δｔ，ｎΔｔにおける単位等価回路２０の端子間電圧Ｖ1n-1/2，Ｖ1n-3/2，同じく電流Ｉ1n-1およびステップＳ１０７２で設定された等価回路２を流れる電流Ｉn＝Ｈnにより表される。したがって、時刻（ｎ＋１／２）Δｔにおいては、これらの電圧値および電流値は既知となるため、Ｖ1n+1/2も演算可能となる。また、Ｉ1nは、数４８に示されるように、既知となっているＶ1n-1/2，Ｉnおよび数４７で演算されたＶ1n+1/2により表される。したがって、Ｉ1nも演算可能となる。
【０１６５】
上述の通り、数４７および数４８で表される第２式は、微分方程式を含まない数式であるため、本導体表面解析プログラムをＣＰＵ７０等に過大な負荷をかけることなく実行することができる。
【０１６６】
続いて、ＣＰＵ７０が電界演算手段７２２として機能し、ステップＳ１０７３において演算されたＶ1n+1/2をもとのＶn+1/2に加算して新たなＶn+1/2とする（ステップＳ１０７４）。すなわち、以下の数８０のごとく演算する。
【数８０】
2005233841.tif 000084
【０１６７】
今、ステップＳ１０７１より、もとのＶn+1/2は、初期値０であるので、これにステップＳ１０７３で演算されたＶ1n+1/2を加算して、新たなＶn+1/2＝Ｖ1n+1/2が得られることとなる。
【０１６８】
プログラムは、ｉ＝ｉ＋１を指令し（ステップＳ１０７５）、伝送線路等価回路２の段数ｐを超えるか否かを判断する（ステップＳ１０７６）。そして、ステップＳ１０７６においてｉが伝送線路等価回路２の有する単位等価回路２０の段数ｐを超えるまで（例えば、ｐ＝５の場合、ステップＳ１０７６においてｉ＝６＞ｐとなるまで）ステップＳ１０７３〜Ｓ１０７６を繰り返す（ステップＳ１０７６においてＮｏ）。
【０１６９】
ステップＳ１０７３〜Ｓ１０７６を所定回数（ｐ回）繰り返す間、ステップＳ１０７４において、Ｖn+1/2は、以下の数８１に示すように、その回におけるステップＳ１０７４の演算結果Ｖin-1/2が繰り返し足し合わされることとなる。
【数８１】
2005233841.tif 000085
【０１７０】
ステップＳ１０７５の結果、ｉ＞ｐとなった場合（ステップＳ１０７６においてＹｅｓ）、上記数８１の結果を以下の数８２に示すように、時刻（ｎ＋１／２）Δｔにおける電界Ｅn+1/2の値を得る。この値は、図１３のステップＳ１０７における磁界計算に用いられることとなる。
【数８２】
2005233841.tif 000086
【０１７１】
以上より、本実施形態に係る導体表面解析プログラムは、以下の効果を奏することとなる。すなわち、１）演算結果として出力される電界Ｅn+1/2は、ＦＤＴＤのセルサイズおよび時間ステップΔｔが含まれていないため、本実施形態の導体表面解析を用いることによる解析精度は、上記数３６で示され、上述の通り十分に小さい近似誤差α（ω）のみに影響されることとなり、非常に高い解析精度を達成することができる。また、２）上記数４７および数４８で表される第２式は、微分方程式を含まない数式であるため、本導体表面解析プログラムをＣＰＵ７０等に過大な負荷をかけることなく実行することができる。
【０１７２】
なお、上記の回路定数解析プログラムと導体表面解析プログラムとを組み合わせて１つのプログラムとすることも可能である。もちろん、本発明の導体表面解析プログラムは、回路定数解析プログラムとは独立して実行可能なプログラムであり、ＲＬ定数化にあたっては、上記回路定数解析プログラムと異なる方法により演算することとしてもよい。
【実施例１】
【０１７３】
以下に、本発明のＦＤＴＤ定式化の前提となった伝送線路等価回路２の入力インピーダンス（Ｚin(ω)）と導体表面インピーダンス（ＺS(ω)）との比較を実施例として示す。なお、以下に本発明の実施例を示すが、本発明はもとより下記実施例によって制限を受けるものではなく、前後記の趣旨に適合し得る範囲で適当に変更を加えて実施することも勿論可能であり、それらはいずれも本発明の技術的範囲に含まれる。
【０１７４】
ここでは、単位等価回路の接続段数ｐを５段とし、導体が「銅」である２次元ＦＤＴＤモデルを用いて例示する。つまり、銅の導体表面インピーダンスを近似する。
【０１７５】
上記実施例において、銅の導電率σ＝５．８×１０7Ｓ／ｍ、比誘電率１．０、比透磁率１．０とする。
【０１７６】
まず、ｐ＝５における数３のＺ’in(ω)で銅の表面インピーダンスを近似したときの回路定数Ｒi，Ｌiを以下の表１に示す。
【表１】
2005233841.tif 000087
【０１７７】
上記回路定数は、ω1／２π＝３ＧＨｚ，ω2／２π＝７ＧＨｚ，ω3／２π＝５ＧＨｚの条件で決定した。銅の表面インピーダンスＺS(ω)と表１の回路定数を用いて図２の回路を構成したときの入力インピーダンスＺin(ω)との比較結果を図１６に示す。図１６においては、リアクタンス成分の近似状態を示している。図１６から明らかなように、広帯域に渡って極めて高精度な近似結果が得られた。
【実施例２】
【０１７８】
続いて、以下に、間隔を１０ｍｍに保った２枚の平行導体板間をＴＥＭ波の電磁波が１ｍの長さを伝搬するモデルについて、本発明導体表面解析プログラムを用いたＦＤＴＤ法によって解析した結果を実施例として示す。本実施例においては、導体板として、銅を想定し、上記実施例１で得られた表１の抵抗値とインダクタンス値を用いて伝送損失の周波数特性を調べた。比較例として、畳み込み積分による計算も試みた。この畳み込み積分は、前述の非特許文献８で提案された導体境界処理（近似の影響があるためにProny近似をせず、過去のデータに遡って畳み込み積分を行う方法であり、非特許文献９の式（３），（４）により計算する）を用いた方法であり、導体表面での電磁界の過去の値のすべてを記憶して計算するため、精度が高いと考えられるものである。
【０１７９】
図１７は本発明の実施例２における伝搬損失の周波数特性を示すグラフである。実線が、本実施例の結果を示すものであり、破線が、比較例の結果を示すものである。図１７に示すように、破線と実線とは、１．０ＧＨｚ以下の周波数範囲において少しのずれが見られる他はほとんど重なって区別が困難になっており、両者がよく一致していることが分かる。すなわち、本実施例は、伝搬損失の周波数特性が精度よく計算可能であることが示された。
【０１８０】
したがって、本実施形態に係る導体表面解析プログラムを用いれば、広帯域にわたって近似精度を高く維持可能であることが示された。
【０１８１】
なお、本発明は上記実施形態および実施例に限定されるものではなく、その趣旨を逸脱しない限りにおいて、上述したもの以外に種々の変更を行うことが可能である。
【図面の簡単な説明】
【０１８２】
【図１】本発明の実施形態において導体表面インピーダンス（Ｚs(ω)）を近似するＲＬ等価回路の概略図である。
【図２】本発明の実施形態における導体表面上の電界Ｅtを演算するための伝送線路等価回路の概略図である。
【図３】本実施形態に係る回路定数解析プログラムのメインルーチンのフローチャートである。
【図４】本実施形態に係る回路定数解析プログラムを構成するサブルーチン２１のフローチャートである。
【図５】本実施形態に係る回路定数解析プログラムを構成するサブルーチン２２のフローチャートである。
【図６】本実施形態に係る回路定数解析プログラムを構成するサブルーチン３１のフローチャートである。
【図７】本実施形態に係る回路定数解析プログラムを構成するサブルーチン３２のフローチャートである。
【図８】本実施形態に係る回路定数解析プログラムを構成するサブルーチン４１のフローチャートである。
【図９】本実施形態に係る回路定数解析プログラムを構成するサブルーチン４２のフローチャートである。
【図１０】本実施形態に係る回路定数解析プログラムを構成するサブルーチン５１のフローチャートである。
【図１１】本実施形態に係る回路定数解析プログラムを構成するサブルーチン５２のフローチャートである。
【図１２】本実施形態に係る導体表面解析プログラムを実行処理するコンピュータの概念図である。
【図１３】本実施形態に係る導体表面解析プログラムのメインフローチャートである。
【図１４】本実施形態に係る導体表面解析プログラムのサブフローチャートである。
【図１５】導体表面解析プログラムを適用する場合の一例としての導体表面付近の様子を示す図である。
【図１６】本発明の実施例１についての近似結果を示すグラフである。
【図１７】本発明の実施例２における伝搬損失の周波数特性を示すグラフである。
【符号の説明】
【０１８３】
１    ＲＬ等価回路
２    伝送線路等価回路
１０    ＲＬ単位等価回路
１１    抵抗
１２    コイル
１３    短絡伝送線路
２０    単位等価回路
【図１】
2005233841.tif 000090
【図２】
2005233841.tif 000091
【図３】
2005233841.tif 000092
【図４】
2005233841.tif 000093
【図５】
2005233841.tif 000094
【図６】
2005233841.tif 000095
【図７】
2005233841.tif 000096
【図８】
2005233841.tif 000097
【図９】
2005233841.tif 000098
【図１０】
2005233841.tif 000099
【図１１】
2005233841.tif 000100
【図１２】
2005233841.tif 000101
【図１３】
2005233841.tif 000102
【図１４】
2005233841.tif 000103
【図１５】
2005233841.tif 000104
【図１６】
2005233841.tif 000105
【図１７】
2005233841.tif 000106