restart; with(plots):   #Parallel transport
mydot := (v,w) -> sum(v[ii]*w[ii],ii=1..nops(v));
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]];

a := 1; b := 2.5;
K :=  [cos(u)*(a*cos(v)+b), sin(u)*(a*cos(v)+b), sin(v)];  #This is the surface: A torus, using the a,b defined above

#First fundamental form, normal vector and Christoffel Symbols
Ku := diff(K,u); Kv := diff(K,v);  Nt := mycross(Ku,Kv);  N := expand(1/sqrt(mydot(Nt,Nt)) * Nt):
E :=factor( simplify(mydot(Ku,Ku))); F := simplify(mydot(Ku,Kv)); G := simplify(mydot(Kv,Kv));
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;
g121 := (G*Ev-F*Gu)/dd; g122 :=(E*Gu-F*Ev)/dd;
g221 := (2*G*Fv-G*Gu-F*Gv)/dd; g222 := (E*Gv-2*F*Fv+F*Gu)/dd;


#########Define the curve and solve the differential equations for parallel transport
dc :=[t,t+Pi/4];  #This is the path in the domain of K

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) ));
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) ));
ans :=dsolve({eq1=0,eq2=0,mu(0)=1/sqrt(2),lambda(0)=1/sqrt(2)},{lambda(t),mu(t)});
alpha := subs(ans,lambda(t)): beta := subs(ans,mu(t)):
W := (evalf(expand(subs(u=dc[1],v=dc[2],alpha*Ku+beta*Kv)))); #This is the direction of the stick for parallel transport


#The rest is graphics
soffset := [0,0,4]: #Where to center the tails of the stick so their path can be observed
S := expand([cos(u)*cos(v), sin(u)*cos(v), sin(v)]):
P0 := plot3d(expand(1/4*S+soffset),u=0..2*Pi,v=-Pi/2..Pi/2,color=yellow,style=wireframe): #A little Yellow Sphere for a marker
P1 := plot3d(K,u=0..2*Pi,v=0..2*Pi,color=blue,scaling=constrained,style=wireframe): P1;  #The actual surface plotted

gk := expand(subs(u=dc[1],v=dc[2],K)); #This is the curve to draw on the surface
expand(subs(t=tt,gk+s*W)):  #This is the endpoint of the stick
gn := expand(subs(u=dc[1],v=dc[2],N)): #This is the normal vector
aa := 0; bb := Pi;
fud := 1; #Fudge factor for graphics
C0 := tubeplot(gk,t=aa..bb,radius=.1,numpoints=50,color=green):  #The curve on the surface
C00 := tubeplot(expand(fud*W+soffset),t=aa..bb,radius=.01,numpoints=50,color=red,style=wireframe): #The sticks with their tails together   
n := 7; nnn := 4; for i from 0 to n do  tt := i/n*(bb-aa)+aa ;   #n is the number of frames
C1 := tubeplot(expand(subs(t=tt,gk+s*W)),s=0..1,radius=.1,color=red,numpoints=nnn): C1q := tubeplot(expand(subs(t=tt,gk+s*W)),s=-1..0,radius=.1,color=magenta,numpoints=nnn):
C2 := tubeplot(expand(subs(t=tt,gk+s*gn)),s=0..1,radius=.1,numpoints=nnn,color=blue):
C3 := tubeplot(expand(subs(t=tt,soffset+s*W)),s=0..fud,radius=.01,color=red,numpoints=nnn,style=wireframe):
C4 := tubeplot(expand(subs(t=tt,soffset+fud*W+s*gn)),s=0..1,radius=.01,color=blue,numpoints=nnn,style=wireframe):

PS[i] := display({C1,C2,C3,C4,C0,C00,C1q,P0,P1});  od:
display([PS[j] $j=0..n],insequence=false);