How to calculate the mass of a inhomogenous sphere?

I want to calculated the mass of a sphere based on a threedimensional discret inhomogenous density distribution. Lets say a set of 3x3x3 cubes of different densities is inscribed by a sphere. What is the fastest way to sum up the partitioned masses using Python?

I tried to calculate the volume under the mathematical equation for a sphere: x^2 +y^2 +z^2 = R^2 for the range of one of the cubes using scipy.integrate.dblquad. However, the result is only valid if the boundaries are smaller than the radius of the sphere and repetitive calculation for lets say 50'000 spheres with 27 cubes each w开发者_高级运维ould be quite slow.

On the other hand, the usual equation for CoM caluations could't be used in my opinion, due to the rather coarse and discrete mass distribution.

---

Timing Experiment

You didn't specify your timing constraints, so I've done a little experiment with a nice integration package.

Without optimization, each integral in spherical coordinates can be evaluated in 0.005 secs in a standard laptop if the cubes densities are straightforward functions.

Just as a reference, this is the program in Mathematica:

Clear@f;
(* Define a cuboid as density function *)
iP = IntegerPart;
f[{x_, y_, z_}, {lx_, ly_, lz_}] :=   iP[x - lx] + iP[y - ly] + iP[z - lz] /; 
   lx <= x <= lx + 3 && ly <= y <= ly + 3 && lz <= z <= lz + 3;

f[{x_, y_, z_}, {lx_, ly_, lz_}] := Break[] /; True;

Timing[Table[s = RandomReal[{0, 3}, 3]; (*sphere center random*)
   sphereRadius = Min[Union[s, 3 - s]]; (*max radius inside cuboid *)
   NIntegrate[(f[{x, y, z} - s, -s] /.  (*integrate in spherical coords *)
       {x -> r Cos@th Sin@phi, 
        y -> r Sin@th Sin@phi, 
        z -> r Cos@phi}) r^2 Sin@phi,
       {r, 0, sphereRadius}, {th, 0, 2 Pi}, {phi, 0, Pi}], 
         {10000}]][[1]]  

The result is 52 secs for 10^4 iterations.

So perhaps you don't need to optimize a lot ...

---

I cannot get your exact meaning of inscribed by a sphere. Also I havent tried the scipy.integrate. However, here are some though:

Set a 3x3x3 cube to unit density. Then take the integration for each cube respectively, so you should have the volume cube V_ijk here. Now for each of sphere, you can get the mass of each sphere by summing V_ijk*D_ijk, where the D_ijk is the density of the sphere.

It should be much faster because you do not need to do integration now.

---

You can obtain an analytic formula for the intersecting volume between a cube (or rectangular prism) and a sphere. It won't be easy, but it should be possible. I have done it for an arbitrary triangle and circle in 2D. The basic idea is to decompose the intersection into simpler pieces, like tetrahedra and volumetric spherical triangle sectors, for which relatively simple volume formulas are known. The main difficult is in considering all the possible cases of intersections. Luckily both objects are convex, so you are guaranteed a single convex intersection volume.

An approximate method might be to simply subdivide the cubes until your approximate numerical integration algorithm does work; this should still be relatively fast. Do you know about Pick's Theorem? That only works in 2D, but there are, I believe, 3D generalizations.