% -*-latex-*-
%\nonstopmode
\documentclass[a4paper,11pt]{article}
\usepackage{a4}
\usepackage{amssymb}
\usepackage{amsmath}
% \usepackage{bbm}
\usepackage[latin1]{inputenc}
\usepackage[francais]{babel}
\usepackage[colorlinks,bookmarks=false,linkcolor=blue,urlcolor=blue]{hyperref}
\usepackage{color}
\usepackage{multicol}
\usepackage{graphicx}
\usepackage{caption}

\setlength{\textheight}{240mm}
\setlength{\topmargin}{-2.06cm}
%\setlength{\footskip}{5mm}
\setlength{\textwidth}{16.0cm}
\setlength{\oddsidemargin}{0.2cm}
\setlength{\evensidemargin}{0.2cm}

\pagestyle{plain}
%
\newcommand{\mail}[1]{{\href{mailto:#1}{#1}}}
\newcommand{\ftplink}[1]{{\href{ftp://#1}{ftp://#1}}}
\newcommand{\httplink}[1]{{\href{http://#1}{http://#1}}}
%

% Eigene Enumerate Umgebungen
\newenvironment{enumerate1}{\begin{list} {(\arabic{enumi})}%
{\usecounter{enumi}%
\setlength{\topsep}{0mm}%
\setlength{\partopsep}{0mm}%
\setlength{\itemsep}{0mm}%
\setlength{\labelsep}{2mm}%
\setlength{\labelwidth}{4mm}%
\setlength{\leftmargin}{2mm}%
\addtolength{\leftmargin}{\labelwidth}%
\addtolength{\leftmargin}{\labelsep}%
\setlength{\itemindent}{0mm}%
}}{\end{list}}

\newenvironment{enumeratei}{\begin{list} {(\roman{enumi})}%
{\usecounter{enumi}%
\setlength{\topsep}{0mm}%
\setlength{\partopsep}{0mm}%
\setlength{\itemsep}{0mm}%
\setlength{\labelsep}{2mm}%
\settowidth{\labelwidth}{(viii)}%
\setlength{\leftmargin}{1mm}%
\addtolength{\leftmargin}{\labelwidth}%
\addtolength{\leftmargin}{\labelsep}%
\setlength{\itemindent}{0mm}%
}}{\end{list}}

\newenvironment{enumeratea}{\begin{list}{(\alph{enumi})}%
{\usecounter{enumi}%
\setlength{\topsep}{0mm}%
\setlength{\partopsep}{0mm}%
\setlength{\itemsep}{0mm}%
\setlength{\labelsep}{2mm}%
\setlength{\labelwidth}{5mm}%
\setlength{\leftmargin}{0mm}%
\addtolength{\leftmargin}{\labelwidth}%
\addtolength{\leftmargin}{\labelsep}%
\setlength{\itemindent}{0mm}%
}}{\end{list}}

\def \be {\begin{equation}}
\def \ee {\end{equation}}
\def \dd   {{\rm d}}
\def \dt {\Delta t}
\def \F {\mathbf{F}}
\def \y {\mathbf{y}}
\def \B {\mathbf{B}}
\def \F {\mathbf{F}}
\def \v {\mathbf{v}}
\def \x {\mathbf{x}}
\def \J {\mathbf{J}}
\def \dd   {{\rm d}}
\def \w   {\omega}


\begin{document}
\noindent\textsc{\'Ecole Polytechnique F\'ed\'erale de Lausanne}\\
{\hfill Semestre d'automne 2019\\
Prof. Laurent Villard, Dr Cristian Sommariva
}
\vspace{5mm}
\begin{center}
{\Large\textbf{Physique Num\'erique I -- Exercice 4}}\\
\vspace{5mm}
{\small \`a rendre jusqu'au {\color{red} 18 D\'ecembre 2020 \`a 23h55} sur le site}
%\httplink{https://moodle.epfl.ch/mod/assign/view.php?id=835970}}
\href{https://moodle.epfl.ch/mod/assign/view.php?id=1108038}{https://moodle.epfl.ch/mod/assign/view.php?id=1108038}
\end{center}

\vspace{5mm}
\setcounter{section}{3}
\section{Syst\`eme gravitationel \`a 1,2 et 3 corps. Sonde solaire Parker. Tatooine, planète aux deux soleils}

\begin{figure}[!h]
\begin{center}
\includegraphics[width=7.8cm]{Parker}\hspace{0.2cm}
\includegraphics[width=7.8cm]{TatooineSunset2}
\caption{A gauche: sonde solaire Parker {\small \textcopyright NASA}. A droite: couchers de soleils sur Tatooine \textcopyright Lucas Skywalker.}
\end{center}
\end{figure}

%\begin{tabular}{lr}
%\begin{minipage}{0.5\linewidth}
\noindent 
On s'int\'eresse au mouvement de $N$ ($N\leq 3$) corps soumis aux
forces de gravitation. La force exerc\'ee sur le corps
$i$ par le corps $j$ s'\'ecrit:
\begin{equation}
\vec{F}_{ij}=- \frac{G\ m_i\ m_j}{|\vec{r}_j-\vec{r}_i|^3}\ \left( \vec{r}_i-\vec{r}_j\right)
\end{equation}
On utilise $G = 6.674 \times 10^{-11} {\rm [m^3kg^{-1}s^{-2}]}$ comme valeur de la
constante de gravitation.


%Dans un deuxi\`eme temps, on tiendra compte de la force de tra\^{\i}n\'ee a\'erodynamique caus\'ee par la pr\'esence de l'atmosph\`ere terrestre:
%\be
%\vec{F}_{t, ij}=-\frac{1}{2}\rho(|\vec{r}_i-\vec{r}_j|)SC_x|\vec{v}_i-\vec{v_j}|(\vec{v}_i-\vec{v_j})
%\ee
%o\`u
%\be
%\rho(r)=\rho_0\exp(-(r-R_T)/\lambda)
%\ee
%est la densit\'e de l'air, $\rho_0$ sa valeur au niveau de la mer, $R_T$ le rayon de la Terre, $\lambda$ une \'epaisseur caract\'eristique, $S$ la surface de section de l'objet, $C_x$ le coefficient de tra\^{\i}n\'ee.
%\end{minipage}
%\hspace{2mm}
%\begin{minipage}{0.5\linewidth}
%\includegraphics[width=6cm]{TatooineSunset2}
%\captionof{figure}{Couchers de soleils sur Tatooine}
%\end{minipage}
%\end{tabular}


\subsection{Num\'erique [5 pts]}
\noindent
Le but de cette excercice est d'introduire un sch\'ema num\'erique d'ordre \'elev\'e et un algorithme \`a pas d'int\'egration automatiquement ajustable. On \'ecrira les \'equations du mouvement en coordonn\'ees cart\'esiennes. Vous pouvez vous limiter au cas $N=3$ corps et au mouvement en deux dimensions d'espace $(x,y)$. On mettra le syst\`eme d'\'equations diff\'erentielles du mouvement sous la forme
\begin{equation}
  \label{eqn1}
  \frac{d \mathbf{y} }{d t} = \mathbf{f}(\mathbf{y}),
\end{equation}
avec le vecteur $\mathbf{y}$ qui comprend les positions \emph{et} les vitesses de
\emph{tous les corps}. Vous \'ecrirez un programme \verb;C++; utilisant l'algorithme de
Runge-Kutta du 4\`eme ordre, Eq.(2.123) des notes de cours, et un sch\'ema \`a pas
de temps adaptatif (Notes de cours section 2.5.2) \textbf{vous-m\^emes}.  
{\em On fera en sorte que le dernier pas de temps finisse exactement \`a un temps $t=t_\mathrm{fin}$ donn\'e.} 

\noindent On vous propose quelques
indications de programmation \`a la section \ref{SEC:GUIDE}. Vous pouvez \'egalement r\'eutiliser, en les modifiant, certaines parties des codes des exercices pr\'ec\'edents, ainsi que les \verb;ConfigFile.h; et \verb;ConfigFile.hpp;. On \'ecrira et
programmera aussi les expressions de l'\'energie m\'ecanique et de la quantit\'e
de mouvement du syst\`eme.

\subsection{Orbite de la sonde Parker [14 pts]}
\noindent
La sonde Parker, de masse $m=685$ kg, se trouve \`a une distance de $r_\mathrm{max}=120$ millions de km du centre du soleil. On aimerait qu'elle suive une orbite telle qu'elle s'approche \`a une distance minimale du centre du soleil de $r_\mathrm{min}=13.464378$ millions de km, voir Figure 2. (N.B.: C'est ce qui s'est pass\'e le 27 septembre 2020.)
Dans cette premi\`ere partie, on ne consid\'erera que la masse du Soleil, $m_S=1.989\times 10^{30}$kg, en n\'egligeant les masses des autres plan\`etes. {\em Indications: Pour les simulations num\'eriques, cela revient \`a prendre une masse du 3e corps,  tr\`es
petite, i.e. de l'ordre de $10^{-36}$ fois celle du soleil. On se placera dans le r\'ef\'erentiel h\'eliocentrique. }
\begin{figure}[!h]
\begin{center}
\includegraphics[width=7.5cm]{Fig2}
\caption{Orbite de la sonde Parker.}
\end{center}
\end{figure}
\begin{enumeratea}
\item
{\bf Partie analytique.}
Calculer les composantes du vecteur vitesse de la sonde Parker en $r=r_\mathrm{max}$ pour qu'elle suive l'orbite d\'ecrite ci-dessus.
Calculer la vitesse de la sonde, $v_\mathrm{max}$,  lorqu'elle sera au point le plus proche du soleil, en $r=r_\mathrm{min}$ . 
Calculer la p\'eriode de r\'evolution de la sonde, $T$.
\item
{\bf Pas de temps fixe.} 
R\'esoudre num\'eriquement la trajectoire de la sonde avec un pas de temps $\Delta t$ fixe, sur une dur\'ee de une p\'eriode de r\'evolution, $t_\mathrm{fin}=T$.
Faire des \'etudes de convergence de la position finale de la sonde, de la vitesse finale de la sonde et de la conservation de l'\'energie m\'ecanique $(\max_t |E_\mathrm{mec}(t)-E_\mathrm(0)|)$,
.%de la distance minimale entre la fus\'ee et le Soleil, $\min(|\vec{r}|)$, et de la vitesse maximale de la fus\'ee $\max(|\vec{v}|)$, 
en variant le nombre de pas de temps \verb;nsteps;. 
%{\em Indication: faire des interpolations paraboliques au voisinage de $\min(|\vec{r}|)$ et de $\max(|\vec{v}|)$ pour trouver ces extrema.}
\item
{\bf Pas de temps adaptatif.} 
En utilisant le pas de temps adaptatif, toujours pour $t_\mathrm{fin}=T$,  faire des \'etudes
de convergence de la position finale de la sonde, de la vitesse finale de la sonde et de la conservation de l'\'energie m\'ecanique $(\max_t |E_\mathrm{mec}(t)-E_\mathrm(0)|)$, 
%de la distance minimale entre la fus\'ee et le Soleil, $\min(|\vec{r}|)$, et de la vitesse maximale de la fus\'ee $\max(|\vec{v}|)$, 
en variant $\epsilon$, qui est  la pr\'ecision requise par pas de temps. On repr\'esentera l'erreur en fonction du nombre de pas de temps \verb;nsteps; effectu\'e.
Illuster comment le pas de temps $\Delta t$ change au cours de la simulation.
Comparer les \'etudes de convergence obtenues \`a pas de temps fixe et \`a pas de temps adaptatif. En particulier, comparer le nombre de pas de temps n\'ecessaire pour obtenir un r\'esultat sur la position finale de pr\'ecision
donn\'ee.
\end{enumeratea}

\subsection{Syst\`eme d'\'etoiles doubles [7 pts]} %==========================================================
\noindent 
On consid\`ere maintenant un syst\`eme d'\'etoiles doubles. La premi\`ere \'etoile est identique au soleil,  $m_1=1.989\times 10^{30}$kg. La deuxi\`eme \'etoile a une masse deux fois plus petite que le soleil, $m_2=m_1/2$. On n\'eglige la masse des autres corps c\'elestes.
\begin{enumeratea}
\item
On cherchera analytiquement des trajectoires circulaires uniformes dans
le {\rm r\'ef\'erentiel du centre de masse}, not\'e $\cal{R}_{CM}$, avec une distance entre les deux soleils $d=50$ millions de km. Il s'agit de trouver analytiquement les rayons des trajectoires, ainsi que  la vitesse angulaire $\Omega$ du mouvement.
\item
Effectuer des simulations avec des conditions initiales compatibles avec le calcul analytique ci-dessus et v\'erifier que les propri\'et\'es physiques du mouvement sont bien satisfaites: conservation de l'\'energie m\'ecanique, conservation de la quantit\'e de mouvement et distance entre les deux soleils  $d_{12}(t)=d=50$ millons de km  = constante.
{\it Indication: simuler une dizaine de p\'eriodes de r\'evolution.}
\end{enumeratea}

\subsection{Orbite de Tatooine dans le syst\`eme \`a deux soleils [9 pts]}
\noindent On consid\`ere le m\^eme syst\`eme d'\'etoiles doubles que pr\'ec\'edemment, auquel on rajoute la plan\`ete Tatooine, de masse $m=6\times 10^{24}$kg. On reste dans le r\'ef\'erentiel du centre de masse du syst\`eme double,  $\cal{R}_{CM}$.
\begin{figure}
\begin{center}
\includegraphics[width=7.5cm]{Fig3}
\caption{Condition initiale pour le système de Tatooine.}
\end{center}
\end{figure}

\begin{enumeratea}
\item
Placer la condition initiale des trois corps  comme indiqu\'e sur la Figure 3. On prendra la position initiale de Tatooine sur l'axe des $y$, $y<0$, \`a une distance $r_{3,0}$ du centre de masse des deux soleils. On prendra la vitesse initiale de Tatooine telle qu'elle aurait une trajectoire circulaire si, au lieu des deux soleils de masses $m_1$ et $m_2$, il n'y avait qu'un soleil de masse $M=m_1+m_2$ au centre de masse.
\item
Effectuer une s\'erie de simulations avec le sch\'ema adaptatif, en variant la distance $r_{3,0}$. Faire des simulations d'une centaine de p\'eriodes de r\'evolution du syst\`eme d'\'etoiles doubles. On mesurera la distance de Tatooine au centre de masse, $r_3(t)$. (Pour que la vie sur Tatooine soit possible, il faut que cette distance varie le moins possible). Illuster des cas d'orbites avec peu de variation de $r_3(t)$, des cas avec une variation de l'ordre de 30\%, et des cas o\`u Tatooine est \'eject\'e du syst\'eme solaire double. {\em Indication: le domaine de position initiale int\'eressant \`a explorer est   $r_{3,0}$ entre $0.3d$ et $5d$.}

%\item
%En rajoutant la pr\'esence de l'atmosph\`ere $(\rho_0=1.2\mathrm{kg/m^3})$, ajuster de nouveau la direction
%initiale de Apollo13 pour minimiser l'acc\'el\'eration maximale comme en 4.3(b).
\end{enumeratea}

\subsection{Points de Lagrange [10 pts]}
\noindent
On aimerait savoir s'il est possible de placer Tatooine sur  une orbite telle qu'elle
garde des distances constantes par rapport aux deux soleils.

\begin{enumeratea}
%% \item
%% Calculer les positions d'\'equilibre du 3e corps dans $\mathcal{R}^\prime$ sur l'axe
%% Terre-Lune. Il y en a trois : derri\`ere le Terre, entre la Terre et la
%% Lune et derri\`ere la Lune. Vous utiliserez la fonction \verb+roots+ de Matlab
%% pour r\'esoudre l'\'equation du 5\`eme degr\'e.
%% \item Placer le 3e corps sur l'un ces points avec $v_{\mathcal{R}^\prime}=0$ et v\'erifier qu'il garde bien des
%% distances constantes avec la Terre et la Lune.
%% \item Pour voir si cette position d'\'equilibre est stable, se placer \`a quelque
%% distance de ce point ($\approx1000$km) et observer la trajectoire obtenue.
%% % \vspace{0.1cm}

\begin{tabular}{lr}
\hspace{-0.35cm} \begin{minipage}{0.5\linewidth}
\item On consid\`ere maintenant un syst\`eme d'\'etoiles doubles avec $m_2=m_1/100$.Tatooine est plac\'ee en une position initiale $P_0$ formant un triangle
\'equilat\'eral avec les deux soleils, dit ``point de Lagrange L4'', voir Figure 4. Sa vitesse initiale est
$\vec{v}_3(0)=\Omega r \vec{e}_\theta$, o\`u 
$\Omega =$ vitesse angulaire du syst\`eme des deux soleils calcul\'ee en 4.3 (a), $(r,\theta) =$
coordonn\'ees polaires centr\'ees au $CM$ des deux soleils. 
Simuler une vingtaines de p\'eriodes de r\'evolution.  V\'erifier que la distance entre Tatooine et les deux soleils reste constante au cours du temps.
Placer ensuite Tatooine \`a quelque distance (100'000km ou plus) du point de
Lagrange L4 et repr\'esenter les trajectoires dans le r\'ef\'erentiel $\cal{R}^\prime$ tournant avec le syst\`eme des deux soleils, d'axe de rotation au centre de masse $(CM)$ du syst\`eme. 

\end{minipage} 
\begin{minipage}{0.5\linewidth}
$\mbox{}$\\
\centerline{\includegraphics[width=6cm]{Fig4}}
\captionof{figure}{Tatooine est plac\'ee au voisinage du point de Lagrange L4,  en $P_0$.}
\end{minipage}
\end{tabular}
\item On revient au cas initial du syst\`eme avec $m_2=m_1/2$. 
 En initialisant Tatooine au point de la grange L4, en $P_0$, simuler une dur\'ee de $t_\mathrm{fin}\approx 20$ p\'eriodes de r\'evolution.  Observer comment la distance entre Tatooine et les deux soleils \'evolue au cours du temps. Repr\'esenter les trajectoires dans le r\'ef\'erentiel $\cal{R}^\prime$.

\end{enumeratea}

\subsection{Facultatif  [5pts]}
\noindent
Le nombre de probl\`emes int\'eressants \`a trois corps est immense... faites fonctionner votre imagination, ou cherchez, sur le web ou ailleurs, des id\'ees \`a poursuivre. Ci-dessous quelques pistes possibles.
\begin{enumeratea}
\item
{\bf Plan\`ete X et autres.} Pour continuer la Section 4.5, trouvez les autres points de Lagrange du syst\`eme des deux soleils et \'etudier les trajectoires obtenues aux voisinages de ces points de Lagrange. Observer la stabilit\'e (ou l'instabilit\'e) de ces positions d\'equilibre de $\mathcal{R}^\prime$.
\item
Programmer le sch\'ema de {\bf velocity-Verlet}, Eq.(2.98) des notes de cours, tout en gardant un algorithme \`a pas de temps adaptatif. Comparer les propri\'et\'es de convergence num\'erique et de conservation de l'\'energie avec le sch\'ema Runge-Kutta d'ordre 4. Tenter de {\em tr\`es longues} simulations avec ces deux sch\'emas.
%Bring the boys back home - essayez de trouver la bonne vitesse
%initiale d'Apollo 13 pour les ramener vivants sur Terre, en tenant compte \`a la fois de la Lune et de
%l'atmosph\`ere terrestre (minimiser les valeurs maximales de l'acc\'el\'eration et de la puissance dissip\'ee,
%tout en \'evitant qu'Appollo13 ne rebondisse sur l'atmosph\`ere).
%\item trouver les positions des autres points de Lagrange, et examiner la stabilit\'e des orbites
%\`a leur voisinage.
\end{enumeratea}

\subsection{Soumission du rapport et du code C++ [5 pts de qualit\'e g\'en\'erale]}
\begin{enumeratea}
\item{Pr\'eparer le fichier du rapport en format \verb+pdf+ portant le nom
\verb+RapportExercice4_Nom1_nom2.pdf+, ainsi que le fichier source {\LaTeX}
\verb+RapportExercice4_Nom1_nom2.tex+.
}
\item{Pr\'eparer le fichier source C++ \verb+Exercice4_Nom1_nom2.cpp+.}
\item{Pr\'eparer le fichier script Matlab \verb+Exercice4_Nom1_nom2.m+.}
\item{Le lien de soumission se trouve dans la bo\^ite du premier exercice
sur notre site "Moodle"  et peut directement \^etre atteint en cliquant
\href{https://moodle.epfl.ch/mod/assign/view.php?id=835970}{ici}.}
\end{enumeratea}
%---------------------------------------------------------------------------------------------------------

\subsection{Guide de programmation}
\label{SEC:GUIDE}
Cet exercice vous offre une exp\'erience plus riche que les autres,
parce que vous allez \'ecrire votre code vous-m\^eme. Afin d'en
faciliter l'\'ecriture, veuillez tenir compte des points suivants qui
pourront vous \^etre utiles.

\subsubsection{Conventions physiques}

On \'ecrira les \'equations du mouvement sous la forme
\begin{equation}
  \label{eqn1_annexe}
  \frac{d \mathbf{y} }{d t} = \mathbf{f}(\mathbf{y})
\end{equation}
Ici, le vecteur $\mathbf{y}$ comprend les positions \emph{et} les
vitesses de \emph{tous les corps}.
Vous devez choisir une convention pour $\mathbf{y}$. Une possibilit\'e est donn\'ee par  $\mathbf{y} =
\left(\mathbf{r}_1,\ldots,\mathbf{r}_N,\mathbf{v}_1,\ldots,\mathbf{v}_N\right)$
avec $\mathbf{r}_i$ la position du corps $i$, et $\mathbf{v}_i$ sa vitesse.
Pour le cas \`a deux dimensions, cela devient par exemple ï¿½ d\'efinir
$\mathbf{y}=\left( x_1, y_1 , x_2, y_2,\ldots, v_{1;x}, v_{1;y} , \ldots\right)$.
\begin{center}
\textbf{Choisir une convention pour $\mathbf{y}$ et en d\'eduire la forme correspondante de $\mathbf{f}$!}
\end{center}
On programmera le sch\'ema de Runge-Kutta d'ordre 4
sous la forme donn\'ee par l'Eq.(2.97) des notes de cours.



%\begin{equation}
%\mathbf{y}_{t+1} = \mathbf{y}_t + \frac{\Delta t}{6}\left( \mathbf{k}_1 + 2\mathbf{k}_2 + 2\mathbf{k}_3  + \mathbf{k}_4\right),
%\end{equation}%
%where the $\mathbf{k}$-vectors are given by
%\begin{align}
%\mathbf{k}_1 &= f(\mathbf{y}_t)\\
%\mathbf{k}_2 &= f(\mathbf{y}_t + 1/2 \mathbf{k}_1)\\
%\mathbf{k}_3 &= f(\mathbf{y}_t + 1/2 \mathbf{k}_2)\\
%\mathbf{k}_4 &= f(\mathbf{y}_t + \mathbf{k}_3)
%\end{align}


\subsubsection{La structure g\'en\'erale}

On propose d'organiser l'architecture de l'outil num\'erique autour de trois
\'etapes principales : l'initialisation des conditions initiales, l'\'evolution
du syst\`eme physique et le post-traitement de la sortie.

Pour l'initialisation, vous pouvez reprendre le template \verb+ConfigFile+
de l'exercice 3. Vous pourrez ainsi adapter le script Matlab
\verb+ParameterScan.m+ pour varier un ou plusieurs param\`etres.

\begin{sloppypar}
On propose d'utiliser une structure similaire aux exercices pr\'ec\'edents,
avec une classe \texttt{Exercice4} dont le constructeur appelle \verb+ConfigFile+
pour lire les param\`etres, et dont la m\'ethode \verb+run()+ effectue
l'it\'eration temporelle en appelant la m\'ethode \verb+step()+ et l'\'ecriture
des r\'esultats \verb+printOut()+.
\end{sloppypar}

Dans la boucle temporelle de \verb+run()+, vous devrez inclure des modifications pour pouvoir ex\'ecuter soit avec un pas de temps $\Delta t$ fixe, soit avec le sch\'ema $\Delta t$ adaptatif.

\begin{sloppypar}
La m\'ethode \verb+step()+ pourra faire appel \`a une fonction \texttt{std::valarray acceleration(std::valarray const\& y)} qui calculera les composantes du vecteur $\mathbf{f}(\mathbf{y})$ (voir la sous-section suivante pour l'utilisation des \texttt{valarrays}). L'utilisation des r\'ef\'erences constantes permet d'\'eviter la copie du tableau \`a chaque appel de la fonction.
\end{sloppypar}

La m\'ethode \verb+printOut()+ peut utiliser une fr\'equence d'\'echantillonage pour \'ecrire tous les \texttt{nsampling} pas de temps afin d'\'eviter les fichiers de sortie trop volumineux.

Une fois que vous avez g\'en\'er\'e votre
fichier de sortie, il vous faut le lire au sein de fonctions ou de scripts Matlab afin d'en
faire ressortir l'information pertinente. Commencez avec une fonction simple
comme par exemple le tra\c{c}age de trajectoires, puis soyez cr\'eatifs et essayez
d'avoir des figures simples qui synth\'etisent votre compr\'ehension.

\subsubsection{Quelques propositions techniques et d'approche}

\begin{enumeratea}
\setlength{\itemsep}{0.1cm}
\item \textbf{Utiliser} \verb;valarrays; pour les vecteurs. Cela
  facilitera  le calcul. En effet, pour des vecteurs  $\mathbf{a}$, $\mathbf{b}$,
  $\mathbf{c}$ de la m\^eme dimension cod\'es sous forme de valarrays, il suffit d'\'ecrire
  \begin{center}
  \verb;c=a+b;
  \end{center}
  pour faire la somme sans avoir \`a \'ecrire une boucle sur les indices.

\item Le concept de \verb;slices; est extr\^emement utile pour acc\'eder \`a
  un groupe d'indices d'un vecteur. Par exemple, imaginons un vecteur
  $$\mathbf{m}=(m_0,m_1, \ldots , m_3)$$ dont on veut remplacer les
  composantes $m_2$ et $m_3$ par les composantes $n_0$ et $n_1$ du vecteur
  $\mathbf{n}=(n_0, n_1)$.

  Si $\mathbf{m}$ and $\mathbf{n}$ sont du type \verb;valarray;, en
  utilisant la commande slice, il suffit d'ï¿½crire

  \begin{center}
    \verb;m[std::slice(2,2,1)] = n;
  \end{center}
  L'op\'eration inverse
  \begin{center}
    \verb;n = m[std::slice(2,2,1)];
  \end{center}
  est aussi possible.

  En fait, la fonction \verb;slice(i,j,k); commence avec l'indice $i$ du vecteur et
  s\'electionne $j$ indices r\'eguli\`erement espac\'es de $k$ positions.

\item Il est indispensable que vous structuriez et planifiez votre
  code avant de vous lancer dans l'\'ecriture:
  \begin{itemize}
  \item {Essayer de faire un  dessin avec des liaisons, les t\^aches et
      les variables principales.}
  \item {Discuter en groupe votre mod\`ele}.
  \end{itemize}

\item Utiliser les documentations pour \verb;C++; (par exemple \texttt{www.cplusplus.com})

% \item Mettre un ``flag'', valarray d'entiers qui valent 0 ou 1, avec 0 pour indiquer si on veut forcer le
% corps correspondant ï¿½ rester fixe. On multipliera les termes du vecteur de fonctions
% $\mathbf{f}(\mathbf{y})$ par le flag correspondant.

\end{enumeratea}
\end{document}


\subsection{Guide de programmation} \label{SEC:GUIDE}
Cet exercice vous offre une exp\'erience plus riche que les autres,
parce que vous allez \'ecrire votre code vous-m\^eme. Afin d'en
faciliter l'\'ecriture, veuillez tenir compte des points suivants qui
pourront vous \^etre utiles.

\subsubsection{Conventions physiques}

Le but de cette excercice est d'\'etudier le probl\`eme \`a $N$ corps en
utilisant les \'equations du mouvement sous la forme
\begin{equation}
  \label{eqn1}
  \frac{d \mathbf{y} }{d t} = \mathbf{f}(\mathbf{y})
\end{equation}
Ici, le vecteur $\mathbf{y}$ comprend les positions \emph{et} les
vitesses de \emph{tous les corps}.
Vous devez choisir une convention pour $\mathbf{y}$ ainsi que pour le
vecteur $\mathbf{f}(\mathbf{y})$. Une possibilit\'e est donn\'ee par  $\mathbf{y} =
\left(\mathbf{r}_1,\ldots,\mathbf{r}_N,\mathbf{v}_1,\ldots,\mathbf{v}_N\right)$
avec $\mathbf{r}_i$ la position du corps $i$, et $\mathbf{v}_i$ sa vitesse.
Pour le cas \`a deux dimensions, cela devient par exemple \`a d\'efinir
$\mathbf{y}=\left( x_1, y_1 , x_2, y_2,\ldots, v_{1;x}, v_{1;y} , \ldots\right)$.
\begin{center}
\textbf{Choisir une convention et en d\'eduire la forme correspondante de $\mathbf{f}$!}
\end{center}
En utilisant  \eqref{eqn1},  on peut programmer le sch\'ema de Runge-Kutta
sous la forme donn\'ee par l'Eq.(2.85) des notes de cours.



%\begin{equation}
%\mathbf{y}_{t+1} = \mathbf{y}_t + \frac{\Delta t}{6}\left( \mathbf{k}_1 + 2\mathbf{k}_2 + 2\mathbf{k}_3  + \mathbf{k}_4\right),
%\end{equation}%
%where the $\mathbf{k}$-vectors are given by
%\begin{align}
%\mathbf{k}_1 &= f(\mathbf{y}_t)\\
%\mathbf{k}_2 &= f(\mathbf{y}_t + 1/2 \mathbf{k}_1)\\
%\mathbf{k}_3 &= f(\mathbf{y}_t + 1/2 \mathbf{k}_2)\\
%\mathbf{k}_4 &= f(\mathbf{y}_t + \mathbf{k}_3)
%\end{align}


\subsubsection{La structure g\'en\'erale}

On propose d'organiser l'architecture de l'outil num\'erique autours de trois
\'etapes principales : initialisation des conditions initiales, l'\'evolution
du syst\`eme physique et le traitement post-diagnostic de la sortie.

Pour l'initialisation, on vous propose de cr\'eer un script matlab
qui g\'en\`ere une batterie de conditions initiales. Voici une fonction Matlab permettant
de cr\'eer de tels fichiers
\small
\begin{verbatim}
function printFile(folder,corpsIndex,r,v,m,G,precision,nsampling,dt,dtMax,tFin)

  corpsfile=sprintf('%s/corps_%d.in',folder,corpsIndex);
  fid = fopen(corpsfile,'w');
  nbCorps=size(r,1);
  for iCorps=1:nbCorps
    fprintf(fid,'%g %g %g %g %g\n',r(iCorps,1),r(iCorps,2),
                                   v(iCorps,1),v(iCorps,2),
                                   m(iCorps));
  end
  fclose(fid);

  configfile=sprintf('%s/configuration_%d.in',folder,corpsIndex);
  outputfile=sprintf('%s/output_%d.dat',folder,corpsIndex);
  fid=fopen(configfile,'w');
  fprintf(fid,'%g\n',G);
  fprintf(fid,'%d\n',precision);
  fprintf(fid,'%d\n',nsampling);
  fprintf(fid,'%g\n',dt*3600*24*365);
  fprintf(fid,'%g\n',dtMax*3600*24*365);
  fprintf(fid,'%g\n',tFin*3600*24*365);
  fprintf(fid,'%d\n',nbCorps);
  fprintf(fid,'%s\n',corpsfile);
  fprintf(fid,'%s\n',outputfile);
  fclose(fid);
\end{verbatim}
\normalsize
Chaque syst\`eme de trois corps(Terre, Lune, asteroide)
aura un couple de fichier : le premier "configuration.in" contient les param\`etres
physiques g\'en\'eraux (G, dt, temps de fin,...) et le second "corps.in" contient
les coordonn\'ees et les vitesses au temps z\'ero.

On vous propose d'\'ecrire une classe Exercice4 sur le mod\`ele de l'exercice 3
qui lira le fichier de configuration dans le constructeur puis les conditions initiales
dans la m\'ethode init() dont voici un exemple
\small
\begin{verbatim}
  void init(){
    last=0;
    t=0;
    // Lecture du fichier de positions des corps
    FILE *pf1;
    pf1=fopen(systemPath.c_str(),"r");
    int ix,iy;
    for (int pos=0;pos<nb;pos++){
      ix=2*pos;
      iy=ix+1;
      fscanf(pf1,"%lf %lf %lf %lf %lf\n",&r[ix],&r[iy],&v[ix],&v[iy],&m[pos]);
      cout << m[pos] << " " << r[ix] << " " << r[iy] ;
      cout << " " << v[ix] << " " << v[iy] <<endl;
    }
    fclose(pf1);
  }
\end{verbatim}
\normalsize
Il vous est sugg\'r\'e de stocker les donn\'ees (position, vitesse et masse) au
sein de tableau std::valarray. La dynamique du syst\`eme peut \^etre d\'ecrite \`a
l'aide de trois m\'ethodes: run(), step(...) et acc\'el\'eration(...).
Vous nous n'aurez cependant plus besoin d'utiliser de m\'ethode virtuelle, ce qui
devrait simplifier le code. Concernant la performance, si vous avez besoin de mettre
\`a jour un tableau au sein d'une m\'ethode et que vous passez ce tableau \`a l'appel
de la m\'ethode, pensez \`a utiliser les r\'ef\'erence sur objet par utilisation du
"\&", ceci afin d'\'eviter les recopies qui peuvent \^etre contre productives.

Concernant les donn\'ees de sortie, elles vont devenir tr\`es volumineuses si vous n'en
prenez pas soin. On vous sugg\`ere donc de d\'efinir une fr\'equence
d'\'echantillonage afin de na pas sauvegarder les positions \`a chaque pas de temps mais
d'ajuster au besoin.

Maintenant que vous avez g\'en\'er\'e votre
fichier de sortie, il vous faut le lire au sein de fonctions Matlab afin d'en
faire ressortir l'information pertinante. Commencez avec une fonction simple
comme par exemple le tra\c cage de trajectoires puis soyez cr\'eatifs et essayez
d'avoir des figure simples qui synth\'etisent votre compr\'ehension.

\subsubsection{Quelques propositions techniques et d'approche}

\begin{enumeratea}
\setlength{\itemsep}{0.1cm}
\item \textbf{utiliser} \verb;valarrays; pour les vecteurs. Cela
  facilitera  le calcul. En effet, pour des vecteurs  $\mathbf{a}$, $\mathbf{b}$,
  $\mathbf{c}$ de la m\^eme dimension cod\'es sous forme de valarrays, il suffit d'\'ecrire
  \begin{center}
  \verb;c=a+b;
  \end{center}
  pour faire la somme sans avoir \`a \'ecrire une boucle sur les indices.

\item le concept de \verb;slices; est extr\^emement utile pour acc\'eder \`a
  un groupe d'indices d'un vecteur. Par exemple, imaginons un vecteur
  $$\mathbf{m}=(m_0,m_1, \ldots , m_3)$$ dont on veut remplacer les
  composantes $m_2$ et $m_3$ par les composantes $n_0$ et $n_1$ du vecteur
  $\mathbf{n}=(n_0, n_1)$.

  Si $\mathbf{m}$ and $\mathbf{n}$ sont du type \verb;valarray;, en
  utilisant la commande slice, il suffit d'\'ecrire

  \begin{center}
    \verb;m[std::slice(2,2,1)] = n;
  \end{center}
  L'op\'eration inverse
  \begin{center}
    \verb;n = m[std::slice(2,2,1)];
  \end{center}
  est aussi possible.

  En fait, la fonction \verb;slice(i,j,k); commence avec l'indice $i$ du vecteur et
  s\'electionne $j$ indices r\'eguli\`erement espac\'es de $k$ positions.

\item Il est indispensable que vous structuriez et planifiez votre
  code avant de vous lancer dans l'\'ecriture:
  \begin{itemize}
  \item \textbf{Essayer de faire un  dessin avec des liaisons, les t\^aches et
      les variables principales.}
  \item \textbf{Discuter en groupe votre mod\`ele}.
  \end{itemize}

\item utiliser les documentations pour \verb;C++; (par example \texttt{www.cplusplus.com})

\item Mettre un ``flag'', valarray d'entiers qui valent 0 ou 1, avec 0 pour indiquer si on veut forcer le
corps correspondant \`a rester fixe. On multipliera les termes du vecteur de fonctions
$\mathbf{f}(\mathbf{y})$ par le flag correspondant.

\end{enumeratea}
\end{document}
