(*^
::[	Information =

	"This is a Mathematica Notebook file.  It contains ASCII text, and can be
	transferred by email, ftp, or other text-file transfer utility.  It should
	be read or edited using a copy of Mathematica or MathReader.  If you 
	received this as email, use your mail application or copy/paste to save 
	everything from the line containing (*^ down to the line containing ^*)
	into a plain text file.  On some systems you may have to give the file a 
	name ending with ".ma" to allow Mathematica to recognize it as a Notebook.
	The line below identifies what version of Mathematica created this file,
	but it can be opened using any other version as well.";

	FrontEndVersion = "Macintosh Mathematica Notebook Front End Version 2.2";

	MacintoshStandardFontEncoding; 
	
	fontset = title, inactive, nohscroll, noKeepOnOnePage, preserveAspect, groupLikeTitle, center, M7, bold, e8,  24, "Times"; 
	fontset = subtitle, inactive, nohscroll, noKeepOnOnePage, preserveAspect, groupLikeTitle, center, M7, bold, e6,  18, "Times"; 
	fontset = subsubtitle, inactive, nohscroll, noKeepOnOnePage, preserveAspect, groupLikeTitle, center, M7, italic, e6,  14, "Times"; 
	fontset = section, inactive, nohscroll, noKeepOnOnePage, preserveAspect, groupLikeSection, grayBox, M22, bold, a20,  18, "Times"; 
	fontset = subsection, inactive, nohscroll, noKeepOnOnePage, preserveAspect, groupLikeSection, blackBox, M19, bold, a15,  14, "Times"; 
	fontset = subsubsection, inactive, nohscroll, noKeepOnOnePage, preserveAspect, groupLikeSection, whiteBox, M18, bold, a12,  12, "Times"; 
	fontset = text, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  14, "Times"; 
	fontset = smalltext, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Times"; 
	fontset = input, nowordwrap, preserveAspect, groupLikeInput, M42, N23, bold,  12, "Courier"; 
	fontset = output, output, inactive, nowordwrap, preserveAspect, groupLikeOutput, M42, N23, L-5,  12, "Courier"; 
	fontset = message, inactive, nowordwrap, preserveAspect, groupLikeOutput, M42, N23,  12, "Courier"; 
	fontset = print, inactive, nowordwrap, preserveAspect, groupLikeOutput, M42, N23,  12, "Courier"; 
	fontset = info, inactive, nowordwrap, preserveAspect, groupLikeOutput, M42, N23,  12, "Courier"; 
	fontset = postscript, PostScript, formatAsPostScript, output, inactive, nowordwrap, preserveAspect, groupLikeGraphics, M7, l34, w282, h287,  12, "Courier"; 
	fontset = name, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7, italic,  10, "Times"; 
	fontset = header, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Geneva"; 
	fontset = leftheader, nohscroll, cellOutline, L2,  12, "Palatino"; 
	fontset = footer, inactive, nohscroll, noKeepOnOnePage, preserveAspect, center, M7,  12, "Geneva"; 
	fontset = leftfooter, cellOutline, blackBox, L2,  12, "Palatino"; 
	fontset = help, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  10, "Times"; 
	fontset = clipboard, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Geneva"; 
	fontset = completions, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Courier"; 
	fontset = special1, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Geneva"; 
	fontset = special2, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Geneva"; 
	fontset = special3, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Geneva"; 
	fontset = special4, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Geneva"; 
	fontset = special5, inactive, nohscroll, noKeepOnOnePage, preserveAspect, M7,  12, "Geneva"; 
	paletteColors = 128; showRuler; automaticGrouping; currentKernel; 
]
:[font = title; inactive; preserveAspect; startGroup]

ImplicitSolids.m
:[font = subsubtitle; inactive; preserveAspect]
Duane Storti and Mark Ganter

Department of Mechanical Engineering, FU-10
University of Washington
Seattle, WA 98195

and

Manufacturing Systems Technology Center
Washington Technology Centers
Seattle, WA 98195

storti@u.washington.edu
ganter@u.washington.edu
;[s]
3:0,0;212,1;259,0;260,-1;
2:2,16,12,Times,2,14,0,0,0;1,13,10,Courier,2,12,0,0,0;
:[font = section; inactive; pageBreak; preserveAspect; startGroup]
Implementation
:[font = subsection; inactive; noPageBreak; preserveAspect; startGroup]
Preparation
:[font = subsubsection; inactive; preserveAspect; startGroup]
Create the context for this Notebook:
:[font = input; initialization; preserveAspect]
*)
BeginPackage["ImplicitSolids`"];
(*
:[font = input; initialization; preserveAspect]
*)
Needs[ "Graphics`ParametricPlot3D`" ] ; 

