Chapter Seven:

Advanced Three-Dimensional Graphics

The package ParametricPlot3D.m must be loaded for most of the graphics work in this chapter. However, it is built-in to version 2.0.

7.1 Optimizing a Root Cellar

c = 0.211728;
f[x_, t_] := Cos[2 Pi t - c x] Exp[-c x]
Plot3D[f[x, t], {x, 0, 30}, {t, 0, 2}];

f[x, t]

D[f[x, t], t]

The following use of Solve fails in 2.0. I think it's a bug.

Solve[N[%] == 0, t]

version 2.0 Notes

ParametricPlot3D is built into version 2. And it should be possible to combine plots and space curves and lines as in this example without resorting to ParametricPlot3D for the surface, but simply using Plot3D. And SpaceCurve is obsolete; space curves are generated by ParametricPlot3D.

Next graphic takes a while

<<:Graphics:ParametricPlot3D.m;

Show[

ParametricPlot3D[{x, t, f[x, t]},
	{x, 0, 30, 0.75}, {t, 0., 40/12, 1/12}],

SpaceCurve[{15, t, f[15, t] + .01}, {t, 0, 40/12, 1./12}],
		
Release[Table[SpaceCurve[
	{x, 0.0336 x + n, f[x, 0.0336 x + n]}, {x, 0, 30, 1}],
		{n, 0, 1}]],

Graphics3D[{Thickness[.0005],
		Line[{{30, .5, -1}, {30, .5, 0}}],
		Line[{{15, 40/12, -1}, {15, 40/12, f[15, 40/12]}}]}],

PlotRange->{{0, 30}, {0, 40/12}, {-1, 1}},
BoxRatios->{1, 1, .8}, Boxed->False,
Axes->Automatic, AxesEdge->{{1, -1}, Automatic, Automatic},
Ticks->{{0, 15, 30}, {0, .5 , 1, 2, 3}, {-.5, 0, .5, 1}},
AxesLabel->{"x: Depth, ft.", "t: Years", None},
ViewPoint->{2, 2, 2}]

Color root cellar from Color Plate 2

The following cell contains the code used to generate the colored root cellar surface on Color Plate 2. It uses ModifiedParametricPlot3D as explained in the torus cell below.

exp = .4;
Show[
ModifiedParametricPlot3D[{x, t, tt = f[x, t],  (u=tt+1//N;

Which[u>1,RGBColor[1,(u-1)^exp,0],u>=0,
		RGBColor[u^exp,0,(1-u)^exp]])},
				{x, 0, 30, (0.75)}, {t, 0, 40/12, (1/12)},
				DisplayFunction->Identity],

Graphics3D[GrayLevel[0]],
SpaceCurve[{15, t, f[15, t] + .001}, {t, 0, 40/12, 1./12},
				DisplayFunction->Identity],
		
Release[Table[SpaceCurve[
	{x, 0.0336 x + n, f[x, 0.0336 x + n]}, {x, 0, 30, 1},
				DisplayFunction->Identity],
		{n, 0, 1}]],

Graphics3D[{Thickness[.0005], 
		Line[{{30, .5, -1}, {30, .5, 0}}],
		Line[{{15, 40/12, -1}, {15, 40/12, f[15, 40/12]}}]}],



PlotRange -> {{0, 30}, {0, 40/12}, {-1, 1}},
BoxRatios -> {1, 1, .8},
DisplayFunction->$DisplayFunction,
Lighting -> False, Boxed -> False,
Axes -> Automatic, AxesEdge -> {{1, -1}, Automatic, Automatic},
Ticks -> {{0, 15, 30}, {0, .5 , 1, 2, 3}, {-.5, 0, .5, 1}},
AxesLabel -> {"Depth, feet", "Years", None},
ViewPoint -> {2, 2, 2}];

7.2 A Parametric Trochoid

ParametricPlot3D[{a, t - a Sin[t], 1 - a Cos[t]},
	{a, 0, 2, .1}, {t, 0., 5 Pi, Pi/16},
	BoxRatios->{4, 5 Pi, 3}, Boxed->False,
	ViewPoint->{1.8, 0, 1.7}]

