
\documentclass{article}
\usepackage[a4paper,margin=2cm]{geometry}
\usepackage[T1]{fontenc}
\usepackage[latin1]{inputenc}%
\usepackage[garamond]{mathdesign}
\usepackage{pst-eqdf,pst-node,pst-tools}
\usepackage{array,amsmath}
\newpsstyle{vecteurA}{arrowinset=0.05,arrowsize=0.1,linecolor={[rgb]{1 0.5 0}}}
\newpsstyle{vecteurB}{arrowinset=0.05,arrowsize=0.1,linecolor={[rgb]{0 0.5 1}}}
\newpsstyle{vecteurC}{arrowinset=0.1,arrowsize=0.2,linecolor={[rgb]{1 0 0}}}
%%%%%%%%%%%%%%%%%%
\title{Gravitation : le problème des deux corps avec PSTricks\\ partie 1}
\date{3 juillet 2\,012}
\begin{document}
\maketitle
\section{Présentation}
 Cette première partie aborde uniquement le problème théorique et l'établissement des relations et formules indispensables à la réalisation des figures avec `PSTricks' et à l'animation avec le package `\textsf{animate}'. Elle pourra donc sembler superflue à tous ceux qui connaissent bien ce problème classique. De très bons livres traitent le \textit{problème des deux corps} de façon très claire comme celui de José-Philippe Pérez dans `\textit{Mécanique}' aux éditions Masson ou de façon très complète(et très claire aussi) comme celui publié par le \textsc{cnes} aux éditions \textsc{cepadues} : `\textit{Le mouvement du satellite}'.

 La deuxième partie détaillera la procédure suivie pour réaliser les schémas et les animations. Toutefois, on peut, dès à présent, avoir accès au code des figures  dans le fichier source de ce document.
\section{Étude théorique}
\newcommand\GravAlg{%
  y[2]|y[3]|%
  M1*(y[4]-y[0])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5|%
  M1*(y[5]-y[1])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5|%
  y[6]|y[7]|%
  M2*(y[0]-y[4])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5|%
  M2*(y[1]-y[5])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5}
%%  0  1   2   3  4  5   6   7
%% x1 y1 x'1 y'1 x2 y2 x'2 y'2
On considère un système de deux corps en interaction gravitationnelle $M_1$ de masse $m_1$ et $M_2$ de masse $m_2$ dans le repère galiléen \textit{inertiel} $\mathcal{R}$. Ils sont supposés ponctuels.
\begin{center}
\begin{pspicture}(-3,-1)(5,6)
\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]%
\pstVerb{/x01 -2 def /y01 1 def
         /x02  5 def /y02 4 def
         /xr0 x02 x01 sub def
         /yr0 y02 y01 sub def
         /M1 3 def
         /M2 1 def
         /Mt M1 M2 add def
         /xG0 x01 M1 mul x02 M2 mul add Mt div def
         /yG0 y01 M1 mul y02 M2 mul add Mt div def}%
\pnode(0,0){O}
\pnode(!x01 y01){M1}
\pnode(!x02 y02){M2}
\pnode(!xG0 yG0){C}
\pscircle[fillstyle=solid,fillcolor=red!50](M1){0.21}
\pscircle[fillstyle=solid,fillcolor=blue!50](M2){0.07}
\psdot(C)
\psline[linestyle=dotted](M1)(M2)
\psline[linestyle=dashed](M2)(O)(M1)
\psline{<->}(5,0)(0,0)(0,6)
\uput[dl](O){$O$}
\uput[r](0,5.8){$y$}
\uput[u](4.8,0){$x$}
\uput[l](0,5.8){$\mathcal{R}$}
\uput{0.22}[l](M1){$M_1$}
\uput{0.1}[u](M2){$M_2$}
\uput[u](C){$C$}
\pcline[offset=5pt,linestyle=none](M1)(M2)
\ncput*[nrot=:U]{$r$}
\rput(M1){\psline[style=vecteurC]{->}(!xr0 5 div yr0 5 div)\uput[u](!xr0 10 div yr0 10 div){$\overrightarrow{F}_{2/1}$}}
\rput(M2){\psline[style=vecteurC]{->}(!xr0 5 div neg yr0 5 div neg)\uput[u](!xr0 10 div neg yr0 10 div neg){$\overrightarrow{F}_{1/2}$}}
\end{pspicture}
\end{center}
On note $\overrightarrow{r}=\overrightarrow{M_1M_2}$. $M_2$ subit de la part de $M_1$ une force attractive $\overrightarrow{F}_{1/2}$ et réciproquement $M_1$ subit de la part de $M_2$ une force attractive $\overrightarrow{F}_{2/1}$ telles que :
\[
\overrightarrow{F}_{1/2}=-\mathcal{G}\frac{m_1m_2}{r^3}\overrightarrow{r}\quad \overrightarrow{F}_{2/1}=+\mathcal{G}\frac{m_1m_2}{r^3}\overrightarrow{r}
\]
Bien sûr : $\overrightarrow{F}_{1/2}+\overrightarrow{F}_{2/1}=\overrightarrow{0}$. Par conséquent, le centre de masse $C$ du système $\{M_1,M_2\}$ est, dans le repère $\mathcal{R}$, soit immobile, soit en mouvement de translation rectiligne uniforme suivant les conditions initiales.

Dans $\mathcal{R}$ la position de $C$ est déterminée par la relation :
\[
\overrightarrow{OC}=\frac{m_1\overrightarrow{OM_1}+m_2\overrightarrow{OM_2}}{m_1+m_2}
\]
\begin{equation}
\overrightarrow{v_{C/\mathcal{R}}}=\frac{m_1\overrightarrow{v_1}+m_2\overrightarrow{v_2}}{m_1+m_2}=\overrightarrow{v_0}
\label{vc}
\end{equation}
$\overrightarrow{v_0}$ étant déterminé par les valeurs initiales de $\overrightarrow{v_1}$ et $\overrightarrow{v_2}$.
L'étude peut donc se poursuivre dans le repère lié au centre de masse $\mathcal{R}^*$, lui-même galiléen.

On pose : $\overrightarrow{r_1}=\overrightarrow{CM_1}$ et $\overrightarrow{r_2}=\overrightarrow{CM_2}$. On a toujours : $\overrightarrow{r}=\overrightarrow{M_1M_2}=\overrightarrow{r_2}-\overrightarrow{r_1}$.
\begin{center}
\begin{pspicture}(-3.5,-1)(5.5,8.25)
%\psframe(-3.5,-1)(5.5,9)
\pstVerb{/x01 -2 def /y01 1 def
         /x02  4 def /y02 5 def
         /xr0 x02 x01 sub def
         /yr0 y02 y01 sub def
         /M1 3 def
         /M2 1 def
         /Mt M1 M2 add def
         /xG0 x01 M1 mul x02 M2 mul add Mt div def
         /yG0 y01 M1 mul y02 M2 mul add Mt div def}%
\pnode(0,0){O}
\pnode(!x01 y01){M1}
\pnode(!x02 y02){M2}
\pnode(!xG0 yG0){C}
\pscircle[fillstyle=solid,fillcolor=red!50](M1){0.21}
\pscircle[fillstyle=solid,fillcolor=blue!50](M2){0.07}
\psdot(C)
\psline[linestyle=dotted](M1)(M2)
\psline[linestyle=dashed](M2)(O)(M1)
\psline[linecolor=gray]{<->}(5,0)(0,0)(0,6)
\rput(C){\psline[style=vecteurA]{<->}(5.5,0)(0,0)(0,6)%
         \uput[l](0,5.8){$\mathcal{R^*}$}
         \psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt](-3,-3)(6,6)%
         \uput[r](0,5.8){$y^*$}
         \uput[u](5.5,0){$x^*$}}
\uput[dl](O){$O$}
\uput[r](0,5.8){$y$}
\uput[u](4.8,0){$x$}
\uput[l](0,5.8){$\mathcal{R}$}
\uput{0.22}[l](M1){$M_1$}
\uput{0.1}[u](M2){$M_2$}
\uput[ul](C){$C$}
\pcline[offset=5pt,linestyle=none](M1)(M2)
\ncput*[nrot=:U]{$r$}
\rput(M1){\psline[style=vecteurC]{->}(!xr0 5 div yr0 5 div)\uput[u](!xr0 10 div yr0 10 div){$\overrightarrow{F}_{2/1}$}}
\rput(M2){\psline[style=vecteurC]{->}(!xr0 5 div neg yr0 5 div neg)\uput[u](!xr0 10 div neg yr0 10 div neg){$\overrightarrow{F}_{1/2}$}}
\end{pspicture}
\end{center}