(*
:[font = input; initialization; preserveAspect; endGroup; endGroup]
*)
Needs[ "Geometry`Rotations`" ] ;
(*
:[font = subsection; inactive; preserveAspect; startGroup]
Usage messages
:[font = input; initialization; wordwrap; noKeepOnOnePage; preserveAspect]
*)
ImplicitSolids::usage = "ImplicitSolids is an implicit function solid modeling package.";
(*
:[font = subsubsection; inactive; preserveAspect; startGroup]
Coordinate Names
:[font = input; initialization; noKeepOnOnePage; preserveAspect; endGroup]
*)
x::usage = "Reserved for use as a Cartesian coordinate.";
y::usage = "Reserved for use as a Cartesian coordinate.";
z::usage = "Reserved for use as a Cartesian coordinate.";
(*
:[font = subsubsection; inactive; noPageBreakInGroup; preserveAspect; startGroup]
"Primitive" solids
:[font = input; initialization; noPageBreakInGroup; wordwrap; preserveAspect; endGroup]
*)
Sphere::usage = "Sphere creates a defining expression for a sphere of radius 1 centered at the origin.";

Cuboid::usage = "Cuboid[ (n: 6) ] creates a defining expression for an nth-degree cuboid with edge length ~=2 and centroid at the origin.  Cuboids have rounded corners which become sharper as n is increased.";

Superellipsoid::usage = "Superellipsoid[ e1, e2 ] creates a defining expression for a superellipsoid: \n
(x^(2/e2) + y^(2/e2))^(e2/e1) + y^(2/e1) ";

Extrude::usage = "Extrude[f(x,y), a , (n: 6) ] creates a defining expression for a degree-n extruded solid of length a along the z-axis with mid-section given by f(x,y).";
(*
:[font = subsubsection; inactive; noPageBreakInGroup; preserveAspect; startGroup]
"Primitive" Sections
:[font = input; initialization; noPageBreakInGroup; wordwrap; noKeepOnOnePage; preserveAspect; endGroup]
*)
Circle::usage = "Circle defines a circle of unit radius  centered at the origin of the x-y plane.";

Squaroid::usage = "Squaroid[ (n:6) ] defines a degree n squaroid of edge 2 centered at the origin of the x-y plane. Squaroid (without arguments) defines a squaroid of degree DEFAULTDEGREE. Squaroids have rounded corners which become sharper as n is increased.";
(*
:[font = subsubsection; inactive; preserveAspect; startGroup]
Movement and Scaling operators
:[font = input; initialization; wordwrap; noKeepOnOnePage; preserveAspect]
*)
Translate::usage = "Translate[ f, {a, b, c}] adjusts the defining expression f(x,y,z) corresponding to translating the object with displacement vector {a, b, c}.";

Rotate::usage = "Rotate[f, coordinateAxis, angle] performs a coordinate transformation on the defining expression f(x,y,x) corresponding to rotating the solid by the specified angle (in radians) about the specified coordinate axis (x, y, or z). To specify the angle in degrees, use the built in constant Degree; e.g. Pi/6 and 30 Degree are equivalent values for the third argument." ;

Permute::usage = "Permute[ f ] permutes x, y, z to y, z, x to realign extrusions along other axes.";
(*
:[font = input; initialization; wordwrap; noKeepOnOnePage; preserveAspect]
*)
Stretch::usage = "Stretch[ f, {a, b, c}] adjusts the defining expression to stretch the solid by factors a, b, and c in the x,y and z directions respectively.";
(*
:[font = input; initialization; wordwrap; noKeepOnOnePage; preserveAspect; endGroup]
*)
Magnify::usage = "Magnify[ f, m ] stretches all dimensions by a factor of m.";
(*
:[font = subsubsection; inactive; pageBreak; preserveAspect; startGroup]
 Algebraic Pseudo-Boolean operators
:[font = input; initialization; wordwrap; noKeepOnOnePage; preserveAspect; endGroup]
*)
Intersect::usage = "Intersect[{f1,..., fm}, (n: 6) ] produces a defining expression for an nth-order approximation of the intersection of the objects defined by f1 through fm. Intersect produces a smoothly blended surface which approaches an exact Boolean intersection (with possibly sharp edges) as n is increased.";