7.3 The Polar View

Clear[f]
f[x_, y_] := x y / (x^2 + y^2)

Plot3D[f[x, y], {x, -2, 2}, {y, -2, 2},
 BoxRatios->{1.5, 1, 1.5}, ViewPoint->{-1.5, -2.8, 2.8},
 AxesEdge->{{-1, -1}, {-1, -1}, Automatic},
 AxesLabel->{"x", "y", None}]

f[r Cos[theta], r Sin[theta]]

Simplify[% /. Cos[x_]^2 -> 1 - Sin[x]^2]

ParametricPlot3D[{r Cos[theta], r Sin[theta], Sin[2 theta]/2},
					{r, 0, 1}, {theta, 0, 2 Pi}, 
	BoxRatios->{1.5, 1, 1.5}, ViewPoint->{-1.5, -2.8, 2.8},
	AxesLabel->{"x", "y", None}, Axes->Automatic,
	AxesEdge->{{-1, -1}, {-1, -1}, Automatic}];

ParametricPlot3D[{r Cos[theta], r Sin[theta], Sin[2 theta]/2},
					{r, 0, 1, .1}, {theta, 0, 2 Pi, Pi/30}, 
	Boxed->False, BoxRatios->{1.5, 1, 1.5},
	ViewPoint->{-1.5, -2.8, 2.8}];

Clear[g]
g[x_,y_] := x^2 (1+y)^3 + 7 y^2

ParametricPlot3D[{x = 4 r Cos[theta], y = r Sin[theta], 
					Max[-20, Min[20, g[x, y]]]},
	{r, 0, 3, .07}, {theta, 0, 2 Pi, Pi/25},
PlotRange-> {{-12, 12}, {-3, 3}, {-20, 20}},
Axes->Automatic, BoxRatios->{1, 1, .6},
Boxed->False, ViewPoint->{-1.1, -.4, .7},
AxesEdge->{Automatic, Automatic, {-1, +1}},
AxesLabel->{"x", "y", "z"},
AmbientLight->GrayLevel[.1]]

7.4 Rotating Circles to Produce a Torus or MÜbius Strip

Torus

rotate[t_] :=  {{Cos[t], -Sin[t], 0},
                {Sin[t],  Cos[t], 0},
                {  0,       0,    1}};
circle[t_] := {4, 0, 0} + {Cos[t], 0, Sin[t]};
rotate[theta] . circle[phi]

ParametricPlot3D[
		{(4 + Cos[phi]) Cos[theta],
		 (4 + Cos[phi]) Sin[theta],
		 Sin[phi]},
{phi, 0, 2 Pi, Pi/12}, {theta, 0, 2 Pi, Pi/12}, Boxed->False]

The version of ParametricPlot3D included with version 1.2 does not support user-defined coloring;  this is fixed in version 2.0. A modification that does allow for such coloring can be found in the firsition of the book by Roman Maeder [Mae, p. 223], and its use as follows will yield the shaded torus in Figure 7.7. For the following to work the version-1.2-user will have to copy, rename, and execute the appropriate code from Maeder's book. Alternatively, the user can use 1.2's ParametricPlot3D but delete the coloring instructions from any uses of ModifiedParametricPlot3D in this chapter.

ModifiedParametricPlot3D[
	{(4 + Cos[phi])*Cos[theta],
	(4+Cos[phi])*Sin[theta],
	Sin[phi],
	GrayLevel[Abs@Sin[theta/2]]},
{phi, 0, 2 Pi, Pi/12}, {theta, 0, 2 Pi, Pi/12},
Boxed->False, Lighting->False]

MÜbius Strip

circle[t_] := {4, 0, 0} + s {Cos[t], 0, Sin[t]}
ParametricPlot3D[rotate[theta] . circle[theta/2],
	{s, -1, 1, 2}, {theta, 0, 2 Pi, Pi/30},
	Boxed->False, ViewPoint->{2.7, -1.7, 2.3}]

Cylindrical Coordinates Diagram

ParametricPlot3D[{Cos[theta], Sin[theta], z},
					{theta, 0, 2 Pi, Pi/16}, {z, 0, 2, 2}]