Dans le repère $\mathcal{R}^*$, la quantité de mouvement du système $\{M_1,M_2\}$ est nulle : $\overrightarrow{p}^*=(m_1+m_2)\overrightarrow{v}_C^*=\overrightarrow{0}$.
Pour chacun des deux corps, la quantité de mouvement dans $\mathcal{R}^*$ est :
\[
\overrightarrow{p}_1^*=m_1\overrightarrow{v}_1^*=m_1(\overrightarrow{v}_{1/\mathcal{R}}-\overrightarrow{v}_{C/\mathcal{R}})=-\overrightarrow{p}^*_2
\]
D'après (\ref{vc}) qui donne  $\overrightarrow{v}_{C/\mathcal{R}}$, on obtient :
\[
\overrightarrow{p}_1^*=m_1\left(\overrightarrow{v}_{1/\mathcal{R}}-\frac{m_1\overrightarrow{v}_{1/\mathcal{R}}+m_2\overrightarrow{v}_{2/\mathcal{R}}}{m_1+m_2}\right)=%
 \frac{m_1m_2}{m_1+m_2}(\overrightarrow{v}_{1/\mathcal{R}}-\overrightarrow{v}_{2/\mathcal{R}})=-\overrightarrow{p}^*_2
\]
Puisque, on a défini $\overrightarrow{r}$ par $\overrightarrow{r_2}-\overrightarrow{r_1}$, en posant $\mu\dfrac{m_1m_2}{m_1+m_2}$ :
\[
\overrightarrow{p}^*_1=-\mu\dfrac{\mathrm{d}\overrightarrow{r}}{\mathrm{d}t}\qquad \overrightarrow{p}^*_2=\mu\dfrac{\mathrm{d}\overrightarrow{r}}{\mathrm{d}t}
\]
La quantité de mouvement de $M_2$ est identique à celle d'une particule fictive de masse $\mu$, appelée \textit{masse réduite}, et de vitesse $\overrightarrow{v}_{2/\mathcal{R}}-\overrightarrow{v}_{1/\mathcal{R}}$ : vitesse relative de $M_2$ par rapport à $M_1$, idem pour $M_1$.

On considère cette particule fictive de masse $\mu$ située au point $M$ tel que : $\overrightarrow{CM}=\overrightarrow{r}$.

La relation fondamentale de la dynamique appliquée à chacun des deux corps donne :
\[
\left\{
\begin{array}[m]{l}
 \dfrac{\mathrm{d}\overrightarrow{p}^*_1}{\mathrm{d}t}=\overrightarrow{F}_{2/1}\\[1em]
\dfrac{\mathrm{d}\overrightarrow{p}^*_2}{\mathrm{d}t}=\overrightarrow{F}_{1/2}
\end{array}
\right.
\]
En remplaçant $\overrightarrow{p}^*_1$ (ou $\overrightarrow{p}^*_2$) par son expression en fonction de $\mu$, on a :
\begin{equation}
 \mu\dfrac{\mathrm{d}^2\overrightarrow{r}}{\mathrm{d}t^2}=\overrightarrow{F}_{1/2}=-\mathcal{G}\frac{m_1m_2}{r^3}\overrightarrow{r}
\label{RFD}
\end{equation}
Nous désignons cette force par $\overrightarrow{F}$. Remarquons qu'elle peut s'écrire, en notant la masse totale $(m_1+m_2)$  du système $\{M_1,M_2\}$ par $m_t$~:
\begin{equation}
\overrightarrow{F}=-\mathcal{G}\frac{\mu m_t}{r^3}\overrightarrow{r}
\label{RF}
\end{equation}
En conclusion, cette particule fictive $M$, de masse $\mu$ positionnée en $\overrightarrow{CM}=\overrightarrow{r}$ est attirée par un corps fictif de masse égale à la masse totale du système $(m_1+m_2)$, placé en $C$. Sa position et sa vitesse initiales sont déterminées par les positions et vitesses initiales de $M_1$ et $M_2$.
\[
\overrightarrow{r}_0=\overrightarrow{r}_{2_0}-\overrightarrow{r}_{1_0} \quad \mbox{et} \quad \overrightarrow{v}^*_0=\overrightarrow{v}_{0/R}=\overrightarrow{v}_{2_0}^*-\overrightarrow{v}_{1_0}^*=
\overrightarrow{v}_{{0_2}/R}-\overrightarrow{v}_{{0_1}/R}
\]
Si nous connaissons le mouvement de $M$ nous pourrons en déduire les mouvements respectifs de $M_1$ et $M_2$, puisque d'après la définition du centre de masse :
\begin{equation}
\overrightarrow{r_1}=-\frac{m_2}{m_1+m_2}\overrightarrow{r}\quad\mbox{et}\quad \overrightarrow{r_2}=\frac{m_1}{m_1+m_2}\overrightarrow{r}
\label{r1r2}
\end{equation}
Les trajectoires de $M_1$ et $M_2$ se déduiront de celle de $M$ par une homothétie de centre C.
\begin{center}
\begin{pspicture}(-3.5,-1)(6.5,8.25)
%\psframe(-3.5,-1)(5.5,9)
\pstVerb{/x01 -2 def /y01 1 def
         /x02  4 def /y02 5 def
         /xr0 x02 x01 sub def
         /yr0 y02 y01 sub def
         /M1 3 def
         /M2 1 def
         /Mt M1 M2 add def
         /xG0 x01 M1 mul x02 M2 mul add Mt div def
         /yG0 y01 M1 mul y02 M2 mul add Mt div def}%
\pnode(0,0){O}
\pnode(!x01 y01){M1}
\pnode(!x02 y02){M2}
\pnode(!xG0 yG0){C}
\pscircle[fillstyle=solid,fillcolor=red!50](M1){0.21}
\pscircle[fillstyle=solid,fillcolor=blue!50](M2){0.07}
\pscircle[hatchcolor=gray,fillstyle=vlines,hatchwidth=0.02,hatchsep=0.05](C){0.25}
\psdot(C)
\psline[linestyle=dotted](M1)(M2)
\psline[linestyle=dashed](M2)(O)(M1)
\psline[linecolor=gray]{<->}(5,0)(0,0)(0,6)
\rput(C){\psline[style=vecteurA]{<->}(5.5,0)(0,0)(0,6)%
         \uput[l](0,5.8){$\mathcal{R^*}$}
         \psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt](-3,-3)(7,6)%
         \uput[r](0,5.8){$y^*$}
         \uput[u](5.5,0){$x^*$}
         \pnode(!xr0 yr0){M}
         \pscircle*(M){0.05}
         \pcline[offset=0.3,linewidth=0.01,arrowsize=0.1]{|->}(0,0)(M)
         \ncput*[nrot=:U]{$\overrightarrow{r}$}
         \psline[linestyle=dotted](M)\uput[r](M){$M(\mu)$}
         \rput(M){\psline[style=vecteurC]{->}(!xr0 5 div neg yr0 5 div neg)\uput[dr](!xr0 10 div neg yr0 10 div neg){$\overrightarrow{F}_{1/2}$}}}
\uput[dl](O){$O$}
\uput[r](0,5.8){$y$}
\uput[u](4.8,0){$x$}
\uput[l](0,5.8){$\mathcal{R}$}
\uput{0.22}[l](M1){$M_1$}
\uput{0.1}[dr](M2){$M_2$}
\uput{0.3}[d](C){$C$}
%\pcline[offset=10pt,linestyle=none](M1)(M2)
%\ncput*[nrot=:U]{$r$}
%\rput(M1){\psline[style=vecteurC]{->}(!xr0 5 div yr0 5 div)\uput[u](!xr0 10 div yr0 10 div){$\overrightarrow{F}_{2/1}$}}
%\rput(M2){\psline[style=vecteurC]{->}(!xr0 5 div neg yr0 5 div neg)\uput[u](!xr0 10 div neg yr0 10 div neg){$\overrightarrow{F}_{1/2}$}}
\end{pspicture}
\end{center}
Pour étudier, dans $\mathcal{R}^*$, le mouvement de ce point matériel fictif $M$ de masse $(\mu)$ soumis à la force centrale $\overrightarrow{F}$, il est avantageux de passer en coordonnées polaires.
\begin{center}
\begin{pspicture}(-3.5,-1)(6.5,8.25)
%\psframe(-3.5,-1)(5.5,9)
\pstVerb{/x01 -2 def /y01 1 def
         /x02  4 def /y02 5 def
         /xr0 x02 x01 sub def
         /yr0 y02 y01 sub def
         /theta_0 yr0 xr0 atan def
         /M1 3 def
         /M2 1 def
         /Mt M1 M2 add def
         /xG0 x01 M1 mul x02 M2 mul add Mt div def
         /yG0 y01 M1 mul y02 M2 mul add Mt div def}%
\pnode(0,0){O}
\pnode(!x01 y01){M1}
\pnode(!x02 y02){M2}
\pnode(!xG0 yG0){C}
\pscircle[fillstyle=solid,fillcolor=red!50](M1){0.21}
\pscircle[fillstyle=solid,fillcolor=blue!50](M2){0.07}
\pscircle[hatchcolor=gray,fillstyle=vlines,hatchwidth=0.02,hatchsep=0.05](C){0.25}
\psdot(C)
\psline[linestyle=dotted](M1)(M2)
\psline[linestyle=dashed](M2)(O)(M1)
\psline[linecolor=gray]{<->}(5,0)(0,0)(0,6)
\rput(C){\psline[style=vecteurA]{<->}(5.5,0)(0,0)(0,6)%
         \uput[l](0,5.8){$\mathcal{R^*}$}
         \psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt](-3,-3)(7,6)%
         \uput[r](0,5.8){$y^*$}
         \uput[u](5.5,0){$x^*$}
         \pnode(!xr0 yr0){M}
         \pscircle*(M){0.05}
         \pcline[offset=0.3,linewidth=0.01,arrowsize=0.1]{|->}(0,0)(M)
         \ncput*[nrot=:U]{$\overrightarrow{r}$}
         \psline[linestyle=dotted](M)\uput[r](M){$M(\mu)$}
         \rput(M){\psline[style=vecteurC]{->}(!xr0 5 div neg yr0 5 div neg)\uput[dr](!xr0 10 div neg yr0 10 div neg){$\overrightarrow{F}$}}}
