// 2007年2月2日に作成する
// Excelで学ぶ物理数学(ナツメ社) 245ページを参考にしました。
// 中田　亨様のホームページを参考にしました。
// http://staff.aist.go.jp/toru-nakata/index-j.html
// 素粒子の新謎解き（リーベル出版）仁平群治著の67ページを参考にしました。

//　四元数を用いて、三次元空間で、回転操作を行ってspin回転動作を描く

#　　t=δ, a=α, b = β, c = γ, s=θ, p=ψ

i@;　// 単位変数であることを宣言する
j@;
k@;

//　回転軸(α, β, γ)
α=0;
β=1;
γ=0;

δ = 0.4; // 回転角度(ラジアン)

// 座標(x, y, z)

r=1; // 半径
θ0 = 0;
θ16;
dθ=π()/16;

ψ0 = 0.0;
ψ40;
dψ = 2π()/40;

IF(1) { // Y10 Y20 では必要、極点近くでは細かく刻む
　　ALL(θ) {
　　　　θ(θ) = (1-COS(θ(θ)))/2*π();
　　};
};

X3;

QUATERNION() { // 四元数の演算を宣言する
　　Q(i) = QUARTER(COS(δ/2), αSIN(δ/2), βSIN(δ/2), γSIN(δ/2));
};

∀(ψ,θ) {
　　x=rSIN(θ(θ))COS(ψ(ψ));
　　y=rSIN(θ(θ))SIN(ψ(ψ));
　　z=rCOS(θ(θ));

// 回転角度(δ)反時計回り

　　QUATERNION() { // 四元数の演算を宣言する
　　　　P(i) = QUARTER(0,x,y,z); // 四元数を生成する
R(i) = QUARTER(COS(θ(θ)/2), 0, 0, SIN(θ(θ)/2));
P(i) = ~R(i)P(i)R(i);// θのπ当たりz軸を中心に2πだけ回転させる
// P(i) = ~R(i)P(i)R(i);
　　　　A(i) = ~Q(i)P(i)Q(i); // y軸を中心にδ/2を２回だけ回転させる
　　};
　　∀(X) {
　　　　Z(ψ,θ,X:0) = ELEMENT(A(i), 0); // 要素を取り出す
　　　　Z(ψ,θ,X:1) = ELEMENT(A(i), 1);
　　　　Z(ψ,θ,X:2) = ELEMENT(A(i), 2);
　　　　Z(ψ,θ,X:3) = ELEMENT(A(i), 3);
　　};

};
MPRT(Z(ψ,θ,X)); // csvファイルに出力する