Show[Graphics3D[{GrayLevel[0],
Polygon[{{1.1,1.1,.5}, {1.1,-1.1,.5}, {-1.1,-1.1,.5}, 
										{-1.1, 1.1,.5}}],
GrayLevel[.5],
Polygon[{{1.1,1.1,0}, {-1.1,-1.1,0}, {-1.1,-1.1,2},											{1.1,1.1,2}}],
GrayLevel[0], Line[{{0,0,0}, {0,0,2.5}}],
Line[{{0,-1.2,0}, {0,1.2,0}}], Line[{{1.2,0,0}, {-1.2,0,0}}]}]]

Show[%, %%, Boxed->False, BoxRatios->{1, 1, 1},
	ViewPoint->{2.4, 1.3, 1.3},
	Axes->Automatic, AxesLabel->{"x", "y", "z"},
	AmbientLight->GrayLevel[.15]]

7.5 Map Coloring on a Torus

T = {{a, b}, {c, d}};
Solve[T . {11, 11} == {12, 12} && 
      T . {12, 0}  == T . {0, 1} + {12, 0} &&
      T . {0, 12}  == T . {1, 0} + {0, 12}, {a, b, c, d}]

Read comment about ModifiedParametricPlot in torus section above!!

spiral[start_, end_, color_] := (*Modified*)ParametricPlot3D[
Append[
 rotate[(144theta + 12phi)/143].circle[(12theta + 144phi)/143],
 	GrayLevel[color]],
{theta, start, end, Pi/24}, {phi, -2 Pi/6, -Pi/6, Pi/12},
	DisplayFunction->Identity]

white1 = spiral[10 Pi +   Pi/12,  21 Pi,            1];
white2 = spiral[ 8 Pi +   Pi/12,  10 Pi - 3 Pi/12,  1];
white3 = spiral[ 2 Pi + 3 Pi/12,   4 Pi -   Pi/12,  1];
white4 = spiral[ 2 Pi - 4 Pi/12,   2 Pi -   Pi/12,  1];
white5 = spiral[      - 3 Pi/12,        -   Pi/12,  1];
white6 = spiral[        3 Pi/12,          7 Pi/12,  1];
green1 = spiral[10 Pi - 3 Pi/12,  10 Pi -   Pi/12,  .5];
green2 = spiral[ 4 Pi -   Pi/12,   4 Pi +   Pi/12,  .5];
green3 = spiral[ 2 Pi -   Pi/12,   2 Pi +   Pi/12,  .5];
green4 = spiral[      -   Pi/12,            Pi/12,  .5];
green5 = spiral[12 Pi - 3 Pi/12,  12 Pi -   Pi/12,  .5];
brown1 = spiral[10 Pi -   Pi/12,  10 Pi +   Pi/12,   0];
brown2 = spiral[ 8 Pi -   Pi/12,   8 Pi +   Pi/12,   0];
brown3 = spiral[ 2 Pi +   Pi/12,   2 Pi + 3 Pi/12,   0];
brown4 = spiral[          Pi/12,          3 Pi/12,   0];
yellow = spiral[ 6 Pi -   Pi/6,    6 Pi,            .3];
red    = spiral[ 4 Pi +   Pi/12,   6 Pi -   Pi/6,   .7];
blue   = spiral[ 6 Pi +   Pi/6,    8 Pi -   Pi/12,  .1];
purple = spiral[ 6 Pi,             6 Pi +   Pi/6,   .9];

Show[white1, white2, white3, white4, white5, white6, brown1,
 	 brown2, brown3, brown4, green1, green2, green3, green4,
	 green5, yellow, red, purple, blue,
	 Lighting->False, ViewPoint->{2, 0.2, 1.2},
	 Boxed->False, RenderAll->False,
	 PlotRange->{{-5, 5}, {-5, 5}, {-1, 1}},
	 DisplayFunction->$DisplayFunction]

Color Version of Torus with diagram as on Color Plate 3

Code

<<:Graphics:ParametricPlot3DColor.m  
(* Modified From Maeder's book; 
   see note in earlier torus discussion *)