\uput[dl](O){$O$}
\uput[r](0,5.8){$y$}
\uput[u](4.8,0){$x$}
\uput[l](0,5.8){$\mathcal{R}$}
\uput{0.22}[l](M1){$M_1$}
\uput{0.1}[dr](M2){$M_2$}
\uput{0.3}[d](C){$C$}
\rput(C){\psline[arrowinset=0.1,arrowsize=0.15]{->}(!theta_0 90 add cos theta_0 90 add sin)
         \psline[arrowinset=0.1,arrowsize=0.15]{->}(!theta_0 cos theta_0 sin)
         \uput{.1}[u](!theta_0 90 add cos theta_0 90 add sin){$\overrightarrow{u}_\theta$}
         \uput[d](!theta_0 cos theta_0 sin){$\overrightarrow{u}_r$}
         \psarc{->}(0,0){1.5}{0}{!theta_0}
         \uput{1.6}[!theta_0 2 div](0,0){$\theta$}}
\end{pspicture}
\end{center}
Pour la suite, retenons que dans $\mathcal{R}*$ :
\[
\overrightarrow{u}_r=
\left|
\begin{array}[m]{l}
\cos\theta\\
\sin\theta
\end{array}
\right.
\quad
\overrightarrow{u}_\theta=
\left|
\begin{array}[m]{l}
-\sin\theta\\
\cos\theta
\end{array}
\right.
\]
\[
\frac{\mathrm{d}\overrightarrow{u}_r}{\mathrm{d}t}%=
%\left|
%\begin{array}[m]{l}
%-\dot{\theta}\sin\theta\\
%\dot{\theta}\cos\theta
%\end{array}
%\right.
=\dot{\theta}\overrightarrow{u}_\theta
\quad
\frac{\mathrm{d}\overrightarrow{u}_\theta}{\mathrm{d}t}=
%\left|
%\begin{array}[m]{l}
%-\dot{\theta}\cos\theta\\
%-\dot{\theta}\sin\theta
%\end{array}
%\right.
-\dot{\theta}\overrightarrow{u}_r
\]
Dans ce repère, la position de $M$ et sa vitesse ont pour expressions respectives :
\[
\overrightarrow{r}=r \overrightarrow{u}_r\quad ;\quad \overrightarrow{v}=\frac{\mathrm{d}\overrightarrow{r}}{\mathrm{d}t}=\frac{\mathrm{d}r}{\mathrm{d}t}\overrightarrow{u}_r + r\frac{\mathrm{d}\overrightarrow{u}_r}{\mathrm{d}t}=\frac{\mathrm{d}r}{\mathrm{d}t}\overrightarrow{u}_r + r\dot{\theta}\overrightarrow{u}_\theta
\]
l'accélération s'écrit :
\[
\overrightarrow{\gamma}=\frac{\mathrm{d}\overrightarrow{v}}{\mathrm{d}t}=\ddot{r}\overrightarrow{u}_r +
 \dot{r}\dot{\theta}\overrightarrow{u}_\theta+
 \dot{r}\dot{\theta}\overrightarrow{u}_\theta +
 r\ddot{\theta}\overrightarrow{u}_\theta -r\dot{\theta}^2\overrightarrow{u}_r
\]
\[
=(\ddot{r} -r\dot{\theta}^2)\overrightarrow{u}_r +
 (2\dot{r}\dot{\theta}+r\ddot{\theta})\overrightarrow{u}_\theta
\]
En résumé :
\[
\overrightarrow{v}=
\left|
\begin{array}{l}
\dot{r}\\
r\dot{\theta}
\end{array}
\right.
\quad  \quad
\overrightarrow{\gamma}=
\left|
\begin{array}{l}
\ddot{r} -r\dot{\theta}^2\\
\dfrac{1}{r}\dfrac{\mathrm{d}}{\mathrm{d}t}(r^2\dot{\theta})
\end{array}
\right.
\]
L'équation (\ref{RFD}) nous donne, en divisant par $\mu$, l'accélération :
\[
 \dfrac{\mathrm{d}^2\overrightarrow{r}}{\mathrm{d}t^2}=-\mathcal{G}\frac{m_t}{r^2}\overrightarrow{u}_r
\]
Et puisque l'accélération radiale est nulle :
\[
\dfrac{1}{r}\dfrac{\mathrm{d}}{\mathrm{d}t}(r^2\dot{\theta})=0
\]
D'où :
\[
r^2\dot{\theta}=\mathrm{C^{ste}}
\]
Il est intéressant ici de calculer le moment cinétique de $M$ et d'exprimer la loi des aires.

Le moment cinétique :
\[
\overrightarrow{\sigma}=\overrightarrow{r} \wedge \mu\overrightarrow{v}
\]
$\overrightarrow{F}$ étant dirigée vers $C$, son moment par rapport à $C$ est nul. Le moment cinétique est donc constant. Il en découle que le mouvement est plan, puisque le vecteur-position $\overrightarrow{r}$ reste toujours perpendiculaire à un vecteur constant. Ce plan est le plan perpendiculaire à $\sigma_0$ contenant $C$.

On peut calculer le moment cinétique $\overrightarrow{\sigma}$, à partir des coordonnées polaires :
\[
\overrightarrow{\sigma}=\mu
\left|
\begin{array}{l}
r\\
0\\
0
\end{array}
\right.
\wedge
\left|
\begin{array}{l}
\dot{r}\\
r\dot{\theta}\\
0
\end{array}
\right.
=
\left|
\begin{array}{l}
0\\
0\\
r^2\dot{\theta}
\end{array}
\right.
=\mu r^2\dot{\theta}\overrightarrow{u}_z
\]
La constante des aires(voir plus loin pourquoi elle s'appelle ainsi) est $\mathrm{C}=\dfrac{\sigma}{\mu}=r^2\dot{\theta}$, c'est aussi la valeur du moment cinétique par unité de masse.

On peut aussi calculer le moment cinétique en coordonnées cartésiennes. Les conditions initiales étant fixées, $\overrightarrow{v}_0$ est donné par ses coordonnées, ainsi que $\overrightarrow{r}_0$.
\[
\overrightarrow{\sigma}=\mu
\left|
\begin{array}{l}
x_0\\
y_0\\
0
\end{array}
\right.
\wedge
\left|
\begin{array}{l}
v_{x_0}\\
v_{y_0}\\
0
\end{array}
\right.
=
\left|
\begin{array}{l}
0\\
0\\
x_0 v_{y_0} - y_0 v_{x_0}
\end{array}
\right.
=\mu(x_0 v_{y_0} - y_0 v_{x_0})\overrightarrow{u}_z
\]
Les livres d'astronomie\footnote{page 43 dans \textit{Le mouvement du satellite} : conférences et exercices de mécanique spatiale. Cepadues-éditions 1983.} préfèrent définir la vitesse par un autre paramètre $\gamma$, \textit{angle de l'horizontale locale }$\overrightarrow{u}_z\wedge \overrightarrow{u}_r$ (c'est-à-dire $\overrightarrow{u}_\theta$) \textit{avec la vitesse}, en se plaçant toujours en coordonnées polaires.
\begin{center}
\begin{pspicture}(-3.5,-1)(6.5,8.25)
%\psframe(-3.5,-1)(5.5,9)
\pstVerb{/x01 -2 def /y01 1 def
         /x02  4 def /y02 5 def
         /xr0 x02 x01 sub def
         /yr0 y02 y01 sub def
         /theta_0 yr0 xr0 atan def
         /M1 3 def
         /M2 1 def
         /Mt M1 M2 add def
         /xG0 x01 M1 mul x02 M2 mul add Mt div def
         /yG0 y01 M1 mul y02 M2 mul add Mt div def}%
\pnode(0,0){O}
\pnode(!x01 y01){M1}
\pnode(!x02 y02){M2}
\pnode(!xG0 yG0){C}
%\pscircle[fillstyle=solid,fillcolor=red!50](M1){0.21}
%\pscircle[fillstyle=solid,fillcolor=blue!50](M2){0.07}
\pscircle[hatchcolor=gray,fillstyle=vlines,hatchwidth=0.02,hatchsep=0.05](C){0.25}
\psdot(C)
%\psline[linestyle=dotted](M1)(M2)
%\psline[linestyle=dashed](M2)(O)(M1)
\psline[linecolor=gray]{<->}(5,0)(0,0)(0,6)
\rput(C){\psline[style=vecteurA]{<->}(5.5,0)(0,0)(0,6)%
         \uput[l](0,5.8){$\mathcal{R^*}$}
         \psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt](-3,-3)(7,6)%
         \uput[r](0,5.8){$y^*$}
         \uput[u](5.5,0){$x^*$}
         \pnode(!xr0 yr0){M}
         \pscircle*(M){0.05}
%         \pcline[offset=0.3,linewidth=0.01,arrowsize=0.1]{|->}(0,0)(M)
%         \ncput*[nrot=:U]{$\overrightarrow{r}$}
         \psline[linestyle=dotted](M)\uput[r](M){$M(\mu)$}
         \rput(M){\psline[style=vecteurC]{->}(!xr0 5 div neg yr0 5 div neg)\uput[dr](!xr0 10 div neg yr0 10 div neg){$\overrightarrow{F}$}
         \psline[linestyle=dotted](!theta_0 90 add cos 2 mul theta_0 90 add sin 2 mul)
         \psline[arrowinset=0.1,arrowsize=0.15]{->}(!theta_0 90 add cos theta_0 90 add sin)
         \psline[arrowinset=0.1,arrowsize=0.15,linecolor=blue]{->}(!theta_0 45 add cos 2 mul theta_0 45 add sin 2 mul)
         \uput[r](!theta_0 45 add cos 1.75 mul theta_0 45 add sin 1.55 mul){$\overrightarrow{v}$}
         \psarc{->}(0,0){1.2}{!theta_0 45 add}{!theta_0 90 add}
         \uput{1.3}[!theta_0 67.5 add](0,0){$\gamma$}}
         }
