with(LinearAlgebra):
with(plots):
with(RealDomain):

printf("This loads LinearAlgebra, plots, RealDomain, AngleOfRotation and AngleOfRotation2 of Rot([1,0,0],theta), Rot(u,theta), u in S^2, DotProd(X,Y), X,Y in S^2, Mul(u,v,theta) = uRot(v,theta), dualMul(u, v, theta) is the dual of Mul(u,v,theta),  PlotSphere() = yellow unit sphere, PlotArc(X,Y) where X,Y in S^2, PlotLine(X,Y) (straight line X to Y), PlotSphericalTriangle(X,Y,Z), ConvertToRange(S::{set,list}),ROT(u) is a rotation A such that eA = u where  e = (1,0,0) and u = (u1,u2,u3) \n"): 

ConvertToRange:=proc(S::{set,list})
local T, T1;
T:=map(t->[op(t)],S);
if nops(T)=2 then
T1:=seq(op(remove(s->type(s,name),T[i])), i=1..2);
return min(T1)..max(T1) else
if type(op(T)[1],name) then return -infinity..op(T)[2] else op(T)[1]..infinity fi; fi;
end proc:

AngleOfRotation:=proc(A)
local tr,cosa,u,a;
u:=Vector[row]([0,1,0]).A;
if evalf(u[3]) < 0 then a:=-1; else a:= 1; fi;
tr:=Trace(A);
cosa:=simplify((tr-1)/2);
if a < 0 then return 2*Pi-arccos(cosa); else return arccos(cosa); fi;
end proc:

AngleOfRotation2:=proc(A)
local tr,cosa,u,a;
u:=Vector[row]([0,1,0]).A;
arctan(u[3],u[2]);
end proc:

#The following from http://inside.mines.edu/fs_home/gmurray/ArbitraryAxisRotation/
Rotate:=proc(X, A, V, theta)
local x,y,z,a,b,c,u,v,w;
description "rotation of X about line thru A in direction of UNIT VECTOR V this was taken from  http://inside.mines.edu/fs_home/gmurray/ArbitraryAxisRotation/";
x,y,z:=seq(X[i],i=1..3);
a,b,c:=seq(A[i],i=1..3);
u,v,w:=seq(-V[i],i=1..3); 
[( a*(v^2+w^2)+u*(-b*v -c*w+u*x+v*y+w*z)
+( (x-a)*(v^2+w^2) + u*(b*v + c*w - v*y - w*z) )*cos(theta) 
 +(b*w - c*v - w*y + v*z)*sin(theta))/(1), 
( b*(u^2+w^2)+v*(-a*u-c*w+u*x+v*y+w*z)
 +( (y-b)*(u^2+w^2) + v*(a*u + c*w - u*x - w*z) )*cos(theta)
 +(-a*w + c*u + w*x - u*z)* sin(theta))/(1), 
(c*(u^2+v^2)+w*(-a*u-b*v+u*x+v*y+w*z)+( (z-c)*(u^2+v^2) + 
w*(a*u + b*v - u*x - v*y))*cos(theta)+
( (a*v - b*u - v*x + u*y)*sin(theta)))/
(1)];
map(simplify,%);  
end proc:

#The following returns the matrix of the rotation by theta about the unit vector U

Rot:=proc(U,theta)
local X;
Matrix(map(X->Rotate(X,[0,0,0],U,theta),[[1,0,0],[0,1,0],[0,0,1]]));
Transpose(%);
map(simplify,%);
end proc:

#Mul(u,v,theta) is the "product" of u,v in S^2 relative to x = Rot([1,0,0],theta)

Mul:=proc(u,v,theta)
Vector[row](u).Rot(v,theta);
return simplify(convert(simplify(%),list));
end:

#dualMul(u,v,theta) is the dual of Mul(u,v,theta)

dualMul:=proc(u, v, theta) 
Vector[row](u).Rot(v,-theta);
return simplify(convert(simplify(%),list));
end proc:

DotProd:=proc(X,Y)
local i;
add(X[i]*Y[i],i=1..3);
simplify(%);
end:

#The following procedure ROT(u) gives a rotation U such that eU = u where  e = (1,0,0) and u = (u1,u2,u3)
ROT:=proc(u)
local e,a,w,nw,NW;
e:=[1,0,0];
a:=arccos(DotProd(u,e));
w:=convert(CrossProduct(Vector(u),Vector(e)),list);
nw:=sqrt(DotProd(w,w));
NW:=expand((1/nw)*w);
Rot(NW,a);
end proc:


PlotSphere:=proc(OPT)
local t,p;
plot3d(1,t=0..2*Pi,p=0..Pi,coords=spherical,scaling=constrained, style=surface,transparency = 0.6,OPT):
end proc:

PlotArc:=proc(X,Y)
local Z,t,W,a,theta;
Z:=CrossProduct(CrossProduct(Vector(X),Vector(Y)),Vector(X));
Z:=convert(Z,list);
a:=sqrt(DotProd(Z,Z));
Z:=expand((1/a)*Z);
theta:=arccos(DotProd(X,Y));
plot3d(expand(cos(t)*X+sin(t)*Z),t=0..theta,color="Red");
end proc:

PlotLine:=proc(X,Y)
pointplot3d([X,Y],style=line, color = "Green");
end proc:

PlotSphericalTriangle:=proc(X,Y,Z,s)
display([PlotSphere(color=yellow),PlotArc(X,Y),PlotArc(X,Z),PlotArc(Y,Z)],color=s);
end proc: