Please write or print out all programs and results. Many people find it convinie nt to create a document,
using Word (or any other word processor) and then copy and paste programs, output, figures, and so
forth into the document. Alternative ly I have inserted a document in blackboard in which you can
insert your answers (instructions what to do are in the document). When you are finished you can pr
int out the document and hand it in to the office or upload to the course Blackboard page. For full
marks the program should also be neat, tidy and commented.
The Solar System is a gravitational system consisting of eight planets and the Sun. We will simulat
and animate the evolution of the system, compare the numerical integration to an analytic solution
and empirically determine Kepler’s third law.
The equation of motion for a planet around a central star is given by Newton’s two laws, F = ma and
GM m
r = ma
(1)
r 3
with m the mass of the Sun and G the gravitational constant. The trajectory of a planet can be
calculated from the initial values [r 0 ,v 0 ] by solving the differential equation above. In this project we
assume that the planets rotate in a circular orbit soley considering the x and y coordinates. From the
solutions of the orbits we can measure the orbital periods and fit the relations between the size of ht
eorbit and the orbital period (which should give Kepler’s law). Read the instructions caregully to the
end first, to make sure you ask all the questions you have so that you can finish the exercise yourself.
F grav = −
1. Read in the radius of orbit of the planets: Write a function that reads in the initial parameters
of the planets from the file solarsystem.txt. The columns are the name of the planet, the orbital
radius (in so-called Astronomical Units, AU). Use the numpy loadtxt command and look up its
options. The function should return a numpy array with two columns: planets, radius. The initial
positions can be take as the radius. [HINT: the input file is a combination of text and floats so you
will need to take this into account.]
2. Calculate the initial velocities: Write a function that calculate p the initial velocities of the
planets, assuming they have the Kepler velocity of a vircular orbit (v = GM/r with M the mass of
the Sun). We will use units of AU and year, i.e. express the velocity in AU/yr. [HINT: scale to the
velocity of the Earth, which in these units is 2π.]
3. Convert to x and y values: Write a function that converts the position and velocities with the
angle to x, y, v x and v y values (NOTE that the velocities are perpendicular to the orbit!). Probably
best to make a sketch to make sure you get all the sin and cos right. The function should return 4
values as a list r.
4. Write the derivatives for the differential equations: Write a funciton deriv that has r nd
t as input variables, with r a list with 4 values: x, y, v x and v y . Rewrite the second order differential
equation for the vectors r and a in four first order differential equations for x, y, v x and v y . The
funciton returns a lit with the derivatives of x, y, v x and v y . Make sure you have GM in the right
units (x, y, in AU, v x , v y in AU/yr and a x , a y in AU/yr 2 ).
5. Calculate the orbits: Define a numpy array t that spans 1000 years with 100000 points and
in a for loop over the plants calculate the initial values of x, y, v x and v y using the initial data and
conversion routine. Use [login to view URL] function to calculate the values of x, y, v x and v y
as a funciton of time.
1PHYS 311 – Computational Physics I: Class Project
6. Determine the periods: Write a function that takes the x or y values of the solved orbit and t
as input and returns the period of the orbit and its error, based on the fact that during one period the
x and y values twice cross zero and thus change sign. Use the following steps to to calculate these:
i) use the [login to view URL]() fucntion to get an array with +1 and -1 values for the elements that are
positive and negative
ii) use the [login to view URL]() funtion to calculate the differences between the consecutive elements. This
array has values of 0 unless the original array changes sign
iii) use [login to view URL]() function to get the indices of the array where this happens
iv) use this to mask the times t (t[np.where...)]). Finnally use the [login to view URL]() function again to
calculate the time differences between the zero corssings (i.e. half of the period)
Calculate the period. Alternatively, if we get to this part during the course, you can use [login to view URL]
to calculate the fast Fourier transform. Use [login to view URL] to calculate the frequencies of the
resulting spectrum. The orbital frequency (i.e. 1/P ) should be the frequency corresponding to the
largest value of the spectrum (in the first half of the spectrum). The function should return the period.
7. Determine Kepler’s law: Use the function above to determine the period on the orbits of each
of the planets and plot the period againts the radius of hte orbit a. Use [login to view URL] fit
or [login to view URL] to fit a power law to the data P ∝ a γ (decide on whether you want to fit the actual
values or the log of the values). Compare the fitted value of γ to the expected value of γ = 3/2.
8. Plot the Solar System Plot the orbits of the eight planets in a properly annotated plot.
Optional make an animation of the Solar System.
Table of constants
Solar Mass
Gravitational Constant
Astronomical Unit
Year
M ⊙
G
AU
yr
1.99 x 10 30
6.67384 x 10 −11
1.496 x 10 11
3.15569 x 10 7
kg
m 3 /kg/s 2
m
s
GOOD LUCK AND HAVE SOME FUN!
2
Hi there, I have seen the PDF and the questions listed there. There are 8 questions and I can do it. The submission date is 22nd. I can deliver you before that. Please let me know if you are okay with this. Regards.