Unite::usage = "Unite[{f1,..., fm}, (n: 6) ] produces a defining expression for an nth-degree approximation of the union of the objects defined by f1 through fm. Unite produces a smoothly blended surface which approaches an exact Boolean union (with possibly sharp edges) as n is increased.";

Difference::usage = "Difference[{f1,...fm},{g1,..., gn}, (n: 6) ] produces a defining expression for an nth-degree approximation of the intersection of the objects defined by f1,...,fm with the complement of the union of objects defined by g1,...,gn; i.e. the union of the f's with the union of the g's removed. Difference produces a smoothly blended surface which approaches an exact Boolean difference (with possibly sharp edges) as n is increased.";

Complement::usage = "Complement[f] produces a defining expression for the complement of the object defined by f.";

(*
:[font = subsubsection; inactive; pageBreak; preserveAspect; startGroup]
Visualization functions
:[font = input; initialization; wordwrap; noKeepOnOnePage; preserveAspect; endGroup]
*)
SphericalRender::usage = "SphericalRender[f, max, ndiv, options] produces a SphericalPlot3D plot of the defining expression f(x,y,z) after xyzToPolar transforms f to spherical polar coordinates. Options are passed to Show via SphericalPlot3D. This function should only be used to visualize objects which have a star-shaped surface (i.e., have one intersection with any ray of length max extending from the origin). For other objects, use RayRender.";

RayRender::usage = "RayRender[f, {x, y, z}, \n
	{xMin, xMax, xDivisions}, \n
	{yMin, yMax, yDivisions}, \n
	{zMin, zMax, zDivisions}] uses raycasting to visualize the object defined by the function f which depends on variables x, y, and z. The viewer looks along the x-axis from a great distance so that the viewing frustum is the right rectangular parallelepiped with front face x=xMin, rear face x=xMax, left face y=ymin, right face y=yMax, bottom face z=zMin and top face z=zMax. yDivisions and zDivisions control the number of pixel values which are calculated along each row and column respectively. Since RayRender can take a long time, as each row of pixels is calculated, a message is printed updating the progress of the computation.";
(*
:[font = subsubsection; inactive; preserveAspect; startGroup]
Coordinate transformation 

:[font = input; initialization; wordwrap; noKeepOnOnePage; preserveAspect; endGroup; endGroup]
*)
xyzToPolar::usage = "xyzToPolar[f] converts a function f(x,y,z) to its spherical polar version f(r,theta, phi).";
(*
:[font = subsection; inactive; pageBreak; preserveAspect; startGroup]
Definitions
:[font = input; initialization; preserveAspect]
*)
Begin["`Private`" ];
(*
:[font = input; initialization; preserveAspect]
*)
Off[General::spell1]; (* Turn off spelling warnings *)
(*
:[font = input; initialization; preserveAspect]
*)
protected = Unprotect[ Circle, Complement, Cuboid ];

(*
:[font = input; initialization; preserveAspect]
*)
Clear[x, y, z]; 

(*
:[font = text; inactive; preserveAspect]
x, y, and z are reserved as Cartesian coordinate names. 

Any previous definitions are removed.
:[font = input; initialization; preserveAspect]
*)
DEFAULTDEGREE = 6;
(*
:[font = text; inactive; preserveAspect]
DEFAULTDEGREE  is the default value for the degree of the approximate Boolean operators and some "primitive" solids and sections.  Increasing DEFAULTDEGREE gives sharper corners (the blending surfaces "cling" more tightly to the objects being joined) and longer computing time for visualizations. 

IMPORTANT NOTE: Choosing an odd value of the blending degree can cause problems if defining expressions are allowed to have negative values!
:[font = subsubsection; inactive; preserveAspect; startGroup]
Primitive Solids
:[font = input; initialization; preserveAspect]
*)
Sphere := (x^2 + y^2 + z^2);
(*
:[font = input; initialization; preserveAspect]
*)
Cuboid[n_: DEFAULTDEGREE] := x^n + y^n + z^n;
(*
:[font = input; initialization; preserveAspect]
*)
Extrude[ f_, a_, n_: DEFAULTDEGREE] := f + (2z/a)^n;
(*
:[font = input; initialization; preserveAspect; endGroup]
*)
Superellipsoid[e1_, e2_] := 

	( (x^2)^(1/e2) + (y^2)^(1/e2) )^(e2/e1) + (z^2)^(1/e1) ;
