/*
    Oy 86,597    原点の x,y 座標
    uy 89,220    10pF のスパン
    2V 141,246   2V ...
    4V 192,307
    6V 244,334
    8V 294,350
    10V 347,362
    12V 397,371
*/
oy : 597 $ uy : 220 $
array(dat, 1, 6) $     /* 電圧、容量 (F) */
dat[1,0] : 130 $
dat[1,1] : 246 $
dat[1,2] : 307 $
dat[1,3] : 334 $
dat[1,4] : 350 $
dat[1,5] : 362 $
dat[1,6] : 371 $

for i : 0 while i < 7 do (
    dat[0,i] : 2 * i,
    if i = 0 then dat[1, i] : 18e-12 else    /* 18p 固定 */
    dat[1,i] : 1e-12 * 10**float((oy - dat[1,i])/(oy - uy)),
    print(dat[0,i], dat[1,i])
) $

/* c0v, c1v, c2v はバイアス電圧「値」とは無関係。単に low, mid, hi を表す。 */
c : 1 $    /* mid の index */
low : dat[0,0] $
mid : dat[0,c] $
hi  : dat[0,6] $
c0v : dat[1,0] $
c1v : dat[1,c] $
c2v : dat[1,6] $
Vj : 1.0 $          /* 適当な初期値 */
M : 1.0 $
Cjo : c0v $

/*
    ln(c0v/c1v) / ln(c0v/c2v) = ln((Vj+1)/(Vj+0)) / ln((Vj+2)/(Vj+0)) を解く
*/
lop : log(c0v/c1v) / log(c0v/c2v) $  /* ここの log() は自然対数 */
fn(x) := log((x+mid)/(x+low)) / log((x+hi)/(x+low)) $
/* print( "lop", lop ) $  目標値 */
f : 0 $
for i : 0 while i < 50 do (
    if abs((fn(Vj) - lop)/lop) < 1e-6 then return(f : 1),
    df : (fn(Vj + 1e-6) - fn(Vj))/1e-6,  /* 微係数を数値計算してしまう悪いヤツ */
    Vj : Vj - (fn(Vj) - lop)/(2.0 * df)
    /* print( Vj, fn(Vj) - lop ) */
) $
if f # 1 then error("") $

/* 他のパラメーターを計算 */
m : log(c0v/c1v) / log((Vj+mid)/(Vj+low)) $
print("Cjo", Cjo, "Vj", Vj, "m", m) $
fcap(Vbias) := Cjo / (1 + Vbias/Vj)**m $
plot2d(fcap(x), [x,low,hi]) $

/* おまけ */
fcap6(Vbias) := Cjo / (1 + Vbias/0.12658)**0.33098 $
fcap4(Vbias) := Cjo / (1 + Vbias/0.19803)**0.36646 $
fcap2(Vbias) := Cjo / (1 + Vbias/0.52876)**0.47707 $
plot2d([fcap6(x),fcap4(x),fcap2(x)], [x,low,24 /* hi */]) $