ソフトウェア/OpenFDTD

(519d) 更新


公開メモ

OpenFDTD を使ってプリント基板上の高速信号の伝播を理解したい

例えば4層基板上のマイクロストリップラインの周囲の電場をこんな風に計算するとか。

microstrip1.png microstrip2.png

目次

OpenFDTD について

株式会社 EEM により開発&メンテナンスされていたオープンソースの FDTD シミュレーションソフト。

Wikipedia:FDTD法 は主にマクスウェル方程式を数値的に解いて電磁波のシミュレーションを行う際に使われる数値計算手法の一種。

株式会社 EEM は2022年9月30日に解散したため、現在は大賀さんが個人的に (?) メンテナンスを続けてらっしゃるようです。

OpenFDTD のダウンロード&説明

以前は http://e-em.co.jp/OpenFDTD にあったのですが、

現在は http://emoss.starfree.jp/ で公開されています。

1.1 概要

(1) 本ソフトウェアについて

本ソフトウェアOpenFDTDはFDTD法(Finite Difference Time Domain method)[1][2] を用いた汎用的な電磁界シミュレーターです。 電磁界の広い用途の解析に使用することができます。

(2) オープンソース

本ソフトウェアOpenFDTDはオープンソースのフリーソフトです。 自由にダウンロードしてお使いください。

(3) 動作環境

動作環境はWindowsまたはLinuxです。 Windows環境には簡易GUIと実行プログラムが付属しています。

GPU やクラスターマシン、スーパーコンピュータ用の高速演算コードまで含まれているすばらしいものです。

Windows 用には GUI も提供されていて、モデルの作成や計算結果のグラフ表示が可能です。

計算例

プリント基板 (PCB) のシミュレーションに使うことを想定した注意点

プリント基板 (PCB) のシミュレーションでは各部のサイズを厳密に設定しなければならならず、セルの大きさ分の誤差が計算結果に影響する。

このため、FDTD 計算の基本を押さえておかないと計算結果が合わず困ったことになります。

FDTD 計算の基本

FDTD では計算したい空間を「直方体のセル」に分割して計算する。

ここではセルの数を $N_x\times N_y\times N_z$ とする。

「セルの頂点」が格子点 $(x_i,y_j,z_k)$ にあたる。頂点数は $(N_x+1)\times (N_y+1)\times (N_z+1)$ である。

電場・電流はセルの辺に沿って計算され、その際には「辺の中心」が基準点となる。 たとえば $E_x$ を計算するには $(x_i+(x_{i+1}-x_i)/2,y_j,z_k)=((x_i+x_{i+1})/2,y_j,z_k)$ が基準点で、その数は $N_x\times (N_y+1)\times (N_z+1)$ である。

基準点の周りに「その電場が仮定される範囲」が取られる。

  • $x_i<x<x_{i+1}$
  • $y_j-(y_j-y_{j-1})/2<y<y_j+(y_{j+1}-y_j)/2$
  • $z_k-(z_k-z_{k-1})/2<z<z_k+(z_{k+1}-z_k)/2$

がその範囲になるため、セルの範囲、$E_{x,i+1/2,j,k}$ の範囲、$E_{y,i,j+1/2,k}$ の範囲、$E_{z,i,j,k+1/2}$ の範囲はすべて異なることに注意が必要。さらに、隣り合うセルのサイズが異なる場合には「電場セル」の中心は「辺の中心」とも一致しない。

磁場・磁流はセルの面に垂直に計算され、面の中心が基準点となる。例えば $H_{z}$ は $(x_{i+1/2},y_{j+1/2},z_k)$ が基準点となる。

マテリアルの範囲

誘電体や金属を配置する際、その範囲は「セル単位」ではなく「座標空間内に厳密に設定された範囲」になる。

どういうことかというと、電場や磁場を計算する際、その「基準点」の座標がどこかによって、どのマテリアルのパラメータを使って計算するかが決まる。

FDTD が基礎とする Maxwell 方程式の形は

$$ \frac{\epsilon_r}{c}\partial_t\bm E=\bm\nabla\times\eta\bm H-\underbrace{(\eta\sigma)\bm E}_\text{電流密度項} $$ $$ \frac{\mu_r}c\partial_t\eta\bm H=-\bm\nabla\times\bm E-\underbrace{(\frac{\sigma_m}\eta)\eta\bm H}_\text{磁流密度項} $$

で、特に材質によって大きく変わる誘電率 $\epsilon_r$ と伝導率 $\sigma$ は電場の方程式だけに含まれるので、「格子点」や「セルの中心」ではなくむしろ「『辺の中心位置』にどの材質が設定されているか」が結果に影響してくるという点に注意が必要なのだ。

FDTD における導体は「ジャングルジム」のようなもので、その辺にしか電流は流れない。 物体をおかしなところで切ると、片方の先がどこにもつながっていないような中途半端な辺がたくさん 生えた「毛深い」物体ができてしまう。「毛」はアンテナっぽく働くけれど、毛と毛の間で横方向に 電流が流れないので思ったのと違う動作をしかねない。

また「辺」は「仮想的な太さ」を持つので、全体としての伝導度や抵抗を設計値と合わせたい場合には その「太さ」を正しく考慮しないと値が異なってしまう。

