clear;

// *** 入力パラメータ ***
m = 1;			// 物体の質量 (g)
b = 0.15;		// 速度の2乗に比例する空気抵抗の係数 (g/s)
g = 9.8;		// 重力加速度 (m/s^2)
v0 = 10;		// 初期速度 (m/s)
deg = 30;		// 角度

// 解くべき方程式の定義
function dx = magnus(t,x)
  // dx/dt = vx
  dx(1) = x(2)
  // dy/dt = vy
  dx(3) = x(4)
  // dvx/dt = - (sw/m) * vy - (b/m) * vx * √(vx^2 + vy^2)
  dx(2) = - sw * x(4) / m - b / m * x(2) * sqrt(x(2) ^ 2 + x(4) ^ 2)
  // dvy/dt = (sw/m) * vx - (b/m) * vy * √(vx^2 + vy^2) - g
  dx(4) = sw * x(2) / m - b / m * x(4) * sqrt(x(2) ^ 2 + x(4) ^ 2) - g
endfunction

// 度からラジアンへの変換
function ret = torad(angle)
  ret = %pi * angle / 180;
endfunction

// *** 常微分方程式の計算 ***
// 時間ベクトル
T = linspace(0,5,1000);	// 2(s)後まで計算
th = torad(deg);		// 度からラジアンへの変換
// 常微分方程式の数値解
for k = 1:3 do
  sw = 1.5 * k;			// 比例定数と回転ベクトルの積をパラメータとして変化させる
  vx0 = v0 * cos(th);		// 水平方向(x)の初速度
  vy0 = v0 * sin(th);		// 垂直方向(y)の初速度
  X0 = [0; vx0; 0; vy0];	// 初期条件(位置と速度)ベクトル
  X = ode(X0, 0, T, magnus);	// 常微分方程式ソルバ
  plot2d(X(1,:),X(3,:),k);
end

// グラフの体裁
zoom_rect([0, 0, 10, 5]);	// 地面まで落下したあとは表示しない
xlabel("x position (m)");
ylabel("y position (m)");