ソフトウェア/OpenFDTD の変更点
更新- 追加された行はこの色です。
- 削除された行はこの色です。
- ソフトウェア/OpenFDTD へ行く。
- ソフトウェア/OpenFDTD の差分を削除
[[公開メモ]] &katex();
* OpenFDTD を使ってプリント基板上の高速信号の伝播を理解したい [#v3733737]
例えば４層基板上のマイクロストリップラインの周囲の電場をこんな風に計算するとか。
&ref(microstrip1.png,,50%);
&ref(microstrip2.png,,33%);
** 目次 [#taf6b842]
#contents
** OpenFDTD について [#ta0c33d3]
株式会社 EEM により開発＆メンテナンスされていたオープンソースの FDTD シミュレーションソフト。
[[Wikipedia:FDTD法]] は主にマクスウェル方程式を数値的に解いて電磁波のシミュレーションを行う際に使われる数値計算手法の一種。
株式会社 EEM は2022年9月30日に解散したため、現在は大賀さんが個人的に (?) メンテナンスを続けてらっしゃるようです。
** OpenFDTD のダウンロード＆説明 [#v73424f6]
以前は 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 も提供されていて、モデルの作成や計算結果のグラフ表示が可能です。
* 計算例 [#m6251b49]
- [[ソフトウェア/OpenFDTD/マイクロストリップライン]]
- [[ソフトウェア/OpenFDTD/解放端・固定端]]
- [[ソフトウェア/OpenFDTD/５０Ω終端]]
- [[ソフトウェア/OpenFDTD/パッドを作って素子を載せる]]
* プリント基板 (PCB) のシミュレーションに使うことを想定した注意点 [#k7b22011]
プリント基板 (PCB) のシミュレーションでは各部のサイズを厳密に設定しなければならならず、セルの大きさ分の誤差が計算結果に影響する。
このため、FDTD 計算の基本を押さえておかないと計算結果が合わず困ったことになります。
** FDTD 計算の基本 [#t2a1160b]
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)$ が基準点となる。
** マテリアルの範囲 [#j3fd261d]
誘電体や金属を配置する際、その範囲は「セル単位」ではなく「座標空間内に厳密に設定された範囲」になる。
どういうことかというと、電場や磁場を計算する際、その「基準点」の座標がどこかによって、どのマテリアルのパラメータを使って計算するかが決まる。
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 では、
- 設定したマテリアルの範囲の「境界上の点」は「マテリアルの内部」として考慮されること、
- 複数のマテリアルに含まれる点は「最後に指定したマテリアル」で上書きされること、
にも注意が必要だ。
** 具体例で考える [#kfb4e75b]
例えば厚さ 1.6 mm の両面基板を考え、裏のグランドプレーンを $z=-1.6\ \text{mm}$ に配置したいとする。
それには任意の $k$ に対して $z_k=-1.6\ \text{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_{k+1}-z_{k-1})/2$ となる。この面が系の端にある場合、例えば $z_0$ にある場合は、$z_{-1}$ の格子面が存在しないため、代わりに反対側の格子面までの距離が使われて、$z_1-z_0$ がトータルの厚さになる。これに伝導率を掛けてシート抵抗が得られる。
この電流による磁場は、本来太さを持つ辺を流れる電流が、あたかも辺の上に集中しているかのように考えて計算される。
同様に、パターン面を $z=0\ \text{mm}$ として、そこに厚さゼロの金属を配置すればそこに電流を流すことができる。
同様に、パターン面を $z=0\ \text{mm}$ として、そこに厚さゼロの金属を配置すればそこに電流を流せる。
この間を誘電体、例えば FR-4 で埋めるには、グランドプレーンやパターン面の金属を配置する「前に」 $\epsilon_r=4.3$ の材質を $-1.6\ \text{mm}\le z\le 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 とかしかないので、これを正しく計算できるほどメッシュを細かくするのは現実的ではなく、何かうまくごまかしながら計算していくしかない感じ？？？
影響を減らすにはやはり、パターンの外周のすぐ外の格子面までの距離を十分に小さくとるしかないのだと思う。とはいえ例えばパターンの厚さは 15 um とかしかないので、これを正しく計算できるほどメッシュを細かくするのは現実的ではなく、何かうまくごまかしながら計算していくしかない感じ？？？
* 結果表示用 Igor Pro スクリプト [#xca6dd8d]
計算結果を図にするのが大変なので、Igor Pro で直接読み込んで表示できるようにした。
いろいろ手抜きをしているので、今後も改善していかないと。
FDTD 計算結果は ofd.out というファイルにごそっと書き出されている。
グラフ化などのポストプロセスはこのデータを読み込んで、画像に加工している。
なので、ofd.out を Igor で読み込めば、好きなように表示可能！
- 電場を xy, yz, zx 面で切って表示するところまで書いた。
- Igor6 でも動くように少し直した。int64 に対応していないので、int32 で２回読んでつなぎ合わせる必要があった。あと、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
* スミスチャート表示用ライブラリ [#z2586223]
上記のプログラムから呼び出して使っている。
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
Counter: 8454 (from 2010/06/03),
today: 2,
yesterday: 39