/*
    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 $
fpprec : 6 $
for i : 0 while i < 7 do (
    print(dat[0,i], fcap(dat[0,i]), fcap(dat[0,i])/dat[1,i])
) $
fpprec : 16 $

/* 浮遊容量を考慮した式 */
fcap_cp(Vbias, Cjo, Vj, m, Cp) := Cjo / (1 + Vbias/Vj)**m + Cp $

/* 最小自乗誤差。容量の対数を目安にしている。 */
mmse(Cjo, Vj, m, Cp) := (
    sum : 0,
    for i : 0 while i < 7 do (
        s : fcap_cp(dat[0,i], Cjo, Vj, m, Cp)/dat[1,i] - 1,
        sum : sum + s**2
    ),
    sum
) $

/* Cjo, Vj, m, Cp のフィット */
Cp : 0.2e-12 $
fct : 1e-3 $
fcu : float(1 / 12) $
for i : 0 while i < 50 do (
    u0 : mmse(Cjo, Vj, m, Cp),
    u1 : mmse((1+fct)*Cjo, Vj, m, Cp)/u0,
    u2 : mmse(Cjo, (1+fct)*Vj, m, Cp)/u0,
    u3 : mmse(Cjo, Vj, (1+fct)*m, Cp)/u0,
    u4 : mmse(Cjo, Vj, m, (1+fct)*Cp)/u0,
    /* print(i, "u0=", u0), */
    /* print(u1, u2, u3, u4), */
    Cjo : Cjo / u1**fcu,
    Vj  : Vj  / u2**fcu,
    m   : m   / u3**fcu,
    Cp  : Cp  / u4**fcu,
    0
) $

(
    fpprec : 6,
    print("Cjo", Cjo, "Vj", Vj, "m", m, "Cp", Cp),
    fpprec : 16,
    print("")
) $

/* 一覧表 */
fpprec : 6 $
for i : 0 while i < 7 do (
    print(dat[0,i], fcap(dat[0,i]), fcap_cp(dat[0,i], Cjo, Vj, m, Cp)/dat[1,i])
) $