(*:Name: Graphics`ContourPlot3D` *)

(* :Title: ContourPlot3D *)

(* :Author: Tom Wickham-Jones *)

(* :Summary:
This package introduces ContourPlot3D and ListContourPlot3D,
the three-dimensional analogs of the built-in functions
ContourPlot and ListContourPlot.
*)

(* :Context: Graphics`ContourPlot3D` *)

(* :Package Version: 1.0 *)

(* :Mathematica Version: 2.2 *)

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

(* :History:
	Created summer 1992 by Tom Wickham-Jones.
*)

(* :Keywords:
        Implicit, IsoContour.
*)


BeginPackage[ "Graphics`ContourPlot3D`", "Utilities`FilterOptions`"]


ContourPlot3D::usage = "ContourPlot3D[ fun, {x, xmin, xmax}, {y, ymin, ymax},
	{z, zmin, zmax}] plots the surface implicitly defined by
	fun[ x, y, z] == 0.  Setting Contours -> {val1, val2,...} will
	plot the surfaces for values val1, val2, etc."

ListContourPlot3D::usage = "ListContourPlot3D[ data, opts] takes a 
	three-dimensional data set interpreted as a representation
	of a function fun[ x, y, z], where the ranges of x, y, and z
	are set by the MeshRange option.   It then plots the surface 
	implicitly defined by fun[x, y, z] == 0.   Setting 
	Contours -> {val1, val2,...} will plot the surfaces for 
	values val1, val2, etc."

