ya 威士忌加水的酒精濃度

numeric A[13],B[7],C[6][12];
A[1]=9.982012300E2;
A[2]=-1.929769495E2;
A[3]=3.891238958E2;
A[4]=-1.668103923E3;
A[5]=1.352215441E4;
A[6]=-8.829278388E4;
A[7]=3.062874042E5;
A[8]=-6.138381234E5;
A[9]=7.470172998E5;
A[10]=-5.478461354E5;
A[11]=2.234460334E5;
A[12]=-3.903285426E4;

B[1]=-2.0618513E-1;
B[2]=-5.2682542E-3;
B[3]=3.6130013E-5;
B[4]=-3.8957702E-7;
B[5]=7.1693540E-9;
B[6]=-9.9739231E-11;

C[1][1]=1.693443461530087E-1;
C[1][2]=-1.046914743455169E1;
C[1][3]=7.196353469546523E1;
C[1][4]=-7.047478054272792E2;
C[1][5]=3.924090430035045E3;
C[1][6]=-1.210164659068747E4;
C[1][7]=2.248646550400788E4;
C[1][8]=-2.605562982188164E4;
C[1][9]=1.852373922069467E4;
C[1][10]=-7.420201433430137E3;
C[1][11]=1.285617841998974E3;

C[2][1]=-1.193013005057010E-2;
C[2][2]=2.517399633803461E-1;
C[2][3]=-2.170575700536993;
C[2][4]=1.353034988843029E1;
C[2][5]=-5.029988758547014E1;
C[2][6]=1.096355666577570E2;
C[2][7]=-1.422753946421155E2;
C[2][8]=1.080435942856230E2;
C[2][9]=-4.414153236817392E1;
C[2][10]=7.442971530188783;

C[3][1]=-6.802995733503803E-4;
C[3][2]=1.876837790289664E-2;
C[3][3]=-2.002561813734156E-1;
C[3][4]=1.022992966719220;
C[3][5]=-2.895696483903638;
C[3][6]=4.810060584300675;
C[3][7]=-4.672147440794683;
C[3][8]=2.458043105903461;
C[3][9]=-5.411227621436812E-1;

C[4][1]=4.075376675622027E-6;
C[4][2]=-8.763058573471110E-6;
C[4][3]=6.515031360099368E-6;
C[4][4]=-1.515784836987210E-6;

C[5][1]=-2.788074354782409E-8;
C[5][2]=1.345612883493354E-8;

/* p = 酒精質量 / 總質量
q = vol_alc / vol_all
d = mass / volume
_s:strong, _w:weak
_a:alcohol, _q:aqua */
qs = qs / 100;
ps = alq2p(qs,t,A,B,C);
ds = ald(ps,t,A,B,C);
ms = ds * vs;
msa = ps * ms;
msq = ms - msa;

da = ald(1,t,A,B,C);
dq = ald(0,t,A,B,C);
mq = dq * vq;

mwa = msa;
mwq = msq + mq;
mw = mwa + mwq;
pw = mwa / mw;
dw = ald(pw,t,A,B,C);
vwa = mwa / da;
vw = mw / dw;
qw = vwa / vw;

print(vw,100*qw);
print(-256); /*空白行*/
print(round(vs, 2), round(ms * 1E-3, 2), round(ds * 1E-3, 6));
print(round(vwa, 2), round(mwa * 1E-3, 2), round(da * 1E-3, 6));
print(round(msq/dq,2),round(msq*1E-3,2),round(dq*1E-3,6));
print(round(vw, 2), round(mw * 1E-3, 2), round(dw * 1E-3, 6));

function alq2p(q, t, A, B, C) {
/*
尋找使 q 達到給定數值的 p

牛頓法
f(x) = 0, x = ?
x_n+1 = x_n - f(x_n) / f'(x_n)

數值微分
f'(x) = (f(x+h)-f(x)) / h
h: delta_x

V：體積
M：質量
ρ：密度 [kg/m³]
V_alc = 1/ρ_alc * M_alc
V_sol = 1/ρ_sol * M_sol
ρ_sol = ald(p)
ρ_alc = ald(1)
M_sol = M_alc + M_aqa

p = M_alc / M_sol
M_alc = p * M_sol
M_aqa = M_sol - M_alc = (1-p) * M_sol

q = V_alc / V_sol
= (1/ρ_alc * p * M_sol) / (1/ρ_sol * M_sol)
= 1/ρ_alc * p * ρ_sol
= p * ald(p) / ald(1)

f(p) = p * ald(p) / ald(1) - q = 0, p = ?
p_n+1 = p_n - f(p_n) / f'(p_n)
*/
h = 1 / 2^10;
eps = 1 / 2^20;
if (q > 1) { /* 初始化 q */
q = 1;
p = q;
} elseif (q <= 0) {
q = 0;
p = 0.1; /* 避免除以 0 */
} else {
p = q;
}
cnt = 20;
k = 0;
do {
k = k + 1;
prep = p;
ald1 = ald(1,t,A,B,C);
fp = p * ald(p,t,A,B,C) / ald1 - q;
fph = (p+h) * ald(p+h,t,A,B,C) / ald1 - q; /*f(p+h)*/
dfp = (fph - fp) / h; /*f'(p)*/
p = p - fp / dfp;
if (prep > eps) {
accura = abs((p-prep)/prep); /* 精度檢查 */
} else {
accura = abs(p-prep); /* 鄰域 0 */
}
} while ((accura>eps) and (k<=cnt)); /* 結束判斷 */

return (p);
}


function ald(p, t, numeric A[], numeric B[], numeric C[][]) {
/*
根據國際酒精表 OIML R22 1975
由重量比計算密度
rho：密度 [kg/m^3]
p：酒精濃度 [%mass] 例如 12%：p=0.12
t：溫度 [℃]
*/
rho = A[1];
for (k=2;k<=12;k=k+1) {
rho = rho + A[k]*p^(k-1);
}
for (k=1;k<=6;k=k+1) {
rho = rho + B[k]*(t-20)^k;
}
for (k=1;k<=11;k=k+1) {
rho = rho + C[1][k]*p^k*(t-20)^1;
}
for (k=1;k<=10;k=k+1) {
rho = rho + C[2][k]*p^k*(t-20)^2;
}
for (k=1;k<=9;k=k+1) {
rho = rho + C[3][k]*p^k*(t-20)^3;
}
for (k=1;k<=4;k=k+1) {
rho = rho + C[4][k]*p^k*(t-20)^4;
}
for (k=1; k <= 2; k=k+1) {
rho = rho + C[5][k] * p^k * (t - 20)^5;
}
return(rho);
}