(* :Title: ImplicitPlot3D *) (* :Author: Steven Wilkinson, Northern Kentucky University *) (* :Summary: Plot3D[ ] plots a function of two variables. Many simple surfaces (e.g., ellipsoids, hyperboloids etc.) are not the graphs of functions of two variables. ImplicitPlot3D[ ] allows the user to plot surfaces implicitly defined by an equation of three variables. *) (* :Context: Graphics`ImplicitPlot3D` *) (* :Package Version: 2.0 *) (* :History: V2.0 by Steven Wilkinson, July 1992 *) (* :Keywords: solution set, 3D graphics, implicitly defined surfaces *) (* :Warning: The algorithm used by ImplicitPlot3D is essentially cubic with respect to the desired resolution. Therefore, it can be time and memory intensive. *) (* :Mathematica Version: 2.0 *) (* :Limitations: The package works by taking the box given by the range of the three variables and dividing it into a number of subboxes as specified by the PlotPoints option. Using the intermediate value theorem, the package checks each of the subboxes to see if the surface passes through any of its edges. The algorithm is actually a little faster than cubic because of the following modification. The user can have the package make a number of prepasses. The first pass divides the original box coarsely. The package retains only those subboxes in which it can detect the surface. These subboxes are then repeatedly subdivided in the same manner retaining only those new subboxes containing the surface. The limitations of this package (which are especially emphasized by the modification) is that if the surface intersects an edge of a subbox an even number of times, the package may not be able to discern that the surface has entered that subbox. The coarser the subdivisions, the more likely this is to happen.*) BeginPackage[ "Graphics`ImplicitPlot3D`","Utilities`FilterOptions`" ] ImplicitPlot3D::usage = "ImplicitPlot3D[eqn,{x,a,b},{y,c,d},{z,e,f}] draws a graph of the set of points with coordinates (x,y,z) in the box a2x2b, c2y2d, e2z2f that satisfy the equation eqn. ImplicitPlot3D[func,{x,a,b},{y,c,d},{z,e,f}] draws a graph of the set of points with coordinates (x,y,z) in the box a2x2b, c2y2d, e2z2f that satisfy the equation func = 0."; Passes::usage = "Passes is an option for ImplicitPlot3D. It is used to specify how many times to subdivide the plot range in order to reach the stated value of PlotPoints. This option is ignored if the PlotPoints option as a list of lists." Options[ImplicitPlot3D] = {Compiled -> True, PlotPoints -> {{5,5,5},{3,3,3}}, Passes -> Automatic}; Begin["`Private`"] vector = {{1,0,0},{0,1,0},{0,0,1},{0,0,0}}; vertex = Flatten[Table[{i,j,k},{i,0,1},{j,0,1},{k,0,1}],2]; edges = Select[ Flatten[ Table[{vertex[[i]],vector[[j]]},{i,8},{j,4}] ,1], (Max[#[[1]]+#[[2]]]<2)& ]; prepassEdges = Select[ Flatten[ Table[{vertex[[i]],vector[[j]]},{i,7},{j,3}] ,1], (Max[#[[1]]+#[[2]]]<2)& ]; ImplicitPlot3D[ lhs_ == rhs_, xr:{_,_,_},yr:{_,_,_},zr:{_,_,_},opts___ ] := ImplicitPlot3D[lhs - rhs,xr,yr,zr,opts] ImplicitPlot3D[ func_,{x_,xMin_,xMax_},{y_,yMin_,yMax_},{z_,zMin_,zMax_}, opts___ ] := Module[{cubes,cubePointer,f,boxes}, {cubes,cubePointer,divisions} = ReleaseHold[ Hold[ makeCubeList[PlotPoints,Passes] ]/. {opts} /. Options[ImplicitPlot3D] ]; f = If[Compiled /. {opts} /. Options[ImplicitPlot3D], Compile[{x,y,z},func], Function[{x,y,z},func] ]; boxes = {{{xMin,xMax},{yMin,yMax},{zMin,zMax}}}; If[Length[cubePointer] > 1, boxes = Fold[ prepassBoxes[ f,#1,#2,cubes,divisions,cubePointer ]&, boxes, Range[Length[divisions] - 1] ] ]; boxes = Flatten[ finalBoxes[ f,#,divisions[[-1]], cubes[[ cubePointer[[-1]] ]] ]& /@ boxes ,1]; Show[ Graphics3D[boxes],FilterOptions[Graphics3D,opts] ] ] (* Functions makeCubeList, standardizePlotPoints, checkErr, factorInteger2, adjustLength, and reduceList are used to take the options PlotPoints and Passes and put them in a standard form. *) makeCubeList[plotPoints_,passes_]:= Module[{plotPts,boxPtr,plotPtsReduced}, plotPts=standardizePlotPoints[plotPoints,passes]; If[plotPts===$Failed,Return[$Failed] ]; {plotPtsReduced,boxPtr} = reduceList[plotPts]; { Apply[ Flatten[Array[{#1-1,#2-1,#3-1}&,{#1,#2,#3}],2]&, plotPtsReduced, 1 ], boxPtr, plotPts } ] standardizePlotPoints[plotPts_,pss_] := checkErr /@ plotPts /; Depth[plotPts] == 3 standardizePlotPoints[plotPts_,pss_] := standardizePlotPoints[{plotPts},pss] /; Depth[plotPts] == 1 standardizePlotPoints[plotPts_,pss_] := Module[{}, Message[ImplicitPlot3D::ppta,plotPoints]; $Failed ] /; Depth[plotPts] != 2 ImplicitPlot3D::ppta = "`1` is not a list of numbers or a list of lists of numbers."; standardizePlotPoints[plotPts_,passes_] := Module[{plotPtsFactored,pss}, plotPtsFactored = factorInteger2 /@ checkErr[plotPts]; pss = Which[ passes===Automatic,Max[Length /@ plotPtsFactored ], passes > 0, Ceiling[passes], True,1 ]; Reverse[ Select[ Transpose[ adjustLength[#,pss]& /@ plotPtsFactored ], # != {1,1,1}& ] ] ] checkErr[{p_Integer}] := {p,p,p} /; p > 0 checkErr[{p_Integer,q_Integer,r_Integer}] := {p,q,r} /; p > 0 && q > 0 && r > 0 checkErr[p___] := Message[ImplicitPlot3D::pptb,p] ImplicitPlot3D::pptb = "`1` is not a list of either 1 or 3 positive integers."; factorInteger2[n_] := Flatten[ Table[#[[1]],{#[[2]]}]& /@ FactorInteger[n] ] adjustLength[factoredList_,n_] := factoredList /; Length[factoredList] == n adjustLength[factoredList_,n_] := Join[ Table[1,{n-Length[factoredList]} ], factoredList ] /; Length[factoredList] < n adjustLength[factoredList_,n_] := Nest[ Sort[ Rest[ ReplacePart[#,#[[1]]*#[[2]],2] ] ]&, factoredList, Length[factoredList] - n ] reduceList[x_List] := Module[{reducedx}, reducedx = Union[x]; {reducedx, Flatten[Position[reducedx,#]& /@ x]} ] (* Functions subtractAltered, sameFaceQ, orderVertices, and reorderedVertices are used to find the planar convex hull of the given points. *) subtractAltered[x_,y_] := x-y /; IntegerQ[x]&&IntegerQ[y] subtractAltered[x_,y_] := z SetAttributes[subtractAltered,Listable]; sameFaceQ[{v1_,e1_},{v2_,e2_}]:= Module[{z}, Apply[Times,subtractAltered[v1+z*e1,v2+z*e2] ] == 0 ] orderVertices[{ordered_,unordered_}] := { Join[ordered,unordered],{} } /; Length[unordered] < 2 orderVertices[{ordered_,unordered_}] := { Join[ordered,{unordered[[1]]} ], Rest[unordered] } /; sameFaceQ[ordered[[-1,2]],unordered[[1,2]] ] orderVertices[{ordered_,unordered_}] := {ordered,RotateLeft[unordered]} reorderedVertices[vertices_] := {} /; Length[vertices] < 2 reorderedVertices[vertices_] := Transpose[vertices][[1]] /; Length[vertices] < 4 reorderedVertices[vertices_] := Transpose[ First[ FixedPoint[ orderVertices, {{vertices[[1]]},Rest[vertices]}, Length[vertices](Length[vertices]-1)/2 ] ] ][[1]] (* Functions locusQ, prepassCubeQ, prepassLocusQ, prepassBoxes, ppBox, finalBoxes are used to determine if and where the surface intersects a given box. *) locusQ[h_,{i1_,j1_,k1_},{{i2_,j2_,k2_},{0,0,0}}] := Module[ {value = Abs[h[i1+i2,j1+j2,k1+k2]]}, NumberQ[value] && value == 0. ] locusQ[h_,{i1_,j1_,k1_},{{i2_,j2_,k2_},{i3_,j3_,k3_}}] := Module[ {value = h[i1+i2,j1+j2,k1+k2] h[i1+i2+i3,j1+j2+j3,k1+k2+k3] }, NumberQ[value] && value < 0. ] prepassCubeQ[g__] := Apply[Or,Map[prepassLocusQ[g,#]&,prepassEdges] ] prepassLocusQ[h_,{i_,j_,k_},{{i1_,j1_,k1_},{i2_,j2_,k2_}}] := Module[ {value = h[i+i1,j+j1,k+k1]*h[i+i1+i2,j+j1+j2,k+k1+k2] }, NumberQ[value] && value <= 0. ] prepassBoxes[func_,boxes_,n_,cb_,divis_,ptr_] := Flatten[ Map[ ppBox[func,#,n,cb,divis,ptr]&, boxes ], 1 ] ppBox[ func_,{{u1_,u2_},{v1_,v2_},{w1_,w2_}},n_,cb_,divis_,ptr_ ] := Module[{f,du,dv,dw,nu,nv,nw,cubes}, {nu,nv,nw} = divis[[n]]; cubes = cb[[ptr[[n]] ]]; du=(u2-u1)/nu; dv=(v2-v1)/nv; dw=(w2-w1)/nw; f[i_,j_,k_] := f[i,j,k] = func[i*du+u1,j*dv+v1,k*dw+w1]; Return[ Apply[{{#1*du+u1,(1+#1)*du+u1}, {#2*dv+v1,(1+#2)*dv+v1}, {#3*dw+w1,(1+#3)*dw+w1}}&, Select[cubes,prepassCubeQ[f,#]&],{1} ] ] ] finalBoxes[ func_,{{u1_,u2_},{v1_,v2_},{w1_,w2_}},{nu_,nv_,nw_},cb_ ] := Module[{f,du,dv,dw,coord,plate}, du=(u2-u1)/nu; dv=(v2-v1)/nv; dw=(w2-w1)/nw; f[i_,j_,k_] := f[i,j,k] = func[i*du+u1,j*dv+v1,k*dw+w1]; coord[{i1_,j1_,k1_},{{i2_,j2_,k2_},{0,0,0}}] := {u1+(i1+i2)*du,v1+(j1+j2)*dv,w1+(k1+k2)*dw}; coord[{i1_,j1_,k1_},{{i2_,j2_,k2_},{i3_,j3_,k3_}}] := Module[{f1=f[i1+i2,j1+j2,k1+k2],f2}, f2=f1/(f[i1+i2+i3,j1+j2+j3,k1+k2+k3]-f1); {u1 + du(i1+i2-i3*f2), v1 + dv(j1+j2-j3*f2), w1 + dw(k1+k2-k3*f2)} ]; plate[i_,j_,k_] := Map[{coord[{i,j,k},#],#}&, Select[edges,locusQ[f,{i,j,k},#]&] ]; Return[ Polygon /@ reorderedVertices /@ Select[Apply[plate,cb,{1}], (# != {})&] ] ] End[] (* "`Private`" *) EndPackage[] (* "Graphics`ImplicitPlot3D`" *)