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

Contour3DSubject::usage=
"Contour3DSubject[x,y,z] is the compiled function in the MathLink code."

ContourPlot3D::usage=
"ContourPlot3D[expr,{x,x1,x2},{y,y1,y1},{z,z1,z2}] draws the polygons of the function expr=0.\n
 ContourPlot3D[expr,{x,x1,x2},{y,y1,y1},{z,z1,z2},Contours->{c1,c2,..}] draws the polygons for expr=c1, expr=c2, ..."

Global`FindContour3D::usage=
"FindContour3D[level,x1,x2,y1,y2,z1,z2,nx,ny,nz,nrec] returns the polygons
 in the volume (x1,y1,z1), (x2,y2,z2) with Contour3DSubject[x,y,z]==level."
  
Global`ContourPolygons3D::usage=
"ContourPolygons3D[funcname,level,x1,x2,y1,y2,z1,z2,nx,ny,nz,nrec] returns
the polygons in the volume (x1,y1,z1), (x1,y2,z2) with\n
funcname[x,y,z]==level."

POVPolySurface::usage=
"POVPolySurface[polynom,options] is a new graphics 3d object. 
 Graphics3D[POVPolySurface[..]] will draw the 3d surface poly==0.
 The options for drawing are BoundingBox,VolumeSubdivisions, 
 VolumeRecursion, and PolySurfaceColor."

BoundingBox::usage="BoundingBox is a option for POVPolySurface, giving the volume for searching and drawing the surface."

VolumeSubdivisions::usage="VolumeSubdivisions is a option for POVPolysurface giving the nunber of subintervals in every space direction for drawing."

VolumeRecursion::usage="VolumeRecursion is a option for POVPolysurface giving the number of subintervals if a intersection is found."

PolySurfaceColor::usage="PolySurfaceColor is a option for POVPolysurface giving the surface color of the polynomial surface."

SphericalEigensystem::usage=
  "First argument of HydrogenDensity3D for eigensystem in spherical coordinates.";
ParaboloidalEigensystem::usage=
  "First argument of HydrogenDensity3D for eigensystem in paraboloidal coordinates.";  
NuclearCharge::usage=
  "NuclearCharge is a option for HydorgenDensity3D."

HydrogenDensity3D::usage=
"HydrogenDensity3D[SphericalEigensystem,n,l,m,{x1,x2},{y1,y2},{z1,z2},opts]
   finds the surface of constant density for the sperical eigenfunctions
   of hydrogen.\n
 HydrogenDensity3D[ParaboloidalEigensystem,n1,n2,m,{x1,x2},{y1,y2},{z1,z2},opts]
   finds the surface of constant density for the eigenfunctions in
   paraboloidal coordinates. The level of the density can be selected by the
   Contours option of the function."