また OpenFDTD では、

  • 設定したマテリアルの範囲の「境界上の点」は「マテリアルの内部」として考慮されること、
  • 複数のマテリアルに含まれる点は「最後に指定したマテリアル」で上書きされること、

にも注意が必要だ。

具体例で考える

例えば厚さ 1.6 mm の両面基板を考え、裏のグランドプレーンを $z=-1.6\ \text{mm}$ に配置したいとする。

そこである $k$ に対して $z_k=-1.6\ \text{mm}$ となるよう格子点を配置して、$z=-1.6\ \text{mm}$ に「厚さゼロ」の完全導体を配置する。

マテリアルを設定した領域の境界と一致する点にはそのマテリアルが適用されるので、これにより、この面内にメッシュ状に配置された長方形の「辺」が完全導体の特性を持つようになり、この面内で電流が自由に流れるようになる。

辺は「太さ」を持つため、ここに完全導体ではなく有限伝導度の金属を置くなら、その厚さにも注意が必要だ。 この時、導体の厚さは次の格子面の中点までと定義されるから、その値は $(z_{k+1}-z_{k-1})/2$ となる。この面が系の端にある場合、例えば $z_0$ にある場合は、$z_{-1}$ の格子面が存在しないため、代わりに反対側の格子面までの距離が使われて、$z_1-z_0$ がトータルの厚さになる。これに伝導率を掛けてシート抵抗が得られる。

この電流による磁場は、本来太さを持つ辺を流れる電流が、あたかも辺の上に集中しているかのように考えて計算される。

同様に、パターン面を $z=0\ \text{mm}$ として、そこに厚さゼロの金属を配置すればそこに電流を流せる。

この間を誘電体、例えば FR-4 で埋めるには、グランドプレーンやパターン面の金属を配置する「前に」 $\epsilon_r=4.3$ の材質を $-1.6\ \text{mm}\le z\le 0\ \text{mm}$ に配置すると良い。後で配置すると先に置いた金属が上書きされてしまう。

グランドプレーンやパターンには「辺の太さ」分だけの仮想的な厚みがあるが、誘電体内で $E_z$ を計算する際にはその金属板の厚みは考慮されない。例えば $z_{k+1/2}$ にある $z$ 方向に向く辺は正しく $z_k$ から $z_{k+1}$ までの長さを持つ。つまりグランドプレーンの仮想的な厚さを無視して $z_k$ まで突き刺さっていて、その計算に誘電体の誘電率と伝導率が使われるのだ。

だから、$z_k<z<z_{k+1/2}$ の範囲は面内の電場を計算するときには金属、面直の電場を計算するときには誘電体であるように扱われる不思議領域になる。そこが問題になるなら、金属面のすぐ隣の面までの面間隔を小さくすることで影響を最小化できるはず。

パターンの太さを考えるときも十分注意しなければならない。格子点から格子点までの長さをピッタリパターンの太さに合わせ、両端の格子点を含む形でパターンに対応する金属領域を設定するのが、金属の形状を定義してそのまわりの誘電体内の電場を計算するには最適と思われるが、そのときパターンの幅は辺の仮想的な太さ分だけ設定領域をはみ出すので、パターンに抵抗率を設定する場合にはその分だけ抵抗が下がってしまう。

影響を減らすにはやはり、パターンの外周のすぐ外の格子面までの距離を十分に小さくとるしかないのだと思う。とはいえ例えばパターンの厚さは 15 um とかしかないので、これを正しく計算できるほどメッシュを細かくするのは現実的ではなく、何かうまくごまかしながら計算していくしかない感じ???

結果表示用 Igor Pro スクリプト

計算結果を図にするのが大変なので、Igor Pro で直接読み込んで表示できるようにした。

いろいろ手抜きをしているので、今後も改善していかないと。

FDTD 計算結果は ofd.out というファイルにごそっと書き出されている。

グラフ化などのポストプロセスはこのデータを読み込んで、画像に加工している。

なので、ofd.out を Igor で読み込めば、好きなように表示可能!

  • 電場を xy, yz, zx 面で切って表示するところまで書いた。
  • Igor6 でも動くように少し直した。int64 に対応していないので、int32 で2回読んでつなぎ合わせる必要があった。あと、wave1[inf] のようなのがエラーになるので wave1[numpnts(wave1)-1] のように書き換えた。
  • NPoint = 0 でおかしくならないよう書き換えた。
  • Zin と Reflection を表示できるようにした。
  • cEpx などの計算が間違ってた(cHpx は合ってるのか?)
  • NFreq2 = 0 でも大丈夫にしようとしたのだけれど、これはゼロにはならないのかも。
  • スミスチャートを表示できるようにした
#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3		// Use modern global access method and strict wave access.
#include "SmithChart"

// load_ofd_out("C:Users:osamu:documents:OpenFDTD:ofd.out")
// のようにして読み込む。

Structure ofd_out_header
	char	title[256]
	int32	Nx, Ny, Nz
	int32	NExl, NExh, NEyl, NEyh, NEzl, NEzh
	int32	NHxl, NHxh, NHyl, NHyh, NHzl, NHzh
	int32	NFreq1, NFreq2, NFeed, NPoint, Niter, Ntime
	int32	SolverMaxIter, SolverNOut
	double	Dt
	int32	NGline
EndStructure

