Cone Gravity

Consider an arbitrary point in the cone:

$\vec{R} = r\cos{\theta} \hat{i} + r\sin{\theta} \hat{j} + z \hat{k}$

The position vector of the point at which the force is to be computed is:

$\vec{p} = 2 \hat{i} + 0 \hat{j} + 0 \hat{k}$

The calculation is done in cylindrical coordinates. The elementary force vector is:

$d\vec{F} = \frac{G(\rho dV)\left(\vec{R}-\vec{p}\right)}{\lvert\vec{R}-\vec{p}\rvert ^3} = \frac{G\left(\rho \left(r d\theta dr dz\right)\right)\left(\vec{R}-\vec{p}\right)}{\lvert\vec{R}-\vec{p}\rvert ^3}$

Plugging in expressions, simplifying and recognising the pattern of summation yields:

$F_x = \int_{0}^{2} \int_{0}^{1 - \frac{z}{2}} \int_{0}^{2\pi} \frac{r\,\left(r\,\cos\left(\mathrm{\theta}\right)-2\right)}{{\left(r^2-4\,\cos\left(\mathrm{\theta}\right)\,r+z^2+4\right)}^{3/2}} d\theta dr dz$

$F_y = \int_{0}^{2} \int_{0}^{1 - \frac{z}{2}} \int_{0}^{2\pi} \frac{r^2\,\sin\left(\mathrm{\theta}\right)}{{\left(r^2-4\,\cos\left(\mathrm{\theta}\right)\,r+z^2+4\right)}^{3/2}}d\theta dr dz$

$F_z = \int_{0}^{2} \int_{0}^{1 - \frac{z}{2}} \int_{0}^{2\pi} \frac{r\,z}{{\left(r^2-4\,\cos\left(\mathrm{\theta}\right)\,r+z^2+4\right)}^{3/2}} d\theta dr dz$

Each of these monsters is solved using Wolfram-Alpha to obtain a net magnitude of $\boxed{\mid \vec{F} \mid \approx 0.50253}$

A code-based solution is as follows:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

clear all
clc

G    = 1;
rho  = 1;

% Point of interest:
rp   = [2;0;0];

% Initialisation of resultant force vector:
F = [0;0;0];

% Discretization of continuous space:
dz     = 1e-2;
dtheta = 1e-2;
dr     = 1e-2;

for z = 0:dz:2
    for r = 0:dr:(1 - 0.5*z)
        for theta = 0:dtheta:2*pi  

            % Element volume:
            dV = r*dr*dtheta*dz;

            % Arbitrary position vector (Cylindrical coordinates):
            R  = [r*cos(theta);r*sin(theta);z];

            % Newton's law of gravitation:
            dF = (G*rho*dV*(R-rp))/(norm(R-rp)^3);

            % Numerical integration
            F  = F + dF;
        end
    end
end

% Resultant force:
Answer = norm(F) % Standard vector 2-norm

% Answer = 0.511477965814259 (approximate)