\documentstyle[a4wide,12pt]{article}

\begin{document}
\section*{Introduction to numerical projects}


Here follows a brief recipe and recommendation on how to write a report for each
project.
\begin{itemize}
\item Give a short description of the nature of the problem and the eventual 
numerical methods you have used.
\item Describe the algorithm you have used and/or developed. Here you may find it convenient
to use pseudocoding. In many cases you can describe the algorithm
in the program itself.

\item Include the source code of your program. Comment your program properly.
\item If possible, try to find analytic solutions, or known limits
in order to test your program when developing the code.
\item Include your results either in figure form or in a table. Remember to
       label your results. All tables and figures should have relevant captions
       and labels on the axes.
\item Try to evaluate the reliabilty and numerical stability/precision
of your results. If possible, include a qualitative and/or quantitative
discussion of the numerical stability, eventual loss of precision etc. 

\item Try to give an interpretation of you results in your answers to 
the problems.
\item Critique: if possible include your comments and reflections about the 
exercise, whether you felt you learnt something, ideas for improvements and 
other thoughts you've made when solving the exercise.
We wish to keep this course at the interactive level and your comments can help
us improve it.
\item Try to establish a practice where you log your work at the 
computerlab. You may find such a logbook very handy at later stages
in your work, especially when you don't properly remember 
what a previous test version 
of your program did. Here you could also record 
the time spent on solving the exercise, various algorithms you may have tested
or other topics which you feel worthy of mentioning.
\end{itemize}



\section*{Format for electronic delivery of report and programs}
%
The preferred format for the report is a PDF file. You can also
use DOC or postscript formats. 
As programming language we prefer that you choose between C/C++, Fortran or Python.
The following prescription should be followed when preparing the report:
\begin{itemize}
\item Use Fronter to hand in your projects, log in  at 
blyant.uio.no and choose 'fellesrom fys3150'.
Thereafter you will see an icon to the left with 'hand in' or 'innlevering'.
Click on that icon and go to the given project. 
There you can load up the files within the deadline.
\item Upload {\bf only} the report file and the source code file(s) you have developed.
The report file should include all of your discussions and a list of the codes you have developed. 
Do not include library files which are available at the course homepage, unless you have
made specific changes to them.
\item Comments  from us on your projects, approval or not, corrections to be made 
etc can be found under
your Fronter domain and are only visible to you and the teachers of the course.

\end{itemize}
Finally, 
we do prefer that you work two and two together. Optimal working groups consist of 
2-3 students. You can then hand in a common report. 





\section*{Project 5, The solar system, deadline   November 30  12am (midnight)}


We study first a hypothetical solar system with one planet, say Earth, which orbits around the Sun.
The only force in the problem is gravity. Newton's law of gravitation  is given by a force $F_G$
\[
F_G=\frac{GM_{\mathrm{sun}}M_{\mathrm{Earth}}}{r^2},
\]
where $M_{\mathrm{sun}}$ is the mass of the Sun and $M_{\mathrm{Earth}}$ is the mass of Earth. The gravitational constant is $G$ and $r$ is the distance between Earth and the Sun.
We assume that the sun has a mass which is much larger 
than that of Earth. We can therefore safely neglect the 
motion of the sun in this problem.  In the first part of this project, your aim is to compute the motion
of the Earth using different methods for solving ordinary differential equations.

We assume that the orbit of Earth around the Sun 
is co-planar, and we take this to be the $xy$-plane.
Using Newton's second law of motion we get the following equations
\[
\frac{d^2x}{dt^2}=\frac{F_{G,x}}{M_{\mathrm{Earth}}},
\]
and 
\[
\frac{d^2y}{dt^2}=\frac{F_{G,y}}{M_{\mathrm{Earth}}},
\]
where $F_{G,x}$ and $F_{G,y}$ are the $x$ and $y$ components of the gravitational force. 
\begin{enumerate}
\item[a)]  Rewrite the above second-order ordinary differential equations as a set of coupled first order
differential equations. Write also these equations in terms of dimensionless variables.
As an alternative to the usage of dimensionless variables, you could also
use so-called  astronomical units (AU as abbreviation). If you choose the latter set of units, 
one astronomical unit of length, known as 1 AU, is the average distance between the Sun and Earth, that is
$1$ AU = $1.5\times 10^{11}$ m.  It can also be convenient to use years instead of seconds since years match
better the solar system. The mass of the Sun is $M_{\mathrm{sun}}=M_{\odot}=2\times 10^{30}$ kg. The mass of Earth is
$M_{\mathrm{Earth}}=6\times 10^{24}$ kg. The mass of other planets like Jupiter is 
$M_{\mathrm{Jupiter}}=1.9\times 10^{27}$ kg and its distance to the Sun is 5.20 AU. Similar numbers for Mars
are $M_{\mathrm{Mars}}=6.6\times 10^{23}$ kg and  1.52 AU, for Venus $M_{\mathrm{Venus}}=4.9\times 10^{24}$ kg and  0.72 AU and for example Saturn are $M_{\mathrm{Saturn}}=5.5\times 10^{26}$ kg and  9.54 AU.