Structure odf_feed
	char	dir, dummy1, dummy2, dummy3 // 4 byte alignment
	int32	i, j, k
	double	dx, dy, dz
	double	volt
	double	delay
	double	z0
EndStructure

function/t odf_feed2str(feed)
STRUCT odf_feed &feed
	string s = ""
	sprintf s, "dir:%d;i:%d;j:%d;k:%d", feed.dir, feed.i, feed.j, feed.k
	sprintf s, "%s;dx:%g;dy:%g;dz:%g", s, feed.dx, feed.dy, feed.dz
	sprintf s, "%s;volt:%g;delay:%g;z0:%g", s, feed.volt, feed.delay, feed.z0
	return s
end

function load_ofd_out(path)
string path
	variable f
	open/r f as path

	STRUCT ofd_out_header header
	fbinread f, header
	string/g title = header.title
	variable/g Nx = header.Nx, Ny = header.Ny, Nz = header.Nz
	variable/g NEx = header.NExl + header.NExh * 2^32
	variable/g NEy = header.NEyl + header.NEyh * 2^32
	variable/g NEz = header.NEzl + header.NEzh * 2^32
	variable/g NHx = header.NHxl + header.NHxh * 2^32
	variable/g NHy = header.NHyl + header.NHyh * 2^32
	variable/g NHz = header.NHzl + header.NHzh * 2^32
	variable/g NFreq1 = header.NFreq1, NFreq2 = header.NFreq2, NFeed = header.NFeed
	variable/g NPoint = header.NPoint, Niter = header.Niter, Ntime = header.Ntime
	variable/g SolverMaxIter = header.SolverMaxIter, SolverNOut = header.SolverNOut
	variable/g Dt = header.Dt, NGline = header.NGline
	
	if(nex!=nx*(ny+1)*(nz+1))
		abort "unexpected NEx"
	endif
	if(ney!=(nx+1)*ny*(nz+1))
		abort "unexpected NEx"
	endif
	if(nez!=(nx+1)*(ny+1)*nz)
		abort "unexpected NEx"
	endif
	if(nhx!=(nx+1)*(ny+2)*(nz+2))
		abort "unexpected NHx"
	endif	
	if(nhy!=(nx+2)*(ny+1)*(nz+2))
		abort "unexpected NHx"
	endif	
	if(nhz!=(nx+2)*(ny+2)*(nz+1))
		abort "unexpected NHx"
	endif	
	
	make/O/N=(Nx+1) Xn;	fbinread /B=3 /F=5 f, Xn
	make/O/N=(Ny+1) Yn;	fbinread /B=3 /F=5 f, Yn
	make/O/N=(Nz+1) Zn;	fbinread /B=3 /F=5 f, Zn
	make/O/N=(Nx) Xc;	fbinread /B=3 /F=5 f, Xc
	make/O/N=(Ny) Yc;	fbinread /B=3 /F=5 f, Yc
	make/O/N=(Nz) Zc;	fbinread /B=3 /F=5 f, Zc
	make/O/N=(Niter) Eiter;	fbinread /B=3 /F=5 f, Eiter
	make/O/N=(Niter) Hiter;	fbinread /B=3 /F=5 f, Hiter
	make/O/N=(SolverMaxIter+1, NFeed) VFeed
	make/O/N=(SolverMaxIter+1, NFeed) IFeed
	if (NFeed > 0) 
		fbinread /B=3 /F=5 f, VFeed
		fbinread /B=3 /F=5 f, IFeed
	endif
	make/O/N=(SolverMaxIter+1, NPoint) VPoint; 
	if (NPoint > 0)
		fbinread /B=3 /F=5 f, VPoint
	endif
	make/O/N=(NFreq1) Freq1; fbinread /B=3 /F=5 f, Freq1
	make/O/N=(NFreq2) Freq2; 
	if (NFreq2>0)
		fbinread /B=3 /F=5 f, Freq2
	endif
	setscale d,0,1,"Hz",freq1,freq2
	
	variable i
	STRUCT odf_feed feedn
	make/O/N=(NFeed)/T Feed
	for(i=0; i<NFeed; i+=1)
		fbinread/B=3 f, feedn
		Feed[i] = odf_feed2str(feedn)
	endfor

	make/O/N=(NFeed, NFreq1)/C Zin
	make/O/N=(NFeed, NFreq1) Ref
	make/O/N=(NFeed, NFreq2) Pin0
	make/O/N=(NFeed, NFreq2) Pin1
	setscale d,0,1,"Ohm",Zin
	if (NFeed > 0)
		fbinread /B=3 /F=5 f, Zin
		fbinread /B=3 /F=5 f, Ref
		if (NFreq2>0)
			fbinread /B=3 /F=5 f, Pin0
			fbinread /B=3 /F=5 f, Pin1
		endif
	endif
	make/O/N=(NFreq1, NPoint)/C Spara
	if (NPoint > 0)
		fbinread /B=3 /F=5 f, Spara
	endif
	make/O/N=(3, 2, NGline) Gline
	if (NGline > 0)
		fbinread /B=3 /F=5 f, Gline
	endif

