(*********************************************************************

	Copyright 1989 by Roman E. Maeder
	
	Adapted from
	Roman E. Maeder: Programming in Mathematica, Addison-Wesley, 1989.

	Permission is hereby granted to make copies of this file for
	any purpose other than direct profit, or as part of a
	commercial product, provided this copyright notice is left
	intact.  Sale, other than for the cost of media, is prohibited.
	
	Permission is hereby granted to reproduce part or all of
	this file, provided that the source is acknowledged.

 *********************************************************************)


BeginPackage["RMPackages`NewParametricPlot3D`"]

ParametricPlot3D::usage =
	"ParametricPlot3D[{x,y,z,(style)}, {u,u0,u1,(du)}, {v,v0,v1,(dv)}, (options..)]
	plots a 3D parametric surface. Options are passed to Show[]"

PointParametricPlot3D::usage =
	"PointParametricPlot3D[{x,y,z}, {u,u0,u1,(du)}, {v,v0,v1,(dv)}, (options..)]
	plots a two-parameter set of points in space. Options are passed to Show[]."

SpaceCurve::usage = "SpaceCurve[{x,y,z}, {u,u0,u1,(du)}, (options..)]
	plots a 3D parametric curve. Options are passed to Show[]."

PointSpaceCurve::usage = "PointSpaceCurve[{x,y,z}, {u,u0,u1,(du)}, (options..)]
	plots a one-parameter set of points in space. Options are passed to Show[]"

SphericalPlot3D::usage = "SphericalPlot3D[r, {theta-range}, {phi-range}, (options...)]
	plots r as a function of the angles theta and phi.
	SphericalPlot3D[{r, style}, ...] uses style to render each surface patch"

CylindricalPlot3D::usage = "CylindricalPlot3D[z, {r-range}, {phi-range}, (options...)]
	plots z as a function of r and phi.
	CylindricalPlot3D[{z, style},  ...] uses style to render each surface patch"

Begin["`Private`"]

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] ]
	]  /; TensorRank[vl] == 3 && Dimensions[vl][[3]] == 3