\uput[dl](O){$O$}
\uput[r](0,5.8){$y$}
\uput[u](4.8,0){$x$}
\uput[l](0,5.8){$\mathcal{R}$}
%\uput{0.22}[l](M1){$M_1$}
%\uput{0.1}[dr](M2){$M_2$}
\uput{0.3}[d](C){$C$}
\rput(C){\psline[arrowinset=0.1,arrowsize=0.15]{->}(!theta_0 90 add cos theta_0 90 add sin)
         \psline[arrowinset=0.1,arrowsize=0.15]{->}(!theta_0 cos theta_0 sin)
         \uput{.1}[u](!theta_0 90 add cos theta_0 90 add sin){$\overrightarrow{u}_\theta$}
         \uput[d](!theta_0 cos theta_0 sin){$\overrightarrow{u}_r$}
         \psarc{->}(0,0){1.5}{0}{!theta_0}
         \uput{1.6}[!theta_0 2 div](0,0){$\theta$}}
\end{pspicture}
\end{center}
Dans ce cas, on peut exprimer la constante des aires par :
$
\mathrm{C}=rv\cos\gamma
$
\begin{center}
\begin{pspicture}(-2,-2)(7,5)
\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]
\pstVerb{/par 2 def
         /exc 0.7 def
         /Phi 30 def
         /radius {par 1 exc t Phi sub cos mul add div} def
         /xE {radius t cos mul neg} def
         /yE {radius t sin mul neg} def}%
\parametricplot[plotpoints=360]{0}{360}{
    xE
    yE}%
\pscustom[fillstyle=vlines,hatchwidth=0.02,hatchcolor=gray]{%
    \psline(0,0)(! /t 200 def xE yE)
\parametricplot[plotpoints=360]{200}{210}{
    xE
    yE}%
   \psline(! /t 210 def xE yE)(0,0)}
\parametricplot[plotpoints=360,arrows=->]{200}{210}{
   /radius {2.25 1 exc t Phi sub cos mul add div} def
    xE
    yE}%
\psarc{->}(0,0){3}{20}{30}
\uput{3.2}[25](0,0){\psframebox[linestyle=none,fillstyle=solid,fillcolor=white,framesep=0pt]{d$\theta$}}
\rput(5,2.3){\psframebox[linestyle=none,fillstyle=solid,fillcolor=white,framesep=0pt]{d$\mathcal{A}$}}
%\psline[linearc=0.25]{->}(2.2,1.8)(3.5,1.8)(3.5,1.5)
\pcline[offset=-0.3,linestyle=none](0,0)(! /t 200 def xE yE)
\ncput*[nrot=:U]{$r$}
\uput[l](0,4.8){$\mathcal{R^*}$}
\uput[r](0,4.8){$y^*$}
\uput[u](7,0){$x^*$}
\uput[dl](0,0){$C$}
\psline[style=vecteurA]{<->}(7,0)(0,0)(0,5)
\end{pspicture}
\end{center}
Géométriquement lorsque $\theta$ varie de $\mathrm{d}\theta$, la surface balayée par le rayon vecteur vaut :
\[
\mathrm{d}\mathcal{A}=\frac{1}{2}r^2\mathrm{d}\theta
\]
On définit ainsi la vitesse \textit{aérolaire}, qui est l'aire balayée par unité de temps, en fonction de C :
\[
\frac{\mathrm{d}A}{\mathrm{d}t}=\frac{\mathrm{C}}{2}
\]
On retrouve ainsi la deuxième loi de Képler.

Pour déterminer l'équation de la trajectoire, on utilise la méthode de Binet, qui consiste en un changement de variable en posant $u=\dfrac{1}{r}$.
\[
\dfrac{\mathrm{d}r}{\mathrm{d}t}=-\dfrac{1}{u^2}\dfrac{\mathrm{d}u}{\mathrm{d}\theta}\dfrac{\mathrm{d}\theta}{\mathrm{d}t}
\]
Nous avons la constante des aires $\mathrm{C}=r^2\dot{\theta}=\dfrac{\dot{\theta}}{u^2}$, par conséquent :
\[
\dfrac{\mathrm{d}r}{\mathrm{d}t}=-\mathrm{C}\dfrac{\mathrm{d}u}{\mathrm{d}\theta}
\]
\[
\dfrac{\mathrm{d}^2r}{\mathrm{d}t^2}=-\mathrm{C}\dfrac{\mathrm{d}^2u}{\mathrm{d}\theta^2}\dfrac{\mathrm{d}\theta}{\mathrm{d}t}
\]
On remplace $\dot{\theta}$ par $\mathrm{C}u^2$ :
\[
\dfrac{\mathrm{d}^2r}{\mathrm{d}t^2}=-\mathrm{C^2}u^2\dfrac{\mathrm{d}^2u}{\mathrm{d}\theta^2}
\]
Nous avons vu qu'en appliquant la loi de Newton, la composante de l'accélération $\gamma$ suivant $\sqrt[n]{u}_r$ vaut :
\[
\ddot{r} -r\dot{\theta}^2=-\mathcal{G}\frac{m_t}{r^2}
\]
Ce qui s'écrit avec la variable $u$ :
\[
-\mathrm{C^2}u^2\dfrac{\mathrm{d}^2u}{\mathrm{d}\theta^2}-\mathrm{C^2}u^3=-\mathcal{G}m_tu^2
\]
Ce qui, après simplification, donne :
\begin{equation}
\psframebox[linestyle=none,fillstyle=solid,fillcolor={[rgb]{1 1 0.5}}]{%
\dfrac{\mathrm{d}^2u}{\mathrm{d}\theta^2}+u=\dfrac{\mathcal{G}m_t}{\mathrm{C^2}}}
\label{equ}
\end{equation}
La solution générale de cette équation peut s'écrire :
\[
u=A\cos(\theta-\varphi)+\dfrac{\mathcal{G}m_t}{\mathrm{C^2}}
\]
$A$ et $\varphi$ sont fixés par les conditions initiales. Le rayon-vecteur a pour expression :
\[
r=\frac{1}{A\cos(\theta-\varphi)+\dfrac{\mathcal{G}m_t}{\mathrm{C^2}}}
\]
Il s'écrit encore ainsi :
\[
r=\dfrac{\dfrac{\mathrm{C^2}}{\mathcal{G}m_t}}{A\dfrac{\mathrm{C^2}}{\mathcal{G}m_t}\cos(\theta-\varphi)+1}
\]
Si on pose $p=\dfrac{\mathrm{C^2}}{\mathcal{G}m_t}$ et $\mathrm{e}=Ap$, on reconnaît l'équation d'une conique dont l'un des foyers est $C$, de paramètre $p$ et d'excentricité $\mathrm{e}$.
\begin{equation}
\psframebox[linestyle=none,fillstyle=solid,fillcolor={[rgb]{1 1 0.5}}]{%
r=\dfrac{p}{1+e\cos(\theta-\varphi)}}
\label{equr}
\end{equation}
Déterminons les constantes en fonction des conditions initiales.

Nous avons vu que la constante des aires C pouvait se calculer de deux façons suivant le choix fait pour définir la vitesse initiale, soit $
\mathrm{C}=r_0v_0\cos\gamma_0$, soit $C=x_0 v_{y_0} - y_0 v_{x_0}$.  Par conséquent $p$ est bien déterminé.

Lorsque $\theta=\varphi$ alors $r=\dfrac{p}{1+\mathrm{e}}$, ce qui est la plus petite valeur possible de $r$. On est donc au périgée P de la conique, on note, en ce point, le rayon-vecteur par $\overrightarrow{r}_\mathrm{P}$. En ce point la vitesse $\overrightarrow{v}_\mathrm{P}$ est perpendiculaire à $\overrightarrow{r}_\mathrm{P}$ et la constante des aires peut s'écrire : $\mathrm{C}=r_\mathrm{P}v_\mathrm{P}$.

L'angle $\varphi$ est l'inclinaison du grand axe de la conique avec $Ox$, plus précisément l'angle du rayon-vecteur au périgée, $\overrightarrow{r}_P$ avec $Ox$.

