The Chase

The Klein quartic is the zero locus of the homogeneous polynomial $x^3 y + y^3 z + z^3 x$ inside $\mathbb{CP}^2$. With real eyes, it is a genus-3 surface with the largest possible automorphism group. Echoing Henry Segerman, I would love to find a geometric immersion in $\mathbb{R}^3$ by applying some nice map $\mathbb{CP}^2 \to \mathbb{R}^3$, but in the absence of that, I took the real part, thought of as living in the real projective plane $\mathbb{RP}^2$.

Of course, I'm actually taking the curve as sitting in the sphere $S^2$, which double-covers $\mathbb{RP}^2$, so you see two curves (green and pink) which are antipodal images of each other. Then I let the sphere rotate around the $y$-axis and take the stereographic projection to the plane.

Here's the code, which required a slightly fancy bit of business with RegionFunction to make the two curves different colors:

Stereo[{x_, y_, z_}] := {x/(1 - z), y/(1 - z)};
InverseStereo[{x_, y_}] := {2 x/(1 + x^2 + y^2), 
   2 y/(1 + x^2 + y^2), (x^2 + y^2 - 1)/(1 + x^2 + y^2)};

Module[{cols},
 cols = RGBColor /@ {"#FF358B", "#AEEE00", "#333333"};
 ParallelTable[
  Show[ContourPlot[
      Evaluate[(Quiet[
           x[[1]]^3 x[[2]] + x[[2]]^3 x[[3]] + x[[3]]^3 x[[1]]] 
               /. {x -> (RotationMatrix[t, {0, 1, 0}].InverseStereo[{u, v}])}) == 0], 
        {u, -#, #}, {v, -#, #}, PlotPoints -> 20, PlotRange -> #, 
      ContourStyle -> Directive[cols[[1]], Thickness[.015]], 
      Frame -> False, ImageSize -> 500, 
      RegionFunction -> 
       Function[{u, v}, (-2 #[[1]] - 2 #[[2]] + 1 - #[[1]]^2 - #[[2]]^2 
           &[Stereo[RotationMatrix[t, {0, 1, 0}].InverseStereo[{u, v}]]]) > 0], 
      Background -> cols[[3]]] &[3], 
   ContourPlot[
      Evaluate[(Quiet[
           x[[1]]^3 x[[2]] + x[[2]]^3 x[[3]] + x[[3]]^3 x[[1]]] 
               /. {x -> (RotationMatrix[t, {0, 1, 0}].InverseStereo[{u, v}])}) == 0], 
        {u, -#, #}, {v, -#, #}, PlotPoints -> 20, PlotRange -> #, 
      ContourStyle -> Directive[cols[[2]], Thickness[.015]], 
      RegionFunction -> 
       Function[{u, v}, (-2 #[[1]] - 2 #[[2]] + 1 - #[[1]]^2 - #[[2]]^2 
           &[Stereo[RotationMatrix[t, {0, 1, 0}].InverseStereo[{u, v}]]]) < 0]] &[3]], 
   {t, 0, 2 π}]
]