#Maple proof that W_s defined by (26), (38), (46) is equivalent
#we restrict ourselves to 2D trajectories, the proof for 3D is similar

restart;

with(plots):
with(LinearAlgebra):

# we write the matrix H and Hdot in (31) in first order Taylor approximation. Note that 
# H and Hdot defined by (28)-(31) have a vanishing last row.
H := Matrix([
[h11,h12,h13], 
[h21,h22,h23],
[0,0,0]])
+ t*Matrix([
[h11_t,h12_t,h13_t], 
[h21_t,h22_t,h23_t],
[0,0,0]]);
Hdot := Matrix([
[h11_t,h12_t,h13_t], 
[h21_t,h22_t,h23_t],
[0,0,0]]);


I_z := Matrix([
[1,0,0], 
[0,1,0]]);

#32)
J := Multiply(I_z,Multiply(H, Transpose(I_z)));
Jdot := Multiply(I_z,Multiply(Hdot, Transpose(I_z)));

#(25),(33)
S := simplify((J + Transpose(J))/2);
Sdot := simplify((Jdot + Transpose(Jdot))/2);
W := simplify((J - Transpose(J))/2);

#we compute the Matrix E used from (34) on
EV := Eigenvectors(S);
HE := EV[2];
E := Matrix([
[HE[1,1]/sqrt(HE[1,1]^2+HE[2,1]^2), HE[1,2]/sqrt(HE[1,2]^2+HE[2,2]^2)], 
[HE[2,1]/sqrt(HE[1,1]^2+HE[2,1]^2), HE[2,2]/sqrt(HE[1,2]^2+HE[2,2]^2)]]);

#compute time partial of E
E_t := Matrix([
[diff(E[1,1],t),diff(E[1,2],t)], 
[diff(E[2,1],t),diff(E[2,2],t)]]):

#since H was given in Taylor expansion around t=0, and all derivatives 
#values are computed, it is sufficient to continue the proof for t=0
t := 0;

#(34),(35)
Sbar := Multiply(Multiply( Transpose(E), S),E);
Sdotbar := Multiply(Multiply( Transpose(E), Sdot),E);

#(37)
u3 := Sdotbar[2,1] / (Sbar[1,1]-Sbar[2,2]);

#(36)
Wbar_s := Matrix([
[0,-u3], 
[u3,0]]);

#(38)
W_s := Multiply(Multiply( E, Wbar_s), Transpose(E));

#(26) computation of W_s by considering strain rotation rate tensor
W_s_26 := Multiply( -E, Transpose(E_t));


#####now computation W_s by unsteadiness minimization

#(23),(41)
VS := Multiply(H , Vector([x,y,1]));
VTS := Multiply(Hdot , Vector([x,y,1]));
V := simplify(Vector( [VS[1],VS[2]] ));
V_t := simplify(Vector( [VTS[1],VTS[2]] ));

#compute M (under (42))
Xp := Vector([-y,x]);
Vp := Vector( [  -V[2],  V[1]]); 
MS := simplify(Multiply( -J,Xp) + Vp);
M := Matrix([
[ MS[1], J[1,1],J[1,2], Xp[1], 1,0 ], 
[ MS[2], J[2,1],J[2,2], Xp[2], 0,1 ]]);

#Vector U of unknowns in (43)
U := Vector(6,symbol=u);

#(45): err is the integral in (38) to be minimised
ERR := V_t - Multiply( M,U);
err := int(int( Multiply(Transpose(ERR),ERR) , x=x0-d..x0+d),y=y0-d..y0+d);

#minimize (45) for unknown U
sol1 := solve({
diff(err,U[1])=0,
diff(err,U[2])=0,
diff(err,U[3])=0,
diff(err,U[4])=0,
diff(err,U[5])=0,
diff(err,U[6])=0},
{U[1],U[2],U[3],U[4],U[5],U[6]});

solve({
diff(err,U[1])=0,
diff(err,U[2])=0,
diff(err,U[3])=0,
diff(err,U[4])=0,
diff(err,U[5])=0,
diff(err,U[6])=0,
U[5]=0,
U[6]=0},
{U[1],U[2],U[3],U[4],U[5],U[6]});


#(46)
W_s_46 := -Matrix([
[ 0, eval(U[1],sol1)  ], 
[ -eval(U[1],sol1),0]]);

##

