// **********************************************
// 
//                  共通部分
// 
// **********************************************
clear;
// 表示領域のサイズ
//xymax = 2000;		// 表示領域を2000*2000に設定
xymax = 400;		// 表示領域を400*400に設定
// グラフの縦横比を等しくする
h = scf(); 		// ウィンドウを作成
ha = h.children(1);	// Axes（座標軸）オブジェクトへのハンドルを取得
ha.isoview = "on";	// 座標軸の縦横比を等しくする
ha.data_bounds = [- xymax,- xymax; xymax, xymax]; // 座標軸表示範囲の設定 

// **********************************************
// 
//                α粒子の運動
// 
// **********************************************
// *** 解くべき常微分方程式 ***
function dx = ruther(t,x)
  r = sqrt(x(1) ^ 2 + x(3) ^ 2);
  dx(1) = x(2);
  dx(2) = 0.0614 * x(1) / r ^ 3;
  dx(3) = x(4);
  dx(4) = 0.0614 * x(3) / r ^ 3;
endfunction

// *** 時間ベクトル ***
T = linspace(0,100000,1000);

// *** α粒子の初速度 ***
v0 = 0.053;
vx0 = v0;
vy0 = 0;

// *** α線源の位置 ***
// x座標
x0 = - 2000;

// y座標を変えながらα粒子の運動を計算
for y0 = -95:10:95
  X0 = [x0; vx0; y0; vy0];
  X = ode(X0,0,T,ruther);
  plot(X(1,:),X(3,:),'-r');
end

// 描画範囲の設定
zoom_rect([- xymax,- xymax, xymax, xymax]);

// **********************************************
// 
//              α粒子の受ける力
// 
// **********************************************
// *** 力を計算する位置(空間グリッド)の設定 ***
nxy = 20;	// グリッドの分割数
x = linspace(- xymax, xymax, nxy);
y = linspace(- xymax, xymax, nxy);
[X,Y] = ndgrid(x,y);

// *** 受ける力を計算するための関数 ***
// x方向の力
function fx = fx(x,y)
  r = sqrt(x .^ 2 + y .^ 2);
  fx = 2 * 79 * 200 / 137 .* x ./ (r .^ 3 + %eps);
endfunction
// y方向の力
function fy = fy(x,y)
  r = sqrt(x .^ 2 + y .^ 2);
  fy = 2 * 79 * 200 / 137 .* y ./ (r .^ 3 + %eps);
endfunction

// *** 力の計算とプロット ***
// α粒子の受ける力の計算
FX = fx(X,Y);
FY = fy(X,Y);
// α粒子の受ける力のベクトルをプロット
champ(x,y,FX,FY);