// CExyz, CHxyz のメモリ構造 from ofd.h
//
//#define NEX(i,j,k) ((i)*(Ny+1)*(Nz+1)+(j)*(Nz+1)+(k))
//#define NEY(i,j,k) ((j)*(Nz+1)*(Nx+1)+(k)*(Nx+1)+(i))
//#define NEZ(i,j,k) ((k)*(Nx+1)*(Ny+1)+(i)*(Ny+1)+(j))
//#define NHX(i,j,k) ((i)*(Ny+2)*(Nz+2)+(j+1)*(Nz+2)+(k+1))
//#define NHY(i,j,k) ((j)*(Nz+2)*(Nx+2)+(k+1)*(Nx+2)+(i+1))
//#define NHZ(i,j,k) ((k)*(Nx+2)*(Ny+2)+(i+1)*(Ny+2)+(j+1))
//
//#define CEX_r(f,i,j,k) cEx_r[(f*NEx)+NEX(i,j,k)]
//#define CEY_r(f,i,j,k) cEy_r[(f*NEy)+NEY(i,j,k)]
//#define CEZ_r(f,i,j,k) cEz_r[(f*NEz)+NEZ(i,j,k)]
//#define CHX_r(f,i,j,k) cHx_r[(f*NHx)+NHX(i,j,k)]
//#define CHY_r(f,i,j,k) cHy_r[(f*NHy)+NHY(i,j,k)]
//#define CHZ_r(f,i,j,k) cHz_r[(f*NHz)+NHZ(i,j,k)]
//#define CEX_i(f,i,j,k) cEx_i[(f*NEx)+NEX(i,j,k)]
//#define CEY_i(f,i,j,k) cEy_i[(f*NEy)+NEY(i,j,k)]
//#define CEZ_i(f,i,j,k) cEz_i[(f*NEz)+NEZ(i,j,k)]
//#define CHX_i(f,i,j,k) cHx_i[(f*NHx)+NHX(i,j,k)]
//#define CHY_i(f,i,j,k) cHy_i[(f*NHy)+NHY(i,j,k)]
//#define CHZ_i(f,i,j,k) cHz_i[(f*NHz)+NHZ(i,j,k)]

	if (NFreq2 > 0)
		make/O/N=(nx, ny, nz, NFreq2)/C cEx, cEy, cEz
	
		make/O/N=(nz+1, ny+1, nx+0, NFreq2) cEx_r; fbinread /B=3 /F=4 f, cEx_r
		make/O/N=(nz+1, ny+1, nx+0, NFreq2) cEx_i; fbinread /B=3 /F=4 f, cEx_i
		cEx = (cmplx(cEx_r[r][q][p][s], cEx_i[r][q][p][s])+cmplx(cEx_r[r+1][q][p][s], cEx_i[r+1][q][p][s])+cmplx(cEx_r[r][q+1][p][s], cEx_i[r][q+1][p][s])+cmplx(cEx_r[r+1][q+1][p][s], cEx_i[r+1][q+1][p][s]))/4
		killwaves cEx_r, cEx_i
		
		make/O/N=(nx+1, nz+1, ny+0, NFreq2) cEy_r; fbinread /B=3 /F=4 f, cEy_r
		make/O/N=(nx+1, nz+1, ny+0, NFreq2) cEy_i; fbinread /B=3 /F=4 f, cEy_i
		cEy = (cmplx(cEy_r[p][r][q][s], cEy_i[p][r][q][s])+cmplx(cEy_r[p+1][r][q][s], cEy_i[p+1][r][q][s])+cmplx(cEy_r[p][r+1][q][s], cEy_i[p][r+1][q][s])+cmplx(cEy_r[p+1][r+1][q][s], cEy_i[p+1][r+1][q][s]))/4
		killwaves cEy_r, cEy_i
	
		make/O/N=(ny+1, nx+1, nz+0, NFreq2) cEz_r; fbinread /B=3 /F=4 f, cEz_r
		make/O/N=(ny+1, nx+1, nz+0, NFreq2) cEz_i; fbinread /B=3 /F=4 f, cEz_i
		cEz = (cmplx(cEz_r[q][p][r][s], cEz_i[q][p][r][s])+cmplx(cEz_r[q+1][p][r][s], cEz_i[q+1][p][r][s])+cmplx(cEz_r[q][p+1][r][s], cEz_i[q][p+1][r][s])+cmplx(cEz_r[q+1][p+1][r][s], cEz_i[q+1][p+1][r][s]))/4
		killwaves cEz_r, cEz_i
	
		make/O/N=(nx+1, ny+1, nz+1, NFreq2)/C cHx, cHy, cHz
	
		make/O/N=(nz+2, ny+2, nx+1, NFreq2) cHx_r; fbinread /B=3 /F=4 f, cHx_r
		make/O/N=(nz+2, ny+2, nx+1, NFreq2) cHx_i; fbinread /B=3 /F=4 f, cHx_i
		cHx = (cmplx(cHx_r[r][q][p][s], cHx_i[r][q][p][s])+cmplx(cHx_r[r+1][q][p][s], cHx_i[r+1][q][p][s])+cmplx(cHx_r[r][q+1][p][s], cHx_i[r][q+1][p][s])+cmplx(cHx_r[r+1][q+1][p][s], cHx_i[r+1][q+1][p][s]))/4
		killwaves cHx_r, cHx_i
	
		make/O/N=(nx+2, nz+2, ny+1, NFreq2) cHy_r; fbinread /B=3 /F=4 f, cHy_r
		make/O/N=(nx+2, nz+2, ny+1, NFreq2) cHy_i; fbinread /B=3 /F=4 f, cHy_i
		cHy = (cmplx(cHy_r[p][r][q][s], cHy_i[p][r][q][s])+cmplx(cHy_r[p+1][r][q][s], cHy_i[p+1][r][q][s])+cmplx(cHy_r[p][r+1][q][s], cHy_i[p][r+1][q][s])+cmplx(cHy_r[p+1][r+1][q][s], cHy_i[p+1][r+1][q][s]))/4
		killwaves cHy_r, cHy_i
	
		make/O/N=(ny+2, nx+2, nz+1, NFreq2) cHz_r; fbinread /B=3 /F=4 f, cHz_r
		make/O/N=(ny+2, nx+2, nz+1, NFreq2) cHz_i; fbinread /B=3 /F=4 f, cHz_i
		cHz = (cmplx(cHz_r[q][p][r][s], cHz_i[q][p][r][s])+cmplx(cHz_r[q+1][p][r][s], cHz_i[q+1][p][r][s])+cmplx(cHz_r[q][p+1][r][s], cHz_i[q][p+1][r][s])+cmplx(cHz_r[q+1][p+1][r][s], cHz_i[q+1][p+1][r][s]))/4
		killwaves cHz_r, cHz_i
	endif
	
	variable/g num
	fbinread /B=3 /F=5 f, num
	
	make/O/N=(nx+1) cEpx;	cEpx = xn
	make/O/N=(ny+1) cEpy;	cEpy = yn
	make/O/N=(nz+1) cEpz;	cEpz = zn
	setscale d,0,1,"m", cEpx, cEpy, cEpz
	
	make/O/N=(nx+2) cHpx;	cHpx = p == 0 ? xn[0]-(xn[1]-Xn[0])/2 : (p < nx ? (xn[p-1]+xn[p])/2 : xn[p-1] + (xn[p-1]-xn[p-2])/2)
	make/O/N=(ny+2) cHpy;	cHpy = p == 0 ? yn[0]-(yn[1]-yn[0])/2 : (p < ny ? (yn[p-1]+yn[p])/2 : yn[p-1] + (yn[p-1]-yn[p-2])/2)
	make/O/N=(nz+2) cHpz;	cHpz = p == 0 ? zn[0]-(zn[1]-zn[0])/2 : (p < nz ? (zn[p-1]+zn[p])/2 : zn[p-1] + (zn[p-1]-zn[p-2])/2)
	setscale d,0,1,"m", cHpx, cHpy, cHpz
			
	FStatus f
	if(V_logEOF != V_filePos)
		print "unexpected file length"
		printf "logEOF = %f, filepos = %f\n", V_logEOF, V_filePos
	endif

	close f
	
	show_ofd()