Global`ShowCharges::usage=
"ShowCharges is a option for ElectroStaticPotential3D. It can be True, False, 
or a number. If ShowCharges is True the position of the charges are drawn as 
points, if it is a number the charges are drawn as points with PointSize of 
the value of ShowCharges.";

ElectroStaticPotential3D::usage=
"ElectroStaticPotential3D[{q1,{rx1,ry1,rz1},..},{x1,x2},{y1,y2},{z1,z2},opts] searches
for the equipotential surfaces of the potential produced by 
the point charges q1, ... at the space positions r={rx1,ry1,rz1} in the
volume x in [x1,x2], y in [y1,y2] and z in [z1,z2]. The potential values
for the surfaces can be given by the Contours option. The surface style
can be given in the ContourStlye option. The volume subdivisions as well 
as the optional recursion can be given by the PlotPoints option."

Begin["`Private`"]

ContourPlot3D::nolink=
  "No MathLink supported on this system, uses Graphics`ContourPlot3D`."

If[$LinkSupported, (* then *)
  lnk=CallLinkProg["MLBin","contour3"],
(* else *)
  lnk=$Failed;
  Message[ContourPlot3D::nolink]
 ]

PlotPoints::pltp=
"PlotPoints->`1` must be a list of two, three or four integers.";

convertPlotP[a_]:=
  Module[{nx,ny,nz,rec},
      Switch[a,
       (*case*) {_Integer,_Integer},
           (*:*) nx=ny=nz=First[a]; rec=Last[a],
       (*case*) {_Integer,_Integer,_Integer}, 
           (*:*) {nx,ny,nz}=a; rec=0,
       (*case*) {_Integer,_Integer,_Integer,_Integer}, 
           (*:*) {nx,ny,nz,rec}=a,
       _, Message[PlotPoints::pltp,a];
          Return[$Failed]
     ];
     {nx,ny,nz,rec}
    ]

Options[ContourPlot3D]:=
  {
    Compiled->True,
    PlotPoints->{5,3},
    Contours->{0.0},
    ContourStyle->{}
  };

ContourPlot3D::nointpnt=
"PlotPoints must be an integer or a list of one or two integers."
ContourPlot3D::norealval=
"Argument `1` does not evaluate to real numbers at the positions `2`."
ContourPlot3D::norlevel=
"Contours must be a list of numbers."

ContourPlot3D[expr_,
  {x_Symbol,x1_?NumberQ,x2_?NumberQ},
  {y_Symbol,y1_?NumberQ,y2_?NumberQ},
  {z_Symbol,z1_?NumberQ,z2_?NumberQ},opt___]:=
  Block[{fun,comp,ppoints,nx,ny,nz,r=0,box,
         levels,lstyle,poly,gropt,iterlimit,tvec},
    comp = Compiled /.{opt} /. Options[ContourPlot3D];
    ppoints= PlotPoints /. {opt} /. Options[ContourPlot3D];
    ppoints=convertPlotP[ppoints];
    If[ppoints!=$Failed || TrueQ[Length[ppoints]==4],
      {nx,ny,nz,r}=ppoints,     
      Message[ContourPlot3D::nointpnt];
      Return[$Failed]
     ];
    levels=Contours /. {opt} /. Options[ContourPlot3D];
    levels=N[levels] /. 0->0.0;
    If[!VectorQ[levels,NumberQ],
      Message[ContourPlot3D::norlevel];
      Return[$Failed]
     ];
    lstyle=ContourStyle /. {opt} /. Options[ContourPlot3D];
    If[lstyle=!={}, (* then *)
      lstyle=Map[Flatten[List[#]] &,lstyle];
      (* bring contour style to the length of levels *)
      While[Length[lstyle]<Length[levels],
        lstyle=Join[lstyle,lstyle];
       ];
      lstyle=Take[lstyle,Length[levels]];
     ];
    If[!MatchQ[expr,_[_Symbol,_Symbol,_Symbol]],(* then, a user expression *)
      If[comp, (* then *)
        fun=Compile[{{x,_Real},{y,_Real},{z,_Real}},expr],
      (* else *)
        fun=Function[{x,y,z},expr]
       ],
      (* else, a defined function *)
       fun=Head[expr]
     ];
    (* test, if the function evaluates to real numbers *)
    tvec=Map[
           Apply[fun,#] &,
           N[{
              {x1,y1,z1},
              0.5*{x1+x2,y1+y2,z1+z2},
              {x2,y2,z2}} /. 0->0.0
            ]
          ];
    If[!VectorQ[tvec,(Head[#]==Real) &],
      Message[ContourPlot3D::norealval,
              expr,{{x1,y1,z1},{x1+x2,y1+y2,z1+z2}/2,{x2,y2,z2}}];
      Return[$Failed]
     ];
    iterlimit=$IterationLimit;
    $IterationLimit=Infinity;
    box=N[{x1,x2,y1,y2,z1,z2}] /. 0-> 0.0;
    If[fun =!=Contour3DSubject,(* then, Mathematica function *)
      poly=Global`ContourPolygons3D[fun,#, Sequence @@ box,
                    nx,ny,nz,r] & /@ levels,
      (* else, the internal function of the MathLink program *)
      poly=Global`FindContour3D[#,Sequence @@ box,
                    nx,ny,nz,r] & /@ levels
     ];
    If[poly===$Failed, poly];
    $IterationLimit=iterlimit;
    If[lstyle=!={} && Length[levels]==Length[lstyle],
      poly=MapThread[Join,{lstyle,poly},1]
     ];
    Show[
      Graphics3D[
        poly,
        FilterOptions[Graphics3D,opt]
       ]
     ]    
   ] 

(* Polynom surface objects *)

(* recreate the polynom *)
POVCoeffs[order_Integer, {u_Symbol,v_Symbol,w_Symbol}]:=
  Module[{i,j,k},
    Reverse[
      Flatten[
        Table[
          If[i+j+k<=order,
            u^i*v^j*w^k],
          {i,0,order},
          {j,0,order},
          {k,0,order}
         ]
       ]
     ] //Select [#, (# =!=Null) &] &
    ]

Options[POVPolySurface]:=
  {BoundingBox->Automatic,
   VolumeSubdivisions->8,
   VolumeRecursion->2,
   PolySurfaceColor->SurfaceColor[RGBColor[1,1,1]]
  }
   
Format[POVPolySurface[c_List,n_Integer,var_List,opt___],OutputForm]:=
  SequenceForm["POVPolySurface[",c.POVCoeffs[n,var]," == 0]",var]

If[$VersionNumber>=3,
  Format[POVPolySurface[c_List,n_Integer,var_List,opt___],StandardForm]:=
    SequenceForm["POVPolySurface[",c.POVCoeffs[n,var]," == 0][",var,"]"]
 ]

Clear[mselect]

mselect[poly_?NumericQ,_]:=poly
mselect[poly_Plus,v_Symbol] := Select[poly,FreeQ[#,v] &]
mselect[prod_Times,v_Symbol]:=0
mselcct[a_,{}]:=a
mselect[poly_,vs_List]:=mselect[mselect[poly,First[vs]] , Rest[vs]]

numCoeff[poly_,what_,vars_]:=
  Module[{c,cc},
    cc=c=Coefficient[poly,what];
    c=mselect[c,vars];       
    If[NumericQ[c], N[c], 0]
   ]

POVPolySurface[poly_?PolynomialQ,opt___Rule]:=
  Module[{var,ord,any,coeff,c},
    var=Variables[poly];
    ord=Max[
          Append[
            Exponent[poly,#] & /@ var,
            Exponent[poly //. (#-> any & /@ var),any]
           ]
         ];
    c=POVCoeffs[ord,var];
    coeff=numCoeff[poly,#,var] & /@ Drop[c,-1];
    AppendTo[coeff,poly /. Thread[var ->Table[0,{Length[var]}]]];
    POVPolySurface[coeff,ord,var,{FilterOptions[POVPolySurface,opt]}]
   ]

POVPolySurface::ibound=
"BoundingBox must be a list of three pairs of real numbers.";
POVPolySurface::isdiv=
"VolumeSubdivisions must be an integer or a list of three integers.";
POVPolySurface::irec=
"VolumeRecursion must be an integer.";

POVPolySurface /:
  Graphics3D[POVPolySurface[coeff_,ord_,_List,opt_],nopt___]:=
    Module[{poly,bbox,sdiv,rec,scolor},
       bbox=BoundingBox /. opt /. Options[POVPolySurface];
       If[Automatic==bbox,
         bbox = PlotRange /. {nopt} /. Options[Graphics3D];
         If[Automatic == bbox, bbox={{-1,1},{-1,1},{-1,1}}];
        ];
       bbox=N[Flatten[bbox]] /. 0 ->0.0;
       If[Length[bbox] =!= 6,
         Message[POVPolySurface::ibound];
         Return[Graphics3D[{},nopt]]
        ];
       sdiv=VolumeSubdivisions /. opt /. Options[POVPolySurface];
       If[IntegerQ[sdiv],sdiv=Table[sdiv,{3}]];
       If[!VectorQ[sdiv,IntegerQ] && Length[sdiv]!=3,
         Message[POVPolySurface::isdiv];
         Return[Graphics3D[{},nopt]]
        ];       
       rec=VolumeRecursion /. opt /. Options[POVPolySurface];
       If[!IntegerQ[rec],
         Message[POVPolySurface::irec];
         Return[Graphics3D[{},nopt]]
        ];
       poly=Global`PolynomSurface3D[
              N[coeff] /. 0-> 0.0,ord,
              Sequence @@ bbox,Sequence @@ sdiv,rec];
       If[poly==$Failed,
         Return[poly]
        ];
       scolor=PolySurfaceColor /. opt /. Options[POVPolySurface];
       Graphics3D[{scolor,poly},FilterOptions[Graphics3D,nopt]]
      ]

