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);