(*
:[font = subsubsection; inactive; preserveAspect; startGroup]
Sections
:[font = input; initialization; preserveAspect]
*)
Circle := (x^2 +y^2);
(*
:[font = input; initialization; preserveAspect; endGroup]
*)
Squaroid[n_:DEFAULTDEGREE] := (x^n + y^n);
(*
:[font = subsubsection; inactive; pageBreak; preserveAspect; startGroup]
Movement operators
:[font = input; initialization; preserveAspect]
*)
Translate[expr_, {xcentroid_, ycentroid_, zcentroid_}] := 

    Module[{tmp},
		MakeFun[expr, 

	    tmp[x, y, z]][x-xcentroid,y-ycentroid,z-zcentroid]
		];
(*
:[font = input; initialization; preserveAspect]
*)
Rotate[expr_, axis_, angle_] :=
    Module[{coordList, vec, transformationMatrix,rules},
		coordList = {x,y,z};
		While[ Not[SameQ[First[coordList],axis]], 

			coordList = RotateRight[coordList] ];
		vec = Drop[ coordList, 1];
		transformationMatrix = N[RotationMatrix2D[angle]];
		rules = Table[
			vec[[i]]->(transformationMatrix . vec)[[i]],
			{i,2}];
		expr /. rules] /; MemberQ[{x,y,z}, axis] ;
(*
:[font = input; initialization; preserveAspect]
*)
Rotate[expr_, axis_, angle_] :=
	( Message[ Rotate::badAxis, axis] ;
		/; Not[ MemberQ[{x,y,z}, axis] ] );
(*
:[font = input; initialization; preserveAspect]
*)
Rotate::badAxis = "`1` is not an allowed rotation axis. \n 

The second argument must be  x, y, or z.";
(*
:[font = input; initialization; preserveAspect]
*)
Permute[f_] := 

	Module[{tmp},
		MakeFun[f, tmp[x, y, z]][y, z, x] ];
(*
:[font = input; initialization; preserveAspect]
*)
Stretch[f_, {a_, b_, c_}] := 

	Module[{tmp},
		MakeFun[f, tmp[x, y, z]][x/a, y/b, z/c] ];
(*
:[font = input; initialization; preserveAspect; endGroup]
*)
Magnify[f_, n_] := Stretch[f, {n, n, n} ];
(*
:[font = subsubsection; inactive; preserveAspect; startGroup]
Coordinate transformation
:[font = input; initialization; preserveAspect; endGroup]
*)
xyzToPolar[f_] := 

    Module[{tmp},
		MakeFun[f, tmp[x, y, z]][r Sin[theta] Cos[phi], 

			r Sin[theta] Sin[phi], r Cos[theta]]
		];
