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

>

>

>