/*
    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
    dat[1,i] : 1e-12 * 10**float((oy - dat[1,i])/(oy - uy)),
    print(dat[0,i], dat[1,i])
) $

c0v : dat[1,0] $    /* 18p 固定 */
c6v : dat[1,3] $    /* あとで必要なら手動修正 */
c12v : dat[1,6] $

Vj : 1.0 $          /* 適当な初期値 */
M : 1.0 $
Cjo : c0v $

/*
    ln(c0v/c6v) / ln(c0v/c12v) = ln((Vj+6)/(Vj+0)) / ln((Vj+12)/(Vj+0)) を解く
*/
lop : log(c0v/c6v) / log(c0v/c12v) $  /* ここの log() は自然対数 */
fn(x) := log((x+6)/(x+0)) / log((x+12)/(x+0)) $
/* 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/c6v) / log((Vj+6)/(Vj+0)) $
print("Cjo", Cjo, "Vj", Vj, "m", m) $
fcap(Vbias) := Cjo / (1 + Vbias/Vj)**m $
plot2d(fcap(x), [x,0,12]) $