2.1    Newton’s law of universal gravitation

Newton’s law of universal gravitation states that the gravitational force between any two objects is proportional to their masses, and inversely proportional to the square of the     distance between them:

In SI units, the gravitational constant G = 6 .674 x 10− 11  N m2  kg−2 . This is very small: gravity is a very weak force compared to electromagnetic and nuclear forces.  In this       assignment we will use a different system of units in which G = 1.

The image below illustrates Newton’s law of universal gravitation. A point mass m1          attracts another point mass m2  by a force F2  pointing along the line intersecting both   points.  In the same way, m2  attracts m1  with an equal force F1  pointing in the opposite direction.

Consider two masses m1  and m2 , as shown in the figure. The gravitational force exerted on m1  due to m2  will be

where r21  is the position vector from m2  to m1 , i.e. r1 - r2

Similarly, the force exerted on m2  due to m1  will be

and clearly F12  = -F21 : the forces are equal and opposite, as they should be according to Newton’s third law.

From Newton’s second law, we also have F = ma, so therefore we can see that the accelerations experienced by m1  due to m2  and vice versa will be

These accelerations will not be equal in magnitude unless the objects have equal masses.

2.2    Converting Newton’s laws to a numerical model

Let’s consider a system with two objects: a star with mass M and a planet with mass m . In a solar system like ours M > m and so the acceleration experienced by the star due to the planet will be small compared to the acceleration experienced by the planet due to the star. The gravitational force on the planet due to the star is then

To keep things simple, we will assume the acceleration on the star due to the planet is  small enough that we can neglect it, and put the star in a fixed position at the origin of our co-ordinate system. Thus rMm  = rm  - rM  becomes rm . We’ll drop the now              unnecessary subscript on rm , and replace the unit vector with r/|r|, to give a               simpler-looking expression:

We can work out the motion of the planet straight from this, just by using the laws of motion:

In order to calculate the motion of the planet around the star, we need to take the initial position and velocity of the planet as the starting positions, and integrate in order to       calculate the force.  Here we will not do this exactly but use a numerical method to get   an approximate result. This approach is particularly useful in cases where it is not            possible to find an exact solution, e.g. when we are taking general relativity into account or dealing with more than two bodies.  Here we use Euler’s method, as in module 9. This is the simplest, most naive method and definitely not the most accurate, as you will learn in later computing courses.

We replace the differentials to small finite differences and rewrite the force as

and then rearrange this to give the change in velocity of the planet in the time δt as

where δv and δt are very small finite quantities. We can define δt (usually known as the timestep) in our code as a parameter.

Similarly, we can write the change in the planet’s position over the same δt as δr = vδt

Given our starting point of the initial position and velocity of the planet, and a fixed timestep δt, we can calculate how r and v change:

We can use these equations to calculate the path of the planet, and create an animation of this, just as we simulated the path of a projectile with air resistance in module 9.

3.   Task instructions

3.1    Part A (40% of marks)

●  Download the template notebook provided.

● Complete and edit this template to create your own self-contained notebook, completing the tasks described.

● The template includes outline functions for you to complete to

1) calculate the gravitational force between two objects;

2)  update the position of the planet in each time step;

3) display an animation of a planet orbiting a star.

● You will also need to complete the introduction text cell, and add further code and text cells in the ”investigation” section.

● The template also includes some additional code to help you ensure each function is working correctly. You should leave the code in these cells unchanged. You may add your own tests if you wish, but these should be in a separate cell.

3.2    Part B

1) Two planets (20% of marks)

a) Copy the animate planet function from part A into a new code cell, and    adapt it so that instead of a single planet it shows the paths of two separate planets in the gravitational field of the same star. You should not take into  account the gravitational attraction between the planets in this section, only the force between each planet and the star.

Instead of a single position and velocity, your function will need to take as arguments two positions and velocities (a position and velocity for each    planet) and update two positions and velocities.

Your function should call move planet from part A, and hence make use of the existing force function. You do not need to repeat or change these      functions.

b) Test your code using a star with mass 1000 and the following initial conditions:

● first planet: position (1 .85, 0, 0), velocity (0 , 23 .25, 0)

● second planet: position (1 .96, 0 .4, 0), velocity (-4 .44, 21 .9, 0)

Describe and comment on the behaviour in a text cell.

2) Two planets with gravitational interactions (30% of marks)

a) Create a new function, similar to move planet, but this time updating the  positions of two planets and taking into account the gravitational attraction between the planets as well as between each planet and the star.

In addition to the positions and velocities of the two planets, this function will need to take two more arguments to represent the masses of the planets.

Your function should call the force function from part A: you do not need to repeat or change this function.

b) Create a new function, similar to animate planet but this time dealing with two planets and using the new function you have written to move the planets taking into account the gravitational attraction between the planets.

c) Test your code using the same initial conditions as for the previous case, and with each planet having a mass of 2 .0 in the units we are using for this         assignment.

Describe and comment on the behaviour in a text cell.

3)  Further investigation (10% of marks)

Try varying some of the parameters of your simulation: the masses of the star and the planets, and the starting positions and velocities of the planets. There are       many aspects you could investigate, but you should pick one or two. Add a few     relevant animations and text cells to describe and explain what you observe, along with a conclusion.

For example, you could

● calculate the required velocity to achieve a stable circular orbit with a given radius, and see how close or how massive the planets have to be to have a significant effect on each other’s orbit;

● turn one of the planets into a moon orbiting the other, by careful choice of masses and initial conditions. (This is hard.)