L'énergie mécanique se conserve au cours du mouvement, elle est la somme de l'énergie cinétique et de l'énergie potentielle. En fonction des conditions initiales elle s'écrit(on rappelle que $\mu$ est la masse(\textit{réduite})) de la particule fictive $M$ étudiée :
\[
\mathcal{E}=\frac{1}{2}\mu v_0^2-\mathcal{G}\frac{\mu(m_1+m_2)}{r_0}\quad \text{ou}\quad \frac{1}{2}\mu v_0^2-\mathcal{G}\frac{m_1m_2}{r_0}
\]
Au périgée, elle s'exprimera par :
\[
\mathcal{E}=\frac{1}{2}\mu v_\mathrm{P}^2-\mathcal{G}\frac{m_1m_2}{r_\mathrm{P}}
\]
Pour exprimer l'énergie dans le cas général, calculons la vitesse en fonction de $u$. Rappelons que :
\[
\overrightarrow{v}=
\left|
\begin{array}{l}
\dot{r}\\
r\dot{\theta}
\end{array}
\right.
\]
Calculons séparément $\dot{r}$ et $r\dot{\theta}$ :
\[
\dot{r}=\dfrac{\mathrm{d}r}{\mathrm{d}u}\dfrac{\mathrm{d}u}{\mathrm{d}\theta}\dfrac{\mathrm{d}\theta}{\mathrm{d}t}=
   -\dfrac{1}{u^2}\dfrac{\mathrm{d}u}{\mathrm{d}\theta}\dot{\theta}=
   -\mathrm{C}\dfrac{\mathrm{d}u}{\mathrm{d}\theta}
\]
\[
r\dot{\theta}=\dfrac{1}{u}\mathrm{C}u^2=\mathrm{C}u
\]
\[
v^2=\mathrm{C^2}\left[\left(\frac{\mathrm{d}u}{\mathrm{d}\theta}\right)^2+u^2\right]
\]
\[
v^2=\mathrm{C^2}\left[
 \frac{\mathrm{e^2}}{p^2}\sin^2(\theta-\varphi)+\frac{1}{p^2}\left[1+2\mathrm{e}\cos(\theta-\varphi)+\mathrm{e^2}\cos^2(\theta-\varphi)\right]
  \right]
\]
\[
v^2=\frac{\mathrm{C^2}}{p^2}\left[1+2\mathrm{e}\cos(\theta-\varphi)+\mathrm{e^2} \right]
\]
L'énergie cinétique s'exprime par :
\[
E_C=\frac{1}{2}\mu\frac{\mathrm{C^2}}{p^2}\left[ 1+2\mathrm{e}\cos(\theta-\varphi)+\mathrm{e^2} \right]
\]
On se souvient que $p=\dfrac{\mathrm{C^2}}{\mathcal{G}m_t}$, il vient :
\[
E_C=\frac{1}{2}\frac{\mathcal{G}\mu m_t}{p}\left[ 1+2\mathrm{e}\cos(\theta-\varphi)+\mathrm{e^2} \right]
\]
Pour l'énergie potentielle nous avons :
\[
E_p=-\frac{\mathcal{G}\mu m_t}{r}=-\frac{\mathcal{G}\mu m_t}{p}\left[1+e\cos(\theta-\varphi)\right]
\]
et pour l'énergie mécanique, nous obtenons après simplification :
\[
\mathcal{E}=E_C+E_p=\frac{\mathcal{G}\mu m_t}{2p}(\mathrm{e^2}-1)
\]
On en déduit l'excentricité :
\[
\mathrm{e^2}=\frac{2p\mathcal{E}}{\mathcal{G}\mu m_t}+1\Longrightarrow \mathrm{e}=\sqrt{\frac{2p\mathcal{E}}{\mathcal{G}\mu m_t}+1}
\]
On peut l'exprimer en fonction de la constante des aires en posant $K=\mathcal{G}m_t$ et en remplaçant $p$ par $\dfrac{\mathrm{C^2}}{K}$ :
\begin{equation}
\psframebox[linestyle=none,fillstyle=solid,fillcolor={[rgb]{1 1 0.5}}]{%
\mathrm{e}=\sqrt{\frac{2\mathrm{C^2}\mathcal{E}}{\mu K^2}+1}}
\label{exc}
\end{equation}
Il reste à déterminer $\varphi$, c'est le problème le plus simple à résoudre.
\[
r_0=\frac{p}{1+\mathrm{e}\cos(\theta_0-\varphi)}
\]
On en déduit :
\begin{equation}
\psframebox[linestyle=none,fillstyle=solid,fillcolor={[rgb]{1 1 0.5}}]{%
\varphi=\theta_0-\arccos\left[\left(\frac{p}{r_0}-1\right)\frac{1}{\mathrm{e}}\right]}
\label{phi}
\end{equation}
On distingue trois cas :
\begin{itemize}
\item
  $\mathrm{e}<1$ ou $\mathcal{E}<0$ : ellipse
\item
  $\mathrm{e}=1$ ou $\mathcal{E}=0$ : parabole
\item
 $\mathrm{e}=>$ ou $\mathcal{E}>0$ : hyperbole
\end{itemize}
Et on s'intéresse maintenant au mouvement elliptique. Les conditions initiales $(\overrightarrow{r}_0,\overrightarrow{v}_0)$ doivent vérifier :
\[
\mathcal{E}=\frac{1}{2}\mu v_0^2-\mathcal{G} \frac{m_1 m_2}{r_0}<0
\]
\begin{center}
\begin{pspicture}(-2,-3)(5,7)
\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]
\pstVerb{
/arccos {
   dup
   dup mul neg 1 add sqrt
   exch
   atan
} def
        /G 1 def
        /x01 2 def
        /y01 2 def
        /v0x1 .05 def
        /v0y1 0.25 def
        /x02 -2 def
        /y02 -0.5 def
        /v0x2 -0.25 def
        /v0y2 -0.5 def
        /M1 3 def
        /M2 1 def
        /Mt M1 M2 add def
        /mu M1 M2 mul M1 M2 add div def % masse réduite
        /K G Mt mul def
        /xG0 M1 x01 mul M2 x02 mul add M1 M2 add div def
        /yG0 M1 y01 mul M2 y02 mul add M1 M2 add div def
% conditions initiales pour le point réduit
        /xr0 x01 x02 sub def
        /yr0 y01 y02 sub def
        /theta_0 yr0 xr0 atan def
        /r0 xr0 dup mul yr0 dup mul add sqrt def
        /v0xr v0x1 v0x2 sub def
        /v0yr v0y1 v0y2 sub def
        /v0r_2 v0xr dup mul v0yr dup mul add def
% constante des aires
        /Cste xr0 v0yr mul yr0 v0xr mul sub def % Cste des aires
        /Energie 0.5 mu mul v0r_2 mul % 1/2 mv^2
                 G M1 M2 mul mul r0 div sub
         def
         /par Cste dup mul K div def % p
         /exc 2 Cste dup mul mul Energie mul mu div K dup mul div 1 add sqrt def % e
         /a_2 par 1 exc dup mul sub div def % demi-grand axe
         /b_2 a_2 1 exc dup mul sub sqrt mul def % demi-petit axe
         /c_2 a_2 exc mul def % distance focale
         % phase
         /Phi theta_0 par r0 div 1 sub exc div arccos sub def
         /rP par 1 exc add div def
         /rA par 1 exc sub div def
% vitesses à l'apogée et au périgée
         /vA G Mt mul 2 rA div 1 a_2 div sub mul sqrt def
         /vP G Mt mul 2 rP div 1 a_2 div sub mul sqrt def
% positions du périgée et de l'apogée
         /xP rP Phi cos mul def
         /yP rP Phi sin mul def
         /xA rA Phi cos mul neg def
         /yA rA Phi sin mul neg def
         /xW xA xP add 2 div def
         /yW yA yP add 2 div def
% periode
       /periode 6.28 a_2 3 exp G div Mt div sqrt mul def
      /radius {par 1 exc t Phi sub cos mul add div} def
      /xE {radius t cos mul} def
      /yE {radius t sin mul}def
        }%
\parametricplot[plotpoints=360]{0}{360}{xE yE}%
\uput[l](0,6.8){$\mathcal{R^*}$}
\uput[r](0,6.8){$y^*$}
\uput[u](5,0){$x^*$}
\uput[l](0,0){$C$}
\psdot(0,0)
\psline[style=vecteurA]{<->}(5,0)(0,0)(0,7)
\pnode(!xr0 yr0){M0}\psdot(M0)
\rput(M0){\psline[unit=2,style=vecteurA]{->}(!v0xr v0yr)\uput[r](0,0){$M_0$}\uput[r](!v0xr 2 mul v0yr 2 mul){$\overrightarrow{v}_0$}}
\psline(M0)
\pcline[offset=0.25,linestyle=none](0,0)(M0)
\ncput*[nrot=:U]{$r_0$}
\psarc[style=vecteurB]{->}(0,0){1.5}{0}{!theta_0}
\uput{1.6}[!theta_0 2 div](0,0){$\theta_0$}
\pnode(!rA Phi cos mul neg rA Phi sin mul neg){A}
\psline[linestyle=dotted](A)
\uput[ur](A){$A$}
\pnode(!rP Phi cos mul rP Phi sin mul){P}
\psline[linestyle=dotted](P)
\uput[dl](P){$P$}
\rput(-2,0){\psPrintValue[decimals=4]{exc}}
\pnode(!xW yW){W}\psdot(W)\uput[r](W){$\Omega$}
\rput{!Phi}(W){\pnode(!0 b_2){B1}\pnode(!0 b_2 neg){B2}\psline[linestyle=dotted](B1)(B2)}
\uput[ul](B2){$B$}
\pcline[offset=-0.25,linestyle=none](W)(B2)
\ncput*{$b$}
\pcline[offset=0.25,linestyle=none](W)(A)
\ncput*{$a$}
\pcline[offset=0.25,linestyle=none](0,0)(W)
\ncput*{$c$}
\psline[linecolor=blue](W)(B2)
\psline[linecolor=red](W)(A)
\psline[linecolor=green](W)(0,0)
\rput{! 180 Phi add}(W){\psdot(!c_2 0)}
\psarcn{->}(0,0){0.5}{0}{!Phi}\uput{0.55}[!Phi 2 div](0,0){$\varphi$}
\rput{!Phi}(P){\pnode(! 0 vP){vP}\psline[style=vecteurA]{->}(vP)}
\uput[r](vP){$\overrightarrow{v}_P$}
\rput{!Phi}(A){\pnode(! 0 vA neg){vA}\psline[style=vecteurA]{->}(vA)}
\uput[u](vA){$\overrightarrow{v}_A$}
\end{pspicture}
\end{center}
L'un des foyers de l'ellipse est le centre attracteur $C$. L'autre est le symétrique de $C$ par rapport à $\Omega$.

Au périgée, le rayon-vecteur vaut : $r_P=\dfrac{p}{1+\mathrm{e}}$, à l'apogée $r_A=\dfrac{p}{1-\mathrm{e}}$.

Le demi-grand axe est égal à : $a=\dfrac{r_P+r_A}{2}=\dfrac{p}{1-\mathrm{e^2}}$.

