sticks.mw

 > restart; with(plots):

 > mydot := (v,w) -> sum(v[ii]*w[ii],ii=1..nops(v)); # Define Dot Product

 > mycross := (v,w) -> [v[2]*w[3]-v[3]*w[2], v[3]*w[1] - v[1]*w[3], v[1]*w[2]-v[2]*w[1]]; #Define Cross Product

 >

 > a := 1; b := 2.5;  #Set parameters for the torus: a is inner tube radius, b is the outer radius

 > S :=  [cos(u)*(a*cos(v)+b), sin(u)*(a*cos(v)+b), sin(v)];  #The Torus

 > Su := diff(S,u); Sv := diff(S,v):  Nt := mycross(Su,Sv);  N := expand(1/sqrt(mydot(Nt,Nt)) * Nt):  #Partials and unit normal

 > E :=factor( simplify(mydot(Su,Su))); F := simplify(mydot(Su,Sv)); G := simplify(mydot(Sv,Sv));

 > Eu := diff(E,u); Ev := diff(E,v); Fu := diff(F,u); Fv := diff(F,v); Gu := diff(G,u); Gv := diff(G,v);

 > dd :=2*(E*G-F^2);

 > g111 := (G*Eu-2*F*Fu+F*Ev)/dd; g112 := (2*E*Fu-E*Ev-F*Eu)/dd;  #Christoffel Symbols

 > g121 := (G*Ev-F*Gu)/dd; g122 :=(E*Gu-F*Ev)/dd; #Christoffel Symbols

 > g221 := (2*G*Fv-G*Gu-F*Gv)/dd; g222 := (E*Gv-2*F*Fv+F*Gu)/dd; #Christoffel Symbols

 >

 > ######### Pick a curve on the surface ######

 > cc := 1;

 > dc :=[cc+t,t];  #equation of the curve in the domain

 > #ut := diff(dc[1],t); vt := diff(dc[2],t);

 > ;

 >

 > #Solve paralel transport equations for a vector field lambda(t)*Su + mu(t)*Sv

 > eq1 := eval(subs(v=dc[2],u=dc[1],ut=diff(dc[1],t),vt=diff(dc[2],t),diff(lambda(t),t) +lambda(t)*(g111*ut + g121*vt) + mu(t)*(g121*ut+g221*vt) )); #Equations for parallel transpot

 > eq2 := eval(subs(v=dc[2],u=dc[1],ut=diff(dc[1],t),vt=diff(dc[2],t),diff(mu(t),t) + lambda(t)*(g112*ut+g122*vt) + mu(t)*(g122*ut+g222*vt) )); #Equations for parallel transport

 > ans :=dsolve({eq1=0,eq2=0,mu(0)=1/sqrt(2),lambda(0)=1/sqrt(2)},{lambda(t),mu(t)}); #Solve for the equations. The starting conditions are in here with mu(0) and lambda(0)

 > alpha := subs(ans,lambda(t)); beta := subs(ans,mu(t)); #Extract the solutions

 > ;W := simplify(evalf(expand(subs(u=dc[1],v=dc[2],alpha*Su+beta*Sv)))):  #Substitute to get the vector field for parallel transport

 >

 >

 >

 > #Now do the graphics

 > soffset := [0,0,4];  #Location for the sphere

 > Sph := expand([cos(u)*cos(v), sin(u)*cos(v), sin(v)]): #Sphere

 > P0 := plot3d(expand(1/4*Sph+soffset),u=0..2*Pi,v=-Pi/2..Pi/2,color=yellow,style=wireframe): #Plot of the sphere

 > P1 := plot3d(S,u=0..2*Pi,v=0..2*Pi,color=blue,scaling=constrained,style=wireframe): #Plot of the surface

 >

 > gk := expand(subs(u=dc[1],v=dc[2],S));  #Equation for the curve on S

 >

 >

 >

 > gn := expand(subs(u=dc[1],v=dc[2],N)):  #Equation of the normal vector

 > aa := 0; bb := Pi; #Start and end parameters for the curve

 > fud := 1; #Scale for the length of the parallel vector field in the plot

 > C0 := tubeplot(gk,t=aa..bb,radius=.1,numpoints=50,color=green): #Plot of curve on S

 > C00 := tubeplot(expand(fud*W+soffset),t=aa..bb,radius=.01,numpoints=50,color=red,style=wireframe):  # The plot of the vector field with all tails at the center of the sphere

 > n := 20; nnn := 4; for i from 0 to n do  tt := i/n*(bb-aa)+aa ;  # Loop to move along the curve

 > C1 := tubeplot(expand(subs(t=tt,gk+s*W)),s=0..1,radius=.1,color=red,numpoints=nnn):  #The stick on the surface

 > C2 := tubeplot(expand(subs(t=tt,gk+s*gn)),s=0..1,radius=.1,numpoints=nnn,color=blue):  #The normal on the surface

 > C3 := tubeplot(expand(subs(t=tt,soffset+s*W)),s=0..fud,radius=.01,color=red,numpoints=nnn,style=wireframe): #The stick on the sphere

 >

 > PS[i] := display({C1,C2,C3,C0,C00,P0,P1});  od:

 > display([PS[j] \$j=0..n],insequence=false);

 >

 >

 >