C++

# Optional prize problem

If you finish the rest of the lab early and would like to enter for the Part 1B computing prize then you may attempt this exercise. There are no marks for this exercise. In order to be considered for the Prize, your solution must be received by Prof Roberto Cipolla by one week after the end of the Lent Full Term 2010, in the form specified at the end of this document.

Numerically integrate the orbits of celestial bodies in order to predict when solar eclipses will occur. Given below are the positions, velocities, masses, and radii of the Moon, Earth, and Sun at midnight on the beginning of 7 January, 2008. The Cartesian coordinates are given relative to the barycentre of the solar system, with axial directions determined by the orbital plane of the earth. Positions are given in astronomical units (AU), and velocities are given in AU per day:

 Body x (AU) y (AU) z (AU) vx (AU/day) vy (AU/day) vz (AU/day) mass (kg) radius (m) Moon -2.698040026379963×10-1 9.478530342547004×10-1 -2.904143638119555×10-4 -1.626731661042607×10-2 -4.779136494151825×10-3 2.734493101618888×10-5 7.349×1022 1.7374×106 Earth -2.697539791224478×10-1 9.505128628199724×10-1 -8.468828035093709×10-5 -1.683798264597689×10-2 -4.785904833082302×10-3 -9.942252208381242×10-8 5.9736×1024 6.3781×106 Sun 1.177670680164471×10-4 4.956736129503275×10-3 -7.039430413584526×10-5 -6.401003620910089×10-6 5.939607737084592×10-7 1.106864330607326×10-7 1.9891×1030 6.95500×108

The formula for gravitational attractive force between bodies A and B is given by

 FAB = G⋅mA⋅mB dAB2

where G is the gravitational constant, mX is the mass of body X in kilograms, and dAB is the distance between the two bodies in meters. Use these values for G and the astronomical unit:

 G = 6.67300×10-11 m3⋅kg-1⋅s-2 1 AU = 1.49598×1011 m

## RUNGE-KUTTA INTEGRATION

The fourth-order Runge-Kutta method (RK4) for numerical integration evaluates the derivative at 4 different points at each step. Consider the initial value problem:

 y' = f(y) y(t0) = y0

Choose a step size h. Then the RK4 update step is given by:

k1=f(yn)
k2=f(yn + ½⋅h⋅k1)
k3=f(yn + ½⋅h⋅k2)
k4=f(yn + h⋅k3)
yn+1=
 yn + h (k1 + 2k2 + 2k3 + k4) 6

The per-step error of the RK4 method is on the order of h5. In contrast, the Euler method has a per-step error of h2, requiring much smaller h for the same error bound.

## EVALUATION

Try the Euler method, the RK4 method, and perhaps other numerical integration methods (e.g. Gaussian Quadrature) to compute the trajectories of the bodies, and calculate the expected times for solar eclipses. Try varying the step size of integration and the angular alignment tolerance for what should be considered an eclipse. A list of predicted eclipse dates is given by Wikipedia. How far into the future can your program correctly predict eclipses relative to the list? How does this vary with step size and integration method?

Your answer must be submitted entirely on paper, and sent to Prof Roberto Cipolla through the internal mail (hand it in at the Baker Building Enquiry office). It should contain the following parts:

• Description of the technique you have used to solve the problem.
• Source code for the programs you have written. Please include appropriate comments.
• The estimated eclipse times and any other computed data.