La distance focale est égale à : $c=a-r_P=\dfrac{p\mathrm{e}}{1-\mathrm{e^2}}=a\mathrm{e}$.

Pour calculer le demi-petit axe de l'ellipse $b$, on peut rappeler que l'ellipse est lieu des points dont la somme des distances aux foyers est égale à $2a$. Au point $B$ cette somme est égale à $2a$, par conséquent le théorème de Pythagore appliqué dans le triangle rectangle $C\Omega B$ donne la valeur du demi-axe :  $b=\sqrt{a^2-c^2}=a\sqrt{1-\mathrm{e^2}}$.

Rappelons que $\mathrm{e^2}=\dfrac{2p\mathcal{E}}{\mathcal{G}\mu m_t}+1$. On peut exprimer $\mathcal{E}$ en fonction de $\mathrm{e}$ :

\[
\mathcal{E}=\frac{\mathcal{G}\mu m_t}{2p}(\mathrm{e^2}-1)=-\frac{\mathcal{G}\mu m_t}{2a}\quad\text{ou}\quad\mathcal{E}=-\frac{\mathcal{G}m_1m_2}{2a}
\]
Pour le système $\{M_1,M_2\}$,  son énergie s'exprime très simplement en fonction du grand axe de l'ellipse $2a$.

Calculons la vitesse tout au long de l'ellipse :

\[
\mathcal{E}=\frac{1}{2}\mu v^2-\frac{\mathcal{G}\mu m_t}{r}=-\frac{\mathcal{G}\mu m_t}{2a}
\]
\[
v^2=\mathcal{G}m_t\left(\frac{2}{r}-\frac{1}{a}\right)
\]
On en déduit les vitesses au périgée et à l'apogée :
\[
v_P=\sqrt{\frac{\mathcal{G}m_t}{p}}(1+\mathrm{e})\qquad v_A=\sqrt{\frac{\mathcal{G}m_t}{p}}(1-\mathrm{e})
\]
La vitesse est maximale au périgée et minimale à l'apogée.

La période de révolution se  détermine à partir de la vitesse \textit{aérolaire}. Cette vitesse vaut :
\[
\mathcal{V}=\frac{\mathrm{C}}{2}=\frac{\sqrt{p\mathcal{G}m_t}}{2}
\]
Sachant que surface de l'ellipse est $S=\pi a b$ :
\[
T=\frac{2\pi a b}{\sqrt{p\mathcal{G}m_t}}
\]
On remplace $b$ par $a\sqrt{(1-\mathrm{e^2})}$ et $p$ par $a(1-\mathrm{e^2})$ :
\[
T=\frac{2\pi a^2 \sqrt{(1-\mathrm{e^2})}}{\sqrt{a(1-\mathrm{e^2})\mathcal{G}m_t}}=\frac{2\pi a^2}{\sqrt{a\mathcal{G}m_t}}
\]
En élevant les deux membres au carré on retrouve la troisième loi de Képler :
\begin{equation}
\psframebox[linestyle=none,fillstyle=solid,fillcolor={[rgb]{1 1 0.5}}]{%
\frac{T^2}{a^3}=\frac{4\pi^2}{\mathcal{G}m_t}}
\label{Kepler3}
\end{equation}
\section{Le mouvement des 2 corps dans le repère du centre de masse}
Nous avons vu (\ref{r1r2}), que les trajectoires de $M_1$ et $M_2$ se déduisent de celle de $M$ par une homothétie de centre~C.
\[
\overrightarrow{r_1}=-\frac{m_2}{m_1+m_2}\overrightarrow{r}\quad\mbox{et}\quad \overrightarrow{r_2}=\frac{m_1}{m_1+m_2}\overrightarrow{r}
\]
Les trajectoires sont donc deux ellipses dont les caractéristiques, paramètre, demi-grand axe et demi-petit axe, se déduisent de l'ellipse de la particule de masse réduite dans le même rapport.

\begin{itemize}
  \item Paramètre : $p_1=p\dfrac{m_2}{m_1+m_2}$.
  \item rayon-vecteur à l'apogée  : $r_{A_1}=r_A\dfrac{m_2}{m_1+m_2}$
  \item rayon-vecteur au périgée  : $r_{P_1}=r_P\dfrac{m_2}{m_1+m_2}$
  \item demi-grand axe : $a_1=\dfrac{r_{A_1}+r_{P_1}}{2}=\dfrac{m_2}{m_1+m_2}\dfrac{r_{A}+r_{P}}{2}=a\dfrac{m_2}{m_1+m_2}$
  \item excentricité, on vérifie facilement qu'elle est inchangée : $\mathrm{e_1}=\dfrac{r_{A_1}-r_{P_1}}{r_{A_1}+r_{P_1}}=\mathrm{e}$
\end{itemize}
pour la deuxième :
\begin{itemize}
  \item Paramètre : $p_2=p\dfrac{m_1}{m_1+m_2}$.
  \item rayon-vecteur à l'apogée  : $r_{A_2}=r_A\dfrac{m_1}{m_1+m_2}$
  \item rayon-vecteur au périgée  : $r_{P_2}=r_P\dfrac{m_1}{m_1+m_2}$
  \item demi-grand axe : $a_1=\dfrac{r_{A_2}+r_{P_2}}{2}=\dfrac{m_1}{m_1+m_2}\dfrac{r_{A}+r_{P}}{2}=a\dfrac{m_1}{m_1+m_2}$
  \item excentricité : $\mathrm{e_2}=\dfrac{r_{A_2}-r_{P_2}}{r_{A_2}+r_{P_2}}=\mathrm{e}$
\end{itemize}
Leurs équations respectives s'écrivent :
\[
r_1=\dfrac{p_1}{1+\mathrm{e}\cos(\theta-\varphi)} \quad r_2=\dfrac{p_2}{1+\mathrm{e}\cos(\theta-\varphi)}
\]
Dans l'exemple suivant $\{m_1=3,m_2=1\}$. Les vitesses initiales dans $\mathcal{R}^*$ sont indiquées par une flèche. La trajectoire de $M_1$ est en bleu, celle de $M_2$ en rouge et celle du point fictif est en pointillés.
\begin{center}
\begin{pspicture}(-8,-10)(3,3)
\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]
\pstVerb{
/arccos {
   dup
   dup mul neg 1 add sqrt
   exch
   atan
} def
        /G 1 def
        /x01 2 def
        /y01 2 def
        /v0x1 .2 def
        /v0y1 0.25 def
        /x02 -3 def
        /y02 0 def
        /v0x2 -0.25 def
        /v0y2 -0.5 def
        /M1 3 def
        /M2 1 def
        /Mt M1 M2 add def
        /mu M1 M2 mul M1 M2 add div def % masse réduite
        /K G Mt mul def
        /xG0 M1 x01 mul M2 x02 mul add M1 M2 add div def
        /yG0 M1 y01 mul M2 y02 mul add M1 M2 add div def
        /vG0x M1 v0x1 mul M2 v0x2 mul add M1 M2 add div def
        /vG0y M1 v0y1 mul M2 v0y2 mul add M1 M2 add div def
        /K1 M2 Mt div neg def
        /K2 M1 Mt div def
% conditions initiales pour le point réduit
        /xr0 x02 x01 sub def
        /yr0 y02 y01 sub def
        /theta_0 yr0 xr0 atan def
        /r0 xr0 dup mul yr0 dup mul add sqrt def
        /v0xr v0x2 v0x1 sub def
        /v0yr v0y2 v0y1 sub def
        /v0r_2 v0xr dup mul v0yr dup mul add def
% constante des aires
        /Cste xr0 v0yr mul yr0 v0xr mul sub def % Cste des aires
        /Energie 0.5 mu mul v0r_2 mul % 1/2 mv^2
                 G M1 M2 mul mul r0 div sub
         def
         /par Cste dup mul K div def % p
         /exc 2 Cste dup mul mul Energie mul mu div K dup mul div 1 add sqrt def % e
         /a_2 par 1 exc dup mul sub div def % demi-grand axe
         /b_2 a_2 1 exc dup mul sub sqrt mul def % demi-petit axe
         /c_2 a_2 exc mul def % distance focale
         % phase
         /Phi theta_0 par r0 div 1 sub exc div arccos sub def
         /rP par 1 exc add div def
         /rA par 1 exc sub div def
% vitesses à l'apogée et au périgée
         /vA G Mt mul 2 rA div 1 a_2 div sub mul sqrt def
         /vP G Mt mul 2 rP div 1 a_2 div sub mul sqrt def
% positions du périgée et de l'apogée
         /xP rP Phi cos mul def
         /yP rP Phi sin mul def
         /xA rA Phi cos mul neg def
         /yA rA Phi sin mul neg def
         /xW xA xP add 2 div def
         /yW yA yP add 2 div def
% periode
       /periode 6.28 a_2 3 exp G div Mt div sqrt mul def
      /radius {par 1 exc t Phi sub cos mul add div} def
      /xE {radius t cos mul} def
      /yE {radius t sin mul}def
        }%