(*
:[font = subsubsection; inactive; pageBreak; preserveAspect; startGroup]
Algebraic Pseudo-Boolean Operators
:[font = input; initialization; preserveAspect]
*)
Intersect[ fList_List, n_: DEFAULTDEGREE ] := 

	Apply[Plus, Map[ (#^n)& , fList] ] ;
(*
:[font = input; initialization; preserveAspect]
*)
Difference[ fList_List, gList_List, n_: DEFAULTDEGREE] :=
	( 1 / Apply[Plus, Map[(#^-n)&, fList] ]+ 

		Apply[Plus, Map[(#^-n)&, gList] ] ) ;
(*
:[font = input; initialization; preserveAspect]
*)
Complement[ f_ ] := 1/f;
(*
:[font = input; initialization; preserveAspect; endGroup]
*)
Unite[ fList_List, n_: DEFAULTDEGREE ] := 

	1 / Apply[Plus, Map[ (#^-n)& , fList] ];
(*
:[font = subsubsection; inactive; preserveAspect; startGroup]
Visualization by 3D Spherical Plot
:[font = input; initialization; preserveAspect; endGroup]
*)
SphericalRender[f_, max_, ndiv_, opts___] := 

  Module[{}, 

	Off[FindRoot::precwarn];
	SphericalPlot3D[
		FindRoot[ xyzToPolar[f] == 1,
			{r, {0.05,max}}, AccuracyGoal -> 8, MaxIterations->25,
			Compiled -> False][[1,2]],
		{theta, 0, Pi, Pi/ndiv},
		{phi, 0, 2*Pi, Pi/ndiv}, opts] ;
	On[FindRoot::precwarn];
	Null];
(*
:[font = subsubsection; inactive; preserveAspect; startGroup]
Visualization by RayCasting
:[font = text; inactive; preserveAspect]
(1) Make function from expression
:[font = input; initialization; preserveAspect]
*)
Off[RuleDelayed::rhs]; 

(* turn off message which the next expression produces *)

MakeFun[expr_, 

    name_Symbol[xPar_Symbol, yPar_Symbol, zPar_Symbol]] := 

    	(name[xPar_, yPar_, zPar_] = expr; name);

On[RuleDelayed::rhs]; (* turn the message back on *)
(*
:[font = text; inactive; pageBreak; preserveAspect]
(2) Fix Root
:[font = text; inactive; preserveAspect]
Evaluate definining function and determine where it first decreases past unity  on [0,1]. This is achieved by simply stepping along and evaluating until a function value less than one is obtained, then interpolating. Note that objects thinner than (xMax-xMin)/xDivisions may be missed.	(The built-in rootfinder FindRoot cannot guarantee not missing even larger objects.) 

Also add new values to range of rootfinder to account for anomalous conditions, i.e. if no crossing is found  return BackGround and if the function is initially less than one, we are inside the object and a message is returned along with the value ForeGround.
:[font = input; initialization; preserveAspect]
*)
StepRoot[CompiledFunName_, y_, z_, {xMin_, xMax_, xDivisions_}] :=
  Module[{dx = (xMax-xMin)/xDivisions, 

  		fc = CompiledFunName, 

  		fold, fnew, xold, xnew=xMin},
	fnew = fc[xMin,y,z];
	Which[fnew<1, 

	    (Message[RayRender::badxMin, y, z]; ForeGround),
      	    fnew == 1., xMin,
      	    True,
       		(While[fnew>1,
        		(If[xnew > xMax, Return[BackGround] ];
        		xold = xnew; fold = fnew;
       			xnew += dx; fnew = fc[xnew,y,z])];
		N[xold + (xnew-xold)(1-fold)/(fnew-fold)]
			) 

		] 

	];
(*
:[font = input; initialization; preserveAspect]
*)
RayRender::badxMin = 

  "Object in front of viewing frustum at y = `1`, z = `2`.";
(*
:[font = text; inactive; pageBreak; preserveAspect]
(3) Calculate GrayLevel 
:[font = text; inactive; preserveAspect]
   (a) Direction Cosine function
:[font = input; initialization; preserveAspect]
*)
MakeDirCosFun[expr_, 

	name_Symbol[xPar_Symbol, yPar_Symbol, zPar_Symbol]] := 

    	MakeFun[N[ Abs[D[expr,xPar]] / 

		Sqrt[D[expr,xPar]^2 + D[expr,yPar]^2+ 

			D[expr,zPar]^2] ],
		name[xPar, yPar, zPar]];
(*
:[font = text; inactive; preserveAspect]
Create a compiled version of the direction cosine function for faster evaluation.
:[font = input; initialization; preserveAspect]
*)
MakeCompiledDirCosFun[expr_, 

	name_Symbol,xPar_Symbol, yPar_Symbol, zPar_Symbol] := 


	Module[{tmp},
		MakeDirCosFun[expr, tmp[x, y, z]];
		name = Compile[{xPar, yPar, zPar}, tmp[x,y,z] ] ]
(*
:[font = input; initialization; preserveAspect]
*)
MakeCompiledDirCosFun[expr_, 

	name_CompiledFunction,
	xPar_Symbol, yPar_Symbol, zPar_Symbol] := 


	Message[MakeCompiledDirCosFun::recompile];
(*
:[font = input; initialization; preserveAspect]
*)
MakeCompiledDirCosFun::recompile =
	"Second argument must be a symbol, not a CompiledFunction.";
(*
:[font = text; inactive; pageBreak; preserveAspect]
 	(b) Gray level calculation
:[font = input; initialization; preserveAspect]
*)
BackgroundGrayLevel = 0.; 

(* GrayLevel displayed where background is visible *)
(*
:[font = input; initialization; preserveAspect]
*)
ForegroundGrayLevel = 1.0;
(* GrayLevel displayed where the object lies *)
(* in front of viewing frustum *)
(*
:[font = input; initialization; preserveAspect]
*)
CalculateGrayLevel[CompiledDirCosFunName_, 

	BackGround, y_, z_] = BackgroundGrayLevel;
(*
:[font = input; initialization; preserveAspect]
*)
CalculateGrayLevel[CompiledDirCosFunName_, 

	ForeGround, y_, z_] := ForegroundGrayLevel;
(*
:[font = input; initialization; preserveAspect]
*)
CalculateGrayLevel[CompiledDirCosFunName_, 

	root_, y_, z_] := CompiledDirCosFunName[root, y, z]
     

(*
:[font = input; initialization; preserveAspect]
*)


DetermineGrayLevel[CompiledFunName_, 

		   CompiledDirCosFunName_, 

           y_, z_, xMin_, xMax_, xDivisions_] :=
	CalculateGrayLevel[CompiledDirCosFunName,
		StepRoot[CompiledFunName, y, z,
			{xMin, xMax, xDivisions}], y, z]
(*
:[font = text; inactive; preserveAspect]
(4) Construct GrayLevel Table 
:[font = input; initialization; preserveAspect]
*)
MakeTable[CompiledFunName_, CompiledDirCosFunName_, 

   xMin_, xMax_, xDivisions_,
   yMin_, yMax_, yDivisions_, 

   zMin_, zMax_, zDivisions_] :=
    

Module[{i = 0, j = 1,
		dy = (yMax - yMin)/yDivisions,
		dz = (zMax - zMin)/zDivisions},
    Table[
		 (	If[Mod[i++, yDivisions + 1]==0,
		 		Print["Computing row ", j, " of  
",zDivisions +1];
		 		j++;
		 		i = 1
		 		]; 

		 	DetermineGrayLevel[CompiledFunName, 

		 		CompiledDirCosFunName,  

		 		y, z, xMin, xMax, xDivisions]
		 	),
		{z, zMin, zMax, dz}, {y, yMin, yMax, dy}] ];
(*
:[font = text; inactive; pageBreak; preserveAspect]
(5) Driver
:[font = input; initialization; preserveAspect]
*)
RayCast[expr_,
        {xPar_, yPar_, zPar_},
        xMin_, xMax_, xDivisions_,
        yMin_, yMax_, yDivisions_, 

        zMin_, zMax_, zDivisions_] :=
  

(Clear[CompiledDefFun, CompiledDirCosFun, Intensity];
CompiledDefFun = 

	Compile[ {xPar, yPar, zPar}, expr];
MakeDirCosFun[expr, Intensity[xPar, yPar, zPar] ];	
CompiledDirCosFun = 

	Compile[ {xPar, yPar, zPar}, 

		Intensity[xPar, yPar, zPar]];

MakeTable[
	CompiledDefFun,
    CompiledDirCosFun,
   	xMin, xMax, xDivisions,
	yMin, yMax, yDivisions, 

	zMin, zMax, zDivisions] )
(*
:[font = input; initialization; preserveAspect]
*)
RayRender[Expression_,
          {xParName_, yParName_, zParName_},
          {xMin_, xMax_, xDivisions_}, 

          {yMin_, yMax_, yDivisions_}, 

          {zMin_, zMax_, zDivisions_}] :=
          

Show[
  Graphics[
	Raster[
	    RayCast[Expression,
			{xParName, yParName, zParName},
			xMin, xMax, xDivisions,
			yMin, yMax, yDivisions, 

			zMin, zMax, zDivisions],
	    {{0,0}, {1, 1}}]
	],
    AspectRatio -> (zMax - zMin) / (yMax - yMin)];   

(*
:[font = input; initialization; preserveAspect]
*)
Protect[ Evaluate[ protected ] ]; 

(* restore protection of system symbols *)
(*
:[font = input; initialization; preserveAspect; endGroup; endGroup]
*)
End[ ]; (* end of private context *)
(*
:[font = subsection; inactive; pageBreak; preserveAspect; startGroup]
Finish
:[font = input; initialization; preserveAspect]
*)
(* protect exported symbols *)


(* package name *)
Protect[ImplicitSolids]; 


(* primitives *)
Protect[Sphere, Cuboid, Extrude, Superellipsoid]; 


(*sections *)
Protect[Circle, Squaroid]; 


(* operators *)
Protect[Translate, Rotate, Permute, Stretch, Magnify];

(* blended pseudo-Booleans *) 

Protect[Intersect, Unite, Difference]; 


(* coordinate transformation *)
Protect[xyzToPolar]; 


(*visualization functions *)
Protect[SphericalRender, RayRender]; 

(*
:[font = input; initialization; preserveAspect]
*)
On[General::spell1]; (* Turn spelling warnings back on *)
(*
:[font = input; initialization; preserveAspect; endGroup; endGroup; endGroup]
*)
EndPackage[ ]; (* end of package context *)
(*
^*)
