(* (c) 1990, Andrew J. Hanson, Indiana University.
   Modified by Stewart Dickson, The Post Group.
   Solution of Fermat's equation z1^n + z2^n = 1.   *)

MakePolygons[vl_List] :=
    Block[{l = vl,
           l1 = Map[RotateLeft, vl],
           mesh},
        mesh = {l, RotateLeft[l], RotateLeft[l1], l1};
        mesh = Map[Drop[#, -1]&, mesh, {1}];
        mesh = Map[Drop[#, -1]&, mesh, {2}];
        Polygon /@ Transpose[ Map[Flatten[#, 1]&, mesh] ]
    ]

u1[a_,b_] := .5 (E^(a + I*b) + E^(- a - I*b))
u2[a_,b_] := .5 (E^(a + I*b) - E^(- a - I*b))

z1k[a_,b_,n_,k_] := E^(k*2*Pi*I/n)*u1[a,b]^(2.0/n)
z2k[a_,b_,n_,k_] := E^(k*2*Pi*I/n)*u2[a,b]^(2.0/n)

(* Then make n^2 meshes with something like: *)

n=5; alpha = 0.25 Pi;
Do[ Do[ surface[k1+1,k2+1] =
      MakePolygons[Table[
        N[{Re[z1k[a,b,n,k1]], Re[z2k[a,b,n,k2]],
             Cos[alpha]*Im[z1k[a,b,n,k1]] + Sin[alpha]*Im[z2k[a,b,n,k2]]}],
         {a,-1.0,1.0,0.2}, {b,0, .5 Pi, .05 Pi}]],
    {k2,0, n - 1}],
  {k1,0, n - 1}]

(* where MakePolygons can be concocted from the source code for
   parametric plot 3D, or just use ParametricPlot3D directly.  
   But, as I mentioned, I have never succeeded in plotting the
   whole surface for n>2 because I run out of memory.  I can
   only plot a couple of meshes at a time, not all n^2 of them.

    -Andy *)

(* Plotting notes:
    n = even:     {a,-1.0,1.0,0.4}, {b,0, .5 Pi, .1 Pi}],
    n = odd:      {a,-1.0,1.0,0.2}, {b,0, .5 Pi, .05 Pi}],
*)

(* To Show:
   Show[Table[surface[k1, k2], {k1, 1, n}, {k2, 1, n}]]

   n = 3, use ViewPoint->{1.1771, 2.2947, -2.1905}
   n = 4, use ViewPoint->{1.0812, -0.9385, -3.0659} - SPD *)

VecLen[{x_, y_, z_}] := Sqrt[N[x ^2 + y ^2 + z ^2]];
Cross[u_List, v_List] := Det[{{{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}, u, v}];

centerview = Which[n == 2, {1.3, -2.4, 2.0},
	n == 3, {1.1771, 2.2947, -2.1905},
	n == 4, {1.0812, -0.9385, -3.0659},
	n == 5, {1.3, -2.4, 2.0},
	n == 6, {1.3, -2.4, 2.0} ];

centerline = {0.0, 0.0, -1.0};
eyeline = Cross[centerline, centerview];
ueyeline = eyeline / VecLen[eyeline];
viewdist = VecLen[centerview];
eyedist = viewdist / 24;
leftview = centerview + eyedist * ueyeline;
rightview = centerview - eyedist * ueyeline + {0.0, 0.0, 0.0001};
leftplot = Graphics3D[Table[surface[k1, k2], {k1, 1, n}, {k2, 1, n}],
	Background->RGBColor[0, 0, 0],
	AspectRatio->.47,
        ViewPoint->leftview, Boxed->False, RenderAll->False];
filename = "Fermat_left.ps";
Display[filename, leftplot];
Clear[leftplot];

rightplot = Graphics3D[Table[surface[k1, k2], {k1, 1, n}, {k2, 1, n}],
	Background->RGBColor[0, 0, 0],
	AspectRatio->.47,
        ViewPoint->rightview, Boxed->False, RenderAll->False];
filename = "Fermat_right.ps";
Display[filename, rightplot];
