// 微分幾何 (朝倉書店) 細野忍著 4000円+税

// 曲面x(u,v)の接平面が決まる様子 14ﾍﾟｰｼﾞ図1.6

# t=δ, a=α, b = β, c = γ, e=θ

//　回転軸(α, β, γ)、 α^2+β^2+γ^2=1である必要がある
α=0;
β=0;
γ=1;

//δ = 0.2; // 回転角度(ラジアン)
δ = 0.5; // 回転角度(ラジアン)
//δ = 0; // 回転角度(ラジアン)

u0 = 0;
u25;
du = π()/25;

v0 = 0;
v25;
dv = 2π()/25;

x4;

QUATERNION() { // 四元数の演算を宣言する
　　Q() = QUARTER(COS(δ/2), αSIN(δ/2), βSIN(δ/2), γSIN(δ/2));
};
_X(u,v) { // 立体射影座標
　　RETURN(2u/(1+u^2+v^2));
};
_Y(u,v) {
　　RETURN(2v/(1+u^2+v^2));
};
_Z(u,v) {
　　RETURN((-1+u^2+v^2)/(1+u^2+v^2));
};
_XX(u,v) { // 極座標
　　RETURN(SIN(u)COS(v));
};
_YY(u,v) {
　　RETURN(SIN(u)SIN(v));
};
_ZZ(u,v) {
　　RETURN(COS(u));
};
∀(x,u,v) {
// 回転角度(δ)反時計回り

　　QUATERNION() { // 四元数の演算を宣言する
　　　　P() = QUARTER(0,_XX(u(u),v(v)),_YY(u(u),v(v)),_ZZ(u(u),v(v))); // 局座標で球面の四元数を生成する
//　　　　P() = QUARTER(0,_X(u(u),v(v)),_Y(u(u),v(v)),_Z(u(u),v(v))); // 南半球を立体射影座標で四元数を生成する
　　　　A() = ~Q()P()Q(); // δ/2を２回だけ回転させる
　　};
　　R(x:0,u,v) = ELEMENT(A(), 0); // 要素を取り出す
　　R(x:1,u,v) = ELEMENT(A(), 1);
　　R(x:2,u,v) = ELEMENT(A(), 2);
　　R(x:3,u,v) = ELEMENT(A(), 3);
　　R(x:4,u,v) = _XX(u(u),v(v))^2+_YY(u(u),v(v))^2+_ZZ(u(u),v(v))^2; // 常に1であることを確認する
};

TITLE(u,"\n%s=%.1f","VAL");
TITLE(v,"%s=%.1f","VAL");
OUT(R(x,v,u),R(x,u,v));
TITLE(u,"%s=%.1f","VAL");
TITLE(v,"\n%s=%.1f","VAL");
SKIPTITLE();
OUT(R(x,u,v),R(x,u,v));
#COMMENT0
p0 = -1;
p5;
dp = 0.4;
q0 = -1;
q5;
dq = 0.4;
u = 0.5;
v = 0.5;
h = 0.001;

T?;
T0 = _XX(u,v); // 接点の位置を求める
T1 = _YY(u,v);
T2 = _ZZ(u,v);
T3 = (_XX(u+h,v)-_XX(u,v))/h; // x/∂u
T4 = (_YY(u+h,v)-_YY(u,v))/h; // y/∂u
T5 = (_ZZ(u+h,v)-_ZZ(u,v))/h; // z/∂u
T6 = (_XX(u,v+h)-_XX(u,v))/h; // x/∂v
T7 = (_YY(u,v+h)-_YY(u,v))/h; // y/∂v
T8 = (_ZZ(u,v+h)-_ZZ(u,v))/h; // z/∂v

∀(x,p,q) {
// 回転角度(δ)反時計回り

　　QUATERNION() { // 四元数の演算を宣言する
　　　　P() = QUARTER(0,T0+T3*p(p)+T6*q(q),T1+T4*p(p)+T7*q(q),T2+T5*p(p)+T8*q(q)); // 四元数を生成する
　　　　A() = ~Q()P()Q(); // δ/2を２回だけ回転させる
　　};
　　S(x:0,p,q) = ELEMENT(A(), 0); // 要素を取り出す
　　S(x:1,p,q) = ELEMENT(A(), 1);
　　S(x:2,p,q) = ELEMENT(A(), 2);
　　S(x:3,p,q) = ELEMENT(A(), 3);
};

SKIPTITLE();
TITLE(p,"\n%s=%.1f","VAL");
TITLE(q,"%s=%.1f","VAL");
OUT(S(x,q,p),S(x,p,q));
TITLE(p,"%s=%.1f","VAL");
TITLE(q,"\n%s=%.1f","VAL");
SKIPTITLE();
OUT(S(x,p,q),S(x,p,q));

W?;

x = SIN(u)COS(v); // 極座標(u,v)から位置(x,y,z)を求める
y = SIN(u)SIN(v);
z = COS(u);

U = x/(1-z); // 南半球上のアドレス（座標、パラメータ）を求める
V = y/(1-z);

W0 = _X(U,V); // 接点の位置を求める
W1 = _Y(U,V);
W2 = _Z(U,V);
W3 = (_X(U+h,V)-_X(U,V))/h; // x/∂U
W4 = (_Y(U+h,V)-_Y(U,V))/h; // y/∂U
W5 = (_Z(U+h,V)-_Z(U,V))/h; // z/∂U
W6 = (_X(U,V+h)-_X(U,V))/h; // x/∂V
W7 = (_Y(U,V+h)-_Y(U,V))/h; // y/∂V
W8 = (_Z(U,V+h)-_Z(U,V))/h; // z/∂V
W9 = W0^2+W1^2+W2^2;
TRACE("接点の値%f", W9);
∀(x,p,q) {
// 回転角度(δ)反時計回り

　　QUATERNION() { // 四元数の演算を宣言する
　　　　P() = QUARTER(0,W0+W3*6p(p)+W6*6q(q),W1+W4*6p(p)+W7*6q(q),W2+W5*6p(p)+W8*6q(q)); // 四元数を生成する
　　　　A() = ~Q()P()Q(); // δ/2を２回だけ回転させる
　　};
　　S(x:0,p,q) = ELEMENT(A(), 0); // 要素を取り出す
　　S(x:1,p,q) = ELEMENT(A(), 1);
　　S(x:2,p,q) = ELEMENT(A(), 2);
　　S(x:3,p,q) = ELEMENT(A(), 3);
};

SKIPTITLE();
TITLE(p,"\n%s=%.1f","VAL");
TITLE(q,"%s=%.1f","VAL");
OUT(S(x,q,p),S(x,p,q));
TITLE(p,"%s=%.1f","VAL");
TITLE(q,"\n%s=%.1f","VAL");
SKIPTITLE();
OUT(S(x,p,q),S(x,p,q));
#END

TRACE("T(%f,%f,%f) W(%f,%f,%f)", T0,T1,T2,W0,W1,W2);