Finally,  mass units can be obtained by using the fact that Earth's orbit is almost circular around the Sun.
For circular motion we know that the force must obey the following relation
\[
F_G= \frac{M_{\mathrm{Earth}}v^2}{r}=\frac{GM_{\odot}M_{\mathrm{Earth}}}{r^2},
\]
where $v$ is the velocity of Earth. 
The latter equation can be used to show that
\[
v^2r=GM_{\odot}=4\pi^2\mathrm{AU}^3/\mathrm{yr}^2.
\]
Discretize the above differential equations and set up an algorithm for solving these equations using the so-called Euler-Cromer method discussed in the lecture notes, chapter 13.
\item[b)]  Write then a program which solves the above differential equations for the Earth-Sun system
using the Euler-Cromer method. Find out which initial value for the velocity that gives a circular orbit
and test the stability of your algorithm as function of different time steps $\Delta t$. Find a possible maximum
value $\Delta t$ for which the Euler-Cromer method does not yield stable results. Make a plot of the results you obtain for the position of Earth (plot the $x$ and $y$ values) orbiting  the Sun.

Check also for the case of a circular orbit that both the kinetic and the potential energies are constants.
Check also that the angular momentum is a constant. Explain why these quantities
are conserved.
\item[c)] Modify your code by implementing the fourth-order Runge-Kutta method and compare the stability of your results by repeating the steps in b). Compare the stability of the two methods, in particular as functions
of the needed step length $\Delta t$. 

\item[d)] Kepler's second law states that the line joining a planet to the Sun sweeps out equal areas in equal times. Modify your code so that you can verify Kepler's second law for the case of an elliptical orbit. 
Compare both the Runge-Kutta method and the Euler-Cromer method and check that the total energy and angular
momentum are conserved. Why are these quantities conserved?  
A convenient choice of starting values are an initial position of 1 AU and an initial velocity of 5 AU/yr.

\item[e)] Consider then a planet which begins at a distance of 1 AU from the sun. Find out by trial and error
what the initial velocity must be in order for the planet to escape from the sun.  Can you find an exact answer?

\item[f)]  We will now study the three-body problem, still with the Sun kept fixed at the center but 
including Jupiter (the most massive planet in the solar system, having a mass that is approximately 1000 times
smaller than that of the Sun) together with Earth. This leads us to a three-body problem. Without Jupiter, Earth's motion is stable and unchanging with time. The aim here is to find out how much Jupiter alters Earth's motion.

The program you have developed can easily be modified by simply adding the magnitude of the force betweem Earth and Jupiter.

This force is given again by 
\[
F_{\mathrm{Earth-Jupiter}}=\frac{GM_{\mathrm{Jupiter}}M_{\mathrm{Earth}}}{r_{\mathrm{Earth-Jupiter}}^2},
\]
where $M_{\mathrm{Jupiter}}$ is the mass of the sun and $M_{\mathrm{Earth}}$ is the mass of Earth. 
The gravitational constant is $G$ and $r_{\mathrm{Earth-Jupiter}}$ is the distance between Earth and Jupiter.

We assume again that the orbits of the two planets are co-planar, and we take this to be the $xy$-plane. 
Modify your first-order differential equations in order to accomodate both the
motion of Earth and Jupiter by taking into account the distance in $x$ and
$y$ between Earth and Jupiter. Set up the algorithm and plot the positions of Earth and Jupiter using the fourth-order Runge-Kutta method.

As you will notice, the influence on Earth from Jupiter is very small. Repeat 
these calculations by increasing the mass of Jupiter by a factor of 10 and 1000
 and plot the position of Earth. Investigate also the effect on Mars (that is replace Earth with Mars).
Comment  your results.




\end{enumerate}


\end{document}