8
$\begingroup$

I'm still pretty new to Mathematica, so I would like to seek advice regarding a geometrical problem.

I am currently trying to define that as an extra condition in the Mathematica code below. 

  reg = ImplicitRegion[x^2/a^2 + y^2/b^2 + z^2/c^2 <= 1 {z, y, x}];
  Volume[reg, Assumptions -> a > 0 && b > 0 && c > 0]

Any one has any idea how to incorporate it into the extra conditions in defining the implicit region?

| improve this question | |
$\endgroup$
  • $\begingroup$ What is the equation of the plane ....? $\endgroup$ – Dr. belisarius Oct 10 '14 at 4:53
  • $\begingroup$ In this case, the equation of the plane (XZ) is y = 0 right? I would need to rotate this plane about the point X and use the rotated plane as a condition in evaluating the volume. Not sure how do I do that though. $\endgroup$ – Corse Oct 10 '14 at 5:01
13
$\begingroup$

Let you have a vector ${\bf p}$, which is perpendicular to the plane and an ellipsoid with axes $(a,b,c)$. The illustration (2D for simplicity):

enter image description here

Mathematica can calculate the numeric value of the clipped volume easily

Nvolume[p_, abc_] := Volume[RegionIntersection[
   ImplicitRegion[{x, y, z}.N[p] > 0, {x, y, z}], 
   Ellipsoid[N[abc] {1, 0, 0}, N[abc]]]]

p = RandomReal[{-1, 1}, 3];
abc = RandomReal[{1, 2}, 3];

Nvolume[p, abc]
(* 16.2584 *)

Mathematica cannot derive the general formula, but it isn't difficult to derive manually. Let us introduce new coordinates

$$ x' = x/a, \quad y' = y/b, \quad z' = z/c. $$

In these coordinates the ellipsoid becomes the unit ball

enter image description here

The Jacobian of this transformation is $J=abc$. In the new coordinates the normalized perpendicular vector is

$$ {\bf n} = \frac{(ap_x,bp_y,cp_z)}{\sqrt{a^2p_x^2+b^2p_y^2+c^2p_z^2}}. $$

Now it is simple to integrate the volume along the axis $\xi$ because the cross section is a circle

$$ V=abc\int_{-n_x}^1\pi(1-\xi^2)d\xi = \pi abc \left(\frac{2}{3} + n_x - \frac{n_x^3}{3}\right) $$

