ya水割りのアルコール度数 実行数: 1833
出来上がった水割りのアルコール度数を計算します。 | ||
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) { /* qが与えられたに値になるようなpを探索 Newton法 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: 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) { /* init q */ q = 1; p = q; } elseif (q <= 0) { q = 0; p = 0.1; /* avoid 0 div */ } 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); /* for neighbour 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] e.g. 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); } |
本ライブラリは会員の方が作成した作品です。 内容について当サイトは一切関知しません。