clear;

// *** 物理定数 ***
rho_h2o = 1.0E6;			// 水の密度 (g/m^3)
rho_air = 1.3E3;			// 空気の密度 (g/m^3)
mu = 1.8E-2;				// 空気の分子粘性係数 (g s/m)
g = 9.8;				// 重力加速度 (m/s^2)

// *** 終端速度の関数 ***
// Re < 1
function v = v1(r)
 v = 2 .* r .^ 2 * (rho_h2o - rho_air) .* g ./ 9 ./ mu
endfunction
// 2 < Re < 500
function v = v2(r)
 v = 2 * r * (4 * (rho_h2o - rho_air) ^ 2 * g ^ 2 / 255 / rho_air / mu) ^ (1 / 3)
endfunction
// 500 < Re < 10E5
function v = v3(r)
 v = sqrt(8 * r * (rho_h2o - rho_air) * g / (3 * 0.44 * rho_air))
endfunction

// *** 雨粒の半径 ***
R1 = logspace(-7,-3);
R2 = logspace(-6,-2);
R3 = logspace(-5,-1);
R = logspace(-7,-1);

// *** 終端速度のプロット ***
xsetech([0,0,1,0.45],[1E-7,1E-6,1E-1,1E3],"ll")
plot(R1,v1(R1),"-r");
plot(R2,v2(R2),"-g");
plot(R3,v3(R3),"-b");
xlabel("Radius (m)")
ylabel("Terminate velocity (m/s)")

// *** レイノルズ数のプロット ***
xsetech([0,0.5,1,0.45],[1E-7,1,1E-1,1E7],"ll")
plot(R,2,"--k");
plot(R,500,"--k");
plot(R,1E5,"--k");
plot(R1,2 .* R1 .* rho_air .* v1(R1) ./ mu,"-r");
plot(R2,2 .* R2 .* rho_air .* v2(R2) ./ mu,"-g");
plot(R3,2 .* R3 .* rho_air .* v3(R3) ./ mu,"-b");
xlabel("Radius (m)")
ylabel("Reynolds number")