##This Maple sheet proofs the equivalence of (2) and (4) in the Eurographics 2013 paper
# Global Selection of Stream Surfaces
#
# J. Martinez Esturo, M. Schulze, C. Rössl, H. Theisel
# University of Magdeburg

restart;
with(LinearAlgebra):

#Jacobian matrix
J := Matrix([[u_x,u_y,u_z],
             [v_x,v_y,v_z],
             [w_x,w_y,w_z]]);

#Velocity field
V := <u,v,w>;

#Tangent of seed curve
Sdot := <usdot,vsdot,wsdot>;

#Second derivative vector of seed curve
Sddot := <usddot,vsddot,wsddot>;

#Without loss of generality we assume than V is in the direction of the x-axis, and that Sdot is in the x-y plane.
#This gives a normal in the direction of the z-axis
v := 0;
w := 0;
wsdot := 0;

#We use this assumption to set the orientation of the normal, which simplifies computations later on. Any alternative of
#assume(u>0,vsdot>0);
#assume(u<0,vsdot<0);
#assume(u<0,vsdot>0);
#works equally well for the proof.
assume(u>0,vsdot<0);

#First and second order partials of stream surface
Xs := Sdot;
Xt := V;

Xss := Sddot;
Xst := Multiply(J,Sdot);
Xtt := Multiply(J,V);

#Classic differential geometric coefficients of first fundamental form
E := Multiply(Transpose(Xs),Xs);
F := Multiply(Transpose(Xs),Xt);
G := Multiply(Transpose(Xt),Xt);

#Unnormalized normal
HNormal := CrossProduct(Xs,Xt);

#Unit normal
NNormal := HNormal / sqrt(HNormal[1]^2 + HNormal[2]^2 + HNormal[3]^2);

#More classic differential geometric coefficients of second fundamental form
L := simplify(Multiply(Transpose(NNormal),Xss));
M := simplify(Multiply(Transpose(NNormal),Xst));
N := simplify(Multiply(Transpose(NNormal),Xtt));

#Now we compute the principal directions P1,P2.
K := Matrix([[dt^2,-ds*dt,ds^2],[E,F,G],[L,M,N]]);

dt1 := eval(solve(Determinant(K)=0,dt)[1],ds=ds1);
dt2 := eval(solve(Determinant(K)=0,dt)[2],ds=ds2);

ds1 := 1;
ds2 := 1;

#Unnormalized principal directions
HP1 := ds1*Xs + dt1*Xt;
HP2 := ds2*Xs + dt2*Xt;

#Principal directions
P1 := HP1 / sqrt(HP1[1]^2 + HP1[2]^2 +HP1[3]^2);
P2 := HP2 / sqrt(HP2[1]^2 + HP2[2]^2 +HP2[3]^2);

#Compute the angle alpha
#Normalized velocity
Vn := V / Norm(V,2);

#Squared cosine of angle between between V and P1:
cosp1q := Multiply(Transpose(P1),Vn)^2;

#Squared sine of angle between V and P1:
sinp1q := Multiply(Transpose(P2),Vn)^2;

#Compute the principal curvatures k1,k2
k1 := solve((k*E-L)*(k*G-N) - (k*F-M)^2,k)[1];
k2 := solve((k*E-L)*(k*G-N) - (k*F-M)^2,k)[2];

#Now we have everything set up to compute e_a using (2)
#Our error function
e_a_2 := simplify(cosp1q*sinp1q*(k2-k1)^2);

#We can simplify e_a (2)
e_a_2s := factor(expand(numer(e_a_2))/expand(denom(e_a_2)));

#############################################

#Now we compute e_a using (4)

B := CrossProduct(V, NNormal);

e_a_4 := (Multiply(Multiply(Transpose(NNormal),J),B))^2 / Norm(V,2)^4;

#############################################
#We compare e_a_4 and e_a_3s, it should be zero to complete the proof
simplify(e_a_2s - e_a_4);