Off[General::spell];  
    Off[General::spell1]; 
Off[Solve::ifun];
Off[Reduce::ifun];
Needs["Graphics`ImplicitPlot`"];
Needs["Graphics`ContourPlot3D`"];
disp := DisplayFunction -> $DisplayFunction

AxesOnSurface[
        eq_] := 
Module[{ϕaxis = 
            Insert[ParametricPlot3D[
                Evaluate[
                  Append[eq[0, ϕ], 
                    RGBColor[0, 1, 0]]], {ϕ, -Pi,  π}, 
                Lighting -> False, PlotPoints -> 40, 
                Boxed -> False, Axes -> False, 
                DisplayFunction -> Identity], 
              Thickness[0.015], {1, 1}], θaxis = 
            Insert[ParametricPlot3D[
                Evaluate[
                  Append[eq[θ, 0], 
                    RGBColor[1, 0, 0]]], {θ, -Pi,  π}, 
                Lighting -> False, PlotPoints -> 40, 
                Boxed -> False, Axes -> False, 
                DisplayFunction -> Identity], 
              Thickness[0.015], {1, 1}]}, 
        Show[{ϕaxis, θaxis}, 
          DisplayFunction -> Identity]]

tor[θ_, ϕ_] := {(Cos[θ] + 4) Cos[ϕ], 
        Sin[θ], (Cos[θ] + 4) Sin[ϕ]}

torus[a_ : -π, b_ : π, c_ : -π, d_ : π, 
          plotpoints_ : {40, 20}] := 
        ParametricPlot3D[
          Evaluate[tor[θ, ϕ]], {ϕ, c, d}, {θ, a, b}, 
          PlotPoints -> plotpoints, Boxed -> False, 
          Axes -> False, DisplayFunction -> Identity];

height[θ_, ϕ_] := tor[θ, ϕ][[3]]

T[α_, β_][θ_, ϕ_] := {α + θ, β + ϕ};

trans = 
        Apply[T, 
          Distribute[{{(π/2), -(π/2)}, {(π/2), -(π/2)}}, List], {1}];

height[i_][θ_, ϕ_] := 
        height @@ trans[[
              i]][θ, ϕ];

gradient[f_, {θ__}] := 
      Through[Through[Apply[Derivative, 
                IdentityMatrix[Length[{θ}]], {1}][
              f]][θ]]

