clear;

// *** 入力パラメータ ***
r = 1E-3;			// 雨粒の半径 (m)

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

// *** 運動方程式 ***
m = rho_h2o * 4 * %pi * r ^ 3 / 3;	// 雨粒の重さ (g)
//k = 0.235 * mu * %pi * r * r / nu;	// 空気抵抗の係数 (g/s)
k = 0.22 * %pi * rho_air * r ^ 2;	// 空気抵抗の係数 (g/s)
v0 = 0;					// 初期速度 (m/s)
// 解くべき方程式の定義
deff("vdot = df(t,v)", "vdot = g - (k / m) * v ^ 2");

// *** 常微分方程式の計算 ***
// 時間ベクトル(s)
T = linspace(0, 5, 101);
// 常微分方程式の数値解
V = ode(v0, 0, T, df);

// *** 速度の描画 ***
// 数値計算(赤丸)
plot(T, V, 'or');
// 解析解(緑破線)
plot(T, sqrt(m * g / k) * tanh(sqrt(k * g / m) .* T), '--g');

// *** グラフの体裁 ***
// 軸ラベル
xlabel("Time (s)");
ylabel("Velocity (m/s)")
// 凡例
legend(["Numerical";"Analytical"],4)