# Finding volume of a segment

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?

• What is the equation of the plane ....? Oct 10, 2014 at 4:53
• 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. Oct 10, 2014 at 5:01

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):

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

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 $$A = \frac{\pi abcp (1-n_x^2)}{\sqrt{a^2p_x^2+b^2p_y^2+c^2p_z^2}}.$$

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

area[p, abc]
(* 6.20243 *)

• this is beautiful method and exposition Oct 10, 2014 at 23:33
• 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? Oct 11, 2014 at 6:32
• @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]}. Oct 11, 2014 at 12:51
• @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? Oct 12, 2014 at 2:29
• @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 Oct 12, 2014 at 2:34

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}}]]}]]]  • thats a nice visualization! thank you Oct 13, 2014 at 1:44 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}]


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"}]


• 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. Oct 10, 2014 at 8:43
• @Corse see update in relation to angles Oct 10, 2014 at 9:16
• 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) Oct 10, 2014 at 11:28
• 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? Oct 10, 2014 at 11:29
• @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 Oct 10, 2014 at 12:09