The goal of this project was to simulate the motion of two planets going in orbit around a star. The program is written to simulate the Earth, Venus and the Sun as a system. It will update the position of both the Earth and Venus based on the total gravitational forces exerted on the respective bodies. The time interval of the updates can be chosen by the user.
Newton's law of universal gravitation states that every point mass in the universe attracts every other point mass, with a force equal to:
where G is the gravitational constant (6.674 x 10-11 m3 kg-1 s-2), r is the distance between the two bodies, and m1 and m2 are their respective masses.
This model is set in a two-dimensional coordinate system, with the Sun placed in the middle. The Sun is assumed to be stationary. The planets start out with their average orbital speed in both x- and y- direction.
One must find the total force exerted on each of the planets in order to update their positions. Once the total force acting on a body is known it can be decomposed into x- and y-components, named Fx and Fy, respectively. In order to do this one needs the angle theta, which is given by:
where dy and dx are given by:
The distance between the bodies can then be easily calculated with the following formula:
Fx and Fy is calculated by:
Once the sum of forces in both x- and y-direction is known one can calculate the new velocities of the bodies:
The velocities can then be used to update the planetary positions using these formulas:
The code is available here.
The following parameters were used for the various calculations needed for the simulation:
To simulate the system a class called 'Body' was written. Objects of this class contain all the information needed to compute and orient in the system. When the program starts the main function of main.cpp creates three bodies: earth, venus and sun.
earth gets starting point (AU, 0), venus gets (-distanceFromSun, 0) and sun (0,0). Since orbital speed is perpendicular to the force vector, and points in the direction of motion, the speed in both x- and y- direction have been chosen as shown in figure 2.
When the bodies have been created the program starts calculating the forces. The function 'force' takes two bodies as input, calculates the gravitational force between them in both in x- and y- directions, and returns the results.
The function 'update' takes all the bodies as input and calculates the next step in the simulation. Step size is chosen by 'deltaT'. Here the total forces on both earth and venus are calculated, and velocities and positions are updated as described in the theory-section.
There are three functions for showing results; 'system_gif', 'distance_graph' and 'orbit_trajectory'. All three take the inputs deltaT and n. deltaT (described above) is entered into the update function. n is the number of updates to be run. It also decides the sizes of the datasets that are used for the graphs. Only one of the three functions can be run at one time, as they all write to the same canvas. These functions need to be uncommented in main(), depending on desired output, which seems like a practical solution.
Below are examples of the three different types of results produced by the program.
The system gif:
Weekly update over one year (366/7 steps) (with tracing)
Discussion and conclusion:
The graph showing the trajectories of the planets seem to agree well with expectations: circular trajectories with constant radius. The graph displaying the distance between the planets show a clear periodicity, with little variation from peak to peak. The gif with tracing made from 24-hour update over 366 days (steps) shows irregular orbits with non-constant radiuses. This could be because of some error with the gif setup as I experienced objects sliding while drawing gifs earlier in the project.