BeginPackage["geostuff`",{"Graphics`Colors`","Graphics`Shapes`"}]​​
chris::usage="chris[glower,gupper,{u,v}] computes the Christoffel symbols for a 2D surface, whose coordinates are given by u and v, whose metric is given by glower. The upper indices of the metric (Inverse) are given by gupper.";
equ::usage="Gives the geodesics equations. Current version assumes two coordinates in x and y. equ[Gamma, {x,y}] gives the equations without initial conditions. equ[Gamma, pt,{x,y}][vector] gives the geodesic equations with initial conditions at pt in direction vector.";​​
geosurf::usage="geosurf[f,pt,{x,y}][vector] gives a Table of points (a preLine) lying on the geodesic from point in the direction vector on the surface as the graph of f. The current version assumes that the function f is in the coordinates x and y.";
geoAndsurf::usage="geoAndsurf[f,range,{x,y}][vector] returns {Graphics3D[geo], Plot3D[surface],options} where geo is the geodesic from point {0,0} in the direction vector on the surface as the graph of f. The current version assumes that the function f is in the coordinates x and y, and the starting point of the geodesic is {0,0,f[0,0]}. Use Show[##,PlotRange->Reasonable]&@@ because the plot range is not good.";
metr::usage="metr[f,{x,y}] gives the metric for the embedded surface, the graph of f, in the variables x and y.";​​
AllGeo::usage=" AllGeo[{x, y}, ind] returns a graphic object showing the geodesic on surface ind, that is the exponentitation of the tangent vector )x,y_.(-3, Sphere), (-2, RP2), (-1, Donut), (0, Flat), (1, HyperPlane).";
DonutGeodesic::usage="DonutGeodesic[x, y] returns a Graphics3D object. A geodesic in the direction (x,y) from a point on the torus, as the surface of a donut, starting at (0,0)=(4,0,0).";
FlatGeodesic::usage="FlatGeodesic[x, y] returns a Graphics object. The geodesic in the direction (x,y) from the middle of the flat square torus.";
HyperPlaneGeodesic::usage="HyperPlaneGeodesic[x, y] returns a Graphics object. The geodesic in the hyperbolic plane originating at at (0,1) and in the direction of (x,y).";
RP2Geodesic::usage="RP2Geodesic[x, y] returns a Graphics3D object. The geodesic in the direction of (x,y) on the real projective plane (RP2), displayed as the northern hemisphere.";
SphereGeodesic::usage="SphereGeodesic[x,y] returns a Graphics3D object. The geodesic in the direction of the tangent vector (x,y) from the north pole on the sphere. More precisely it displays the image of the exponential map of the line segment from (0,0) to (x,y).";​​
Begin["Private`"]​​

Special Surfaces

scaleSphere=1.02​​scaleTorus=1.02

functions for sphere

expsphere[{x_,y_}]:=Module[{r=(x^2+y^2)^(1/2)},If[r0,Return[{0,0,1}]];​​{x/rSin[r],y/rSin[r],Cos[r]}]
​​geosphere[x_,y_,n_Integer]:=Module[{r=(x^2+y^2)^(1/2),k},If[r≥2Pi,Return[greatcircle[x/r,y/r,n]]];​​k=Ceiling[rn];​​expsphere/@Table[t{x,y},{t,0,1,1/k}]]
greatcircle[x_,y_,n_Integer]:=expsphere/@Table[t{x,y},{t,0,2Pi,Pi/(3n)}]

functions for flat torus