end

function show_ofd_sub(c1, c2, c3)
string c1, c2, c3
	wave cEpx, cEpy, cEpz
	wave/C cEx, cEy, cEz
	nvar nx, ny, nz

	make/O/N=(nx,ny) exy
	make/O/N=(ny,nz) eyz
	make/O/N=(nz,nx) ezx
	
	NVAR NFreq2
	wave Freq2

	wave cEp1 = $("cEp"+c1)
	wave cEp2 = $("cEp"+c2)
	wave cEp3 = $("cEp"+c3)
	string g="graph_e"+c1+c2+"_"+c3
	dowindow $g
	if(!v_flag)
		variable/g $("Ep"+c3)=round(100*1000*(cEp3[0]+cEp3[inf])/2)/100
		variable/g Ef = Freq2[0]/1e6
		variable/g Efnowi = 0
		variable/g Efnow = Freq2[0]/1e6
		display
		appendimage $("e"+c1+c2) vs {cEp1, cEp2}
		dowindow/c $g
		Label/w=$g bottom "Position "+c1+" / \\U"
		Label/w=$g left "Position "+c2+" / \\U"
		ModifyImage/w=$g $("e"+c1+c2) ctab= {0.2,*,Rainbow,1}, log=1
		ControlBar/w=$g 20
		CheckBox check0 title="log",proc=CheckLogProc,win=$g,value=1,side=1
		SetVariable $("setFreq"+c3) title="freq (MHz):",win=$g,proc=SetSliceCoordinateProc,limits={(Freq2[0]/1e6),(Freq2[inf]/1e6),1},live=1,size={220,18},value=Ef
		SetVariable $("setEp"+c3) title=(c3+" (mm):"),win=$g,proc=SetSliceCoordinateProc,limits={(1000*cEp3[0]),(1000*cEp3[inf]),0.01},live=1,size={200,18},value=$("Ep"+c3)
		
		make/O/N=2 $("e"+c1+c2+"_cur1")
		wave e12_cur1 = $("e"+c1+c2+"_cur1")
		setscale/p x,(cEp1[0]),(cEp1[inf]-cEp1[0]),waveunits(cEp1,-1),e12_cur1
		setscale d,0,1,waveunits(cEp2,-1),e12_cur1
		appendtograph/w=$g e12_cur1
		execute "e"+c1+c2+"_cur1:=ep"+c2+"/1000"

		make/O/N=2 $("e"+c1+c2+"_cur2")
		wave e12_cur2 = $("e"+c1+c2+"_cur2")
		setscale/p x,(cEp2[0]),(cEp2[inf]-cEp2[0]),waveunits(cEp2,-1),e12_cur2
		setscale d,0,1,waveunits(cEp1,-1),e12_cur2
		appendtograph/w=$g/VERT e12_cur2
		execute "e"+c1+c2+"_cur2:=ep"+c1+"/1000"

		ModifyGraph rgb=(0,0,0), lstyle=2, lsize=0.1
		ColorScale/C/N=text1/F=0/Z=1/A=RC/X=0.00/Y=0.00/E image=$("e"+c1+c2), heightPct=120, log=1
		ColorScale/C/N=text1 logLTrip=0.1,minor=1
		
		if(0==cmpstr(c1,"z"))
			ModifyGraph/w=$g swapXY=1
		endif
	endif
	