POVPolySurface /: Display[fname_String,surf_POVPolySurface,opt___]:=
  Display[fname,Graphics3D[surf],opt]

(* Hydrogen densities *)

hydroCoords={SphericalEigensystem,ParaboloidalEigensystem};  

Options[HydrogenDensity3D]:=
  {PlotPoints->{5,2},
   Contours->{0.01},
   NuclearCharge->1.0
  }   

HydrogenDensity3D[
  coordsys_ /; MemberQ[hydroCoords,coordsys],
  nn1_Integer,ln2_Integer,m_Integer,
  {x1_?NumberQ,x2_?NumberQ},
  {y1_?NumberQ,y2_?NumberQ},
  {z1_?NumberQ,z2_?NumberQ},opt___]:=
  Module[{pltpnt,cntr,parab,sdivs,i,rec,zeff,poly},
    pltpnt=PlotPoints /. {opt}  /.  Options[HydrogenDensity3D];
    cntr=N[Contours /. {opt}  /.  Options[HydrogenDensity3D]];
    zeff=N[NuclearCharge /. {opt} /. Options[HydrogenDensity3D]];
    If[SphericalEigensystem===coordsys,
      parab=0,
      parab=1
     ];
    pltpnt=convertPlotP[pltpnt];
    If[VectorQ[pltpnt,IntegerQ] && TrueQ[Length[pltpnt]==4],
      sdivs=Drop[pltpnt,-1];
      rec=Last[pltpnt],
      sdivs={8,8,8},
      rec=3
     ];
    If[!NumberQ[zeff], zeff=1.0];
    poly=Global`HydrogenDensitySurface3D[
          nn1,ln2,m,zeff,parab,#,
          Sequence @@ (N[{x1,x2,y1,y2,z1,z2}] /. 0->0.0),
          Sequence @@ sdivs,rec] & /@ cntr;
    If[FreeQ[poly,$Failed],
      Show[Graphics3D[poly],FilterOptions[Graphics3D,opt]],
      $Failed
     ]       
   ]

(* Electro static potential of point charges *)

Options[ElectroStaticPotential3D]:=
  {PlotPoints->{5,3},
   Contours->{0.0},
   ContourStyle->{},
   Global`ShowCharges->False
  }  

