> #sfb 15 mar 04
> # modeled after the matlab solution
> #plotsetup(x11);
> L:=6;
> Range:=x=-L..L,y=-L..L;
> fun:=x^2-y^2;
> with(plots):with(linalg):
> # geodesic starts at [x0,y0,fun(x0,y0)] and moves int the
> # v = [v1, v2, ?] direction (in the tangent plane
> # h is the step size, each point is `h' away at least in tangent (length tn)
> # jmax is to prevent infinite loops
> # L is the limit the curve to the graph
> # eventually r is the old point, q is the current point and
> # p is the next point.
> # fx, fy are parials, n is the normal at q
> # geodesic returns the curve as a list of points
> #[[x0,y0,z0],[x1,y1,z1],...[xn,yn,zn]]
> geodesic:=proc(x0,y0,v1,v2,h,jmax,limit)
> local r, q, p, fx, fy, n, t, j, v, tn, p2, t2, tmin, points, fast;
> fast:=array(1..jmax+1);
> q:=[x0,y0,eval(fun,{x=x0,y=y0})]; # current point no old
> fx:=diff(fun,x); fy:=diff(fun,y);
> n:=[eval(fx,{x=x0,y=y0}),eval(fy,{x=x0,y=y0}),-1]; # normal at q
> v:=[v1, v2, n[1]*v1+n[2]*v2];
> tn:=h*v/norm(v,2);
> p:=evalf(q+tn); # the new point;
> j:=0;fast[j+1]:=q;
> while ( (j < jmax) and (abs(q[1]) < L) and (abs(q[2])<L) )
> do
>  #project p back to the surface z = fun(x,y)
>  p2:=solve({fun=z,x=p[1]+n[1]*t,y=p[2]+n[2]*t,z=p[3]+n[3]*t},{x,y,z,t});
>  # this would be easy if the solution was unique, if not find
>  # the nearest solution.
>  #print(p2);
>  if (nops([p2])>1) then
>   #more than one solution
>   tmin:=eval(t,p2[1]);
>   p:=eval([x,y,z],p2[1]);
>   for i from 2 to nops([p2]) do
>    t2:=eval(t,p2[i]);
>    if ( evalb(eval(abs(t2)<abs(tmin)))) then
>     # better solution
>     tmin:=t2;
>     p:=eval([x,y,z],p2[i]);
>    fi;
>   od;
>  else
>   # exactly one solution
>   p:=eval([x,y,z],p2);
>  fi;
>  #print(eval(fun,{x=p[1],y=p[2]})-p[3]);
>  fast[j+1]:=q;
>  # rotate points
>  r:=q; q:=evalf(p);
>  n:=[eval(fx,{x=q[1],y=q[2]}),eval(fy,{x=q[1],y=q[2]}),-1];
>  # project r onto the tangent plane, solution is unique this time
>  p2:=solve({n[1]*x+n[2]*y+n[3]*z=n[1]*q[1]+n[2]*q[2]+n[3]*q[3],
>     x=r[1]+n[1]*t,y=r[2]+n[2]*t,z=r[3]+n[3]*t},{x,y,z,t});
>  v:=q-eval([x,y,z],p2);
>  tn:=h*v/norm(v,2);
>  p:=evalf(q+tn);
>  j:=j+1;
>  print(j);
> od;
> points:=[seq(fast[i],i=1..j)];
> RETURN(points);
> end proc;
> 
> curve1:=geodesic(1,0,-1/3,1,0.3,50,L):
> curve2:=geodesic(-1,0,1/3,1,0.3,50,L):
> 
> display(plot3d(fun,Range,style=wireframe,shading=z),
>  spacecurve([t,t,0],t=-L..L,color=black),
>  spacecurve([t,-t,0],t=-L..L,color=black),
>  pointplot3d(curve1,color=red,connect=true,thickness=3),
>  pointplot3d(curve2,color=green,connect=true,thickness=3),
>  pointplot3d(curve1,color=black,symbol=cross,symbolsize=18),
>  pointplot3d(curve2,color=black,symbol=cross,symbolsize=18)
> );
>