MakePolygons[vl_List] :=
	Block[{l, l1, mesh, cols},
		l = Map[Take[#, 3]&, vl, {2}]; (* the coords *)
		l1 = Map[RotateLeft, l];
		cols = Map[#[[4]]&, vl, {2}];  (* the colors *)
		mesh = {l, RotateLeft[l], RotateLeft[l1], l1};
		mesh = Map[Drop[#, -1]&, mesh, {1}];
		mesh = Map[Drop[#, -1]&, mesh, {2}];
		cols = Drop[cols, -1];
		cols = Map[Drop[#, -1]&, cols, {1}];
		mesh = Polygon /@ Transpose[ Map[Flatten[#, 1]&, mesh] ];
		Flatten[Transpose[{Flatten[cols], mesh}]]
	]  /; TensorRank[vl] == 3 && Dimensions[vl][[3]] == 4


FilterOptions[ command_Symbol, opts___ ] :=
	Block[{keywords = First /@ Options[command]},
		Sequence @@ Select[ {opts}, MemberQ[keywords, First[#]]& ]
	]

Attributes[ParametricPlot3D] = {HoldFirst}

ParametricPlot3D[ fun_,
		{u_, u0_, u1_, du_:Automatic}, {v_, v0_, v1_, dv_:Automatic}, opts___ ] :=
	Block[{plotpoints, ndu = N[du], ndv = N[dv]},
		plotpoints = PlotPoints /. {opts} /. Options[Plot3D];
		If[ du === Automatic, ndu = N[(u1-u0)/(plotpoints-1)] ];
		If[ dv === Automatic, ndv = N[(v1-v0)/(plotpoints-1)] ];
		Show[ Graphics3D[MakePolygons[Table[ N[fun],
		                                     {u, u0, u1, ndu}, {v, v0, v1, ndv}] ]],
		      FilterOptions[Graphics3D, opts] ]
	]  /; NumberQ[N[u0]] && NumberQ[N[u1]] && NumberQ[N[v0]] && NumberQ[N[v1]]


Attributes[PointParametricPlot3D] = {HoldFirst}

PointParametricPlot3D[ fun_,
		{u_, u0_, u1_, du_:Automatic}, {v_, v0_, v1_, dv_:Automatic}, opts___ ] :=
	Block[{plotpoints, ndu = N[du], ndv = N[dv]},
		plotpoints = PlotPoints /. {opts} /. Options[Plot3D];
		If[ du === Automatic, ndu = N[(u1-u0)/(plotpoints-1)] ];
		If[ dv === Automatic, ndv = N[(v1-v0)/(plotpoints-1)] ];
		Show[ Graphics3D[Table[ Point[N[fun]], {u, u0, u1, ndu}, {v, v0, v1, ndv} ]],
		      FilterOptions[Graphics3D, opts] ]
	]  /; NumberQ[N[u0]] && NumberQ[N[u1]] && NumberQ[N[v0]] && NumberQ[N[v1]]


Attributes[SpaceCurve] = {HoldFirst}

SpaceCurve[ fun_, ul:{_, u0_, u1_, du_}, opts___ ] :=
	Show[ Graphics3D[Line[Table[ N[fun], ul ]]], FilterOptions[Graphics3D, opts] ] /;
			NumberQ[N[u0]] && NumberQ[N[u1]] && NumberQ[N[du]]

SpaceCurve[ fun_, {u_, u0_, u1_}, opts___ ] :=
    Block[{plotpoints},
    	plotpoints = PlotPoints /. {opts} /. Options[Plot3D];
	SpaceCurve[ fun, {u, u0, u1, (u1-u0)/(plotpoints-1)}, opts ]
    ]


Attributes[PointSpaceCurve] = {HoldFirst}

PointSpaceCurve[ fun_, ul:{_, u0_, u1_, du_}, opts___ ] :=
	Show[ Graphics3D[Table[ Point[N[fun]], ul ]], FilterOptions[Graphics3D, opts] ] /;
			NumberQ[N[u0]] && NumberQ[N[u1]] && NumberQ[N[du]]

PointSpaceCurve[ fun_, {u_, u0_, u1_}, opts___ ] :=
    Block[{plotpoints},
    	plotpoints = PlotPoints /. {opts} /. Options[Plot3D];
	PointSpaceCurve[ fun, {u, u0, u1, (u1-u0)/(plotpoints-1)}, opts ]
    ]


Attributes[SphericalPlot3D] = {HoldFirst}

SphericalPlot3D[ {r_, style_}, tlist:{theta_, __}, plist:{phi_, __}, opts___ ] :=
	Block[{rs},
		ParametricPlot3D[ {(rs = r) Sin[theta] Cos[phi],
		                   rs Sin[theta] Sin[phi],
		                   rs Cos[theta],
		                   style},
		                  tlist, plist, opts ]
	]

SphericalPlot3D[ r_, tlist:{theta_, __}, plist:{phi_, __}, opts___ ] :=
      ParametricPlot3D[ r{Sin[theta] Cos[phi],
                          Sin[theta] Sin[phi],
                          Cos[theta]},
                        tlist, plist, opts ]


Attributes[CylindricalPlot3D] = {HoldFirst}

CylindricalPlot3D[ {z_, style_}, rlist:{r_, __}, plist:{phi_, __}, opts___ ] :=
	ParametricPlot3D[{r Cos[phi], r Sin[phi], z, style}, rlist, plist, opts]

CylindricalPlot3D[ z_, rlist:{r_, __}, plist:{phi_, __}, opts___ ] :=
	ParametricPlot3D[{r Cos[phi], r Sin[phi], z}, rlist, plist, opts]

End[]

Protect[ParametricPlot3D, PointParametricPlot3D, SpaceCurve,
	PointSpaceCurve, SphericalPlot3D, CylindricalPlot3D]

EndPackage[]