rotate[t_] :=  {{Cos[t], -Sin[t], 0},
				{Sin[t],  Cos[t], 0},
				{  0,       0,    1}};

circle[t_] := {4, 0, 0} + {Cos[t], 0, Sin[t]};

spiral[start_, end_, color_] := ParametricPlot3D[Append[
 rotate[(144 theta+ 12 phi)/143].circle[(12 theta+ 144 phi)/143],
 	Apply[RGBColor, color]],
{theta, start, end, Pi/24}, {phi, -2 Pi/6, -Pi/6, Pi/12},
DisplayFunction -> Identity]

white1 =	spiral[10 Pi +   Pi/12, 21 Pi, 			  {.7,1,1}];
white2 = 	spiral[ 8 Pi +   Pi/12, 10 Pi - 3 Pi/12,  {.7,1,1}];
white3 = 	spiral[ 2 Pi + 3 Pi/12,  4 Pi -   Pi/12,  {.7,1,1}];
white4 = 	spiral[ 2 Pi - 4 Pi/12,  2 Pi -   Pi/12,  {.7,1,1}];
white5 = 	spiral[      - 3 Pi/12,       -   Pi/12,  {.7,1,1}];
white6 = 	spiral[      3 Pi/12,       7 Pi/12,  {.7,1,1}];
green1 =	spiral[10 Pi - 3 Pi/12, 10 Pi -   Pi/12, {.25,1,.25}];
green2 =	spiral[ 4 Pi -   Pi/12,  4 Pi +   Pi/12, {.25,1,.25}];
green3 =	spiral[ 2 Pi -   Pi/12,  2 Pi +   Pi/12, {.25,1,.25}];
green4 =	spiral[      -   Pi/12,           Pi/12, {.25,1,.25}];
green5 =	spiral[12 Pi - 3 Pi/12, 12 Pi -   Pi/12, {.25,1,.25}];
brown1 =	spiral[10 Pi -   Pi/12, 10 Pi +   Pi/12,  {1,.4,0}];
brown2 =	spiral[ 8 Pi -   Pi/12,	 8 Pi +   Pi/12,  {1,.4,0}];
brown3 =	spiral[ 2 Pi +   Pi/12,  2 Pi + 3 Pi/12, {1,.4,0}];
brown4 =	spiral[          Pi/12,         3 Pi/12,  {1,.4,0}];
yellow =	spiral[ 6 Pi -   Pi/6, 	 6 Pi, 			 {1,1,0}];
red    =	spiral[ 4 Pi +   Pi/12,  6 Pi -   Pi/6,  {1,0,.1}];
blue   =	spiral[ 6 Pi +   Pi/6, 	 8 Pi -   Pi/12, {.2,.2,1}];
purple =	spiral[ 6 Pi,		     6 Pi +   Pi/6,  {1,0,1}];

Show[white1, white2, white3, white4, white5, white6, brown1,
	 brown2, brown3, brown4, green1, green2, green3, green4,
	 green5, yellow, red, purple, blue,
	 Lighting -> False, ViewPoint -> {2, 0.2, 1.2},
	 Boxed -> False, RenderAll -> False,
	 PlotRange -> {{-5, 5}, {-5, 5}, {-1, 1}},
	 DisplayFunction -> $DisplayFunction,
	 Background->GrayLevel[0]];

Diagram of Torus With Seven Colored Countries

Bonus re. ∕7.5. Not in Book!

The following is a different approach to generating the colored torus. It computes all patches, not only those that are visible from the front, as is done in the book. The advantage is that the torus can then be rotated via SpinShow. Moreover, the code is shorter and much more elegant, as each country is defined in one line rather than in pieces.

The code has been sped up by replacing the call to rotate and circle with the explicit functions. 

This routine is for a color monitor. Best to set monitor to 16 colors when generating and showing the movie. 

The call to ParametricPlot3DColor.m is a call to a package I made
by modifying ParametricPlot3D.m so that it can color, as 
explained in the discussion of ModifiedParametricPlot above. None of this is necessary in version 2.0. Note that in version 2.0 the mesh size of
the grid is given by the PlotPoints option, not by a 4th element of the
iterator list, as is done throughout this chapter.

