Eg,Ep,Ip としておいてください．
単位は，VおよびAで入力します．
ファイルの例を次に示します．
1 Eg,Ep,Ip
2 0,0,0
3 0,50,0.005
4 0,100,0.0118
5 0,150,0.02
6 -2,17,0
7 -2,50,0.0015
8 -2,100,0.0062
9 -2,150,0.013
10 -2,200,0.021
11 -4,43,0
12 -4,100,0.0027
13 -4,150,0.0075
14 -4,200,0.0147
15 -4,250,0.0232
16 -6,70,0
17 -6,100,8.0E-4
18 -6,150,0.0040
19 -6,200,0.0094
20 -6,250,0.017
21 -8,95,0
22 -8,150,0.0017
23 -8,200,0.0053
24 -8,250,0.0117
25 -8,300,0.02
26 -10,125,0
27 -10,150,4.0E-4
28 -10,200,0.0027
29 -10,250,0.0074
30 -10,300,0.0142
31 -12,157,0
32 -12,200,0.0012
33 -12,250,0.0043
34 -12,300,0.0097
35 -12,350,0.0164
36 -14,183,0
37 -14,200,3.0E-4
38 -14,250,0.0024
39 -14,300,0.0061
40 -14,350,0.0115
41 -16,205,0
42 -16,250,0.0010
43 -16,300,0.0034
44 -16,350,0.0077
45 -18,228,0
46 -18,250,3.0E-4
47 -18,300,0.0017
48 -18,350,0.0048
49 -20,265,0
50 -20,300,7.0E-4
51 -20,350,0.0028
この例では，50点のプレート電流を読み取っていますが，
20点程度でも十分な精度のパラメータを得ることができます．
各グリッド電圧に対して，最低でもカットオフに近いプレート電圧と，
それよりもかなり大きな電圧(特性図上に載っている最大のプレート電圧あるいは
プレート損失近辺のプレート電圧)におけるプレート電流を入力しておきます．
CSVファイルは，テキストエディタで作成することができます． または，Microsoft Excelなどの表計算ソフトで作成し， CSV形式で保存することもできます． 表計算ソフトを用いれば，プレート電流をmA単位で入力して， A単位に変換することも容易でしょう．
1 "Ip.cal" <-
2 function(fn, name, gfn, G.lim=NULL, Ig.ratio=NULL, verbose=FALSE, ...)
3 {
4 # モデルのパラメータを求める
5 # fn: データファイル名
6 # name: パラメータオブジェクト名
7 # gfn: キャリブレーションの状況を示すプレート特性のグラフファイル名.
8 # 指定されていなければグラフを画面に表示する.
9 # G.lim: Ep=Eg のパービアンス(省略可)
10 # Ig.ratio: Ep=Eg のときのカソード電流に対するグリッド電流の比率(省略可)
11 # verbose: 途中経過の表示
12 # ...: 電極間容量および他のパラメータ
13 # Cgp, Cgk Cpk, Ci, Co: 電極間容量
14 # ig=(Ep,Ip,Ig): グリッド電流データ
15 # pentode.chara=(Ep,Eg2,Eg,Ip,Ig2,rp): 五極管動作例
16
17 sig <- function(x) signif(x, 8)
18 print.level <- 1 # 診断出力のレベル
19
20 # キャリブレーション
21 x <- scan(fn, sep=",", what="") # データファイルを読む
22 egepip <- matrix(as.numeric(x[-(1:3)]), ncol=3, byrow=TRUE)
23 eg <- egepip[,1] # グリッド電圧
24 ep <- egepip[,2] # プレート電圧
25 ip <- egepip[,3] # プレート電流
26 n <- length(eg) # キャリブレーションデータのポイント数
27
28 if (!missing(gfn))
29 psi(file=gfn, wh=c(103.9, 68.8))
30
31 # グラフを描く範囲を決定する
32 xlim <- c(0, max(ep)*1.04)
33 ylim <- c(0, max(ip)*1.04)
34 xlab <- eunit1(pretty(xlim, 8), "V")
35 ylab <- eunit1(pretty(ylim), "A")
36
37 # キャリブレーションポイントをグラフに描く(ここではまだ描かない)
38 plot(ep, ip, type="n", xlim=xlim, ylim=ylim,
39 xaxs="i", yaxs="i", xaxt="n", yaxt="n",
40 xlab=paste("Ep", xlab$unit),
41 ylab=paste("Ip", ylab$unit))
42
43 # 同じグリッド電圧のポイントを同じ色で描く
44 ueg <- unique(eg) # 一意のグリッド電圧
45 for (i in 1:length(ueg)) {
46 ok <- eg == ueg[i]
47 points(ep[ok], ip[ok], pch="+", col=(i-1)%%4+3)
48 }
49
50 # 軸を描く
51 axis(1, at=xlab$at, lab=xlab$lab)
52 axis(1, at=xlab$at, lab=FALSE, tck=1, lty=2, lwd=0.5)
53 axis(2, at=ylab$at, lab=ylab$lab)
54 axis(2, at=ylab$at, lab=FALSE, tck=1, lty=2, lwd=0.5)
55
56 if (is.null(G.lim)) {
57 # 最適化の変数からモデルの変数への変換関数
58 decode <- function(x) {
59 z <- exp(x)
60 z[3] <- z[3] / (z[3] + 1) * 0.65 + 0.35
61 z[4] <- z[4] / (z[4] + 1) * 2 - 1
62 list(G=z[1], muc=z[2], alpha=z[3], Ego=z[4])
63 }
64 } else {
65 # 最適化の変数からモデルの変数への変換関数
66 decode <- function(x) {
67 z <- exp(x)
68 z[3] <- z[3] / (z[3] + 1) * 0.65 + 0.35
69 z[4] <- z[4] / (z[4] + 1) * 2 - 1
70 list(G=z[1], muc=z[2], alpha=z[3], Ego=z[4], G.lim=G.lim, Ig.ratio=Ig.ratio)
71 }
72 }
73
74 # モデルの変数から最適化の変数への変換関数
75 encode <- function(x) {
76 x[3] <- (x[3] - 0.35) / 0.65
77 x[3] <- x[3] / (1 - x[3])
78 x[4] <- (x[4] + 1) / 2
79 x[4] <- x[4] / (1 - x[4])
80 log(x)
81 }
82
83 # 最適化する関数
84 f <- function(x) {
85 p <- decode(x) # モデルの変数へ変換
86 if (print.level == 2)
87 str(p)
88 ipe <- Ip(p, Eg=eg, Ep=ep) # モデルによるプレート電流
89 if (any(is.na(ipe))) { # モデルからNAが出てきた時の処置
90 print(egepip[, 1:2]) # (ここにくるはずはない?)
91 print(ipe)
92 str(p)
93 stop()
94 }
95 ipd <- ipe^(2/3) - ip^(2/3) # モデルと実際の値の差
96 sqrt(mean(ipd^2))/mean(ip^(2/3)) # プレート電流の平均値に対する
97 # RMS偏差の大きさ
98 }
99
100 # パラメータ初期値の推定
101 # 増幅率の推定
102 mineg <- min(eg) # データで最も深いグリッド電圧
103 if (mineg < 0) {
104 minep <- min(ep[eg == mineg]) # その中で最もプレート電圧が低いもの
105 muhat <- minep / -mineg # それをカットオフとみなして mu_c を推定
106 } else
107 muhat <- 30 # B級、ポジティブグリッド管用
108
109 # パービアンスの推定
110 mu0 <- 5/3 * muhat # Eg=0の mu
111 ep0 <- ifelse(eg == max(eg), ep, 0)
112 idx <- ep0 == max(ep0) # Eg=0で最もプレート電圧が高いものにより
113 Ghat <- 2/3 * ip[idx] / (0.6 + ep[idx]/mu0)^1.5 # パービアンスを推定
114 if (verbose) {
115 cat("muhat=", muhat, "\n", sep=" ")
116 cat("Ghat=", Ghat, "\n", sep=" ")
117 }
118
119 ini <- encode(c(Ghat, muhat, 0.6, 0.6)) # 他のパラメータの初期値は0.6
120
121 # 最適化
122 optim.out <- optim(ini, f,
123 control=list(parscale=ini, maxit=5000, reltol=1e-10))
124 cat("code=", optim.out$convergence, " err=", optim.out$value, "\n", sep="")
125
126 # パラメータオブジェクトの作成と保存
127 est <- decode(optim.out$par) # 推定されたモデルのパラメータ
128 for (i in 1:length(est)) # 有効数字を8桁にする
129 est[[i]] <- sig(est[[i]])
130 est <- c(est, list(...))
131 est$err <- optim.out$value
132 est$code <- optim.out$convergence
133 if (!is.null(est$pentode.chara)) {
134 pen.Ep <- est$pentode.chara[1]
135 pen.Eg2 <- est$pentode.chara[2]
136 pen.Eg <- est$pentode.chara[3]
137 pen.Ip <- est$pentode.chara[4]
138 pen.Ig2 <- est$pentode.chara[5]
139 pen.rp <- est$pentode.chara[6]
140 g2r <- pen.Ig2/(pen.Ip + pen.Ig2)
141 a <- (1-pen.Ep/(pen.Ep+10))^1.5
142 r <- (g2r - a)/(1 - a) # (1.41)
143 est$g2.r <- sig(r)
144 est$Ea <- sig(pen.Eg2 - pen.Ip * pen.rp * 1.5)
145 }
146 if (!missing(name)) # オブジェクト名が指定されていれば
147 assign(name, est, env=.GlobalEnv) # 大域オブジェクトとする
148
149 # モデルによるプレート特性をグラフに書き加える
150 g.plate.cal(est, ueg, xlim, ylim)
151
152 if (!missing(gfn))
153 dev.off()
154 invisible()
155 }
まず fn で指定されたファイルからデータを読み込み， 各列ごとに分解して eg, ep, ip に付値しておきます．
グラフを描く範囲は，最大のプレート電圧やプレート電流の4%増しとして， すべての点がグラフの枠内に収まるようにしています．
最適化を行なう関数 optim は，
パラメータの定義域として - 
から の範囲を探索しますが，
モデルのパラメータは0以上なので，
最適化のパラメータとモデルのパラメータの間で
相互に変換する関数を定義しておきます．
それが関数 decode と encode で，
最適化のパラメータの指数をとったものがモデルのパラメータ(> 0)となり，
モデルのパラメータの対数をとったものが最適化のパラメータとなります．
また，モデルのパラメータの対数をとることにより，
パラメータ間の桁数の違いもかなり吸収され，
パラメータのスケールの調整に気を遣わなくてすみます．
パラメータ alpha の範囲は0.35から1であり，
パラメータ Ego の範囲は -1 から1なので，変換が異なります．
最適化の対象となる関数は f で， 最適化関数から渡されたベクトルをモデルのパラメータに変換し， モデルによりプレート電流の1.5乗根を求め，実際のプレート電流の1.5乗根との差をとります． この差はプラスにもマイナスにもなりうるので， 誤差の2乗の平均値の平方根をとり， プレート電流の1.5乗根の平均に対する比をできるだけ小さくなるようにします． このように標準化したのは， 異なる真空管に対するフィッティングの良さを比較できるようにするためです． プレート電流の1.5乗根をとるのは， 各グリッド電圧に対するプレート電圧とプレート電流の組ができるだけ直線上に並ぶようにするためです． 1.5乗根をとることにより，カットオフ近辺の誤差が小さくなるように最適化されます．
最適化に入る前に，パラメータの初期値を推定します．
や Ego は先験的に値がわかっているので，それぞれ0.6とします．
については，
最もグリッド電圧が低い場合のカットオフとなるプレート電圧から求めます．
またパービアンスについては，Eg = 0 の特性曲線から求めます．
推定されたパラメータは有効数字を8桁に制限しています． これは，この関数を用いてパラメータの推定をせずに， 本書に載っているパラメータの数値を使用しても 同じ結果が再現されるようにするためです．
Ip.cal の使用例を以下に示します． ここでは，12AU7.csv というファイルからデータを読み込み， t12AU7 というオブジェクトにパラメータを作成します． グラフは画面に表示されます．
> Ip.cal("12AU7.csv", "t12AU7", Cgp=1.5e-12, Ci=1.6e-12, Co=0.4e-12)
Read 153 items # 読み込んだデータの総数(1行目を含む)
code=0 err=0.02995393 # 収束した
> str(t12AU7)
List of 9
$ G : num 0.000548
$ muc : num 13.2
$ alpha: num 0.58
$ Ego : num 0.829
$ Cgp : num 1.5e-12
$ Ci : num 1.6e-12
$ Co : num 4e-13
$ err : num 0.0300
$ code : int 0
電極間容量は，次のパラメータ名で与えます．
| パラメータ名 | 意味 |
| Cgp | グリッド-プレート間容量 |
| Cgk | グリッド-カソード間容量 |
| Cpk | プレート-カソード間容量 |
| Ci | 入力容量 |
| Co | 入力容量 |
グリッド電流のデータがある場合は， プレート電圧(=グリッド電圧)，プレート電流，グリッド電流の3つの値をベクトルとして パラメータ名 ig で与えます． たとえば 801 の場合， Ep = Eg = 100 V のとき Ip = 225 mA, Ig = 48 mA ですので，
> Ip.cal("801.csv", "t801", Cgp=6e-12, Cgk=4.5e-12, Cpk=1.5e-12,
ig=c(100, 225e-3, 48e-3))
のようにして与えます．
五極管の場合は，標準的な動作例からプレート電圧，スクリーングリッド電圧，コントロールグリッド電圧，プレート電流，スクリーングリッド電流，プレート内部抵抗の6つの値をベクトルとしてパラメータ名 pentode.chara で与えます． たとえば 6L6 の場合， Ep = 250 V, Eg2 = 250 V, Eg = - 14 V, Ip = 72 mA, Ig2 = 5 mA, rp = 22.5 kΩ ですから，
> Ip.cal("6L6.csv", "t6L6T", Cgp=0.6e-12, Ci=10e-12, Co=6.5e-12,
pentode.chara=c(250, 250, -14, 72e-3, 5e-3, 22.5e3),
ig=c(10, 18e-3, 32e-3))
のようにして与えます．
できるだけ
Ep = Eg2 のデータを与えます．
ただし，プレート内部抵抗の値については，
プレート特性図ができるだけ再現されるように適宜調整します．