\documentclass[a4paper,12pt]{article}

\usepackage[left=25mm,right=25mm,top=20mm,bottom=20mm]{geometry}
\usepackage{hyperref,graphicx}
\pagestyle{empty}
\renewcommand{\refname}{}

\begin{document}

\begin{center}
{\bf SIMPLIFIED DYNAMICS OF ASTEROID\\ACCESSING OUR SOLAR SYSTEM\\}
\vspace{5mm}
{\it Lud\v{e}k Hyn\v{c}\'{\i}k $^1$, V\'aclav Uhl\'{\i}\v{r} $^2$\\}
\vspace{5mm}
$^1$University of West Bohemia, $^2$Franti\v{s}ek Pe\v{s}ta's Observatory
\end{center}

\noindent{\bf Introduction}

\vspace{3mm}\noindent\rm
Mankind explores the space from time out of mind. There are many bodies in the space influencing more or less the life on our Earth. To describe the motion of various spatial bodies may be useful even for future. This paper does not solve the problem in detail. It develops only a very simple model based on multi-body dynamics that describes spatial motion of a asteroid flying inside our Solar system and it is more useful for students entering the dynamics of bodies and mystery of MATLAB \cite{bal,mat}.\\

\noindent{\bf Physical model}

\vspace{3mm}\noindent\rm Considering the dimension of modelled space, the whole Solar system is simulated as a system of mass points. There is the Sun that is fixed and all the other planets inluding our Earth's Moon circulating around along the given trajectories. These trajectories are fixed by boundary conditions based on astronomical literature \cite{kle,roc,van}. For each planet (or the Earth's Moon), semimajor axis $a$ and excentricity $e$ of its trajectory are defined \cite{kle,roc,van}. Furthermore, period of revolution $T$, perihelium length $\varphi$, argument of perihelium width $\omega$ and inclination related to the ecliptic of Earth $i$ are defined \cite{kle,roc,van}. These characteristics define the initial position and motion of all the planets and the Earth's Moon. Moreover, each body has its mass $m$ and also the radius $R$ of each body supposing spherical shape is considered \cite{kle,roc,van} since the collision detection.\\

\noindent The asteroid is given by initial position and velocity in space (three components for position and three components for velocity). Each planet including the Sun and the Earth's Moon acts by gravitational force to the asteroid. No planet (or the Earth's Moon) deflection is considered. Hence, one have defined the three ordinary differential equations of motion of the second order

\begin{equation}
m\ddot{\bf x}(t)={\bf F}_{total}(t)
\end{equation}

\noindent with initial conditions specifying initial position vector and initial velocity vector

\begin{equation}
{\bf x}(0)={\bf x}_0\mbox{,}\quad\dot{\bf x}(0)={\bf v}_0
\end{equation}

\noindent for three spatial positions where $m$ and ${\bf x}$ are mass and vector of position of the asteroid and ${\bf F}_{total}(t)$ is the sum of all force vectors acting on the asteroid at a given time

\begin{equation}\label{force}
{\bf F}_{total}(t)=\sum\limits_{j=1}^{11}{\bf F}_j(t)
\end{equation}

\noindent where parameter $j$ covers the Sun, the Earth's Moon and nine further planets in our Solar system. These vectors depend on relative position of asteroid and particular interacting body that can be simply developed. The gravitational force vector ${\bf F}$ between two bodies is defined as

\begin{equation}
{\bf F}=\kappa\frac{m_1m_2}{\vert{\bf r}\vert^3}{\bf r}
\end{equation}

\noindent where $\kappa$ is the gravitational constant, $m_1$ and $m_2$ are the masses of the two bodies and the vector ${\bf r}$ is the distance vector betweem the two bodies. To know this vector at each solution time step, one has to know the position of each body at a given time. The position of the asteroid is unknown solution of our problem but each other body in the system has a given trajectory. The current position ${\bf x}_j$ of the particular body of the mass $m_j$	 on the trajectory at time $t$ is computed by