end

Function SetSliceCoordinateProc(sva) : SetVariableControl
	STRUCT WMSetVariableAction &sva

	switch( sva.eventCode )
		case 1: // mouse up
		case 2: // Enter key
		case 3: // Live update
			string win = sva.win
//			Variable dval = sva.dval
			string c1 =  win[7,7]
			string c2 = win[8,8]
			string c3 = win[10,10]

			wave e12 = $("e"+c1+c2)
			wave/C cEx, cEy, cEz
			wave freq2
			nvar Efnow, Efnowi, Ef
			variable Efi
			FindLevel/Q freq2, ef * 1e6
			if(v_flag)
				efi = Efnowi
			else
				efi = x2pnt(freq2, V_Levelx)
			endif
			if(Efnowi==efi)
				if((ef < efnow) && (efi>0))
					efi -= 1
				elseif((ef > efnow) && (efi<numpnts(freq2)-1))
					efi += 1
				endif
			endif
			ef = freq2[efi]/1e6
			Efnowi = efi
			efnow = ef
			
			show_e12("graph_exy_z")
			show_e12("graph_eyz_x")
			show_e12("graph_ezx_y")
			
			break
		case -1: // control being killed
			break
	endswitch

	return 0
End

function show_e12(win)
string win
	string c1 =  win[7,7]
	string c2 = win[8,8]
	string c3 = win[10,10]
	wave e12 = $("e"+c1+c2)
	wave/C cEx, cEy, cEz

	nvar ep3 = $("Ep"+c3)
	wave cEp3 = $("cEp"+c3)
	FindLevel/Q/Q cEp3, ep3/1000
	variable ep3i = V_Levelx
	nvar Efi = Efnowi
	
	if(0==cmpstr(c3,"x"))
		e12 = sqrt(cabs(cEx[ep3i][p][q][Efi])^2+cabs(cEy[ep3i][p][q][Efi])^2+cabs(cEz[ep3i][p][q][Efi])^2)
	elseif(0==cmpstr(c3,"y"))
		e12 = sqrt(cabs(cEx[q][ep3i][p][Efi])^2+cabs(cEy[q][ep3i][p][Efi])^2+cabs(cEz[q][ep3i][p][Efi])^2)
	elseif(0==cmpstr(c3,"z"))
		e12 = sqrt(cabs(cEx[p][q][ep3i][Efi])^2+cabs(cEy[p][q][ep3i][Efi])^2+cabs(cEz[p][q][ep3i][Efi])^2)
	endif
	
	wave freq2
	variable ef = freq2[efi]
	TextBox/w=$win/C/N=text0/F=0/X=0/Y=0/E=2 "freq = "+num2str(ef)+" MHz, "+c3+" = "+num2str(ep3)+" mm"
end

function show_ofd()
	NVAR NFreq2
	if(NFreq2>0)
		variable/g Epx, Epy, Epz
		show_ofd_sub("x", "y", "z")
		show_ofd_sub("y", "z", "x")
		show_ofd_sub("z", "x", "y")
		
		show_e12("graph_exy_z")
		show_e12("graph_eyz_x")
		show_e12("graph_ezx_y")
	endif
	
//	STRUCT WMSetVariableAction sva
//	sva.eventCode = 3
//	sva.win = "graph_exy_z"
//	SetSliceCoordinateProc(sva)
//	sva.win = "graph_eyz_x"
//	SetSliceCoordinateProc(sva)
//	sva.win = "graph_ezx_y"
//	SetSliceCoordinateProc(sva)
	
	NVAR NFeed
	if(NFeed>0)
		show_zin()
	endif
end