<<:Graphics:ParametricPlot3DColor.m;

spiral[theta0_, theta1_, phi0_, phi1_, color_] := ParametricPlot3D[

{(4 + Cos[(144*phi + 12*theta)/143])*Cos[(12*phi + 144*theta)/143], 
 (4 + Cos[(144*phi + 12*theta)/143])*Sin[(12*phi + 144*theta)/143], 
 Sin[(144*phi + 12*theta)/143], 
  
  RGBColor @@color},

{theta, theta0, theta1, Pi/12}, {phi, phi0, phi1, Pi/6},
			DisplayFunction->Identity, RenderAll->False];

purple =	spiral[0,          Pi/6,    Pi/6,     Pi/3,     {1,0,1}];
green  =	spiral[-3 Pi/12, -Pi/12,    Pi/2,  13 Pi/6, {.25,1,.25}];
brown  =	spiral[-Pi/12,    Pi/12,    Pi/3,     2 Pi,    {1,.4,0}];
yellow = 	spiral[-Pi/6,         0,    Pi/6,     Pi/3,     {1,1,0}];
blue   = 	spiral[Pi/6,   23 Pi/12,    Pi/6,     Pi/3,   {.2,.2,1}];
red    =	spiral[Pi/12,   11 Pi/6,       0,     Pi/6,    {1,0,.1}];
white  = 	spiral[Pi/12,  21 Pi/12,    Pi/6,     2 Pi,    {.7,1,1}];

theGraphic = 
Show[green, brown, yellow, blue, purple, red, white,
	Lighting -> False, ViewPoint -> {2, 0.2, 1.2},
	 Boxed -> False, 
	 PlotRange -> {{-5, 5}, {-5, 5}, {-1, 1}},
	 DisplayFunction -> $DisplayFunction,
	 Background->GrayLevel[0.3]];

Definition of SpinShow (from Samples 3 folder distributed with Mathematica)

Needs["Geometry`Rotations`"]

Attributes[SpinShow] = {};

Options[SpinShow] = {
	SpinOrigin -> {0,0,1.5},
	SpinTilt -> {0,0},
	SpinDistance -> 2,
	SpinSteps -> 15,
	SpinRange -> {0, 360 Degree},
	RotateLights -> False,
	Axes -> None};

