(* :Title: Surface of Revolution of Curves *)

(* :Context: Graphics`SurfaceOfRevolution` *)

(* :Author: John M. Novak *)

(* :Summary:
This package provides functions for plotting surfaces generated by
revolving a curve about an axis.  The curve may be expressed as a
function, parametrically, or as a list of points.  Rotation can 
be about an arbitrary axis.
*)

(* :Package Version: 2.1 *)

(* :Mathematica Version: 2.0 *)

(* :Copyright: Copyright 1990-1999, Wolfram Research, Inc.*)

(* :History:
	V1.0 by Kevin McIsaac, November 1990.
	Extensively modified by John M. Novak, April 1991.
	V2.0, Complete rewrite, John M. Novak, February 1992.
	V2.1, Changed Axis to RevolutionAxis, E.C. Martin, October 1992.
*)

(* :Keywords:
	Surface of Revolution, Curve
*)

(* :Sources:
	Rogers, David F. and Adams, J. Alan, Mathematical Elements for
		Computer Graphics, McGraw-Hill, 1976.
*)

BeginPackage["Graphics`SurfaceOfRevolution`",
	"Utilities`FilterOptions`"];

SurfaceOfRevolution::usage =
"SurfaceOfRevolution[fun, {u, u0, u1}, (options...)]
plots the surface of revolution of the curve defined by
the function fun (the dependent variable is in the z plane).
SurfaceOfRevolution[{xfun, zfun}, {u, u0, u1}, (options...)]
plots the surface of revolution of the curve defined 
parametrically by {xfun, zfun}.
SurfaceOfRevolution[{xfun, yfun, zfun}, {u, u0, u1}] is also
allowed.  The range of revolution can be specified after the
range of the variable u by {t, tmin, tmax}.  SurfaceOfRevolution
accepts the option RevolutionAxis, along with ParametricPlot3D options.";

ListSurfaceOfRevolution::usage =
"ListSurfaceOfRevolution[pts, ({t, tmin, tmax})] will generate
a rotated surface from the list of points pts (either pairs in
the x-z plane, or triplets in space), with the variable of
rotation being t, ranging from tmin to tmax (if not specified,
the rotation will be all the way around the axis).  ListSurfaceOfRevolution
accepts the option RevolutionAxis, along with Graphics3D options.";

RevolutionAxis::usage =
"RevolutionAxis is an option for the SurfaceOfRevolution and
ListSurfaceOfRevolution functions.  It is either a pair (representing an axis
in the x-z plane) or a triplet (representing an axis in space).  The default is
{0,0,1} (the z-axis).";

Begin["`Private`"]

(* private numberQ - phase out for NumericQ after V3.0 *)
numberQ[n_] := NumberQ[N[n]]

Options[SurfaceOfRevolution] =
	Join[{RevolutionAxis -> {0,0,1}},
	    Options[ParametricPlot3D]];

SurfaceOfRevolution[func_?(Head[#] =!= List &),
		range:{x_Symbol, _, _},rest___] :=
	SurfaceOfRevolution[{x,0,func},range,rest]

SurfaceOfRevolution[{x_, z_}, rest___] :=
	SurfaceOfRevolution[{x, 0, z}, rest]

SurfaceOfRevolution[{x_,y_,z_},
		vrange:{_Symbol, _?numberQ, _?numberQ},
		Optional[trange:{t_Symbol, _?numberQ, _?numberQ},
			{Unique[],0,2 Pi}],
		opts___?OptionQ] :=
	ParametricPlot3D[Evaluate[
		{x,y,z} . rotationmatrix[RevolutionAxis/.Flatten[{opts}]/.
				Options[SurfaceOfRevolution],
			t]],
		vrange, trange,
		Evaluate[FilterOptions @@
		    Flatten[{ParametricPlot3D,opts, 
		        Options[SurfaceOfRevolution]}]]
	]

Options[ListSurfaceOfRevolution] =
	Join[{PlotPoints -> 15,
	RevolutionAxis -> {0,0,1}},
	Options[Graphics3D]];

SetOptions[ListSurfaceOfRevolution, Axes -> True]

ListSurfaceOfRevolution[
		pts:{{_?numberQ, _?numberQ}..},
		rest___] :=
	ListSurfaceOfRevolution[
		Map[{First[#],0,Last[#]}&,pts],
		rest]

ListSurfaceOfRevolution[
		pts:{{_?numberQ, _?numberQ, _?numberQ}..},
		Optional[trange:{t_Symbol, min_?numberQ, max_?numberQ},
			{Unique[],0,2 Pi}],
		opts___?OptionQ] :=
	Module[{pp, axis, rmatrix},
		{pp, axis} = {PlotPoints, RevolutionAxis}/.
		        Flatten[{opts}]/.
			Options[ListSurfaceOfRevolution];
		listsurfaceplot[Table[rmatrix = rotationmatrix[axis,t];
				Map[# . rmatrix &, pts],
				{t, min, max, (max - min)/(pp - 1)}],
			FilterOptions @@ Flatten[{Graphics3D,opts,
			                Options[ListSurfaceOfRevolution]}]
		]
	]
		

rotationmatrix[{a_,b_}, theta_] :=
	rotationmatrix[{a,0,b},theta]

rotationmatrix[axis_,theta_] :=
	Module[{n1,n2,n3},
		{n1,n2,n3} = axis/Sqrt[Plus @@ (axis^2)]//N;
		{{n1^2 + (1 - n1^2) Cos[theta],
			n1 n2 (1 - Cos[theta]) + n3 Sin[theta],
			n1 n3 (1 - Cos[theta]) - n2 Sin[theta]},
		{n1 n2 (1 - Cos[theta]) - n3 Sin[theta],
			n2^2 + (1 - n2^2) Cos[theta],
			n2 n3 (1 - Cos[theta]) + n1 Sin[theta]},
		{n1 n3 (1 - Cos[theta]) + n2 Sin[theta],
			n2 n3 (1 - Cos[theta]) - n1 Sin[theta],
			n3^2 + (1 - n3^2) Cos[theta]}}//N
	]

(* What follows is a modified version of
	ListSurfacePlot3D from Graphics3D.m *)

listsurfaceplot[vl_,opts___] :=
	Module[{l = vl, l1 = Map[RotateLeft, vl], mesh},
		mesh = {l, l1, RotateLeft[l1], RotateLeft[l]};
		mesh = Map[Drop[#, -1]&, mesh, {1}];
		mesh = Map[Drop[#, -1]&, mesh, {2}];
		Show[Graphics3D[
			Polygon /@ Transpose[ Map[Flatten[#, 1]&, mesh] ]],
			opts]
	]

End[]

EndPackage[]