expflat[{x_,y_}]:={x-otherfloor[x],y-otherfloor[y]}
conpair[{{a_,b_},{c_,d_}}]:=Module[{m,n},m=Min[otherfloor@a,otherfloor@c];​​n=Min[otherfloor@b,otherfloor@d];​​{{a-m,b-n},{c-m,d-n}}]
otherfloor[x_]:=(*Thisfunctionisusedtofixround-offerrors.*)Module[{ans},ans=If[NonNegative[x],IntegerPart[x],IntegerPart[x]-1];​​If[x-ans>.999,ans+.999,ans]]​​​​
geoflat[{a_,b_},{x_,y_}]:=Module[{cross1={},cross2={},allcross,segs},If[x≠0,cross1=((#1-a&)/@Range[Ceiling[Min[a,a+x]],Floor[Max[a,a+x]]])/x];​​If[y≠0,cross2=((#1-b&)/@Range[Ceiling[Min[b,b+y]],Floor[Max[b,b+y]]])/y];​​allcross=Union[cross1,cross2];​​If[Length[allcross]>0,segs=Rest@Transpose[{RotateRight[allcross],allcross}];​​PrependTo[segs,{0,First[allcross]}];​​AppendTo[segs,{Last[allcross],1}],segs={{0,1}}];​​conpair/@Map[(Times[#1,{x,y}]+{a,b})&,segs,{2}]]​​

functions for donut=torus(3,1)

tor2[{x_,y_}]:={Cos[x](3+scaleTorusCos[y]),Sin[x](3+scaleTorusCos[y]),scaleTorusSin[y]}
eq1[{a1_,a2_},{b1_,b2_}]:={x''[t]2Sin[y[t]]/(3+Cos[y[t]])x'[t]y'[t],y''[t]-Sin[y[t]](3+Cos[y[t]])x'[t]^2,x[0]a1,y[0]a2,x'[0]b1,y'[0]b2}
perfun[solut_,per_][t_]:=Module[{axe,n,b},axe=First@First[{x[per],y[per]}/.solut];​​n=IntegerPart[t/per];​​b=t-nper;​​{naxe,0}+First@({x[b],y[b]}/.solut)]
ymax[a2_,b1_]:=Module[{ans,k=Abs[b1(3+Cos[a2])^2]},If[k>4,Return[]];​​If[k≤2,Pi,ArcCos[k-3]]]
periodover4[a2_,b1_]:=Module[{y1=0,y2=ymax[a2,b1],k=Abs[b1(3+Cos[a2])^2]},If[y2y1,Return[1]];​​If[y2<0||y2>Pi||k>4,Return[]];​​If[k≥2,NIntegrate[(3+Cos[y])/Sqrt[(3+Cos[y])^2-k^2],{y,y1,.9999y2}]+(3+Cos[y2])Sqrt[.0001y2]/Sqrt[2Sin[y2]],NIntegrate[(3+Cos[y])/Sqrt[(3+Cos[y])^2-k^2],{y,y1,y2}]]]
critgeo[{xp0_,yp0_},r_,ss_]:=Module[{sg1=Sign[xp0],sg2=Sign[yp0],sol3,tmaxi=9.5,ept=4.428},sol3=NDSolve[{x'[t]sg12/(3+Cos[y[t]])^2,y'[t]sg2Sqrt[5+6Cos[y[t]]+Cos[y[t]]^2],x[0]0,y[0]0},{x,y},{t,0,tmaxi}];​​If[r≤tmaxi,Return[Table[tor2@First[{x[t],y[t]}/.sol3],{t,0,r,ss}]]];​​If[r≤tmaxi+4Pi,Return[Join[Table[tor2@First[{x[t],y[t]}/.sol3],{t,0,tmaxi,ss}],Table[tor2@{sg1(ept+t/2),Pi},{t,0,r-tmaxi,ss}]]]];​​Join[Table[tor2@First[{x[t],y[t]}/.sol3],{t,0,tmaxi,ss}],Table[tor2@{sg1(ept+t/2),Pi},{t,0,4Pi+Mod[r-tmaxi,4Pi]}]]]
​​geotor[{xp1_,yp1_}]:=Module[{xp0,yp0,t0,stepsize,r=Sqrt[16xp1^2+yp1^2]},If[r0,Return[]];​​xp0=xp1/r;​​yp0=yp1/r;​​stepsize=Min[.1,r/100];​​If[xp1≠0&&Abs[Abs[yp1/xp1]-4Sqrt[3]]<.002,Return[critgeo[{xp0,yp0},r,stepsize]]];​​t0=periodover4[0,xp0];​​sol=NDSolve[eq1[{0,0},{xp0,yp0}],{x,y},{t,0,4t0}];​​Table[tor2@perfun[sol,4t0][t],{t,0,r,stepsize}]]

functions for RP2

functions for Hyperbolic Plane

Global functions

Surfaces as Graphs