function show_zin()
	NVAR NFreq1, NFeed
	wave Freq1
	wave/c Zin
	make/O/N=(NFreq1)/C Zinp
	setscale d,0,1,"Ohm",Zinp
	DoWindow graph_zin
	if(!V_Flag)
		display Zinp vs Freq1
		DoWindow/C graph_zin
		appendtograph/w=graph_zin Zinp vs Freq1
		appendtograph/w=graph_zin Zinp vs Freq1
		ModifyGraph/w=graph_zin cmplxMode(Zinp)=1, cmplxMode(Zinp#1)=2, grid=1, rgb(Zinp#1)=(0,0,65535)
		ModifyGraph/w=graph_zin  rgb(Zinp#2)=(0,0,0),cmplxMode(Zinp#2)=3
		Label/w=graph_zin left "Zin / \\U";DelayUpdate
		Label/w=graph_zin bottom "Frequency / \\U"
		ControlBar/w=graph_zin 20
		variable/g feedi
		SetVariable setvar0 proc=SetFeedProc,value=feedi,limits={0,0,1},win=graph_zin
		make/N=2 zincsr
		execute "zincsr:=Ef*1e6"
		setscale x,0,100,"Ohm",zincsr
		appendtograph/w=graph_zin/VERT zincsr
		ModifyGraph/w=graph_zin lstyle(zincsr)=2,rgb(zincsr)=(0,0,0)
	endif
	
	DoWindow graph_ref
	if(!V_Flag)
		make/o Refp
		display/N=graph_ref Refp vs Freq1
		modifygraph/w=graph_ref grid=1, zero(left)=1
		SetAxis/w=graph_ref/A/E=1 left
		Label/w=graph_ref left "Reflection / \\U"
		Label/w=graph_ref bottom "Frequency / \\U"
	endif
	
	DoWindow graph_smith
	if(!V_Flag)
		smith_chart()
		DoWindow/C graph_smith
		make/o zsx, zsy
		appendtograph/w=graph_smith zsy vs zsx
		ModifyGraph/w=graph_smith lsize(zsy)=2,rgb(zsy)=(0,0,0)
	endif
	
	show_zin_sub()
end

Function CheckLogProc(cba) : CheckBoxControl
	STRUCT WMCheckboxAction &cba

	switch( cba.eventCode )
		case 2: // mouse up
			string win = cba.win
			Variable checked = cba.checked
			ModifyImage/w=$(win) $(win[6,8]) log=(checked)
			ColorScale/C/N=text1 log=(checked)
			break
		case -1: // control being killed
			break
	endswitch

	return 0
End

function show_zin_sub()
	ControlInfo/w=graph_zin setvar0
	variable dval = V_Value

	wave/t feed
	TextBox/w=graph_zin/C/N=text0/F=0/X=0/Y=0/E=2 "Feed#"+num2str(dval)+" : "+feed[dval]
	wave/c zinp, zin
	Zinp = Zin[dval][p]
	
	NVAR NFreq1
	make/o/N=(NFreq1) Refp
	setscale d,0,1,"dB",Refp
	variable z0 = numberbykey("z0", feed[dval])
	Refp=20*log(max(cabs((Zinp-z0)/(Zinp+z0)),1e-10))
	
	variable zz0 = 50
	make/o/N=(numpnts(zinp)) zsx = smith_x(real(zinp)/zz0, imag(zinp)/zz0)
	make/o/N=(numpnts(zinp)) zsy = smith_y(real(zinp)/zz0, imag(zinp)/zz0)
end

Function SetFeedProc(sva) : SetVariableControl
	STRUCT WMSetVariableAction &sva

	switch( sva.eventCode )
		case 1: // mouse up
		case 2: // Enter key
		case 3: // Live update
			show_zin_sub()
			break
		case -1: // control being killed
			break
	endswitch

	return 0
End

スミスチャート表示用ライブラリ

上記のプログラムから呼び出して使っている。

SmithChart.ipf

#pragma TextEncoding = "UTF-8"
#pragma rtGlobals=3		// Use modern global access method and strict wave access.

// ====================== スミスチャート(イミタンスチャート)

// 使い方:
// wave/C zc // プロットしたい複素ウェーブ
// variable z0 = 50	// 基準インピーダンス
// make/N=(numpnts(zc)) zsx = smith_x(real(zc)/z0, imag(zc)/z0)	// 反射率に直す
// make/N=(numpnts(zc)) zsy = smith_y(real(zc)/z0, imag(zc)/z0)	// 反射率に直す
// smith_chart()	// グリッドを表示
// appendtograph zsy vs zsx	// トレースを追加
//	ModifyGraph lsize(zsy)=2,rgb(zsy)=(0,0,0)
//
// あるいは:
// smith_display(zc, 50)

function smith_display(zc, z0)
wave/c zc
variable z0
	string wname = nameofwave(zc)+"_smith"
	make/N=(numpnts(zc)) $(wname+"x"), $(wname+"y")
	wave wx = $(wname+"x")
	wave wy = $(wname+"y")
	wx = smith_x(real(zc)/z0, imag(zc)/z0)
	wy = smith_y(real(zc)/z0, imag(zc)/z0)
	smith_chart()
	appendtograph wy vs wx
	ModifyGraph lsize($(wname+"y"))=2,rgb($(wname+"y"))=(0,0,0)
end

// インピーダンスから反射率へ

function smith_x(x, y)
variable x, y
	return (x^2+y^2-1)/((x+1)^2+y^2)
end

function smith_y(x, y)
variable x, y
	return 2*y/((x+1)^2+y^2)
end

// 反射率からインピーダンスへ

function smith_ix(x, y)
variable x, y
	return -(x^2+y^2-1)/((x-1)^2+y^2)
end

function smith_iy(x, y)
variable x, y
	return -2*y/((x-1)^2+y^2)
end

function smith_csr_moved(s)
STRUCT WMWinHookStruct &s
	if(s.eventCode!=7)
		return 0
	endif
	if(0==cmpstr(s.traceName, "smith_gridy") || strlen(s.traceName)==0)
		TextBox/C/N=text0/A=LB/X=0/Y=0/E=2/F=0/B=1 ""
		return 0
	endif
	
	string csr = s.cursorName
	variable h = hcsr($csr, s.winName)
	variable v = vcsr($csr, s.winName)
	variable ix = smith_ix(h, v)
	variable iy = smith_iy(h, v)
	string str
	sprintf str, "reflection = %.1g\n%.3f%+.3fi", sqrt(h^2+v^2), ix, iy
	TextBox/C/N=text0/A=LB/X=0/Y=0/E=2/F=0/B=1 str
end

// イミタンスチャートウィンドウの作成
function smith_chart()
	// グリッドウェーブの作成
	smith_chart_grid()
	
	// グリッドの表示
	wave smith_gridy, smith_gridx
	Display smith_gridy vs smith_gridx
	AppendToGraph smith_gridy vs smith_gridx
	ModifyGraph margin=20,height={Plan,1,left,bottom}
	ModifyGraph rgb(smith_gridy)=(65535,32768,32768),rgb(smith_gridy#1)=(32768,40777,65535)
	ModifyGraph muloffset(smith_gridy#1)={-1,0}
	ModifyGraph tick=3, noLabel=2, standoff=0, axRGB=(65535,65535,65535,0)
	setwindow kwTopWin, hook(smith_hook)=smith_csr_moved, hookevents=4

	// テキスト表示
	SetDrawLayer UserFront

	variable i
	make/o smithx={-1.08,-0.6544,-0.3162,0.011,0.283,0.614,1.022,-1.011,-0.6654,0.6508,0.9632,0.9484,0.636,-0.658,-1.0294}
	make/o smithy={0,0,0,0,0,0,0,0.4192,0.8382,0.772,0.3346,-0.397,-0.8346,-0.908,-0.4742}
	make/o/t smiths={"0","0.2","0.5","1","2","5","∞","0.2j","0.5j","2j","5j","-5j","-2j","-0.5j","-0.2j"}
	for(i=0; i<numpnts(smithx); i+=1)
		SetDrawEnv xcoord= bottom,ycoord= left,textrgb= (65535,32768,32768)
		DrawText smithx[i],smithy[i],smiths[i]
	endfor
	SetDrawEnv textrgb= (65535,32768,32768),textxjust= 2
	DrawText 1,0,"Impedance"
	
	make/o smithx={-1.084,-0.6508,-0.3088,0.0074,1.022,0.2426,0.5734,-1.011,-0.6838,0.5698,0.9338,0.9154,0.5734,-0.6874,-1.0184}
	make/o smithy={-0.08,-0.08,-0.08,-0.08,-0.08,-0.08,-0.08,0.342,0.772,0.8382,0.408,-0.478,-0.8934,-0.8492,-0.397}
	make/o/t smiths={"∞","5","2","1","0","0.5","0.2","5j","2j","0.5j","0.2j","-0.2j","-0.5j","-2j","-5j"}
	for(i=0; i<numpnts(smithx); i+=1)
	 	SetDrawEnv xcoord= bottom,ycoord= left,textrgb= (32768,40777,65535)
		DrawText smithx[i],smithy[i],smiths[i]
	endfor
	SetDrawEnv textrgb= (32768,40777,65535)
	DrawText -0.0386,0,"Admittance"
	
	SetDrawEnv xcoord= bottom,ycoord= left,textrgb= (44253,29492,58982)
	DrawText -0.0074,1.0146,"j"
	SetDrawEnv xcoord= bottom,ycoord= left,textrgb= (44253,29492,58982)
	DrawText -0.0184,-1.0698,"-j"
	
	killwaves smithx, smithy, smiths
end

function smith_chart_grid()
	make/o/N=2000 smith_gridx = nan, smith_gridy = nan
	wave gx = smith_gridx
	wave gy = smith_gridy
	variable i = 0, j

	// 等抵抗線
	make/o smith_zr = {0, 0.2, 0.5, 1, 2, 5}
	for(j=0; j<numpnts(smith_zr); j+=1)
		gx = (i<=p && p<=i+101) ? (1-smith_x(smith_zr[j],0))/2*cos(2*pi*p/100)+(1+smith_x(smith_zr[j],0))/2 : gx
		gy = (i<=p && p<=i+101) ? (1-smith_x(smith_zr[j],0))/2*sin(2*pi*p/100) : gy
		i+= 101
		gx[i]=nan; gy[i]=nan ; i+= 1
	endfor
	
	// 等インダクタンス線
	make/o smith_zi = {-5, -2, -1, -0.5, -0.2, 0, 0.2, 0.5, 1, 2, 5}
	make/o/N=100 smith_zr = smith_ix((p-50)/50, 0)
	for(j=0; j<numpnts(smith_zi); j+=1)
		gx = (i<=p && p<i+100) ? smith_x(smith_zr[p-i], smith_zi[j]) : gx
		gy = (i<=p && p<i+100) ? smith_y(smith_zr[p-i], smith_zi[j]) : gy
		i+= 100
		gx[i]=nan; gy[i]=nan ; i+= 1
	endfor
	
	killwaves smith_zr, smith_zi
	redimension/N=(i) smith_gridx, smith_gridy
end

添付ファイル: filemicrostrip2.png 168件 [詳細] filemicrostrip1.png 167件 [詳細]

Counter: 4413 (from 2010/06/03), today: 1, yesterday: 0