ElectroStaticPotential3D::boxval=
"The bounding box for surface search must evaluate to numbers.";
ElectroStaticPotential3D::cntval=
"The Contours->`1` must evaluate to numbers.";

ElectroStaticPotential3D[
  chargepos:{{_?NumberQ,{_?NumberQ,_?NumberQ,_?NumberQ}}..},
  {x1_,x2_},{y1_,y2_},{z1_,z2_},opts___]:=
  Module[{pltpoint,cntr,cntrstyle,showch,poly,nx,ny,nz,rec,box,pnts,c},
    pltpoint=PlotPoints /. {opts} /. Options[ElectroStaticPotential3D];
    cntr=Contours /. {opts} /. Options[ElectroStaticPotential3D];
    cntrstyle=ContourStyle /. {opts} /. Options[ElectroStaticPotential3D];
    showch=Global`ShowCharges /. {opts} /. Options[ElectroStaticPotential3D];
    If[!VectorQ[pltpoint,IntegerQ],
      Message[PlotPoints::pltp,pltpoint];
       Return[$Failed]
     ];
    pltpoint=convertPlotP[pltpoint];
    If[VectorQ[pltpoint,IntegerQ] || TrueQ[Length[pltpoint]==4],
      {nx,ny,nz,rec}=pltpoint,
      Return[$Failed]
     ];
    box=N[{x1,x2,y1,y2,z1,z2}] /. 0->0.0;
    If[!VectorQ[box,Head[#]==Real &],
      Message[ElectroStaticPotential3D::boxval];
      Return[$Failed]
     ];
    If[!VectorQ[cntr,NumberQ],
      Message[ElectroStaticPotential3D::cntval,cntr];
      Return[$Failed]
     ]; 
    If[cntrstyle!={},
      While[Length[cntrstyle]<Length[cntr],
        cntrstyle=Join[cntrstyle,cntrstyle]
       ];
      cntrstyle=Take[cntrstyle,Length[cntr]]
     ]; 
    cntr=N[cntr] /. 0->0.0;
    poly=Global`EStaticPotentialSurface3D[
           #,Sequence @@ box,nx,ny,nz,rec,N[chargepos] /. 0->0.0] & /@ cntr;
    If[TrueQ[showch] ||NumberQ[showch],
      pnts= Point[Last[#]] & /@ chargepos;
      If[NumberQ[showch],
        PrependTo[pnts,PointSize[showch]]
       ],
      pnts={}
     ];
    If[cntrstyle!={},
      poly=Map[Apply[Join,#] &,Transpose[{List /@ cntrstyle,poly}]]
     ]; 
    gr=Join[pnts,poly];
    Show[
      Graphics3D[
        gr,FilterOptions[Graphics3D,opts]
       ]
     ]
   ]
End[]
EndPackage[]
Null