clear;

// *** 入力パラメータ ***
m = 1.0;			// 粒子の質量
q = 1.0;			// 粒子の電荷
// 一様磁界
bx = 0;
by = 0;
bz = 3.0;
// 一様電解
ex = 0;
ey = 1.5;
ez = 0;

// *** 解くべき常微分方程式の定義 ***
function dx = em(t,x)
  // dx/dt = vx
  dx(1) = x(2);
  // dy/dt = vy
  dx(3) = x(4);
  // dz/dt = vz
  dx(5) = x(6);
  // dvx/dt = (q/m) * Ex + (q/m) * (vy*Bz-vz*By)
  dx(2) = q * ex / m + q * (x(4) * bz - x(6) * by) / m;
  // dvy/dt = (q/m) * Ey + (q/m) * (vz*Bx-vx*Bz)
  dx(4) = q * ey / m + q * (x(6) * bx - x(2) * bz) / m;
  // dvz/dt = (q/m) * Ez + (q/m) * (vx*By-vz*Bx)
  dx(6) = q * ez / m + q * (x(2) * by - x(4) * bx) / m;
endfunction

// 時間ベクトル
T = linspace(0,20,201);

// 初期条件
x0 = 0;
vx0 = - 1;
y0 = 0;
vy0 = 0;
z0 = 0;
vz0 = 0.5;
X0 = [x0; vx0; y0; vy0; z0; vz0];

// 常微分方程式ソルバ
X = ode(X0,0,T,em);

// グラフのプロット
param3d(X(1,:),X(3,:),X(5,:),,flag=[5,4],ebox=[-5,5,-5,5,0,10]);
param3d(X(1,:),X(3,:),zeros(X(5,:)),flag=[5,4],ebox=[-5,5,-5,5,0,10]);
xlabel("x");
ylabel("y");
zlabel("z");