(*******************************************************************
This file was generated automatically by the Mathematica front end.
It contains Initialization cells from a Notebook file, which
typically will have the same name as this file except ending in
".nb" instead of ".m".

This file is intended to be loaded into the Mathematica kernel using
the package loading commands Get or Needs.  Doing so is equivalent
to using the Evaluate Initialization Cells menu command in the front
end.

DO NOT EDIT THIS FILE.  This entire file is regenerated
automatically each time the parent Notebook file is saved in the
Mathematica front end.  Any changes you make to this file will be
overwritten.
***********************************************************************)

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`"]




scaleSphere=1.02
scaleTorus=1.02



expsphere[{x_,y_}]:=
  Module[{r=(x^2+y^2)^(1/2)},If[r\[Equal]0,Return[{0,0,1}]];
    {x/r Sin[r],y/r Sin[r],Cos[r]}]


geosphere[x_,y_,n_Integer]:=
  Module[{r=(x^2+y^2)^(1/2),k},
    If[r\[GreaterEqual]2Pi,Return[greatcircle[x/r,y/r,n]]];
    k=Ceiling[r n];
    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)}]



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_]:=(*This function is used to fix round-off errors.*)
    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\[NotEqual]0,
      cross1=((#1-a&)/@Range[Ceiling[Min[a,a+x]],Floor[Max[a,a+x]]])/x];
    If[y\[NotEqual]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}]]




tor2[{x_,y_}]:={Cos[x] (3+scaleTorus Cos[y]),Sin[x](3+scaleTorus Cos[y]),
    scaleTorus Sin[y]}

eq1[{a1_,a2_},{b1_,b2_}]:={x''[t]\[Equal]2Sin[y[t]]/(3+Cos[y[t]])x'[t]y'[t],
    y''[t]\[Equal]-Sin[y[t]](3+Cos[y[t]])x'[t]^2,x[0]\[Equal]a1,
    y[0]\[Equal]a2,x'[0]\[Equal]b1,y'[0]\[Equal]b2}

perfun[solut_,per_][t_]:=
  Module[{axe,n,b},axe=First@First[{x[per],y[per]}/.solut];
    n=IntegerPart[t/per];
    b=t-n per;
    {n axe,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\[LessEqual]2,Pi,ArcCos[k-3]]]

periodover4[a2_,b1_]:=
  Module[{y1=0,y2=ymax[a2,b1],k=Abs[b1(3+Cos[a2])^2]},
    If[y2\[Equal]y1,Return[1]];
    If[y2<0||y2>Pi||k>4,Return[]];
    If[k\[GreaterEqual]2,
      NIntegrate[(3+Cos[y])/Sqrt[(3+Cos[y])^2-k^2],{y,
            y1,.9999y2}]+(3+Cos[y2])Sqrt[.0001y2]/Sqrt[2 Sin[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]\[Equal]sg1 2/(3+Cos[y[t]])^2,
          y'[t]\[Equal]sg2 Sqrt[5+6Cos[y[t]]+Cos[y[t]]^2],x[0]\[Equal]0,
          y[0]\[Equal]0},{x,y},{t,0,tmaxi}];
    If[r\[LessEqual]tmaxi,
      Return[Table[tor2@First[{x[t],y[t]}/.sol3],{t,0,r,ss}]]];
    If[r\[LessEqual]tmaxi+4 Pi,
      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\[Equal]0,Return[]];
    xp0=xp1/r;
    yp0=yp1/r;
    stepsize=Min[.1,r/100];
    If[xp1\[NotEqual]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}]]



polcoord[{x_,y_}]:={Cos[2Pi x] Cos[Pi y/2],Sin[2 Pi x] Cos[Pi y/2],
    Sin[Pi y/2]}


halfsphere[n_Integer,m_Integer]:=
  Module[{grid2,top,triangs},grid2=Table[{i/n,j/m},{j,0,m-1},{i,0,n}];
    top=Rest@Last[grid2];
    triangs=Transpose[{top,RotateRight[top],Table[{0,1},{Length[top]}]}];
    squas=
      Flatten[Transpose[{Rest[Rest/@grid2],Rest[Rest/@RotateRight[grid2]],
            Rest[Rest/@RotateRight[grid2,{1,1}]],
            Rest[Rest/@RotateRight[grid2,{0,1}]]},{3,1,2,4}],1];
    Polygon/@(Map[polcoord,Join[triangs,squas],{2}])]


halfsphere2[n_Integer,m_Integer]:=
  Module[{grid2,top,triangs},grid2=Table[{i/n,j/m},{j,0,m-1},{i,0,n}];
    top=Rest@Last[grid2];
    squas=
      Flatten[Transpose[{Rest[Rest/@grid2],Rest[Rest/@RotateRight[grid2]],
            Rest[Rest/@RotateRight[grid2,{1,1}]],
            Rest[Rest/@RotateRight[grid2,{0,1}]]},{3,1,2,4}],1];
    Polygon/@(Map[polcoord,Join[{top},squas],{2}])]



exphalfsphere[{x_,y_}]:=
  Module[{r=(x^2+y^2)^(1/2),th},If[r\[Equal]0,Return[{0,0,1}]];
    th=Mod[r,Pi,-Pi/2];
    {x/r Sin[th],y/r Sin[th],Cos[th]}]

georp2[x_,y_,n_Integer]:=
  Module[{r=(x^2+y^2)^(1/2),k},
    If[r\[GreaterEqual]Pi,Return[{halfgreatcircle[x/r,y/r,n]}]];
    If[r\[GreaterEqual]Pi/2,k=Ceiling[(r-Pi/2)n];
      {expsphere/@Table[((r-Pi/2)t-Pi/2){x/r,y/r},{t,0,1,1/k}],
        quartergreatcircle[x/r,y/r,n]},k=Ceiling[r n];
      {expsphere/@Table[t{x,y},{t,0,1,1/k}]}]]

halfgreatcircle[x_,y_,n_Integer]:=
  expsphere/@Table[t{x,y},{t,-Pi/2,Pi/2,Pi/(3n)}]

quartergreatcircle[x_,y_,n_Integer]:=
  expsphere/@Table[t{x,y},{t,0,Pi/2,Pi/(4n)}]



cent[{x_,y_}]:=If[x\[NotEqual]0,{y/x,0}]

radius[{x_,y_}]:=Module[{r=(x^2+y^2)^(1/2)},If[x\[NotEqual]0,r/Abs[x]]]

ang[{x_,y_}]:=If[y\[NotEqual]0,Mod[ArcTan[-x/y],Pi]]

ang[{x_,y_},{o1_,o2_}]:=If[x\[NotEqual]o1,Mod[ArcTan[(y-o2)/(x-o1)],Pi]]

trans[th_,{x_,y_}]:=
  Module[{r=(x^2+y^2)^(1/2),s0,s1},s0=Log[Sqrt[1+Cos[th]]/Sqrt[1-Cos[th]]];
    If[x>0,s1=s0+r,s1=s0-r];
    {Tanh[s1],Sech[s1]}]

vertgeo[y_]:=
  Show[Graphics[{Green,Line[{{-1,0},{1,0}}],Black,Thickness[.01],
        Line[{{0,1},{0,Exp[y]}}],Blue,PointSize[.03],Point[{0,Exp[y]}],
        NavyBlue,Point[{0,1}]}],PlotRange\[Rule]{{-1,1},{0,Max[1,Exp[y]]+1}},
    Axes\[Rule]True,AspectRatio\[Rule]Automatic]

horzgeo[x_]:=
  Show[Graphics[{Green,
        Line[Sort[{{1-2UnitStep[x],0},{N[x]+2UnitStep[x]-1,0}}]],Black,
        Thickness[.01],Line[{{0,1},{x,1}}],Blue,PointSize[.03],Point[{x,1}],
        NavyBlue,Point[{0,1}]}],
    PlotRange\[Rule]{Sort[{1-2UnitStep[x],N[x]+2UnitStep[x]-1}],{0,2}},
    Axes\[Rule]True,AspectRatio\[Rule]Automatic]

RP2=halfsphere2[25,9]



AllGeo[{x_,y_},ind_]:=
  Switch[ind,-3,SphereGeodesic[x,y],-2,RP2Geodesic[x,y],-1,DonutGeodesic[x,y],
    0,FlatGeodesic[x,y],1,HyperPlaneGeodesic[x,y]]

DonutGeodesic[x_,y_]:=
  Module[{pts},pts=If[x\[Equal]0&&y\[Equal]0,{tor2@{0,0}},geotor[{x,y}]];
    Graphics3D[{Green,Torus[3,1],Black,Thickness[.01],Line[pts],
        PointSize[.03],NavyBlue,Point[First[pts]],Blue,Point[Last[pts]]},
      Lighting\[Rule]False,Boxed\[Rule]False,
      ViewPoint\[Rule]{2.605,0.015,1.660}]]


FlatGeodesic[x_,y_]:=
  Graphics[{Green,Thickness[.02],Line[{{0,0},{1,0},{1,1},{0,1},{0,0}}],
      Thickness[.01],Black,Line/@geoflat[{.5,.5},{N@x,N@y}],PointSize[.03],
      NavyBlue,Point[{0.5,0.5}],Blue,Point@expflat[{x,y}+{.5,.5}]},
    AspectRatio\[Rule]Automatic]



HyperPlaneGeodesic[x_,y_]:=
  Module[{th,th2,dest,r,deco,c,rax,rany},If[x\[Equal]0,Return[vertgeo[y]]];
    If[y\[Equal]0,Return[horzgeo[x]]];
    th=ang[{x,y}];
    dest=trans[th,{x,y}];
    th2=ang[dest,{0,0}];
    r=radius[{x,y}];
    c=cent[{x,y}];
    deco=r dest+c;
    If[deco[[1]]>0,rax={-1,deco[[1]]+1},rax={deco[[1]]-1,1}];
    If[IntervalMemberQ[Interval[{th,th2}],Pi/2],rany={0,r+1},
      rany={0,Max[1,deco[[2]]]+1}];
    Graphics[{Green,Line[Transpose[{rax,{0,0}}]],Black,Thickness[.01],
        Circle[c,r,Sort@N[{th,th2}]],Blue,PointSize[.03],Point[deco],NavyBlue,
        Point[{0,1}]},PlotRange\[Rule]{rax,rany},Axes\[Rule]True,
      AspectRatio\[Rule]Automatic]]


RP2Geodesic[x_,y_]:=
  Graphics3D[{Green,RP2,Black,Thickness[.015],
      Line/@(scaleSphere georp2[x,y,10]),NavyBlue,PointSize[.03],
      Point[scaleSphere exphalfsphere[{0,0}]],Blue,
      Point[scaleSphere exphalfsphere[{x,y}]]},Lighting\[Rule]False,
    PlotRange\[Rule]1.1{{-1,1},{-1,1},{-1,1}}]


SphereGeodesic[x_,y_]:=
  Graphics3D[{Green,Sphere[],Black,Thickness[.015],
      Line[scaleSphere geosphere[x,y,10]],NavyBlue,PointSize[.03],
      Point[scaleSphere expsphere[{0,0}]],Blue,
      Point[scaleSphere expsphere[{x,y}]]},Lighting\[Rule]False,
    PlotRange\[Rule]1.1{{-1,1},{-1,1},{-1,1}}]





chris[glo_,gup_,vars_List]:=
  Table[.5 Tr[
        Table[gup[[i,j]](D[glo[[k,l]],vars[[i]]] -D[glo[[i,k]],vars[[l]]]-
                D[glo[[i,l]],vars[[k]]]),{i,2}]],{j,2},{k,2},{l,2}]

equ[gam_List,{x_,y_}]:=Module[{cs=gam/.{x\[Rule]x[t],y\[Rule]y[t]}},
    {x''[t]\[Equal]First[{x'[t],y'[t]}.cs[[1]].{{x'[t]},{y'[t]}}],
      y''[t]\[Equal]First[{x'[t],y'[t]}.cs[[2]].{{x'[t]},{y'[t]}}]}]

equ[gam_List,{a1_,a2_},{b1_,b2_},{x_,y_}]:=
  Join[equ[gam,{x,y}],{x[0]\[Equal]a1,y[0]\[Equal]a2,x'[0]\[Equal]b1,
      y'[0]\[Equal]b2}]



geosurf[f_,{a1_,a2_},{x_,y_}][b1_,b2_]:=Module[{gup,glo,cs,sol1,r,r2},
    glo=metr[f,{x,y}];
    r2=First[{b1,b2}.(glo/.{x\[Rule]a1,y\[Rule]a2}).{{b1},{b2}}];
    If[r2\[LessEqual]0, Return[{{a1,a2}}]];
    r=Sqrt[r2];
    ss=Min[.05,r/100];
    gup=Inverse[glo];
    cs=chris[glo,gup,{x,y}];
    sol1=NDSolve[equ[cs,{a1,a2},{b1/r,b2/r},{x,y}],{x,y},{t,0,r}];
    Table[First[{x[t],y[t]}/.sol1],{t,0,r,ss}]]



metr[f_,{x_,y_}]:=(Outer[Dot,#1,#1,1]&)[{{1,0,D[f,x]},{0,1,D[f,y]}}]

grf[f_,unitnormal_,d_,{x_,y_}][{a_,b_}]:={a,b,(f/.{x\[Rule]a,y\[Rule]b})}+
    d (unitnormal/.{x\[Rule]a,y\[Rule]b})


geoAndsurf[f_,ra_,{x_,y_}][b1_,b2_]:=
  Module[{pts=geosurf[f,{0,0},{x,y}][b1,b2],gg,un},
    un={-D[f,x],-D[f,y],1};
    un=un/(un.un);
    gg=Plot3D[{f/.{x\[Rule]s,y\[Rule]t},RGBColor[1,0,0]},{s,-ra,ra},{t,-ra,
          ra},PlotPoints\[Rule]30,DisplayFunction\[Rule]Identity];
    {Graphics3D[{Thickness[.015],GrayLevel[0],Line[grf[f,un,.08,{x,y}]/@pts]},
        Boxed\[Rule]False,ViewPoint\[Rule]{0,-ra,ra/2},
        PlotRange\[Rule]ra{{-1,1},{-1,1},{-1,1}}],gg,
      DisplayFunction\[Rule]$DisplayFunction}]




End[]

EndPackage[]