(a) Use an iterative approach to answering the questions. In other words, write a program that calculates the gravitational force on the comet, its momentum, and its position in small time steps (i.e. small time intervals). The steps to doing the calculations are listed in Matter \& Interactions, Chapter 3.
 Define the values of constants such as G to use in the program.
 Specify the masses, initial positions, and initial momenta of the interacting objects.
 Specify and appropriate value for , small enough that the objects don't move very far during one update.
 Create a "loop" structure for repetitive calculations to calculate force on each object, the momentum of each object, the position of each object, and the clock reading.
Inside the loop, the program should do the following for each object in the simulation:
(1) Calculate the (vector) forces acting on the object and sum the forces to calculate the net force on the object.
(2) Update the momentum of the object: .
(3) Update the position of the object , where the average velocity can be approximated as . (Note: there are more accurate ways to approximate the average velocity.)
(4) Update the clock reading, .
The loop repeats, performing each of the calculations above for each time step, .
You can use a condition in the while loop or an if statement to stop the loop and print the position and momentum of the comet. Or, you can simply print the position and momentum of the comet and after one orbit, stop the simulation and examine the data in order to answer the questions.
Here is an example program in Python to simulate the motion of the comet. Note that the time interval is converted to seconds.
from visual import *
sun=sphere(radius=1e10, pos=(0,0,0), color=color.yellow)
comet=sphere(radius=1e10, pos=(8.79e10,0,0), color=color.white)
sun.m=2.0e30
comet.m=2.6e14
G=6.7e11
comet.v=vector(0,5.4e4,0)
comet.p=comet.m*comet.v
t=0
dt=24.0*3600
print "time", "\t", "Position", "\t", "Velocity"
while t<=30.0*24*3600:
rate(100)
print t, "\t", comet.pos, "\t", comet.v, comet.p
r=comet.possun.pos
rmag=mag(r)
rhat=r/rmag
F=(G*comet.m*sun.m/rmag**2)*rhat
comet.p=comet.p+F*dt
comet.v = comet.p/comet.m
comet.pos=comet.pos+comet.v*dt
t=t+dt
The program will stop at t = 30 days, which is seconds. The position and velocity printed by the program are:
You can easily print the momentum or calculate it. The momentum of the comet at this instant is
Since the initial conditions and constants are given to two significant figures, the position and momentum at t = 30 days can only be reported to two significant figures even though many more significant figures are used in the calculations.
(b) To find the period of the comet, to the nearest hour, use an infinite while loop (i.e. use while 1: and watch the simulation. You will see that the comet travels very slowly when far from Sun. As it approaches Sun, it dramatically speeds up. Stop the simulation just after it the comet completes one period. Examine the data to see when the xvelocity changes from a positive value to a negative value. That's the time interval during which it reached perihelion.
A sketch of the path of the comet is shown below.
Figure: The comet at perihelion.
The program above, with an infinite while loop, will print the following time, position, and velocity data near perihelion.
Note that between t = 556243200.0 s and t = 556329600.0 s, the comet's yposition changes from a negative value to a positive value, and the comet's xvelocity changes from a positive value to a negative value. Also, the maximum xposition of the comet occurs at t = 556243200.0 s. It is during this time interval that the comet crosses the +x axis and therefore completes one period. We don't exactly the instant when it reaches perihelion. However, our best guess for the period, using this simulation, is . Let's convert t his to years.
This comet will have a period of 17.6 y. Note that this data does not correspond to an actual, observed periodic comet. However, Halley's comet, which you have probably heard of, has a period of 76 y, its mass is between and kg, its distance from Sun at perihelion is m, and its perihelion speed is 54 km/s. Thus, the comet just simulated does have similar orbital parameters as Halley's comet.
The biggest numerical error in the simulation is in the approximation . There are various ways of improving this approximation and therefore finding a more accurate determination of the period of the comet.
(c) Aphelion occurs in this case when the comet crosses the –x axis as shown below.
Figure: The comet at aphelion.
Examine the data to find where the xvelocity changes from a negative value to a positive value. During the same time interval, the yposition should change from a positive value to a negative value. The time interval for the program above where the comet's yposition changes from positive to negative is:
Therefore, judging by the yposition, the comet reaches aphelion between t = 255225600.0 s and t = 255312000.0 s. This is t = 8.09 y. You will notice something quite peculiar about the data shown above. First, when the comet crosses the –x axis, its xvelocity did not change signs but stayed negative. In fact, by looking through the data, we can find the time interval when the xvelocity changes signs as shown below.
During this interval, the xvelocity goes from a negative value to a positive value. Thus, using xvelocity as an indicator, aphelion occurs at t = 276220800.0 s = 8.76 y. Also, during this interval, the xposition is most negative, as expected. The time for the comet to reach aphelion should be half the period, and 17.6 / 2 = 8.8 y which is consistent. This suggests that we should use the xvelocity as an indicator of aphelion rather than the yposition which seems to be much more uncertain. But either way, we see that the xposition, to three significant figures, is 1.94e+12 m. Rounding to two significant figures, the comet's aphelion distance from Sun is , since we can only be accurate to two significant figures anyway. Regardless of some peculiarity in the data, we come to the same conclusion for the aphelion distance of the comet from Sun, to two significant figures.
Because of our approximation that and because numerical errors add up (and we are adding error during thousands of time steps), then our answers will have numerical error. But to only two significant figures, the simulation gives fairly good estimates of the comet's actual period and actual aphelion distance.
