restart; with(plots):  
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]];
#Define the torus
a := 1: b := 2.5:
K := [cos(u)*(a*cos(v)+b), sin(u)*(a*cos(v)+b), sin(v)]; 
#Compute the Normal
Ku := diff(K,u); Kv := diff(K,v):  Nt := mycross(Ku,Kv);  N := expand(1/sqrt(mydot(Nt,Nt)) * Nt):
dc :=[cos(t),sin(t)+s]; #Define the domain curve to plot
gk := expand(subs(u=dc[1],v=dc[2],K)): #The curve on the surface
soffset := [0,0,4]: #Set the location for the sphere to show the image of the Gauss map
gs := expand(subs(u=dc[1],v=dc[2],N+soffset)): #The image of the Gauss map applied to the curve 

#Now for graphics

S := expand([cos(u)*cos(v), sin(u)*cos(v), sin(v)]+soffset): 
P0 := plot3d(S,u=0..2*Pi,v=-Pi/2..Pi/2,color=green,style=wireframe):  #Assign the plot ofthe sphere in green 
P1 := plot3d(K,u=0..2*Pi,v=0..2*Pi,color=blue,scaling=constrained): P1;  #The original torus

n := 10: for i from 0 to n do  ss := i/n*2*Pi;
C1 := tubeplot(expand(subs(s=ss,gk)),t=0..2*Pi,radius=.1,numpoints=150):
C2 := tubeplot(expand(subs(s=ss,gs)),t=0..2*Pi,radius=.1,numpoints=150):
PS[i] := display({C1,C2,P1,P0});
od:
display([PS[j] $j=0..n],insequence=true);