Newton's Law of gravitation tells us the gravitational force on a mass, M, due to another mass, m, is
where G is the gravitational constant (big "G"). The direction of the force on M is towards m. We can introduce the acceleration field due to the mass, m, to be
where r is the distance from the mass, m, and the unit r-vector is outwards from m. The direction of the acceleration at each point is inwards towards the mass (hence the appearance of the negative sign). The acceleration field gives the force per unit "test mass" placed in the field.
The following simulation shows such a gravitational field. Note that we use the acceleration as the field rather than the force, as the force on a test object placed in the field will depend on the mass of the test object.
The simulation initially shows a central mass which is fixed in place. The size of the central mass can be changed using the slider. The strength of the field can also be changed using a slider. The third slider sets the mass of a test object that can be inserted into the field. After setting the mass using the slider, clicking anywhere in the field to insert a mass (but not inside an existing mass). Initially start with interactions on (this means masses collide like hard spheres). When interaction is removed, the test sphere will pass through the central mass. The strength of the field, size of the central mass and number of test masses can be modified while the simulation is running. Modify these variables to observe how a changing gravitational field impacts the test masses.
Some questions to consider while running the simulation are
bigG = 1.0; % Gravitational constant
mass1 = 1.0; % Mass of centre object
mass2 = 1.0; % Mass of falling object
radius1 = 1.0; % Radius of centre object
radius2 = 0.5; % Radius of falling object
dphi0 = 0.2; % Initial tangential velocity
dr0 = 0.0; % Initial radial velocity
phi0 = 0.0; % Initial tangential position
r0 = radius1+radius2+1.0; % Initial radial position
dt = 1.0e-3; % Simulation step size
outputRate = 100; % Number of steps between output cycles
outputLength = 400; % Length of simulation (# output files)
% Formulate our problem as a system of ODES: y' = odefun(t, x)
odefun = @(t, x) [x(2); -bigG*mass1/x(1).^2; dphi0];
xx = [r0; dr0; phi0];
% Functions to calculate the mass position
px = @(xx) xx(1)*sin(xx(3));
py = @(xx) xx(1)*cos(xx(3));
pxc = @(xx) px(xx) - radius2;
pyc = @(xx) py(xx) - radius2;
figure;
d2 = radius2*2;
x = [px(xx)]; y = [py(xx)];
h = plot(x, y);
ch = rectangle('Position',[pxc(xx) pyc(xx) d2 d2],'Curvature',[1,1]);
bh = rectangle('Position',[-radius1 -radius1 2*radius1 2*radius1],...
'Curvature',[1,1]);
axis([-4 4 -4 4]);
while ishandle(h)
% Solve the ODE using a first order method (might be inacurate)
for n = 1:outputRate
xx = xx + odefun(0.0, xx)*dt;
if xx(1) < radius1 + radius2
xx(1) = 2.0*(radius1+radius2) - xx(1);
xx(2) = -xx(2);
end
xx(3) = mod(xx(3), 2.0*pi);
end
% Update the graph
x = [x px(xx)]; set(h, 'xdata', x);
y = [y py(xx)]; set(h, 'ydata', y);
set(ch, 'Position',[pxc(xx) pyc(xx) d2 d2]);
pause(0.05);
drawnow;
end
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
bigG = 1.0; # Gravitational constant
mass1 = 1.0; # Mass of centre object
mass2 = 1.0; # Mass of falling object
radius1 = 1.0; # Radius of centre object
radius2 = 0.5; # Radius of falling object
dphi0 = 0.2; # Initial tangential velocity
dr0 = 0.0; # Initial radial velocity
phi0 = 0.0; # Initial tangential position
r0 = radius1+radius2+1.0; # Initial radial position
dt = 1.0e-3; # Simulation step size
outputRate = 100; # Number of steps between output cycles
outputLength = 400; # Length of simulation (# output files)
odefun = lambda x: np.array([x[1], -bigG*mass1/x[0]**2, dphi0]);
xx = np.array([r0, dr0, phi0]);
px = lambda x: xx[0]*np.sin(xx[2])
py = lambda x: xx[0]*np.cos(xx[2])
pxc = lambda x: px(x) - radius2
pyc = lambda x: py(x) - radius2
fig, ax = plt.subplots()
x = [px(xx)]
y = [py(xx)]
line, = ax.plot(x, y, lw=1)
ax.set_xlim(-4, 4)
ax.set_ylim(-4, 4)
circle1 = plt.Circle((0.0, 0.0), radius1)
circle2 = plt.Circle((px(xx), py(xx)), radius2, fill=False)
ax.add_artist(circle1)
ax.add_artist(circle2)
def animate(step_number):
global xx
# Calculate the new position (1st order, might be inaccurate)
for i in xrange(0, outputRate):
xx = xx + odefun(xx)*dt
if xx[0] < radius1+radius2:
xx[0] = 2.0*(radius1+radius2) - xx[0]
xx[1] = -xx[1]
xx[2] = np.fmod(xx[2], 2.0*np.pi)
# Update line/circle
x.append(px(xx))
y.append(py(xx))
line.set_data(x, y)
circle2.center = px(xx), py(xx)
return line,circle2
ani = animation.FuncAnimation(fig, animate, xrange(1, outputLength),
interval=25, blit=True, repeat=False)
#ani.save('output.mp4')
plt.show()
#include <fstream>
#include <cmath>
const double bigG = 1.0; // Gravitational constant
const double mass1 = 1.0; // Mass of centre object
const double mass2 = 1.0; // Mass of falling object
const double radius1 = 1.0; // Radius of centre object
const double radius2 = 0.5; // Radius of falling object
const double dphi0 = 0.2; // Initial tangential velocity
const double dr0 = 0.0; // Initial radial velocity
const double phi0 = 0.0; // Initial tangential position
const double r0 = radius1+radius2+1.0; // Initial radial position
const double dt = 1.0e-3; // Simulation step size
const unsigned outputRate = 100; // Number of steps between output cycles
const unsigned outputLength = 400; // Length of simulation (# output files)
int main(int argc, char** argv) {
double r = r0, dr = dr0;
double phi = phi0, dphi = dphi0;
double t = 0.0;
std::ofstream fp("output.dat");
fp << "# t\tr\tphi\n";
fp << t << '\t' << r << '\t' << phi << '\n';
for (unsigned i = 0; i < outputLength; ++i) {
for (unsigned j = 0; j < outputRate; ++j) {
// Use a first order approximation (might be inaccurate)
double oldr = r;
r += dt*dr;
dr -= dt*bigG*mass1/oldr/oldr;
// Test bounce condition
if (r < radius1 + radius2) {
r = 2.0*(radius1 + radius2) - r;
dr = -dr;
}
phi = fmod(phi + dt*dphi, 2.0*M_PI);
}
// Write output
fp << t << '\t' << r << '\t' << phi << '\n';
}
return 0;
}
set terminal gif animate 10
set output 'output.gif'
set xrange [-4:4]
set yrange [-4:4]
set size ratio -1
x(r, phi) = r*sin(phi)
y(r, phi) = r*cos(phi)
radius2 = 0.5
set object 1 circle at 0,0 size 1.0 fc rgb "navy"
stats "output.dat" nooutput
do for [i=1:int(STATS_records)] {
plot "output.dat" u (x($2, $3)):(y($2, $3)) every ::::i w lines notitle, \
"output.dat" u (x($2, $3)):(y($2, $3)):(radius2) \
every ::i::i with circles notitle
}