\parametricplot[plotpoints=360,linestyle=dotted]{0}{360}{xE yE}%
\parametricplot[plotpoints=360,linecolor=blue]{0}{360}{xE K1 mul yE K1 mul}%
\parametricplot[plotpoints=360,linecolor=red]{0}{360}{xE K2 mul yE K2 mul}%
\uput[l](0,2.8){$\mathcal{R^*}$}
\uput[r](0,2.8){$y^*$}
\uput[u](3,0){$x^*$}
\uput[l](0,0){$C$}
\psdot(0,0)
\psline[style=vecteurA]{<->}(3,0)(0,0)(0,3)
\pnode(!xr0 yr0){M0}\psdot(M0)
\rput(M0){\psline[unit=2,style=vecteurA]{->}(!v0xr v0yr)\uput[l](0,0){$M_0$}\uput[l](!v0xr 2 mul v0yr 2 mul){$\overrightarrow{v}_0$}}
\pcline[offset=0.25,linestyle=none](M0)(0,0)
\ncput*[nrot=:U]{$r_0$}
\psarcn[style=vecteurB]{->}(0,0){0.75}{0}{!theta_0}
\uput{0.75}[!theta_0 2 div 180 add](0,0){$\theta_0$}
\pnode(!rA Phi cos mul neg rA Phi sin mul neg){A}
\psline[linestyle=dotted](A)
%\uput[ur](A){$A$}
\pnode(!rP Phi cos mul rP Phi sin mul){P}
\psline[linestyle=dotted](P)
%\uput[dl](P){$P$}
%\rput(-2,0){\psPrintValue[decimals=4]{exc}}
\pnode(!xW yW){W}\psdot(W)\uput[r](W){$\Omega$}
%\psline[linecolor=blue](W)(B2)
%\psline[linecolor=red](W)(A)
%\psline[linecolor=green](W)(0,0)
%\rput{! 180 Phi add}(W){\psdot(!c_2 0)}
%\psarc{->}(0,0){0.5}{0}{!Phi}\uput{0.55}[!Phi 2 div](0,0){$\varphi$}
%\rput{!Phi}(P){\pnode(! 0 vP){vP}\psline[style=vecteurA]{->}(vP)}
%\uput[r](vP){$\overrightarrow{v}_P$}
%\rput{!Phi}(A){\pnode(! 0 vA neg){vA}\psline[style=vecteurA]{->}(vA)}
%\uput[u](vA){$\overrightarrow{v}_A$}
\pnode(!x01 xG0 sub y01 yG0 sub){M01}\psdot(M01)\uput[r](M01){$M_{0_1}$}
\pnode(!x02 xG0 sub y02 yG0 sub){M02}\psdot(M02)\uput[dr](M02){$M_{0_2}$}
\psline[linestyle=dashed](M0)(M02)(M01)
\rput(M01){\psline[unit=2,linecolor=blue]{->}(!v0x1 vG0x sub v0y1 vG0y sub)}
\rput(M02){\psline[unit=2,linecolor=red]{->}(!v0x2 vG0x sub v0y2 vG0y sub)}
\end{pspicture}
\end{center}
\section{Exemples à développer, réalisés avec pst-eqdf}
\[
\mu\frac{\mathrm{d}^2\overrightarrow{r}}{\mathrm{d}t^2}=-G\frac{m_1m_2}{r^3}\overrightarrow{r}
\]
\[
\frac{\mathrm{d}^2\overrightarrow{r}}{\mathrm{d}t^2}=-G\frac{(m_1+m_2)}{r^3}\overrightarrow{r}
\]

\[
\left\{
\begin{array}[m]{l}
 \ddot{x}=-G\displaystyle\frac{m_1+m_2}{r^3}x\\[1em]
 \ddot{y}=-G\displaystyle\frac{m_1+m_2}{r^3}y
\end{array}
\right.
\]
%  x    y    x'   y'
% y[0] y[1] y[2] y[3]
\def\FictifAlg{%
  y[2]|y[3]|%
  -(M1+M2)*y[0]/(y[0]^2+y[1]^2)^1.5|%
  -(M1+M2)*y[1]/(y[0]^2+y[1]^2)^1.5}
\begin{center}
\begin{pspicture}(-7,-4)(7,8)
\pstVerb{/G 1 def
        /x01 1 def
        /y01 2 def
        /v0x1 .1 def
        /v0y1 0.1 def
        /x02 -1 def
        /y02 -2 def
        /v0x2 -2 def
        /v0y2 0 def
        /M1 2 def
        /M2 20 def
        /mu M1 M2 mul M1 M2 add div def % masse réduite
        /xG M1 x01 mul M2 x02 mul add M1 M2 add div def
        /yG M1 y01 mul M2 y02 mul add M1 M2 add div def
        /x0 x01 x02 sub def
        /y0 y01 y02 sub def
        /v0x v0x1 v0x2 sub def
        /v0y v0y1 v0y2 sub def
        }%
%%  0  1   2   3  4  5   6   7
%% x1 y1 x'1 y'1 x2 y2 x'2 y'2
\def\InitCondred{ x0  y0  v0x  v0y}
\def\InitCond{ x01  y01  v0x1  v0y1 x02 y02  v0x2   v0y2}
\psset{method=rk4}
\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=8pt]%
\psequadiff[whichabs=0,whichord=1,
            plotpoints=1000,algebraic,
            tabname=X1Y1]{0}{50}{\InitCond}{\GravAlg}
\listplot[unit=1,linecolor=red]{X1Y1 aload pop}
\psequadiff[whichabs=4,whichord=5,
            plotpoints=1000,algebraic,
            tabname=X2Y2]{0}{50}{\InitCond}{\GravAlg}
\listplot[unit=1,linecolor=blue]{X2Y2 aload pop}
\psequadiff[whichabs=0,whichord=1,
            plotpoints=1000,algebraic,
            tabname=XrYr]{0}{50}{\InitCondred}{\FictifAlg}
\listplot[unit=1,linecolor=green]{XrYr aload pop}
\psdots(!x01 y01)(!x02 y02)
\psdot[linecolor=red](!xG yG)
\psline{<->}(7,0)(0,0)(0,8)
\end{pspicture}
\end{center}
%
\newpage
%  x    y    x'   y'
% y[0] y[1] y[2] y[3]
\def\FictifAlg{%
  y[2]|y[3]|%
  -(M1+M2)*y[0]/(y[0]^2+y[1]^2)^1.5|%
  -(M1+M2)*y[1]/(y[0]^2+y[1]^2)^1.5}
% on se place au centre de masse
\def\GravAlgIIcorps{%
  y[2]|y[3]|%
  M2*(y[4]-y[0])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5|%
  M2*(y[5]-y[1])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5|%
  y[6]|y[7]|%
  M1*(y[0]-y[4])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5|%
  M1*(y[1]-y[5])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5}
%%  0  1   2   3  4  5   6   7
%% x1 y1 x'1 y'1 x2 y2 x'2 y'2
\newpage
\begin{center}
\def\InitCondred{ xr0 yr0 v0xr  v0yr}
\begin{pspicture}(-7,-4)(7,8)
\pstVerb{/G 1 def
        /x01 2 def
        /y01 2 def
        /v0x1 .05 def
        /v0y1 0.25 def
        /x02 -2 def
        /y02 -0.35 def
        /v0x2 -0.25 def
        /v0y2 -0.5 def
        /M1 3 def
        /M2 1 def
        /Mt M1 M2 add def
        /mu M1 M2 mul M1 M2 add div def % masse réduite
        /xG M1 x01 mul M2 x02 mul add M1 M2 add div def
        /yG M1 y01 mul M2 y02 mul add M1 M2 add div def
        /x0 x01 x02 sub def
        /y0 y01 y02 sub def
        /v0x v0x1 v0x2 sub def
        /v0y v0y1 v0y2 sub def
        /mu M1 M2 mul M1 M2 add div def % masse réduite
        /xG M1 x01 mul M2 x02 mul add M1 M2 add div def
        /yG M1 y01 mul M2 y02 mul add M1 M2 add div def
        /vxG M1 v0x1 mul M2 v0x2 mul add M1 M2 add div def
        /vyG M1 v0y1 mul M2 v0y2 mul add M1 M2 add div def
% conditions initiales pour le point réduit
        /xr0 x01 x02 sub def
        /yr0 y01 y02 sub def
        /r0 xr0 dup mul yr0 dup mul add sqrt def
        /v0xr v0x1 v0x2 sub def
        /v0yr v0y1 v0y2 sub def
        /Lc xr0 v0yr mul yr0 v0xr mul sub mu mul def % moment cinetique
        /K xr0 v0yr mul yr0 v0xr mul sub def
        /Energie 0.5 mu mul v0xr dup mul v0yr dup mul add mul % 1/2 mv^2
%                 Lc dup mul 2 div mu div r0 dup mul div
                 G M1 M2 mul mul r0 div sub
         def
%        /par xr0 v0yr mul yr0 v0xr mul sub dup mul G Mt mul div def % paramètre de l'ellipse
%%%%%%%%%%%%%%
 %      /b_2 par 1 exc dup mul sub sqrt div def % demi-petit axe
        /a_2 G M1 M2 add mul mu mul Energie div 2 div neg def %
 % excentricité
        /exc 1 K 2 exp a_2 G mul Mt mul div sub sqrt def
% paramètre
%        /par2 a_2 1 exc dup mul sub mul def
% periode
       /periode 6.28 a_2 3 exp G div Mt div sqrt mul def
        }%
%%  0  1   2   3  4  5   6   7
%% x1 y1 x'1 y'1 x2 y2 x'2 y'2
\def\InitCondred{ x0  y0  v0x  v0y}
\def\InitCond{ x01  y01  v0x1  v0y1 x02 y02  v0x2   v0y2}
\psset{method=rk4}
\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=8pt]%
\psequadiff[whichabs=0,whichord=1,
            plotpoints=1000,algebraic,
            tabname=X1Y1]{0}{100}{\InitCond}{\GravAlgIIcorps}
\listplot[unit=1,linecolor=red]{X1Y1 aload pop}
\psequadiff[whichabs=4,whichord=5,
            plotpoints=1000,algebraic,
            tabname=X2Y2]{0}{100}{\InitCond}{\GravAlgIIcorps}
