ya Alcohol content of "mizuwari"

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 = mass_alc / mass_all
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); /*blank line*/
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) {
/*
Search for p such that q equals the given value

Newton’s method
f(x) = 0, x = ?
x_{n+1} = x_n - f(x_n) / f'(x_n)

Numerical differentiation
f'(x) = (f(x+h) - f(x)) / h
h: delta_x

V: Volume
M: Mass
rho: density [kg/m^3]
V_alc = 1/rho_alc * M_alc
V_sol = 1/rho_sol * M_sol
rho_sol = ald(p)
rho_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/rho_alc * p * M_sol) / (1/rho_sol * M_sol)
= 1/rho_alc * p * rho_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) { /* initialize q */
q = 1;
p = q;
} elseif (q <= 0) {
q = 0;
p = 0.1; /* avoid division by zero */
} 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); /* Accuracy check */
} else {
accura = abs(p - prep); /* for neighbor 0 */
}
} while ((accura > eps) and (k <= cnt)); /* Termination condition */

return (p);
}


function ald(p, t, numeric A[], numeric B[], numeric C[][]) {
/*
From the International Alcohol Table OIML R22 1975
Calculate density from mass fraction
rho: Density [kg/m^3]
p: Alcohol concentration [% mass] e.g., 12%: p=0.12
t: Temperature [°C]
*/
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;
}