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

(* :Context: Graphics`SurfaceOfRevolution` *)

(* :Author: John M. Novak *)

(* :Summary: Plots surfaces of revolution of curves, specified
	as a function or parametrically, or as a list of points.
	Rotation can be about arbitrary axis. *)

(* :Package Version: 2.0 *)

(* :Mathematica Version: 2.0 *)

(* :History:
	V1.0 by Kevin McIsaac, Nov. 1990
	Extensively modified by John M. Novak, April 1991
	V2.0, Complete rewrite, John M. Novak, Feb. 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[{x-fun,z-fun},{u,u0,u1},(opts...)]
plots the surface of revolution of the curve defined 
parametrically by {x-fun,z-fun}.
SurfaceOfRevolution[{x-fun, y-fun, z-fun}, {u, u0, u1}] is also
allowed. Also, the range of revolution can be specified after the
range of the variable u by {t,tmin,tmax}.  Accepts option 
specifying Axis, 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 used,
the rotation will be all the way around the axis.) Accepts
option specifying Axis, along with Graphics3D options.";

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

Begin["`Private`"]

Options[SurfaceOfRevolution] =
	{Axis -> {0,0,1}};

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[N[#]]&), _?(NumberQ[N[#]]&)},
		Optional[trange:{t_Symbol, _?(NumberQ[N[#]]&), _?(NumberQ[N[#]]&)},
			{Unique[],0,2 Pi}],
		opts:(Rule_ | RuleDelayed_)...] :=
	ParametricPlot3D[Evaluate[
		{x,y,z} . rotationmatrix[Axis/.{opts}/.
				Options[SurfaceOfRevolution],
			t]],
		vrange, trange,
		Evaluate[FilterOptions[ParametricPlot3D,opts]]
	]

Options[ListSurfaceOfRevolution] =
	{PlotPoints -> 15,
	Axis -> {0,0,1}};

ListSurfaceOfRevolution[
		pts_?(MatrixQ[N[#],NumberQ] &)/;Last[Dimensions[pts]] === 2,
		rest___] :=
	ListSurfaceOfRevolution[
		Map[{First[#],0,Last[#]}&,pts],
		rest]

ListSurfaceOfRevolution[
		pts_?(MatrixQ[N[#],NumberQ] &)/;Last[Dimensions[pts]] === 3,
		Optional[trange:{t_Symbol, min_?(NumberQ[N[#]]&), max_?(NumberQ[N[#]]&)},
			{Unique[],0,2 Pi}],
		opts:(Rule_ | RuleDelayed_)...] :=
	Module[{pp, axis, array, rmatrix},
		{pp, axis} = {PlotPoints, Axis}/.{opts}/.
			Options[ListSurfaceOfRevolution];
		listsurfaceplot[Table[rmatrix = rotationmatrix[axis,t];
				Map[# . rmatrix &, pts],
				{t, min, max, (max - min)/(pp - 1)}],
			FilterOptions[Graphics3D,opts],
			Axes -> True
		]
	]
		

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[]