solrules := 
        DeleteCases[(Solve[
                  gradient[height[#1], {θ, ϕ}] == 
                    0, {θ, ϕ}] &) /@ 
            Range[4], {ϕ -> x_, θ -> 
                y_} /; ! (Im[x] == 0 && 
                  Im[y] == 0), ∞];

TorusAndAxes := Show[torus[], AxesOnSurface[tor]]

solutions := 
        Sort[Flatten[(Apply[#1, {θ, ϕ} /. solrules[[1]], {1}] &) /@ trans, 
              1] ⋃ (SameTest -> (tor @@ #1 == 
                      tor @@ #2 &)), 
          height @@ #1 > height @@ #2 &];

Hessian[f_, {x__}] := Flatten[Outer[D, {f}, {x}, {x}], 1]

index[f_, x_List, a_List] := 
      Count[Eigenvalues[
          Hessian[f, x] /. Thread[x -> a]], _?
          Negative]

MorsePolynomial[indices_List, 
        t_] := (∑_(i = Min[indices])^(Max[indices]))Count[indices, 
            i] (t^i)

PiecewiseFlow3D[eq_, height_, l_List, t_List, 
        color_:  RGBColor[1, 0, 0]] := 
      Module[{sols = {}, sol, solrules, s, θ, ϕ}, 
        Fold[(AppendTo[
                sols, {Join @@ (solrules = 
                        NDSolve[
                          Join[Thread[{θ'[s], ϕ'[
                                    s]} == -gradient[
                                    height, {θ[s], ϕ[
                                        s]}]], {θ[
                                  0] == #1[[1]], ϕ[
                                  0] == #1[[2]]}], {θ, ϕ}, {s, 0, #2}]), {s, 0, #2}}]; 
              Flatten[({θ[s], ϕ[s]} /. solrules)] /. 
                s -> #2) &, l, t]; 
        With[{f = eq[θ[s], ϕ[s]][[1]], 
            g = eq[θ[s], ϕ[s]][[2]], 
            h = eq[θ[s], ϕ[
                    s]][[3]]}, ParametricPlot3D[##, 
                    DisplayFunction -> Identity] & @@ {{f, g, h, 
                      color} /. #[[1]], #[[2]]} & /@ 
            Evaluate[sols]]]

LevelPoints[height_, value_, {l1_, l2_, l3_, l4_}, 
        pts_] := Graphics[
          ImplicitPlot[
            Evaluate[height[a, b] == value], {a, l1, l2}, {b, l3, l4}, 
            DisplayFunction -> Identity, 
            PlotPoints -> pts]][[1, 1, 3, 1]]

LevelPoints[height_, 
        value_] := Graphics[
          ImplicitPlot[
            Evaluate[height[a, b] == value], {a, -Pi, 
              Pi}, {b, -Pi, Pi}, DisplayFunction -> Identity, 
            PlotPoints -> 90]][[1, 1, 3, 1]]

LevelCurve[eq_, height_, val_] := 
      Show[Graphics3D[{RGBColor[1, 0, 0], Thickness[0.02], 
                Line[#]} &@(eq @@@ LevelPoints[height, val])], 
        DisplayFunction -> Identity]

LevelCurvesOnTorus := 
      Module[{values = {4.8, 2.8, -2.9, -3.5}, p, q, 
          r}, 
p = 
          Flatten[Map[LevelPoints[height, #] &, values], 1]; 
        q = p /. {x_?NumberQ, y_?NumberQ} -> {x, Pi - y}; 
        r = Join[p, 
            q]; 
Show[{Graphics3D[
              Map[{PointSize[0.02], Point[#]} &, tor @@@ r]], torus[]}, 
          DisplayFunction -> Identity]]

SpecialFlowsOnTorus := 
      With[{startpoints = 
            Take[solutions, {1, 3}]}, PiecewiseFlow3D[tor, 
              height, #, {1, 1, 1, 1, 1}, RGBColor[0, 1, 0]] & /@ 
          N[Join[(Transpose[Take[startpoints, {2, 3}]] + {0, 1/10}) // 
                Transpose, (Transpose[
                      Take[startpoints, {2, 3}]] + {0, -1/10}) // 
                Transpose, (Transpose[startpoints] + {1/10, 0}) // 
                Transpose, (Transpose[startpoints] + {-1/10, 0}) // 
                Transpose]]]

FlowsOnTorus1 := 
      With[{values = {4.8, 
              2.8, -2.9, -3.5}}, 
Show[{PiecewiseFlow3D[tor, height, #, {1, 1, 1}] & /@ 
              Flatten[Map[
                  LevelPoints[height, #, {-Pi, 0, -Pi/2, Pi/2}, 
                      13] &, values], 1]}]]

FlowsOnTorus2 := 
      Module[{w = FlowsOnTorus1, u, v, p}, 
u = 
          w /. Line[{{a_, b_, c_}, {d_, e_, f_}}] -> 
              Line[{{a, -b, c}, {d, -e, f}}]; 
        v = w /. Line[{{a_, b_, c_}, {d_, e_, f_}}] -> 
              Line[{{-a, b, c}, {-d, e, f}}]; 
        p = u /. Line[{{a_, b_, c_}, {d_, e_, f_}}] -> 
              Line[{{-a, b, c}, {-d, e, f}}]; 
        Show[{w, u, v, p}, DisplayFunction -> Identity]]

FlowsOnTorus := 
      Show[{torus[], FlowsOnTorus2, 
          SpecialFlowsOnTorus /. 
            Line[x_] :> {Thickness[0.01], Line[x]}}]

th = Pi/6;

Rt = 
        FunctionExpand[{{1, 0, 0}, {0, Cos[th], Sin[th]}, {0, -Sin[th], 
              Cos[th]}}];

tiltedTor[θ_, ϕ_] := 
        Flatten[Rt . Transpose[{tor[θ, ϕ]}]];

tiltedTorus[a_: - Pi, b_:  Pi, c_: - Pi, d_:  Pi] := 
        ParametricPlot3D[
          Evaluate[tiltedTor[θ, ϕ]], {ϕ, c, d}, {θ, 
            a, b}, PlotPoints -> {40, 20}, Boxed -> False, 
          Axes -> False, DisplayFunction -> Identity];

tiltedHeight[θ_, ϕ_] := tiltedTor[θ, ϕ][[3]]

tiltedHeight[i_][θ_, ϕ_] := 
        tiltedHeight @@ trans[[
              i]][θ, ϕ];

tiltedSolrules := 
        DeleteCases[(Solve[
                  gradient[tiltedHeight[#1], {θ, ϕ}] == 
                    0, {θ, ϕ}] &) /@ 
            Range[4], {ϕ -> x_, θ -> 
                y_} /; ! (Im[x] == 0 && 
                  Im[y] == 0), ∞];

tiltedSolutions := 
        N[Sort[Union[
              Select[Flatten[(Apply[#1, {θ, ϕ} /. 
                            tiltedSolrules[[1]], {1}] &) /@ trans, 1], 
                gradient[tiltedHeight, #] == {0, 0} &], 
              SameTest -> (Simplify[tiltedTor @@ #1] == 
                      Simplify[tiltedTor @@ #2] &)], 
            tiltedHeight @@ #1 > tiltedHeight @@ #2 &]];

tiltedColors := 
        Join[Table[RGBColor[1, 0, 0], {1}], Table[RGBColor[1, 1, 0], {2}], 
          Table[RGBColor[0, 0, 1], {1}]];

tiltedCoords := 
        tiltedTor @@@ N[tiltedSolutions] // Chop;

tiltedColoredCoords := 
        Transpose[{tiltedColors, tiltedCoords}];

tiltedCriticals := 
        Graphics3D[{#1, PointSize[0.04], Point[N[#2]]} & @@@ 
            tiltedColoredCoords];

SpecialFlowsOnTiltedTorus := 
      With[{startpoints = {{-0.42359877559829884`, 
                1.5707963267948966`}, {-0.6235987755982988`, 
                1.5707963267948966`}, {2.4764012244017013`, 
                1.5807963267948966`}, {2.4764012244017013`, 
                1.5607963267948965`}, {-2.7179938779914945`, -1.5707963267948966`}, {-2.5179938779914943`, -1.5707963267948966`}, {-1.3764012244017012`, 
                1.5907963267948966`}, {-1.3764012244017012`, 
                1.5507963267948965`}}}, PiecewiseFlow3D[tiltedTor, 
              tiltedHeight, #, {1, 1, 1, 1, 1}, RGBColor[0, 1, 0]] & /@ 
          startpoints]

FlowsOnTiltedTorus1 := 
      With[{values = {4.3, 
              2.3, -2.5}}, 
Show[{PiecewiseFlow3D[
                  tiltedTor, tiltedHeight, #, {1, 1, 1}] & /@ 
              Flatten[Map[
                  LevelPoints[tiltedHeight, #, 
                      N[{-Pi, Pi, -Pi/2, Pi/2}], 40] &, values], 
                1]}]]

FlowsOnTiltedTorus2 := 
      With[{w = FlowsOnTiltedTorus1}, 
        Show[{w, 
            w /. Line[{{a_, b_, c_}, {d_, e_, f_}}] -> 
                Line[{{-a, b, c}, {-d, e, f}}]}, 
          PlotRange -> All, 
          DisplayFunction -> Identity]]

FlowsOnTiltedTorus := 
      Show[{tiltedTorus[], FlowsOnTiltedTorus2, 
          SpecialFlowsOnTiltedTorus /. 
            Line[x_] :> {Thickness[0.01], Line[x]}}]

vertices[s_] := 
        Union[Flatten[Cases[s, Polygon[l_] :> l, Infinity], 
            1]] ;
SetAttributes[Side, Orderless];
sides[s_] := 
        Union[Cases[
            s /. Polygon[l_] :> 
                Side @@@ Partition[Append[l, First[l]], 2, 1], _Side, 
            Infinity]];
faces[s_] := Cases[s, _Polygon, Infinity];
NumberOfVertices[s_] := Length[vertices[s]]; 
    NumberOfSides[s_] := Length[sides[s]];
NumberOfFaces[s_] := Length[faces[s]];
EulerCharacteristic[s_] := 
        NumberOfVertices[s] - NumberOfSides[s] + 
          NumberOfFaces[s];
star[v_, gr_] := 
        Cases[gr, Polygon[{l___, v, m___}], Infinity];
link[v_, gr_] := 
        Union[Cases[
            star[v, gr] /. 
              Polygon[l_] :> 
                Line /@ Partition[Append[l, First[l]], 2, 
                    1], _?((Head[#] === Line) && (FreeQ[#, 
                        v]) &), Infinity]];
star[v_Side, tr_] := 
        Cases[tr, 
          Polygon[l_] /; 
            MemberQ[Side @@@ Partition[Append[l, First[l]], 2, 1], v], 
          Infinity];
boundary[tr_] := 
      Select[sides[tr], Length[star[#, tr]] == 1 &] /. 
        Side[l__] :> Line[{l}]
OrderedUnion[l_List] := 
        Block[{i}, i[n_] := (i[n] = Sequence[]; n);  
          i /@ l];
closeup[k_, epsilon_:  10^-4] := 
        k /. Polygon[l_] :> 
            Polygon[N[OrderedUnion[Rationalize[l, epsilon]]]];

sphereloop[t_] := 
        ParametricPlot3D[Evaluate[sp[u, t]], {u, 0, 2  Pi}, 
          Boxed -> False, Axes -> False, 
          DisplayFunction -> Identity];

torusloop[s_] := 
      Insert[ParametricPlot3D[
          Evaluate[
            Append[tor[θ, s], 
              RGBColor[1, 0, 0]]], {θ, -Pi,  π}, 
          Lighting -> False, PlotPoints -> 40, Boxed -> False, 
          Axes -> False, DisplayFunction -> Identity], 
        Thickness[0.02], {1, 1}]

TranslateShape[shape_, vec_List] := 
      Block[{tvec = N[vec]}, 
          shape /. {poly : Polygon[_] :> 
                Map[(tvec + #) &, poly, {2}], 
              line : Line[_] :> Map[(tvec + #) &, line, {2}], 
              point : Point[_] :> 
                Map[(tvec + #) &, point, {1}]}] /; 
        Length[vec] == 3

plane[v_] := 
      TranslateShape[
        ParametricPlot3D[{2 (1 - v)*Cos[t], 2  v*Sin[t], 
            s}, {t, -Pi/2, Pi/2}, {s, -1, 1}, Boxed -> False, 
          Axes -> False, 
          DisplayFunction -> Identity], {-1 + 2 v, 0, 
          0}]

sp[u_, v_] := {Cos[u] Cos[v], Sin[u] Cos[v], Sin[v]}

sphere[p_:  25] := 
        ParametricPlot3D[
          Evaluate[sp[u, v]], {u, 0, 2  Pi}, {v, -Pi/2, Pi/2},  
          PlotPoints -> p, Boxed -> False, Axes -> False, 
          DisplayFunction -> Identity];

cyl[pts_:  Automatic] := 
      With[{c1 = 
            ParametricPlot3D[{r Cos[θ], r Sin[θ], 0}, {r, 0, 
                1}, {θ, 0, 2  Pi}, PlotPoints -> pts, 
              DisplayFunction -> Identity], 
          c2 = ParametricPlot3D[{r Cos[θ], r Sin[θ], 1}, {r, 
                0, 1}, {θ, 0, 2  Pi}, PlotPoints -> pts, 
              DisplayFunction -> Identity], 
          c3 = ParametricPlot3D[{Cos[θ], Sin[θ], t}, {θ, 
                0, 2  Pi}, {t, 0, 1}, PlotPoints -> pts, 
              DisplayFunction -> 
                Identity]}, 
Show[{c1, c2, c3}, 
          Boxed -> False, Axes -> False]]

klein[θ_, ϕ_] := {(2 + Cos[ϕ/2] Sin[θ] - 
                Sin[ϕ/2] Sin[2  θ]) Cos[ϕ], (2 + 
                Cos[ϕ/2] Sin[θ] - 
                Sin[ϕ/2] Sin[2  θ]) Sin[ϕ], 
          Sin[ϕ/2] Sin[θ] + 
            Cos[ϕ/2] Sin[2  θ]};

kleinBottle[a_: - Pi, b_:  Pi, c_: - Pi, d_:  Pi, plotpoints_:  25] := 
        ParametricPlot3D[
          Evaluate[klein[θ, ϕ]], {ϕ, c, d}, {θ, a, 
            b}, PlotPoints -> plotpoints, Boxed -> False, 
          Axes -> None, Lighting -> True, 
          DisplayFunction -> Identity];

rectangle := 
        Show[Graphics[{RGBColor[1, 0, 0], 
              Line[{{-π, -π}, {π, -π}}], 
              RGBColor[0, 1, 0], Line[{{π, -π}, {π, π}}], 
              RGBColor[1, 0, 0], Line[{{π, π}, {-π, π}}], 
              RGBColor[0, 1, 0], 
              Line[{{-π, π}, {-π, -π}}]}], 
          Axes -> True, AspectRatio -> Automatic, 
          AxesLabel -> {θ, ϕ}, 
          DisplayFunction -> Identity];

smallrectangle := 
        Graphics[{RGBColor[0, 0, 1], 
            Line[{{-π/2, -π/2}, {π/2, -π/2}}], 
            Line[{{π/2, -π/2}, {π/2, π/2}}], 
            Line[{{π/2, π/2}, {-π/2, π/2}}], 
            Line[{{-π/2, π/2}, {-π/2, -π/
                    2}}]}];

\(patches := 
      With[{t1 = 
            Insert[torus[\(-Pi\) - 0.1, 0.1, \(-Pi\) - 0.1, 0.1], 
              RGBColor[1, 0, 0], {1, 1}], 

          t2 = Insert[torus[\(-Pi\) - 0.1, 0.1, \(-0.1\), Pi + 0.1], 
              RGBColor[0, 0, 1], {1, 1}], 

          t3 = Insert[torus[\(-0.1\), Pi + 0.1, \(-Pi\) - 0.1, 0.1], 
              RGBColor[0, 1, 0], {1, 1}], 

          t4 = Insert[torus[\(-0.1\), Pi + 0.1, \(-0.1\), Pi + 0.1], 
              RGBColor[1, 1, 0], {1, 1}]}, 

        Show[{t1, t2, t3, t4}, Lighting -> False]]\)