\listplot[unit=1,linecolor=blue]{X2Y2 aload pop}
% mouvement de M2 par rapport à M1
\psequadiff[plotpoints=1000,algebraic,
            plotfuncx=y dup 4 get exch 0 get sub ,
            plotfuncy=dup 5 get exch 1 get sub,
            tabname=XY1]{0}{22.7}{\InitCond}{\GravAlgIIcorps}
\listplot[unit=1,linecolor=green]{XY1 aload pop}
% mouvement de M1 par rapport à M2
\psequadiff[plotpoints=1000,algebraic,
            plotfuncx=y dup 0 get exch 4 get sub ,
            plotfuncy=dup 1 get exch 5 get sub,
            tabname=XY2]{0}{22.7}{\InitCond}{\GravAlgIIcorps}
\listplot[unit=1,linecolor=magenta]{XY2 aload pop}
% mouvement de M1 par rapport à G
%\psequadiff[plotpoints=1000,algebraic,
%            plotfuncx=y 0 get
%                      y 4 get M2 mul
%                      y 0 get M1 mul add
%                      Mt div sub ,
%            plotfuncy=1 get
%                      y 1 get
%                      M1 mul
%                      y 5 get M2 mul add
%                      Mt div sub,
%            tabname=XY3]{0}{22.7}{\InitCond}{\GravAlgIIcorps}
%\listplot[unit=1]{XY3 aload pop}
% mouvement de M2 par rapport à G
%\psequadiff[plotpoints=1000,algebraic,
%            plotfuncx=y 4 get
%                      y 4 get M2 mul
%                      y 0 get M1 mul add
%                      Mt div sub ,
%            plotfuncy=5 get
%                      y 1 get M1 mul
%                      y 5 get M2 mul add
%                     Mt div sub,
%            tabname=XY4]{0}{22.7}{\InitCond}{\GravAlgIIcorps}
%\listplot[unit=1,linecolor=red]{XY4 aload pop}
% centre de masse
\psequadiff[plotpoints=1000,algebraic,
            plotfuncx=y dup 0 get M1 mul exch
                      4 get M2 mul add
                      Mt div,
            plotfuncy=dup 1 get M1 mul exch
                      5 get M2 mul add
                      Mt div,
            tabname=XGYG]{0}{100}{\InitCond}{\GravAlgIIcorps}
\listplot[unit=1,linecolor=cyan]{XGYG aload pop}
\psdots(!x01 y01)(!x02 y02)
\psdot[linecolor=red](!xG yG)
\psline{<->}(7,0)(0,0)(0,8)
\rput(-2,0){\psPrintValue[decimals=4]{periode}\hphantom{00000}s}
\end{pspicture}
\end{center}
%
\newpage
\begin{center}
\def\InitCondred{ xr0 yr0 v0xr  v0yr}
\begin{pspicture}(-7,-6)(7,8)
\pstVerb{
/arccos {
   dup
   dup mul neg 1 add sqrt
   exch
   atan
} def
        /G 1 def
        /x01 2 def
        /y01 2 def
        /v0x1 .05 def
        /v0y1 0.25 def
        /x02 -2 def
        /y02 -0.5 def
        /v0x2 -0.25 def
        /v0y2 -0.5 def
        /M1 3 def
        /M2 1 def
        /Mt M1 M2 add def
        /mu M1 M2 mul M1 M2 add div def % masse réduite
        /xG0 M1 x01 mul M2 x02 mul add M1 M2 add div def
        /yG0 M1 y01 mul M2 y02 mul add M1 M2 add div def
        /x0 x01 x02 sub def
        /y0 y01 y02 sub def
        /theta_0 y0 x0 atan def
        /v0x v0x1 v0x2 sub def
        /v0y v0y1 v0y2 sub def
        /mu M1 M2 mul M1 M2 add div def % masse réduite
        /xG M1 x01 mul M2 x02 mul add M1 M2 add div def
        /yG M1 y01 mul M2 y02 mul add M1 M2 add div def
        /vxG M1 v0x1 mul M2 v0x2 mul add M1 M2 add div def
        /vyG M1 v0y1 mul M2 v0y2 mul add M1 M2 add div def
% conditions initiales pour le point réduit
        /xr0 x01 x02 sub def
        /yr0 y01 y02 sub def
        /r0 xr0 dup mul yr0 dup mul add sqrt def
        /v0xr v0x1 v0x2 sub def
        /v0yr v0y1 v0y2 sub def
        /v0r_2 v0xr dup mul v0yr dup mul add def
%         v0r 2 Mt mul r0 div G mul ge { a faire} if
        /Lc xr0 v0yr mul yr0 v0xr mul sub mu mul def % moment cinetique
        /K xr0 v0yr mul yr0 v0xr mul sub def
        /Energie 0.5 mu mul v0xr dup mul v0yr dup mul add mul % 1/2 mv^2
%                 Lc dup mul 2 div mu div r0 dup mul div
                 G M1 M2 mul mul r0 div sub
         def
%        /par xr0 v0yr mul yr0 v0xr mul sub dup mul G Mt mul div def % paramètre de l'ellipse
%%%%%%%%%%%%%%
 %      /b_2 par 1 exc dup mul sub sqrt div def % demi-petit axe
        /a_2 G M1 M2 add mul mu mul Energie div 2 div neg def %
 % excentricité
        /exc 1 K 2 exp a_2 G mul Mt mul div sub sqrt def
% paramètre
        /par2 a_2 1 exc dup mul sub mul def
% periode
       /periode 6.28 a_2 3 exp G div Mt div sqrt mul def
% phase
      /Phi theta_0 par2 r0 div 1 sub exc div arccos sub def
        }%
%%  0  1   2   3  4  5   6   7
%% x1 y1 x'1 y'1 x2 y2 x'2 y'2
\def\InitCondred{ x0  y0  v0x  v0y}
\def\InitCond{ x01  y01  v0x1  v0y1 x02 y02  v0x2   v0y2}
\psset{method=rk4}
%\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=8pt]%
% mouvement de M1 par rapport à G
\psequadiff[plotpoints=1000,algebraic,
            plotfuncx=y 0 get
                      y 4 get M2 mul
                      y 0 get M1 mul add
                      Mt div sub ,
            plotfuncy=1 get
                      y 1 get
                      M1 mul
                      y 5 get M2 mul add
                      Mt div sub,
            tabname=XY3]{0}{23.6}{\InitCond}{\GravAlgIIcorps}
\listplot[unit=1,linecolor=red]{XY3 aload pop}
% mouvement de M2 par rapport à G
\psequadiff[plotpoints=1000,algebraic,
            plotfuncx=y 4 get
                      y 4 get M2 mul
                      y 0 get M1 mul add
                      Mt div sub ,
            plotfuncy=5 get
                      y 1 get M1 mul
                      y 5 get M2 mul add
                     Mt div sub,
            tabname=XY4]{0}{23.6}{\InitCond}{\GravAlgIIcorps}
\listplot[unit=1,linecolor=blue]{XY4 aload pop}
\psline{<->}(7,0)(0,0)(0,8)
\rput(-2,2){\psPrintValue[decimals=4]{periode}\hphantom{0000000}s}
%\parametricplot[plotpoints=360,linestyle=dotted]{0}{360}{%
%    /radius par2 1 exc t theta_0 sub cos mul add div def
%     radius t cos mul
%     radius t sin mul}%
\psdots(!x01 xG0 sub y01 yG0 sub)(!x02 xG0 sub y02 yG0 sub)
\psline(!x01 xG0 sub y01 yG0 sub)(!x02 xG0 sub y02 yG0 sub)
% mouvement de M (masse réduite) par rapport à G
\psequadiff[plotpoints=1000,algebraic,
            plotfuncx=y 4 get
                      y 4 get M2 mul
                      y 0 get M1 mul add
                      Mt div sub
                      Mt mul M1 div,
            plotfuncy=5 get
                      y 1 get M1 mul
                      y 5 get M2 mul add
                      Mt div sub
                      Mt mul M1 div,
            tabname=XYM]{0}{23.6}{\InitCond}{\GravAlgIIcorps}%
\listplot[unit=1,linecolor=gray!80,linewidth=0.1]{XYM aload pop}
\rput(!x01 xG0 sub y01 yG0 sub){\psline[unit=5,linecolor=red!50]{->}(!v0x1 vxG sub v0y1 vyG sub)}
\rput(!x02 xG0 sub y02 yG0 sub){\psline[unit=5,linecolor=blue!50]{->}(!v0x2 vxG sub v0y2 vyG sub)}
% point "réduit"
\pstVerb{/xr0 x01 xG0 sub Mt mul M2 div neg def
         /yr0 y01 yG0 sub Mt mul M2 div neg def}%
\psdot[linecolor=gray](!xr0 yr0)
\psline(!x01 xG0 sub y01 yG0 sub)(!x02 xG0 sub y02 yG0 sub)(!xr0 yr0)
\uput[d](0,0){$C$}
\uput[r](!x01 xG0 sub y01 yG0 sub){$M_1$}
\uput[dr](!x02 xG0 sub y02 yG0 sub){$M_2$}
\uput[l](!xr0 yr0){$M$}
\rput(!xG0 neg yG0 neg){\psline{<->}(7,0)(0,0)(0,8)\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=8pt]}
\parametricplot[plotpoints=360,linecolor=white]{0}{360}{
    /radius par2 1 exc t Phi sub cos mul add div def
     radius t cos mul neg
     radius t sin mul neg}%
\end{pspicture}
\end{center}
\end{document} 