volume[p_, abc_] := π Times @@ abc (2/3 + # - #^3/3) &@@ Normalize[abc p]

volume[p, abc]
(* 16.2584 *)

The result is the same.

Update: OP asks also about the area of the intersection. It is also an interesting question.

Mathematica region functionality is very powerful for numerical computations:

Narea[p_, abc_] := Area[RegionIntersection[ImplicitRegion[{x, y, z}.N[p] == 0, {x, y, z}], 
   Ellipsoid[N[abc] {1, 0, 0}, N[abc]]]]

Narea[p, abc]
(* 6.20243 *)

The analytic formula can be derived using the Dirac $\delta$-function
\begin{multline} A = \int_\text{ellipse} \delta \left({\bf r}\cdot\frac{{\bf p}}{p}\right)d{\bf r} = abcp \int_\text{unit ball} \delta \left(x'ap_x+y'bp_y+z'cp_z\right)d{\bf r}' = \\ \frac{abcp}{\sqrt{a^2p_x^2+b^2p_y^2+c^2p_z^2}}\int_\text{unit ball} \delta \left({\bf r}'\cdot{\bf n}\right)d{\bf r}'. \end{multline} It is the cross section of the unit ball. Hence \begin{equation} A = \frac{\pi abcp (1-n_x^2)}{\sqrt{a^2p_x^2+b^2p_y^2+c^2p_z^2}}. \end{equation}

area[p_, abc_] := π Times @@ abc (1 - #^2) & @@ Normalize[abc p] Norm[p]/Norm[abc p];

area[p, abc]
(* 6.20243 *)
| improve this answer | |
$\endgroup$
  • $\begingroup$ this is beautiful method and exposition $\endgroup$ – ubpdqn Oct 10 '14 at 23:33
  • $\begingroup$ Hi there @ybeltukov, thats a concise analytical method nicely complementing ubpdqn's proposed implementation! I understand that a plane takes {x,y,z}.N[p] = 0, but what does it mean when {x, y, z}.N[p] > 0 ? Also, how could I get a rotation angle of the plane w.r.t horizontal using the normal vector p? $\endgroup$ – Corse Oct 11 '14 at 6:32
  • $\begingroup$ @Corse, {x, y, z}.N[p] > 0 means the half of the space on the one side of the plane. As in ubpdqn's you can use p = {Sin[a], 0, Cos[a]}. $\endgroup$ – ybeltukov Oct 11 '14 at 12:51
  • $\begingroup$ @ybeltukov thanks for the clarification, i realise that the value of 'a' has to be in terms of a negative angle to get the volume under the rotated plane. Sorry to bother, could there be a simple modification to the code to get the corresponding area of the intersecting plane through the ellipsoid? $\endgroup$ – Corse Oct 12 '14 at 2:29
  • $\begingroup$ @Corse, you can change the sign of p (e.g. p = -{Sin[a], 0, Cos[a]}) or change > 0 to < 0 in Nvolume and signs of # and #^3 in volume $\endgroup$ – ybeltukov Oct 12 '14 at 2:34
8
$\begingroup$

Visualization of the @ybeltukov's answer

This is not an answer, just visualization of ybeltukov's answer.

Step 1 direction3D for direction handling.

opt1 = Sequence[Orange, Thick, Arrowheads[.15]];
opt2 = Sequence[PlotRange -> 1, Ticks -> None, Boxed -> False,
   ViewPoint -> {1, 1, 4}, ViewVertical -> {0, 1, 0}, 
   PlotRegion -> {{-.2, 1.05}, {-.2, 1.05}}, ImageSize -> 60];

direction3D[Dynamic[d_]] := DynamicModule[
  {$p = {0, π/2}, dd},
  LocatorPane[Dynamic[$p, ($p = #; dd @@ #) &],
   Graphics3D[
    {{opt1, Dynamic@Arrow[Tube[{{0, 0, 0}, d}, .04]]},
     {Dynamic@Cylinder[{{0, 0, 0}, d/100}]}}, opt2
    ], {{-π, π}, {π, 0}},Appearance -> None, ImageMargins -> 0
   ], Initialization :> (d = {1, 0, 0};
  dd[t_, s_] := (d = {Sin[t] Sin[s], Cos[s], Cos[t] Sin[s]}))
  ]

Step 2 ybeltukov's result.

volume[p_, abc_] := π Times @@ abc (2/3 + # - #^3/3) & @@Normalize[abc p]

step3 visualization

opts = Sequence[Boxed -> False, Axes -> True,
   ViewPoint -> {1, 1, 4}, ViewVertical -> {0, 1, 0},
   AxesOrigin -> {0, 0, 0}, PlotRange -> {{-3, 3}, {-2, 2}, {-2, 2}}];

Manipulate[
 elp = Graphics3D[{Opacity[.7], Orange, 
    Ellipsoid[{0, 0, 0}, {a, b, c}],
    Opacity[1], Black, 
    Text[Style[#, 16, Bold], #2] & @@@ {{"x", {3.1, 0, 0}}, {"y", {0,  2.2, 0}}, {"z", {0, 0, 2.2}}}
    }, opts];
 pln = ContourPlot3D[
   p.{x + a, y, z} == 0, {x, -5, 5}, {y, -5, 5}, {z, -5, 5},
   Mesh -> None,
   ContourStyle -> {Opacity[.5], Darker@Green}];
 Show[elp, pln, Graphics3D[{
    Text[Style[
      StringJoin["V= ", ToString[N@volume[p, {a, b, c}]/Pi], "\[Pi]"],
       20], {2, 2, 0}]}]],
 Pane[Row[{
    Pane[direction3D[Dynamic[p]], ImageMargins -> {{0, 40}, {0, 0}}], 
    Pane[Column[{Control@{{a, 2}, 1, 2}, Control@{b, 1, 2}, 
       Control@{c, 1, 2}}]]}]]]

enter image description here

| improve this answer | |
$\endgroup$
  • $\begingroup$ thats a nice visualization! thank you $\endgroup$ – Corse Oct 13 '14 at 1:44
6
$\begingroup$

Functions for generating ellipsoid $x^2/a^2+y^2/b^2+z^2/c^2=1$ and plane through point $\vec{p}$ with normal $\vec{n}$, i.e $\vec{n}\cdot(\vec{r}-\vec{p})=0$

el[a_, b_, c_, x_, y_, z_] := x^2/a^2 + y^2/b^2 + z^2/c^2
pl[n_, p_, x_, y_, z_] := z /. First@Solve[n.({x, y, z} - p) == 0, z]

Using as an example: a=1, b=2,c=3 and normal {1,0,par}:

Manipulate[
 Column[{Show[
    Plot3D[Evaluate[pl[{1, 0, par}, {-1, 0, 0}, x, y, z]], {x, -3, 
      3}, {y, -3, 3}, Mesh -> None, PlotStyle -> Opacity[0.5]], 
    RegionPlot3D[
     Evaluate[
      reg = (el[1, 2, 3, x, y, z] < 1 && 
         pl[{1, 0, par}, {-1, 0, 0}, x, y, z] > z)], {x, -3, 
      3}, {y, -3, 3}, {z, -3, 3}, BoxRatios -> {1, 1, 1}, 
     Mesh -> False, PlotStyle -> Red, PerformanceGoal -> "Quality"], 
    RegionPlot3D[
     Evaluate[el[1, 2, 3, x, y, z] < 0.9], {x, -3, 3}, {y, -3, 
      3}, {z, -3, 3}, PlotStyle -> Directive[Blue, Opacity[0.2]], 
     BoxRatios -> {1, 1, 1}, Mesh -> None], 
    PlotRange -> Table[{-3, 3}, {3}], ImageSize -> 400],
   StringForm["Volume of red region:`1`", 
    NumberForm[RegionMeasure[ImplicitRegion[reg, {x, y, z}]], 3]],
   StringForm["Volume of ellipsoid:`1`", NumberForm[N@4 Pi 1 2 3/3, 3]]
   }]
 , {par, 0.5, 4}]

enter image description here

UPDATE

In relation to comment:

rot[a_, p_, x_, y_, z_] := pl[{Sin[a], 0, Cos[a]}, p, x, y, z]

This interactive graphic shows relation of normal to plane and what I understand is the desired angle from horizontal plane:

Manipulate[
 Show[RegionPlot3D[
   Evaluate[el[1, 2, 3, x, y, z] < 1], {x, -4, 4}, {y, -4, 4}, {z, -4,
     4}, Mesh -> False, PlotStyle -> Opacity[0.3]], 
  Plot3D[rot[a Degree, {-1, 0, 0}, x, y, z], {x, -4, 4}, {y, -4, 4}, 
   Mesh -> False, PlotStyle -> {Blue, Opacity[0.5]}], 
  Graphics3D[{{Arrow[2 {{0, 0, 0}, {0, 0, 1}}]}, {Red, 
     Arrow[2 {{0, 0, 0}, {Sin[a Degree], 0, Cos[a Degree]}}]
     },
    {Arrow[{{-1, 0, 0}, {1, 0, 0}}]},
    {Arrow[{{-1, 0, 0}, {-1 + 2 Cos[a Degree], 0, -2 Sin[a Degree]}}]}
    }]], {a, 0, 90, Appearance -> "Labeled"}]

enter image description here

| improve this answer | |
$\endgroup$
  • $\begingroup$ hi there ubpdqn, thank you for the great implementation! Just wondering, how would I be able to obtain the volume of the red region in a generalized equation based on a specified rotation angle? I would still need to go through and fully understand your code there. I'm still not too clear as to how you managed to defined the plane (using vector product) and rotation axis. $\endgroup$ – Corse Oct 10 '14 at 8:43
  • $\begingroup$ @Corse see update in relation to angles $\endgroup$ – ubpdqn Oct 10 '14 at 9:16
  • $\begingroup$ you have addressed my exact problem!! but I have a couple of queries I hope you could advise: 1) What is the term 'par' supposed to represent and how do I convert it to the rotation angle 'a' ? 2) When I completely set the value of 'par' to the extreme right, which gives a value of par=4, the corresponding volume of the red region is 11 which is less than half of the volume of ellipsoid. How do I edit par to obtain the correct limits? (i.e. volume for a horizontal and a vertical plane) $\endgroup$ – Corse Oct 10 '14 at 11:28
  • $\begingroup$ Also, to edit the geometry of the ellipse, I would have to edit the values of a_,b_,c_ for each occurence in the el[] function right? $\endgroup$ – Corse Oct 10 '14 at 11:29
  • $\begingroup$ @Corse par just varies the normal vector that defines the plane (1,0,par). You could use the second or updated code rot to define the inclined plane and the relevant inequalities. Yes, I merely chose a=1,b=2,c=3 for illustrative purposes. You can choose what you like. Finally please also see my answer to your question re: rotation about arbitrary axis:mathematica.stackexchange.com/a/61782/1997 $\endgroup$ – ubpdqn Oct 10 '14 at 12:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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