I am having difficulty getting the level sets of a 3D parametric plot. I have a parametrized 3D surface that I suspect is self-intersecting and one way to see this visually is by plotting the level sets (intersection with planes say z=-4,-2,0,4.. etc) and these isocurves will show me visually that there is self-intersection. This is the code I use (apologies for the long equations!):

wspace[t2_, t3_] := 
 Module[{l1 = 1/2, l2 = 1/3, d2 = 2, a2 = 3, v2, v3}, v2 = Tan[t2/2];
  v3 = Tan[t3/2];
  R = (-l1*v2 + l2*v2 - (l1 + l2)*v3)^2 + (v2 + 
       l1*l2*v2 + (1 - l1*l2) v3)^2 + (l1 + 
       l2 + (-l1*v2 + l2*v2) v3)^2 + (1 - 
       l1*l2 - (v2 + l1*l2*v2) v3)^2;
  {(1 + l1^2) (1 + l2^2) (-1 - a2 - v2^2 + a2*v2^2) (1 + v3^2), 
    2*a2*l1^2*l2^2*v2*v3^2 + 2*a2*l1^2*l2^2*v2 - 
     a2*l1^2*l2*v2^2*v3^2 + 2*a2*l1^2*l2*v2*v3 - a2*l1^2*l2 + 
     a2*l1^2*v2^2*v3 + a2*l1^2*v2*v3^2 + a2*l1^2*v2 + a2*l1^2*v3 + 
     a2*l1*l2^2*v2^2*v3^2 - a2*l1*l2^2 - 2*a2*l1*l2*v2^2*v3 + 
     2*a2*l1*l2*v3 - a2*l1*v2^2*v3^2 + a2*l1 + a2*l2^2*v2^2*v3 - 
     a2*l2^2*v2*v3^2 - a2*l2^2*v2 + a2*l2^2*v3 + a2*l2*v2^2*v3^2 + 
     2*a2*l2*v2*v3 + a2*l2 - 2*a2*v2*v3^2 - 2*a2*v2 + 
     d2*l1^2*l2*v2^2*v3^2 - d2*l1^2*l2 + d2*l1^2*v2^2*v3 + 
     d2*l1^2*v2*v3^2 - d2*l1^2*v2 - d2*l1^2*v3 + 
     d2*l1*l2^2*v2^2*v3^2 + 2*d2*l1*l2^2*v2^2 - 2*d2*l1*l2^2*v2*v3 + 
     2*d2*l1*l2^2*v3^2 + d2*l1*l2^2 - 2*d2*l1*l2*v2*v3^2 - 
     2*d2*l1*l2*v2 + d2*l1*v2^2*v3^2 + 2*d2*l1*v2^2 + 2*d2*l1*v2*v3 + 
     2*d2*l1*v3^2 + d2*l1 - d2*l2^2*v2^2*v3 + d2*l2^2*v2*v3^2 - 
     d2*l2^2*v2 + d2*l2^2*v3 + d2*l2*v2^2*v3^2 - d2*l2 + 
     l1^2*l2*v2^2*v3^2 - l1^2*l2 - l1^2*v2^2*v3 - l1^2*v2*v3^2 + 
     l1^2*v2 + l1^2*v3 - l1*l2^2*v2^2*v3^2 - 2*l1*l2^2*v2*v3 - 
     l1*l2^2 + 2*l1*l2*v2^2*v3 + 2*l1*l2*v3 + l1*v2^2*v3^2 - 
     2*l1*v2*v3 + l1 - l2^2*v2^2*v3 + l2^2*v2*v3^2 - l2^2*v2 + 
     l2^2*v3 - l2*v2^2*v3^2 + l2, 
    2*a2*l1^2*l2*v2^2*v3 + 2*a2*l1^2*l2*v2*v3^2 - 2*a2*l1^2*l2*v2 - 
     2*a2*l1^2*l2*v3 - 2*a2*l1*l2^2*v2^2*v3 - 2*a2*l1*l2^2*v2*v3^2 - 
     2*a2*l1*l2^2*v2 - 2*a2*l1*l2^2*v3 + 2*a2*l1*v2^2*v3 - 
     2*a2*l1*v2*v3^2 - 2*a2*l1*v2 + 2*a2*l1*v3 - 2*a2*l2*v2^2*v3 + 
     2*a2*l2*v2*v3^2 - 2*a2*l2*v2 + 2*a2*l2*v3 + 
     d2*l1^2*l2^2*v2^2*v3^2 + d2*l1^2*l2^2*v2^2 + d2*l1^2*l2^2*v3^2 + 
     d2*l1^2*l2^2 + d2*l1^2*v2^2*v3^2 - d2*l1^2*v2^2 - 
     4*d2*l1^2*v2*v3 - d2*l1^2*v3^2 + d2*l1^2 - d2*l2^2*v2^2*v3^2 + 
     d2*l2^2*v2^2 - 4*d2*l2^2*v2*v3 + d2*l2^2*v3^2 - d2*l2^2 - 
     d2*v2^2*v3^2 - d2*v2^2 - d2*v3^2 - d2 - 2*l1^2*l2*v2^2*v3 - 
     2*l1^2*l2*v2*v3^2 - 2*l1^2*l2*v2 - 2*l1^2*l2*v3 + 
     2*l1*l2^2*v2^2*v3 - 2*l1*l2^2*v2*v3^2 + 2*l1*l2^2*v2 - 
     2*l1*l2^2*v3 - 2*l1*v2^2*v3 - 2*l1*v2*v3^2 + 2*l1*v2 + 2*l1*v3 + 
     2*l2*v2^2*v3 - 2*l2*v2*v3^2 - 2*l2*v2 + 2*l2*v3}/R]

