Intersection of surface with parallel planes

Question

Consider the code (adapted from here)

h = x^2 + y^2/9 + z^2/4 - 1;
g = z;
ContourPlot3D[
   {h == 0, g == 0}, {x, -1, 1}, {y, -3, 3}, {z, -2, 2}, 
   MeshFunctions -> {Function[{x, y, z, f}, h - g]}, 
   MeshStyle -> {{Thick, Blue}}, Mesh -> {{0}}, 
   ContourStyle -> Directive[Orange, Opacity[0.5], Specularity[White, 30]]]

Now I have few ideas on the effect of MeshFunctions. Anyway, the result is very nice for me.

I'd like to do the same but with a parallel plane with $z=k$ for other values of $k$ (for example, $k=1$).

So I tried the code

h = x^2 + y^2/9 + z^2/4 - 1;
g = z;
k := 1;
ContourPlot3D[
   {h == 0, g == k}, {x, -1, 1}, {y, -3, 3}, {z, -2, 2}, 
   MeshFunctions -> {Function[{x, y, z, f}, h - g]}, 
   MeshStyle -> {{Thick, Blue}}, Mesh -> {{0}}, 
   ContourStyle -> Directive[Orange, Opacity[0.5], Specularity[White, 30]]]

and the result was

and (after reading some comments) I discovered that I can plot those two planes together simply using {h == 0, g == k, g == 0} obtaining

Question: How to get the sphere (to be true, an ellipsoid) together with two or three planes corresponding to different values of $k$ and their intersection (the blue curves) all on the same figure?

---

ps: I edited the post to help future readings and to follow the rich comments.

Answer 1

It appears that you are interested in showing only the intersections for an arbitrary set of cutting planes parallel to the xy-plane. That can be achieved by making some small modifications to PatoCriollo's answer. Like so:

h = x^2 + y^2/9 + z^2/4 - 1;
With[{cuts = Range[-5/2, 5/2, 1/2]}, 
  ContourPlot3D[h == 0, {x, -1, 1}, {y, -3, 3}, {z, -2, 2}, 
    MeshFunctions -> {Function[{x, y, z, f}, z]}, 
    MeshStyle -> {{Thick, Blue}}, Mesh -> {cuts}, 
    ContourStyle -> Directive[Opacity[0]]]]

Edit

On second thought, there is no need for g at all. The code above has been edited to eliminate g. This is much faster.

Answer 2

 h = x^2 + y^2/9 + z^2/4 - 1;
g = z;
ContourPlot3D[{h == 0, g == 0, g == k}, {x, -1, 1}, {y, -3, 
  3}, {z, -2, 2}, MeshFunctions -> {Function[{x, y, z, f}, z]}, 
 MeshStyle -> {{Thick, Blue}}, Mesh -> {{0, k }}, 
 ContourStyle -> 
  Directive[Orange, Opacity[0.5], Specularity[White, 30]]]

Answer 3

Using BoundaryStyle only:

You can use BoundaryStyle to mark the intersections of the contour surfaces:

h = x^2 + y^2/9 + z^2/4 - 1;
ContourPlot3D[{h == 0, z == 0}, {x, -1, 1}, {y, -3, 3}, {z, -2, 2}, 
 Mesh -> None, ContourStyle -> Directive[Orange, Opacity[0.5], Specularity[White, 30]], 
 BoundaryStyle -> {1 -> None, 2 -> None, {1, _} -> Directive[Thick, Blue]}]

This also works with multiple contour planes:

ContourPlot3D[Evaluate@Prepend[Thread[z == Range[-2, 2, .4]], h == 0],
 {x, -1,  1}, {y, -3, 3}, {z, -2, 2}, Mesh -> None,  ContourStyle -> None,
 BoundaryStyle -> Flatten[{Thread[Range[13] -> None], {1, _} -> Directive[Thick, Red]}]]

Or use

  ContourStyle -> Thread[Directive[Opacity[.2], RandomColor[13]]]
  (*  Hue /@ RandomReal[1, 13] instead of RandomColor[13] in version 9 *)

to get