\begin{equation}
{\bf x}_j(t)=\left[\begin{array}{ccc}\cos{i} & 0 & -\sin{i} \\
                                           0 & 1 &        0 \\
                                     \sin{i} & 0 &  \cos{i} \end{array}\right]
             \left[\begin{array}{ccc}\cos{\omega} & -\sin{\omega} & 0 \\
                                     \sin{\omega} &  \cos{\omega} & 0 \\
                                                0 &             0 & 1 \end{array}\right]
             \left[\begin{array}{c}ae\cos{\frac{2\pi t}{T}+\varphi}                \\
                                   \sqrt{a^2-a^2e^2}\sin{\frac{2\pi t}{T}+\varphi} \\
                                                                                 0 \end{array}\right]
\end{equation}

\noindent where all the variables used for a particular body are described above. The Sun is static. Rewritting it to form

\begin{equation}
{\bf x}_j(t)={\bf R}_j^i{\bf R}_j^{\omega}\bar{\bf x}_j(t)\mbox{,}
\end{equation}

\noindent vector $\bar{\bf x}_j(t)$ indicates the position at the allipse at the given time, matrix ${\bf R}_j^{\omega}$ means the rotation of the ellipse in its plane around its focal point by angle $\omega$ and matrix ${\bf R}_j^i$ means the spatial rotation of the elipse by angle $i$ in relation to the ellipse of the Earth. Hence, the particular force vector ${\bf F}_j$ in equation \ref{force} can be written in form

\begin{equation}
{\bf F}_j(t)=\kappa\frac{m_{j}m}{\vert{\bf x}_j(t)-{\bf x}(t)\vert^3}\left[{\bf x}_j(t)-{\bf x}(t)\right]
\mbox{.}
\end{equation}

\noindent{\bf Model in MATLAB}

\vspace{3mm}\noindent\rm The whole situation can be very shortly implemented into MATLAB. The basic geometrical and dynamical values are saved in a matrix with twelve rows for the Sun, the Earth's Moon, nine particular planets and the asteroid. Particular columns define mass, radius, semimajor axis, excentricity, period of revolution, perihelium length, argument of perihelium width and inclination respectively. Where the particular position inside matrix is not applicable, there is zero (e.g. the Sun has no inclination) or some further specific values (e.g. initial position and velocity for the asteroid).\\

\noindent The system of three ordinary differential equation of the second order is splitted into the system of six ordinary differential equation of the first order using substitution as

\begin{equation}\label{ode}
\begin{array}{ccc}\dot{\bf x}(t)  & = & {\bf v}(t)\\
                  m\dot{\bf v}(t) & = & {\bf F}_{total}(t)\end{array}\quad\Leftrightarrow\quad
\begin{array}{c}\dot{\bf X}(t)={\bf R}(t)\end{array}
\end{equation}

\noindent where

\begin{equation}
{\bf X}(t)=\left[\begin{array}{c}{\bf x}(t)\\
                                 {\bf v}(t)\end{array}\right]\quad\mbox{and}\quad
{\bf R}(t)=\left[\begin{array}{c}{\bf v}(t)\\
                                 \frac{1}{m}{\bf F}_{total}(t)\end{array}\right]
\end{equation}

\noindent with initial condition

\begin{equation}
\dot{\bf X}(0)=\left[\begin{array}{c}{\bf x}_0\\
                                     {\bf v}_0\end{array}\right]\mbox{.}
\end{equation}

\noindent Consider that both vectors ${\bf x}(t)$ and  ${\bf v}(t)$ are defined in space and hence, the both final vectors ${\bf X}(t)$ and ${\bf R}(t)$ have six components. To the system of the first order ordinary differential equations \ref{ode}, {\bf ODE23} inbuilt MATLAB function is applied. Furthermore, the {\bf EventFunction} is used in order to detect possible collision between asteroid and some other body.\\