Contours::usage = 
	StringJoin[ Contours::usage,
	"   Contours is an option of ContourPlot3D and ListContourPlot3D
	which sets the value of the isocontour to be plotted."]

ContourStyle::usage = 
	StringJoin[ ContourStyle::usage,
	"   ContourStyle is an option of ContourPlot3D and ListContourPlot3D
	which specifies the styles of the isocontour to be drawn."]

MaxRecursion::usage = 
	StringJoin[ MaxRecursion::usage,
	"   MaxRecursion is also an option of ContourPlot3D which specifies
	the number of levels of recursion used in each cube."]


Begin[ "`Private`"]

PlotPoints::contpts = "Value of option PlotPoints -> `1` is not
	of the form x, {x, y}, {x, y, z} or {{x1, x2}, {y1, y2}, {z1, z2}}
	where each entry is a positive integer >= 2.  Using default PlotPoints."

ListContourPlot3D::gval = ContourPlot3D::gval = 
	"Function value `1` at grid point x = `2`, y = `3`,
	z = `4` is not a real number."

ContourPlot3D::bdval = 
ListContourPlot3D::bdval =
	"Value of option Contours -> `1` is not a list of real numbers.  Using default Contours."

ListContourPlot3D::badrnk = "The data is not a three-dimensional list of numbers."

Options[ ContourPlot3D] =
	Join[ {Contours -> {0.}, ContourStyle -> {},
	       MaxRecursion -> 1, PlotPoints -> {3,5},
	       Options[ ParametricPlot3D]}]

Options[ ListContourPlot3D] =
	Join[ {Contours -> 0., ContourStyle -> {},
	       MeshRange -> Automatic,
	       Options[ Graphics3D]}]

SetAttributes[ ContourPlot3D, HoldFirst]

ContourPlot3D[ f_, {u_, u0_, u1_},
                   {v_, v0_, v1_},
                   {w_, w0_, w1_}, opts___] :=
    Block[{vals, poly, rec, com, cstyle, un, vn, wn},
	{un,vals,cstyle,com,rec} = {PlotPoints, Contours, ContourStyle, Compiled, 
		MaxRecursion} /.{opts} /. Options[ContourPlot3D] ;
	{un, vn, wn} = CheckPnts[ un] ;
	
	vals = FixVals[ N[ vals], ContourPlot3D] ;
       	cstyle = FixStyle[ cstyle, Length[ vals]] ;
	poly = Map[ conbld[ f, {u, u0, u1, un},
                               {v, v0, v1, vn},
                               {w, w0, w1, wn},
                               rec, com, #]&, vals] ;
	poly = Flatten[ Transpose[ { cstyle, poly}]] ;
	Show[ Graphics3D[ poly, FilterOptions[ Graphics3D, opts]]]
    ]

conbld[ f_, {u_, u0_, u1_, un_},
            {v_, v0_, v1_, vn_},
            {w_, w0_, w1_, wn_}, recur_, comp_, val_] :=
    Block[{func},
         func = If[ comp, Compile[{u,v,w}, f-val], 
                          Function[{u,v,w}, N[ f-val]]];
         CubeRecur[ recur, func, {u0, u1, un},
                                 {v0, v1, vn},
                                 {w0, w1, wn}]
    ]


ListContourPlot3D[ data_, opts___] :=
    Block[{vals, cstyle, mesh, dims, ndata, nu, nv, nw, du, dv, dw},
	If[ TensorRank[ data] =!= 3, 
            Message[ ListContourPlot3D::badrnk];
            Return[ $Failed]] ;
	{vals, cstyle, mesh} = {Contours, ContourStyle, MeshRange} /.
		{opts} /. Options[ContourPlot3D] ;
	 
	cstyle = FixStyle[ cstyle, Length[ vals]] ;
	nu = Reverse[ Dimensions[ data]] ;       (* to follow ListPlot3D *)
	mesh = GetMesh[ mesh, nu] ;
	vals = FixVals[ N[ vals], ListContourPlot3D] ;
	{{u,u1}, {v, v1}, {w,w1}} = mesh;
	{nu, nv, nw} = nu-1 ; 
	ndata = N[ data] ;
	poly = Map[ 
                   listbld[ ndata-#, u, u1, nu, v, v1, nv, w, w1, nw]&,
                   vals] ;
	poly = Flatten[ Transpose[ { cstyle, poly}]] ;
	Show[ Graphics3D[ poly, FilterOptions[ Graphics3D, opts]]]
    ]


listbld[ data_, u_, u1_, nu_, v_, v1_, nv_, w_, w1_, nw_] :=
    Block[{du, dv, dw},
        du = (u1-u)/nu ;
        dv = (v1-v)/nv ;
        dw = (w1-w)/nw ;
	Table[ 
	      MakeCube[ data, u+(i-1) du,
                              v+(j-1) dv,
                              w+(k-1) dw,du,dv,dw,i,j,k],
              {i,nu}, {j,nv}, {k,nw}] 
    ]

FixVals[ vals_, head_] :=
    If[ VectorQ[ vals, (Head[ #] === Real)&],
	 	vals,
		Message[ head::bdval, vals];
                {0.}
      ]
    

FixStyle[ style_, n_] :=
    Block[{len = Length[ style], i},
	If[ ! ListQ[ style] || VectorQ[ style], 
	    Table[ style, {n}],
            Table[ Part[ style, Mod[ i, len]+1], {i,n}]]
        ]


MakeCube[ndata_, u_, v_, w_, du_, dv_, dw_, i_, j_, k_ ] :=
    Block[{cube, xpt, ypt, zpt},
	cube = Flatten[ 
                   Table[ 
                       Part[ ndata, k1, j1, i1],    (* to follow ListPlot3D *)
                           {i1, i, i+1},
                           {j1, j, j+1},
                           {k1, k, k+1}]] ;
	xpt = Transpose[
                 Flatten[
		   Table[ {i1, j1, k1}, 
			{i1, u, u+du, du},
			{j1, v, v+dv, dv},
			{k1, w, w+dw, dw}], 2]] ;
	{xpt, ypt, zpt} = xpt ;
	Which[
            !GoodArray[ cube], BadPoint[ {cube, xpt, ypt, zpt}];{},
            Apply[ SameQ, Map[ # >= 0. &, cube]], {},
            True, NCube1[]]
     ]



GetMesh[ Automatic, dims_] := Transpose[ {{1,1,1}, dims}]

GetMesh[ mesh_, dims_] :=
    If[ Length[ mesh] === 3 && MatrixQ[ mesh, NumberQ],
        mesh, Transpose[ {{1,1,1}, dims}]]


CheckPnts[ xn_ /; And[IntegerQ[ xn],xn>=2]] := {{xn, xn}, {xn, xn}, {xn, xn}}

CheckPnts[ {xn_, yn_} /; CheckNumber[ {xn, yn}]] :=
	{{xn, yn}, {xn, yn}, {xn, yn}}

CheckPnts[ {xn_, yn_, zn_} /; CheckNumber[ {xn, yn, zn}]] :=
	{{xn, xn}, {yn, yn}, {zn, zn}}

CheckPnts[ {xn_, yn_, zn_} /; And @@ (CheckNumber /@  {xn, yn, zn})] :=
	{xn, yn, zn}

CheckPnts[ x_] := (Message[ PlotPoints::contpts,x]; {{3,5},{3,5},{3,5}})


CheckNumber[ x_List] :=
	Apply[ And, Map[ And[ IntegerQ[#], # >= 2]&, x]]

AutoPoints[ u0_, u1_, n_] :=
        N[ (u1 - u0)/(n-1)]

CubeRecur[ rec_, func_, {u0_, u1_, un_},
                        {v0_, v1_, vn_},
                        {w0_, w1_, wn_}] :=
    Block[{x,y,z,du,dv,dw},
        du = AutoPoints[u0, u1, Part[ un,1]] ;
        dv = AutoPoints[v0, v1, Part[ vn,1]] ;
        dw = AutoPoints[w0, w1, Part[ wn,1]] ;
	Flatten[ Table[ NCube[ rec, func,
                               Part[un,2],
                               Part[vn,2], 
                               Part[wn,2],du,dv,dw,x,y,z], 
                                     {x, u0, u1-du, du},
                                     {y, v0, v1-dv, dv},
                                     {z, w0, w1-dw, dw}]]
    ]


NCube[ rec_,func_,un_,vn_,wn_,du_,dv_,dw_,u_,v_,w_] :=
    Block[{vals, x, y, z, cube, xpt, ypt, zpt},
	vals = 
            Transpose[
                Flatten[
                  Table[ {func[ x,y,z], x, y, z}, {x,u,u+du,du},
                                                  {y,v,v+dv,dv},
                                                  {z,w,w+dw,dw}],2]] ;
	cube = Part[ vals, 1];
	Which[
            !GoodArray[ cube], BadPoint[ vals];{},
            Apply[ SameQ, Map[ # >= 0. &, cube]], {},
            True, If[ rec < 1, 
                      xpt =  Part[ vals, 2];
                      ypt =  Part[ vals, 3];
                      zpt =  Part[ vals, 4];
                      NCube1[],
		      CubeRecur[ rec-1, 
                                 func, 
				 {Part[ vals,2,1], Part[ vals,2,8], {un,un}},
				 {Part[ vals,3,1], Part[ vals,3,8], {vn,vn}},
				 {Part[ vals,4,1], Part[ vals,4,8], {wn,wn}}]
     		  ]
	]
    ]

GoodArray[ x_List] :=
    VectorQ[x, MemberQ[{Real, Integer}, Head[#]]&]


GoodArray[_] := False

BadPoint[ vals_List] :=
    Block[{ cube, xpt, ypt, zpt, res},
        {cube, xpt, ypt, zpt} = vals ;
        res = Flatten[ Position[ 
	    Map[ MemberQ[{Real, Integer}, Head[#]]&, cube], False]] ;
	Map[ Message[ ContourPlot3D::gval, 
                      Part[ cube, #],
                      Part[ xpt, #],
                      Part[ ypt, #],
                      Part[ zpt, #]]&, res] ;
    ]

NCube1[ ] :=
    Block[ {res, cross, built},
        built = Table[ False, {12}] ;
        cross = Table[ CrossVert[ Part[ edgelist, i]], {i,12}] ;
	If[ Apply[ Or, cross],
	    Map[ PolyMake, Map[ Build1, trylist]], Print[bb];
	    {}
        ]
    ]

PolyMake[ g_] := Polygon[ g]

PolyMake[ {}] := {}

Build1[ pos_] :=
    Block[{res, nbuilt},
	If[ Part[ built, pos] || Not[ Part[ cross, pos]],
		{},
		nbuilt = Table[ False, {12}] ;
		res = Build[ pos, {}, {}];
		res]
    ]

Build[  pos_, prev_, face_] :=
    Block[{nface, next, res, npos,
           cnt = 1, flag = True},
	nface = Append[ face, Inter[ pos]] ;
	Part[ built, pos] = True ;
	Part[ nbuilt, pos] = True ;
	next = Part[ nextlist, pos] ;
	While[ flag,
	    npos = (Part[ next, cnt]+{-1,-1,-1}).{4,2,1}+1 ;
	    If[ npos =!= prev,
	        Which[ 
                    Part[ nbuilt, npos], 
                        res = nface ; 
                        flag = False,
	            Part[ cross, npos] && Not[ Part[ built, npos]], 
                        res = Build[ npos, pos, nface] ; 
                        flag = False]] ;
	    cnt++ ; 
	    If[ cnt > 6 && flag, res = {}; flag = False] ;
	] ;
	res
    ]
	
Inter[ pos_] :=
    Block[{p1, p2, f, v},
	{p1, p2} = EdgeToVertPos[ Part[ edgelist, pos]] ;
	v = Part[ cube, p1] ;
        f = v/ (v - Part[ cube, p2]) ;
        {
          v = Part[ xpt, p1]; v + f*(Part[ xpt, p2] - v),
          v = Part[ ypt, p1]; v + f*(Part[ ypt, p2] - v),
          v = Part[ zpt, p1]; v + f*(Part[ zpt, p2] - v)
        }
    ]

CrossVert[ edge_] :=
    Block[ {p1, p2},
	{p1, p2} = EdgeToVertPos[ edge] ;
	If[ Part[ cube, p1] >= 0 && Part[ cube, p2] < 0 ||
	    Part[ cube, p1] < 0 && Part[ cube, p2] >= 0,
	    True, False]
    ]


Nexts[ edge_] :=
    Block[{v1, v2, crd, rest},
	{v1, v2} = EdgeToVerts[ edge] ;
	crd = First[ edge] ;
	rest = Rest[ edge] ;
	Join[ {
              MakeEdge[ v1, crd+1],
              MakeEdge[ v1, crd+2],
              MakeEdge[ v2, crd+1],
              MakeEdge[ v2, crd+2]},
              If[ Apply[ SameQ, rest],
	           {Add1[ edge, 2], Add1[ edge, 3]},
	           {Add1[ edge, 3], Add1[ edge, 2]}]]
    ]

MakeEdge[ v_, n_] :=
    Block[{vn, crd, res},
	vn = Add1[ v, Mod[ n-1, 3] + 1] ;
	crd = Part[ Position[ Abs[v - vn], 1], 1, 1] ;
	Prepend[ Drop[ vn, {crd}], crd]
    ]

  
Add1[ l_, n_] :=
    Block[ {tmp},
	tmp = l ;
        Part[tmp, n] = Mod[ Part[ l, n], 2] + 1 ;
        tmp
    ]

EdgeToVerts[ edge_] :=
    If[ Part[ edge, 2] === Part[ edge, 3], 
        {Insert[ Rest[ edge], 1, First[ edge]],
         Insert[ Rest[ edge], 2, First[ edge]]},
        {Insert[ Rest[ edge], 2, First[ edge]],
         Insert[ Rest[ edge], 1, First[ edge]]}]


EdgeToVertPos[ edge_] := 
  {Insert[ Rest[ edge], 1, First[ edge]] + {-1,-1,-1},
   Insert[ Rest[ edge], 2, First[ edge]] + {-1, -1, -1}}.{4,2,1}+1


(* 
  These expresions should be calculated once and stored.
  They should be calculated last.
*)

edgelist = Flatten[ Table[ {i, j, k}, {i,3}, {j,2}, {k,2}], 2]

trylist = {4, 6, 8, 9,10,11,12}	(* don't need to try every edge *)

nextlist = Map[ Nexts, edgelist]

End[]

EndPackage[]

(*:Examples:

ContourPlot3D[Sqrt[ x^2 + y^2 + z^2], 
		{x, 0, 1}, {y, 0, 1},{z, 0, 1}, 
		Contours -> {.5}]


ContourPlot3D[ x y z, {x,-1,1}, {y,-1,1}, {z,-1,1}, Contours -> {.1}]

ContourPlot3D[ Sqrt[ x^2 + y^2 + z^2], {x,-2,2}, {y,-2,2}, {z,-2,2},
				Contours -> {1.}]

ContourPlot3D[ Sin[ Sqrt[ x^2 + y^2 + z^2]], 
		{x,-2,2}, {y,0,2}, {z,-2,2}, 
		Contours -> {.5}]


data = Table[ x^2 + 2*y^2 + 3*z^2,
                   {z, -1, 1,.25},
                   {y, -1, 1, .25},
                   {x, -1, 1, .25}];

ListContourPlot3D[ data, 
	MeshRange -> {{-1,1}, {-1,1}, {-1,1}},
        Contours -> {3.}]
*)

