clear;

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

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

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

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

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

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