\begin{figure}
\centerline{\includegraphics[width=150mm]{Hyncik1.png}}
\vspace{-10mm}\caption{Asteroid's trajectory during 2005 years}\label{Hyncik1}
\end{figure}

\begin{figure}
\centerline{\includegraphics[width=150mm]{Hyncik2.png}\vspace{-10mm}}
\caption{Asteroid crashed into the Sun after 54 years}\label{Hyncik2}
\end{figure}

\noindent\bf{Simulation and results}

\vspace{3mm}\noindent\rm M-file called {\it MainRun.m} defines all the values, calls the {\it SystemODE.m} function that is a right-hand side of the system of the ordinary differential equations and provides figures including results and animation. There is the {\bf EventFunction} called {\it Crash.m} set as the parameter of the {\bf ODE23} solver to detect possible collision and stop the simulation in case of a collision. Supplementary function {\it GetPosition.m} returns the vector of relative positions between the asteroid and particular body at the given time.\\

\noindent Applying sutable initial conditions especially for position and velocity of the asteroid, one can simulate several cases. Case on Figure \ref{Hyncik1} (\href{run:Hyncik1.avi}{see movie}) describes the situation where the asteroids is caught by the Sun and it starts to revolute with decreasing semimajor axis. Let it run further, the asteroid would collide to the Sun. The axis are written in astronomical units where $1AU\approx 1.5\cdot10^8km$. Case on Figure \ref{Hyncik2} (\href{run:Hyncik2.avi}{see movie}) describes the situation where the asteroids crashes into the Sun. The average time of simulation describing two thousand years on regular PC is about ten till twenty seconds seconds.

\vspace{3mm}\noindent{\bf Conclusion}

\vspace{3mm}\noindent\rm Using MATLAB and its inbuilt functions, we have simply numerically implemented the equation of motion of a mass point in the force field. The force filed is the gravity field generated by the bodies of our Solar system with prescribed motion. Changing the initial conditions such as position in space and velocity, one can simulate the trajectory including collision detection for asteroid of a given mass.\\

\noindent{\bf References}\vspace{-15mm}\begin{thebibliography}{99}
\bibitem{bal}Balda, M.: \'Uvod do MATLABu, vydavatelstv\'{\i} Z\v{C}U, 1996
\bibitem{kle}Kleczek, J.: Velk\'a Encyklopedie vesm\'{\i}ru, Academia, 2002, ISBN 80-200-0906-X
\bibitem{mat}MATLAB 7 User's Guide
\bibitem{roc}P\v{r}\'{\i}hoda, P. et al.: Hv\v{e}zd\'a\v{r}sk\'a ro\v{e}enka 20XX, Hv\v{e}zd\'arna a planet\'arium hl.m. Prahy a A\'U AV \v{C}R, Praha 20XX, ISBN 80-902487-8-0
\bibitem{van}Van\'ysek, V.: Z\'aklady astronomie a astrofyziky, Academia, Praha, 1980
\end{thebibliography}

\vspace{2mm}

\noindent\it\textbf{Contacts:} Ing. Lud\v{e}k Hyn\v{c}{\'\i}k, Ph.D., University of West Bohemia in Pilsen, Departments of Mechanics, Univerzitn{\'\i} 8, 30614 Plze\v{n}, Czech Republic, tel.:~\mbox{+420377634709}, fax:~\mbox{+420377634702}, e-mail: \href{mailto:hyncik3@ntc.zcu.cz}{hyncik3@ntc.zcu.cz}, Mgr. V\'aclav Uhl\'{\i}\v{r}, Hv\v{e}zd\'arna Franti\v{s}ka Pe\v{s}ty, P.O. Box 48, 391 02 Sezimovo �st�, tel.:~\mbox{+420605557731}, fax: ~\mbox{+420381261037}, e-mail: \href{mailto:nosuh@seznam.cz}{nosuh@seznam.cz}\\

\noindent\it\textbf{Acknowledgement:} The work is supported by the Ministry of Education of the Czech Republic, project MSM 4977751303.\\

\end{document}