#to show the equivalence (26), (38), (46), we have to show
#W_s = W_s_26 = W_s_46
simplify( W_s - W_s_26);
simplify( W_s - W_s_46);
#The last two expressions give zero matrices, proving the equivalence of 
#(26), (38), (46)

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

#Maple proof that the 3 trajectories in (50)-(54) have a closed-form solution
#for trv written in (55)


restart;
with(LinearAlgebra):

#(50)-(52)
Y1f := (1/10)*Vector([3*cos(q*t),2*sin(q*t)]);
Y2f := (2/25)*Vector([3*cos(q*t - (2/3)*Pi),2*sin(q*t - (2/3)*Pi)]);
Y3f := (3/25)*Vector([3*cos(q*t + (2/3)*Pi),2*sin(q*t + (2/3)*Pi)]);

#(54)
Qf := Matrix([
[cos(p*t) , -sin(p*t)], 
[sin(p*t) ,  cos(p*t)]]);
Bf := Vector([cos(t),sin(t)]);

#(53)
X1f := Multiply( Qf, Y1f) + Bf;
X2f := Multiply( Qf, Y2f) + Bf;
X3f := Multiply( Qf, Y3f) + Bf;




#derivatives of 3 trajectories
X1f_t := Vector([diff(X1f[1],t) , diff(X1f[2],t)]);
X2f_t := Vector([diff(X2f[1],t) , diff(X2f[2],t)]);
X3f_t := Vector([diff(X3f[1],t) , diff(X3f[2],t)]);

X1f_tt := Vector([diff(X1f_t[1],t) , diff(X1f_t[2],t)]);
X2f_tt := Vector([diff(X2f_t[1],t) , diff(X2f_t[2],t)]);
X3f_tt := Vector([diff(X3f_t[1],t) , diff(X3f_t[2],t)]);


#(28)
X := Matrix([
[X1f[1], X2f[1], X3f[1]], 
[X1f[2], X2f[2], X3f[2]],
[1     , 1     , 1    ]]);

#(29)
Xdot := Matrix([
[X1f_t[1], X2f_t[1], X3f_t[1]], 
[X1f_t[2], X2f_t[2], X3f_t[2]],
[0     , 0     , 0    ]]);

#(30)
Xdotdot := Matrix([
[X1f_tt[1], X2f_tt[1], X3f_tt[1]], 
[X1f_tt[2], X2f_tt[2], X3f_tt[2]],
[0     , 0     , 0    ]]);

#inverse of X
XI := simplify(MatrixInverse(X));

#(31)
H := simplify(Multiply(Xdot, XI));

Hdot := simplify(
Multiply(Xdotdot , XI)
-
Multiply(
Multiply(Xdot , XI),
Multiply(Xdot , XI))
);


I_z := Matrix([
[1,0,0], 
[0,1,0]]);

#(32)
J := Multiply(I_z,Multiply(H, Transpose(I_z)));
Jdot := Multiply(I_z,Multiply(Hdot, Transpose(I_z)));

#(25),(33)
S := simplify((J + Transpose(J))/2);
Sdot := simplify((Jdot + Transpose(Jdot))/2);
W := simplify((J - Transpose(J))/2);

#we compute the Matrix E from (34) 
EV := Eigenvectors(S);
HE := EV[2];
E := Matrix([
[HE[1,1]/sqrt(HE[1,1]^2+HE[2,1]^2), HE[1,2]/sqrt(HE[1,2]^2+HE[2,2]^2)], 
[HE[2,1]/sqrt(HE[1,1]^2+HE[2,1]^2), HE[2,2]/sqrt(HE[1,2]^2+HE[2,2]^2)]]);

#compute time partial of E
E_t := Matrix([
[diff(E[1,1],t),diff(E[1,2],t)], 
[diff(E[2,1],t),diff(E[2,2],t)]]):


#(34),(35)
Sbar := Multiply(Multiply( Transpose(E), S),E);
Sdotbar := Multiply(Multiply( Transpose(E), Sdot),E);

#(37)
u3 := simplify(Sdotbar[2,1] / (Sbar[1,1]-Sbar[2,2]));

#(29)
Wbar_s := Matrix([
[0,-u3], 
[u3,0]]);

#(36)
W_s := simplify(Multiply(Multiply( E, Wbar_s), Transpose(E)));

#(27)
trv := simplify(W-W_s);
#this shows (55)