SpinShow[theGraphic_, options___] := Block[
	{localOptions, spinOrigin, spinTilt, spinDistance,
	spinSteps, spinRange, lightSources, rotateLights,
	thetaFactor, phiFactor, rhoFactor,
	thetaIndex, thetaStep, thetaStart},
	
	localOptions = Join[{options}, Options[SpinShow]];
	spinOrigin = SpinOrigin /. localOptions;
	spinTilt = SpinTilt /. localOptions;
	spinDistance = SpinDistance /. localOptions;
	spinSteps = SpinSteps /. localOptions;
	spinRange = SpinRange //. localOptions;
	lightSources = LightSources /. localOptions /. 
					If[Length[theGraphic] == 2, 
						theGraphic[[2]], {}, {}] /.
					Options[Head[theGraphic]];
	rotateLights = RotateLights /. localOptions;
	
	thetaStart = spinRange[[1]];
	thetaStep = (spinRange[[2]] - spinRange[[1]]) /
					spinSteps;
	
	thetaFactor = Switch[rotateLights, 
						 False, 0, 
						 True, 1,
						 Opposite, -1, _, 0];
	phiFactor = If[rotateLights, 1, 0, 0];
	rhoFactor = If[rotateLights, 1, 0, 0];
	
	Do[Show[theGraphic, 
			ViewPoint->(spinOrigin +
				spinDistance Rotate3D[{1,0,0},
				thetaStart + thetaIndex thetaStep,
				spinTilt[[1]],spinTilt[[2]]]),
			LightSources -> 
				Map[{Rotate3D[#[[1]], 
					thetaFactor (thetaStart + 
								thetaIndex thetaStep),
					phiFactor spinTilt[[1]],
					rhoFactor spinTilt[[2]]], #[[2]]}&, 
					lightSources],
			Select[localOptions, FreeQ[#, SpinOrigin] &&
							FreeQ[#, SpinTilt] &&
							FreeQ[#, SpinDistance] &&
							FreeQ[#, SpinSteps] &&
							FreeQ[#, SpinRange] &&
							FreeQ[#, RotateLights]&]
			],
		{thetaIndex, 0, spinSteps-1}
      ];
	]

SpinShow[theGraphic]

7.6 A Fractalized Tetrahedron Is a Cube

Improved cross product code

The following cross product code, due to AndrÎ Deprit (see the letters in vol. 1, no. 2 of the Mathematica journal) is faster than that given in the routine that follows; substituting it will speed up the fractalization a little.

unitnormal[v1_List, v2_List] := (
		cross = RotateLeft[v1 RotateLeft[v2] - RotateLeft[v1] v2]) / 
				norm[cross]

norm[v_] := Sqrt[v . v]
unitnormal[v1_List, v2_List] :=
	(cross = Det[{IdentityMatrix[3], v1, v2}]) / norm[cross]


fractalize[{gray_, {p_, q_, r_}}] := Block[

{pq = (p+q)/2, qr = (r+q)/2, pr = (p+r)/2,
 new = (p+q+r)/3 + unitnormal[q-p, r-p]*norm[p-q]/Sqrt[6.]},

{Map[{gray, Polygon[#]} &, {{p,pq,pr}, {pq,q,qr}, {pr,qr,r}}],
 Map[{GrayLevel[newgray], 
	Polygon[#]} &, {{pr,pq,new}, {pq,qr,new}, {qr,pr,new}}]}]

Koch3D[n_] := (newgray = 1.;
vertices =
 {{Sqrt[3]/3,0,0}, {-Sqrt[3]/6,.5,0}, {-Sqrt[3]/6,-.5,0}}//N;
triangle = {GrayLevel[newgray], Polygon[vertices]};
tetrahedron = {Thickness[.0001],
			   Line[{{0, 0, 1/Sqrt[6.]}, #}] & /@ vertices};
Show[Graphics3D[{tetrahedron,
	Nest[(newgray -= .3; 
	 # /. {gray_, Polygon[v_]} -> fractalize[{gray, v}]) &,
		triangle, n]}],
		Lighting->False, Boxed->False,
		BoxRatios->{Sqrt[3]/2, 1, 1/Sqrt[6]},
		ViewPoint->({2, .3, 1.3}),
		PlotRange->{{-0.3, 0.6}, {-0.52, 0.52}, {-.02, 0.43}}])

Koch3D[2];

cube = {Thickness[.0001], Line /@ Join[
	Transpose[{
		Flatten[Outer[List, {-1,1}, {-1,1}, {-1}], 2],
		Flatten[Outer[List, {-1,1}, {-1,1}, { 1}], 2]}],
	Transpose[{
		Flatten[Outer[List, {-1,1}, {-1}, {-1,1}], 2],
		Flatten[Outer[List, {-1,1}, {+1}, {-1,1}], 2]}],
	Transpose[{
		Flatten[Outer[List, {-1}, {-1,1}, {-1,1}], 2],
		Flatten[Outer[List, {1},  {-1,1}, {-1,1}], 2]}]] //N};

KochPlanet[n_] := (newgray=.95;
tetrahedron = Map[{GrayLevel[newgray], Polygon[#]} &, {
	{{-1,-1,-1}, {1,-1,1}, {-1,1,1}},
	{{-1,-1,-1}, {1,1,-1}, {1,-1,1}},
	{{1,-1,1},   {1,1,-1}, {-1,1,1}},
	{{-1,-1,-1}, {-1,1,1}, {1,1,-1}} } //N];

Show[Graphics3D[{cube,
	Nest[(newgray -= .3; 
	 	# /. {gray_, Polygon[v_]} -> fractalize[{gray, v}]) &,
	  tetrahedron, n]}],
Lighting->False, Boxed->False,
ViewPoint->{-1.65, 1.4, .9},
BoxRatios->{1, 1, 1}])

Attributes[KochPlanet] = Listable;

KochPlanet[2];

ColorVersion of Fractal Tetrahedron as illustrated on cover of book  (animate the four frames!)

code

colors = {{1,.95,.1},{.3,.3,1},{.7,.3,1},{1,.1,.33}};
Show[Graphics[Table[
{RGBColor@@colors[[j]],
 Polygon[{{j,0},{j,1},{j+1,1},{j+1,0}}]}, {j,4}]]];

norm[v_] := Sqrt[v . v]
unitnormal[v1_List, v2_List] :=
	(cross = Det[{IdentityMatrix[3], v1, v2}]) / norm[cross]
	
fractalize[{gray_, {p_, q_, r_}}] := Block[
	{pq = (p+q)/2, qr = (r+q)/2, pr = (p+r)/2,
	 new = (p+q+r)/3 + unitnormal[q-p, r-p]*norm[p-q]/Sqrt[6.]},
{
Map[{gray, Polygon[#]} &, {{p,pq,pr}, {pq,q,qr}, {pr,qr,r}}],
Map[{RGBColor@@colors[[i]], 
	Polygon[#]} &, {{pr,pq,new}, {pq,qr,new}, {qr,pr,new}}]}]
cube = {Thickness[.0001], Line /@ Join[
	Transpose[{
		Flatten[Outer[List, {-1,1}, {-1,1}, {-1}], 2],
		Flatten[Outer[List, {-1,1}, {-1,1}, { 1}], 2]}],
	Transpose[{
		Flatten[Outer[List, {-1,1}, {-1}, {-1,1}], 2],
		Flatten[Outer[List, {-1,1}, {+1}, {-1,1}], 2]}],
	Transpose[{
		Flatten[Outer[List, {-1}, {-1,1}, {-1,1}], 2],
		Flatten[Outer[List, {1},  {-1,1}, {-1,1}], 2]}]] //N};


colors = {{1,.95,.1},{.3,.3,1},{.7,.3,1},{1,.1,.33}};


KochPlanet[n_] := (i=1;
tetrahedron = Map[{RGBColor@@colors[[1]], Polygon[#]} &, {
	{{-1,-1,-1}, {1,-1,1}, {-1,1,1}},
	{{-1,-1,-1}, {1,1,-1}, {1,-1,1}},
	{{1,-1,1},   {1,1,-1}, {-1,1,1}},
	{{-1,-1,-1}, {-1,1,1}, {1,1,-1}} } //N];

Show[Graphics3D[{cube,
		Nest[(i++; 
	 		# /. {gray_, Polygon[v_]} -> fractalize[{gray, v}]) &,
		tetrahedron, n]}],
Lighting -> False, Boxed -> False,
ViewPoint -> {-1.65, 1.4, .9},
BoxRatios -> {1, 1, 1}])
Attributes[KochPlanet] = {};

7.7 A Pretty Example

$Display = OpenWrite["psfile"];
n = 360;
ParametricPlot3D[
 {mesh (r = 2+Sin[8 phi]*.7) Cos[theta = phi + Sin[16 phi]/4], 
  mesh r Sin[theta],
  mesh * (2 - r/5),
  GrayLevel[(3.3 - r)/2]},
	{phi, 0, 2 Pi, 2 Pi / n}, {mesh, 0, 1, 1/5},
		Boxed->False,
		Lighting->False,
		BoxRatios->{1, 1, 1},
		ViewPoint->{0, 0, 6}];

Close["psfile"]; $Display = "stdout";

code for color version in Color Plate 4

(* uses ordinary ParametricPlot3D with True Lighting *)

n = 540;
ParametricPlot3D[
 {mesh (r = 2 + Sin[8 phi]*.7) Cos[theta = phi + Sin[16 phi]/4], 
  mesh r Sin[theta],
  mesh * (2 - r/5)},

	{phi, 0., 2 Pi, 2 Pi/n//N}, {mesh, 0, 1, 1},
		Boxed -> False,
		Lighting -> True,
		BoxRatios -> {1, 1, 1},
		Background->GrayLevel[0],
		ViewPoint -> ({0, 0, 6})];