zfcn = Function[{t2,t3},
  l1=1/2; l2=1/3; d2=2; a2=3;
  v2=Tan[t2/2];
  v3=Tan[t3/2];
  R=(-l1*v2+l2*v2-(l1+l2)*v3)^2+(v2+l1*l2*v2+(1-l1*l2)v3)^2+(l1+l2+(-l1*v2+l2*v2)v3)^2+(1-l1*l2-(v2+l1*l2*v2)v3)^2;
  (2*a2*l1^2*l2*v2^2*v3+2*a2*l1^2*l2*v2*v3^2-2*a2*l1^2*l2*v2-2*a2*l1^2*l2*v3-2*a2*l1*l2^2*v2^2*v3-2*a2*l1*l2^2*v2*v3^2-2*a2*l1*l2^2*v2-2*a2*l1*l2^2*v3+2*a2*l1*v2^2*v3-2*a2*l1*v2*v3^2-2*a2*l1*v2+2*a2*l1*v3-2*a2*l2*v2^2*v3+2*a2*l2*v2*v3^2-2*a2*l2*v2+2*a2*l2*v3+d2*l1^2*l2^2*v2^2*v3^2+d2*l1^2*l2^2*v2^2+d2*l1^2*l2^2*v3^2+d2*l1^2*l2^2+d2*l1^2*v2^2*v3^2-d2*l1^2*v2^2-4*d2*l1^2*v2*v3-d2*l1^2*v3^2+d2*l1^2-d2*l2^2*v2^2*v3^2+d2*l2^2*v2^2-4*d2*l2^2*v2*v3+d2*l2^2*v3^2-d2*l2^2-d2*v2^2*v3^2-d2*v2^2-d2*v3^2-d2-2*l1^2*l2*v2^2*v3-2*l1^2*l2*v2*v3^2-2*l1^2*l2*v2-2*l1^2*l2*v3+2*l1*l2^2*v2^2*v3-2*l1*l2^2*v2*v3^2+2*l1*l2^2*v2-2*l1*l2^2*v3-2*l1*v2^2*v3-2*l1*v2*v3^2+2*l1*v2+2*l1*v3+2*l2*v2^2*v3-2*l2*v2*v3^2-2*l2*v2+2*l2*v3)/R
];

ParametricPlot3D[wspace[t2,t3],{t2,-Pi,Pi},{t3,-Pi,Pi},
  MeshFunctions->Function[{x,y,t2,t3},zfcn[t2,t3]], Mesh->{{-4,-2,0,2}}, MeshStyle -> {{Thick, Blue}},
  PlotPoints->60, BoxRatios->{1,1,1},AxesLabel->{x,y,z},
  PlotRange->{{-4,2},{-2,5},{-4,3}},PlotStyle->Opacity[0.5]]

This yields the image below , but the curve does not look like an intersection with a plane parallel to the xy-plane. Am I missing something? I tried simply using Mesh->{{-2}} and I still don't get a curve that looks like an intersection with a plane.

Attempt to plot level set contours

up vote 5 down vote accepted

Your problem comes from MeshFunction. Why would your nonlinear MeshFunction produce planar curves? Try instead MeshFunctions -> {#3 &}:

ParametricPlot3D[wspace[t2,t3],{t2,-Pi,Pi},{t3,-Pi,Pi},
  MeshFunctions->{#3 &}, Mesh->{{-4,-2,0,2}}, MeshStyle -> {{Thick, Blue}},
  PlotPoints->60, BoxRatios->{1,1,1},AxesLabel->{x,y,z},
  PlotRange->{{-4,2},{-2,5},{-4,3}},PlotStyle->Opacity[0.5]]

enter image description here

  • Wonderful. Thanks. I clearly did not know how to use MeshFunctions. Thanks for illustrating it for me. – quantum May 16 at 15:12

You can also use the option Exclusions:

ParametricPlot3D[wspace[t2, t3], {t2, -Pi, Pi}, {t3, -Pi, Pi}, 
  PlotStyle -> Opacity[.5],  PlotPoints -> 60, Mesh -> None, BoundaryStyle -> None, 
  BoxRatios -> 1, AxesLabel -> {x, y, z}, 
  PlotRange -> {{-4, 2}, {-2, 5}, {-4, 3}},  
  Exclusions -> {zfcn[t2, t3] == -2, zfcn[t2, t3] == 0, zfcn[t2, t3] == 2}, 
  ExclusionsStyle -> Directive[Thick, Blue]]

enter image description here

  • Any recommendation on whether to use Mesh or Exclusions? – anderstood May 17 at 14:53
  • @anderstood, i personally prefer Mesh (because i am more familiar with it). – kglr May 17 at 16:55

Try this:

Manipulate[
 Show[{
   ParametricPlot3D[wspace[t2, t3], {t2, -Pi, Pi}, {t3, -Pi, Pi}, 
    PlotPoints -> ControlActive[10, 60], BoxRatios -> {1, 1, 1}, 
    AxesLabel -> {x, y, z}, PlotRange -> {{-4, 2}, {-2, 5}, {-4, 3}}, 
    PlotStyle -> Opacity[0.5]],
   Graphics3D[{Blue, Opacity[0.8], 
     InfinitePlane[{{0, 1, z}, {1, 0, z}, {0, 0, z}}]}]
   }], {z, -4, 4, Appearance -> "Labeled"}]

enter image description here

Hope it helps. Have fun!

Your Answer

 
discard

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.