
\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}
\usepackage{animate}
\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}}}
%timeline
\begin{filecontents}{pb2corps.dat}
::0x0
::1
::2
::3
::4
::5
::6
::7
::8
::9
::10
::11
::12
::13
::14
::15
::16
::17
::18
::19
::20
::21
::22
::23
::24
::25
::26
::27
::28
::29
::30
::31
::32
::33
::34
::35
::36
::37
::38
::39
::40
::41
::42
::43
::44
::45
::46
::47
::48
::49
::50
::51
::52
::53
::54
::55
::56
::57
::58
::59
::60
::61
::62
::63
::64
::65
::66
::67
::68
::69
::70
::71
::72
::73
::74
::75
::76
::77
::78
::79
::80
::81
::82
::83
::84
::85
::86
::87
::88
::89
::90
::91
::92
::93
::94
::95
::96
::97
::98
::99
::100
::101
::102
::103
::104
::105
::106
::107
::108
::109
::110
::111
::112
::113
::114
::115
::116
::117
::118
::119
::120
::121
::122
::123
::124
::125
::126
::127
::128
::129
::130
::131
::132
::133
::134
::135
::136
::137
::138
::139
::140
::141
::142
::143
::144
::145
::146
::147
::148
::149
::150
::151
::152
::153
::154
::155
::156
::157
::158
::159
::160
::161
::162
::163
::164
::165
::166
::167
::168
::169
::170
::171
::172
::173
::174
::175
::176
::177
::178
::179
::180
::181
::182
::183
::184
::185
::186
::187
::188
::189
::190
::191
::192
::193
::194
::195
::196
::197
::198
::199
::200
\end{filecontents}
%%%%%%%%%%%%%%%%%%
\title{Gravitation : le problème des trois corps avec PSTricks\\ partie 1}
\date{7 juillet 2\,012}
\begin{document}
\maketitle
\section{Système de trois corps, dont l'un est beaucoup plus massif que les deux autres}
\subsection{Les équations}
On considère un système de trois corps en interaction gravitationnelle $M$ de masse $m$, $M_1$ de masse $m_1$ et $M_2$ de masse $m_2$. On suppose que la masse de $M$ est très grande par rapport aux deux autres et que le repère $\mathcal{R}$ lié à $M$ est galiléen. Les deux corps $M_1$ et $M_2$ sont considérés comme des satellites de $M$. Pour l'étude des mouvements, on suppose que les corps sont ponctuels.
\begin{center}
\begin{pspicture}(-3,-1)(6,7)
\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]%
\pstVerb{/x01 -3 def /y01 4 def
         /x02  6 def /y02 6 def
         /xr012 x02 x01 sub def
         /yr012 y02 y01 sub def
         /M0 10 def
         /M1 3 def
         /M2 1 def
         /Mt M1 M2 add def
         /xG012 x01 M1 mul x02 M2 mul add Mt div def
         /yG012 y01 M1 mul y02 M2 mul add Mt div def
         }%
\pnode(0,0){O}
\pnode(!x01 y01){M1}
\pnode(!x02 y02){M2}
\pnode(!xG012 yG012){C12} % centre de masse de M1 et M2
\pscircle[fillstyle=solid,fillcolor=yellow!50](0,0){0.5}
\pscircle[fillstyle=solid,fillcolor=red!50](M1){0.21}
\pscircle[fillstyle=solid,fillcolor=blue!50](M2){0.07}
\psline[linestyle=dotted](M1)(M2)
\psline[linestyle=dashed](M2)(O)(M1)
\psline{<->}(6,0)(0,0)(0,7)
\uput{0.6}[dl](O){$M$}
\uput[r](0,6.8){$y$}
\uput[u](5.8,0){$x$}
\uput[l](0,6.8){$\mathcal{R}$}
\uput{0.22}[l](M1){$M_1$}
\uput{0.1}[u](M2){$M_2$}
\rput(M1){\psline[style=vecteurC]{->}(!xr012 10 div yr012 10 div)\uput[u](!xr012 10 div yr012 10 div){$\overrightarrow{F}_{2/1}$}
          \psline[style=vecteurC]{->}(!x01 2.5 div neg y01 2.5 div neg)\uput[r](!x01 2.5 div neg y01 2.5 div neg){$\overrightarrow{F}_{M/1}$}}
\rput(M2){\psline[style=vecteurC]{->}(!xr012 10 div neg yr012 10 div neg)\uput[u](!xr012 10 div neg yr012 10 div neg){$\overrightarrow{F}_{1/2}$}
          \psline[style=vecteurC]{->}(!x02 5 div neg y02 5 div neg)\uput[r](!x02 5 div neg y02 5 div neg){$\overrightarrow{F}_{M/2}$}}
\end{pspicture}
\end{center}
On note $\overrightarrow{r_{12}}=\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_{12}^3}\overrightarrow{r_{12}}\quad \overrightarrow{F}_{2/1}=+\mathcal{G}\frac{m_1m_2}{r_{12}^3}\overrightarrow{r_{12}}
\]
On note $\overrightarrow{r_{1}}=\overrightarrow{MM_1}$ et $\overrightarrow{r_{2}}=\overrightarrow{MM_2}$. $M_1$ subit de la part de $M$ une force attractive $\overrightarrow{F}_{M/1}$ et $M_2$ subit de la part de $M$ une force attractive $\overrightarrow{F}_{M/2}$ :
\[
\overrightarrow{F}_{M/1}=-\mathcal{G}\frac{mm_1}{r^3_1}\overrightarrow{r_1}\quad \overrightarrow{F}_{M/2}=-\mathcal{G}\frac{mm_2}{r^3_2}\overrightarrow{r_2}
\]
La loi de Newton appliquée à chacune des particules $M_1$ et $M_2$ s'écrit :
\[
m_1\frac{\mathrm{d}^2\overrightarrow{r_1}}{\mathrm{d}t^2}=-\mathcal{G}\frac{m_1m_2}{r_{12}^3}\overrightarrow{r_{12}}-
 \mathcal{G}\frac{mm_1}{r^3_1}\overrightarrow{r_1}
\]
\[
m_2\frac{\mathrm{d}^2\overrightarrow{r_2}}{\mathrm{d}t^2}=\mathcal{G}\frac{m_1m_2}{r_{12}^3}\overrightarrow{r_{12}}-
 \mathcal{G}\frac{mm_2}{r^3_2}\overrightarrow{r_2}
\]
Sachant que $r_{12}=\sqrt{(x_{2}-x_{1})^2+(y_{2}-y_{1})^2}$, que $r_1=\sqrt{x_1^2+y_1^2}$ et $r_1=\sqrt{x_2^2+y_2^2}$, nous obtenons un système de quatre équations différentielles :
\[
\left\{
\begin{array}[m]{l}
 \ddot{x_1}=\hphantom{-}\mathcal{G}\displaystyle\frac{m_2}{r^3_{12}}(x_2-x_1)-\mathcal{G}\frac{m}{r^3_1}x_1\\[1em]
 \ddot{y_1}=\hphantom{-}\mathcal{G}\displaystyle\frac{m_2}{r^3_{12}}(y_2-y_1)-\mathcal{G}\frac{m}{r^3_1}y_1\\[1em]
 \ddot{x_2}=-\mathcal{G}\displaystyle\frac{m_1}{r^3_{12}}(x_2-x_1)-\mathcal{G}\frac{m}{r^3_2}x_2\\[1em]
 \ddot{y_2}=-\mathcal{G}\displaystyle\frac{m_1}{r^3_{12}}(y_2-y_1)-\mathcal{G}\frac{m}{r^3_2}y_2
\end{array}
\right.
\]
\subsection{La résolution numérique avec pst-eqdf}
Voici comment placer sur la pile les différentes variables, ainsi que le système d'équations en notation algébrique. Il faut retenir la notation utilisée pour représenter les variables en fonction de leurs positions respectives sur la pile :
\begin{verbatim}
%   0    1    2    3    4    5    6    7
% y[0] y[1] y[2] y[3] y[4] y[5] y[6] y[7]
% x1    y1  x'1  y'1  x2   y2   x'2  y'2
\def\GravAlgIIcorps{%
  y[2]|y[3]|%
  M2*(y[4]-y[0])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5-M*y[0]/(y[0]^2+y[1]^2)^1.5|%
  M2*(y[5]-y[1])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5-M*y[1]/(y[0]^2+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-M*y[4]/(y[4]^2+y[5]^2)^1.5|%
  M1*(y[1]-y[5])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5-M*y[5]/(y[4]^2+y[5]^2)^1.5}
\end{verbatim}
Les conditions initiales, positions et vitesses respectives de $M_1$ et $M_2$, déterminent l'évolution du système.
\begin{verbatim}
\def\InitCond{x01 y01 v0x1 v0y1 x02 y02 v0x2 v0y2}
\end{verbatim}
Elles sont définies dans l'environnement \verb+\begin{pspicture}+ :
\begin{verbatim}
\pstVerb{
        /x01 -1 def /y01 0 def
        /x02  5 def /y02 0 def
        /v0x1 0 def
        /v0y1 25 def
        /v0x2 0 def
        /v0y2 -10 def
        /M 500 def
        /M1 10 def
        /M2 1 def
        }%
\end{verbatim}
%   0    1    2    3    4    5    6    7
% y[0] y[1] y[2] y[3] y[4] y[5] y[6] y[7]
% x1    y1  x'1  y'1  x2   y2   x'2  y'2
\def\GravAlgIIIcorps{%
  y[2]|y[3]|%
  M2*(y[4]-y[0])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5-M0*y[0]/(y[0]^2+y[1]^2)^1.5|%
  M2*(y[5]-y[1])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5-M0*y[1]/(y[0]^2+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-M0*y[4]/(y[4]^2+y[5]^2)^1.5|%
  M1*(y[1]-y[5])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5-M0*y[5]/(y[4]^2+y[5]^2)^1.5}

Quelle durée choisir pour les calculs ? On suppose que les interactions entre les deux satellites sont négligeables, et on calcule la période de chacun d'eux.
\begin{verbatim}
\rput(-2,0){\psPrintValue[decimals=4]{periode1}\hphantom{000000}s}
\rput(-2,1){\psPrintValue[decimals=4]{periode2}\hphantom{000000}s}
\end{verbatim}
La durée choisie qui doit être identique pour les deux, est la plus grande des périodes.
\begin{verbatim}
\psequadiff[whichabs=0,whichord=1,
            plotpoints=1000,algebraic,
            tabname=X1Y1
            ,saveData,filename=IIIXYM1.dat
            ]{0}{periode2}{\InitCond}{\GravAlgIIIcorps}
\listplot[unit=1,linecolor=red]{X1Y1 aload pop}
\psequadiff[whichabs=4,whichord=5,
            plotpoints=1000,algebraic,
            tabname=X2Y2
            ,saveData,filename=IIIXYM2.dat
            ]{0}{periode2}{\InitCond}{\GravAlgIIIcorps}
\listplot[unit=1,linecolor=blue]{X2Y2 aload pop}
\end{verbatim}

%\newpage
\begin{center}
\def\InitCond{ x01  y01  v0x1  v0y1 x02 y02 v0x2  v0y2}
\psset{method=rk4}
\begin{pspicture}(-6,-6)(6,6)
\pstVerb{
/arccos {
   dup
   dup mul neg 1 add sqrt
   exch
   atan
} def
        /G 1 def
        /x01 -1 def /y01 0 def
        /x02  5 def /y02 0 def
        /v0x1 0 def
        /v0y1 25 def
        /v0x2 0 def
        /v0y2 -10 def
        /M0 500 def
        /M1 10 def
        /M2 1 def
        /xG012 M1 x01 mul M2 x02 mul add M1 M2 add div def
        /yG012 M1 y01 mul M2 y02 mul add M1 M2 add div def
% calcul des périodes en supposant les interactions
% entre les 2 satellites négligeables
% satellite 1
        /Cste1 x01 v0y1 mul def % Cste des aires pour 1
        /Energie1 0.5 M1 mul v0y1 2 exp mul % 1/2 mv^2
                 G M1 M0 mul mul x01 abs div sub def
         /par1 Cste1 dup mul M0 div def %
         /exc1 2 Cste1 dup mul mul Energie1 mul M1 div M0 dup mul div 1 add sqrt def
         /a_1 par1 1 exc1 dup mul sub div def % demi-grand axe
         /periode1 6.283185 a_1 3 exp G div M0 div sqrt mul def
% satellite 2
        /Cste2 x02 v0y2 mul def % Cste des aires pour 2
        /Energie2 0.5 M2 mul v0y2 2 exp mul % 1/2 mv^2
                 G M2 M0 mul mul x02 abs div sub def
         /par2 Cste2 dup mul M0 div def %
         /exc2 2 Cste2 dup mul mul Energie2 mul M2 div M0 dup mul div 1 add sqrt def
         /a_2 par2 1 exc2 dup mul sub div def % demi-grand axe
         /periode2 6.283185 a_2 3 exp G div M0 div sqrt mul def
         }%
\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=8pt](-6,-6)(6,6)%
\psequadiff[whichabs=0,whichord=1,
            plotpoints=1000,algebraic,
            tabname=X1Y1
            ,saveData,filename=IIIXYM1.dat
            ]{0}{periode2}{\InitCond}{\GravAlgIIIcorps}
\listplot[unit=1,linecolor=red]{X1Y1 aload pop}
\psequadiff[whichabs=4,whichord=5,
            plotpoints=1000,algebraic,
            tabname=X2Y2
            ,saveData,filename=IIIXYM2.dat
            ]{0}{periode2}{\InitCond}{\GravAlgIIIcorps}
\listplot[unit=1,linecolor=blue]{X2Y2 aload pop}
\pnode(!x01 y01){M01}
\pnode(!x02 y02){M02}
\pscircle[fillstyle=solid,fillcolor=red!50](0,0){0.25}
\pscircle[fillstyle=solid,fillcolor=gray!50](M01){0.15}
\pscircle[fillstyle=solid,fillcolor=blue!50](M02){0.07}
\rput(M01){\psline[unit=0.2,linecolor=red]{->}(!v0x1 v0y1)}
\rput(M02){\psline[unit=0.2,linecolor=blue]{->}(!v0x2 v0y2)}
\rput(-2,0){\psPrintValue[decimals=4]{periode1}\hphantom{000000}s}
\rput(-2,1){\psPrintValue[decimals=4]{periode2}\hphantom{000000}s}
\end{pspicture}
\end{center}
%\newpage
\subsection{L'animation}
On peut voir sur l'animation, que les périodes précédemment calculées ne sont qu'approximatives, ainsi le satellite le plus éloigné fait un peu plus d'un tour.

On remarquera d'autre part que, sur la durée choisie, les trajectoires sont bouclées : le système des 3 corps est stable avec les conditions initiales choisies.

\begin{center}
%
\def\nFrames{200}% 200 images
\begin{animateinline}[controls,timeline=pb2corps.dat,%loop,
                     begin={\begin{pspicture}(-6,-6)(6,6)},
                     end={\end{pspicture}}]{10}% 10 images/s
\pstVerb{/XY1 [(IIIXYM1.dat) run] def
         /XY2 [(IIIXYM2.dat) run] def}
\psframe*[linecolor={[cmyk]{1 1 0 0.7}}](-6,-6)(6,6)
\psgrid[subgriddiv=0,gridcolor=white,griddots=10,gridlabels=0pt]%
\listplot[linecolor=gray]{XY1 aload pop}
\listplot[linecolor=gray]{XY2 aload pop}
\pscircle[fillstyle=solid,fillcolor=red!50](0,0){0.4}
\newframe
\multiframe{\nFrames}{i=0+10}{% 1 point sur 10
\pstVerb{/X1 XY1 \i\space get def
         /Y1 XY1 \i\space 1 add get def
         /X2 XY2 \i\space get def
         /Y2 XY2 \i\space 1 add get def
        }%
\pscircle[fillstyle=solid,fillcolor=gray!80](!X1 Y1){0.2}
\pscircle[fillstyle=solid,fillcolor=yellow](!X2 Y2){0.1}
}
\end{animateinline}
\end{center}
\newpage
\section{Mouvements relatifs}
\subsection{Mouvements vus de $M_1$}
Les coordonnées de $M_2$ rapportées à $M_1$  sont :
\[
\left\{
\begin{array}[m]{l}
 x_{2/1}=x_2-x_1\\[1em]
 y_{2/1}=y_2-y_1\\[1em]
\end{array}
\right.
\]
Celles de $M$ sont $(-x_1,-y_1)$.
\begin{verbatim}
%   0    1    2    3    4    5    6    7
% y[0] y[1] y[2] y[3] y[4] y[5] y[6] y[7]
% x1    y1  x'1  y'1  x2   y2   x'2  y'2
\psequadiff[plotpoints=1000,algebraic,
            plotfuncx=y dup 0 get exch
                      4 get sub,
            plotfuncy=dup 1 get exch
                      5 get sub,
            tabname=XY21,
%            ,saveData,filename=IIIM21.dat
            ]{0}{periode2}{\InitCond}{\GravAlgIIIcorps}
\psequadiff[plotpoints=1000,algebraic,
            plotfuncx=y 0 get neg,
            plotfuncy=1 get neg,
            tabname=XY01
%            ,saveData,filename=IIIXY01.dat
            ]{0}{periode2}{\InitCond}{\GravAlgIIIcorps}
\end{verbatim}
\begin{center}
\def\InitCond{ x01  y01  v0x1  v0y1 x02 y02 v0x2  v0y2}
\psset{method=rk4}
\begin{pspicture}(-7,-6)(6,6)
\pstVerb{
        /G 1 def
        /x01 -1 def /y01 0 def
        /x02  5 def /y02 0 def
        /v0x1 0 def
        /v0y1 25 def
        /v0x2 0 def
        /v0y2 -10 def
        /M0 500 def
        /M1 10 def
        /M2 1 def
        /xG012 M1 x01 mul M2 x02 mul add M1 M2 add div def
        /yG012 M1 y01 mul M2 y02 mul add M1 M2 add div def
% calcul des périodes en supposant les interactions
% entre les 2 satellites négligeables
% satellite 1
        /Cste1 x01 v0y1 mul def % Cste des aires pour 1
        /Energie1 0.5 M1 mul v0y1 2 exp mul % 1/2 mv^2
                 G M1 M0 mul mul x01 abs div sub def
         /par1 Cste1 dup mul M0 div def %
         /exc1 2 Cste1 dup mul mul Energie1 mul M1 div M0 dup mul div 1 add sqrt def
         /a_1 par1 1 exc1 dup mul sub div def % demi-grand axe
         /periode1 6.283185 a_1 3 exp G div M0 div sqrt mul def
% satellite 2
        /Cste2 x02 v0y2 mul def % Cste des aires pour 2
        /Energie2 0.5 M2 mul v0y2 2 exp mul % 1/2 mv^2
                 G M2 M0 mul mul x02 abs div sub def
         /par2 Cste2 dup mul M0 div def %
         /exc2 2 Cste2 dup mul mul Energie2 mul M2 div M0 dup mul div 1 add sqrt def
         /a_2 par2 1 exc2 dup mul sub div def % demi-grand axe
         /periode2 6.283185 a_2 3 exp G div M0 div sqrt mul def
         }%
\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt](-7,-6)(6,6)%
\psequadiff[plotpoints=1000,algebraic,
            plotfuncx=y dup 4 get exch
                      0 get sub,
            plotfuncy=dup 5 get exch
                      1 get sub,
            tabname=XY21
%            ,saveData,filename=IIIM21.dat
            ]{0}{periode2}{\InitCond}{\GravAlgIIIcorps}
\listplot[unit=1,linecolor=blue]{XY21 aload pop}
\psequadiff[plotpoints=1000,algebraic,
            plotfuncx=y 0 get neg,
            plotfuncy=1 get neg,
            tabname=XYM01
%            ,saveData,filename=IIIXYM01.dat
            ]{0}{periode2}{\InitCond}{\GravAlgIIIcorps}
\listplot[unit=1,linecolor=red]{XYM01 aload pop}
\pnode(!x02 x01 sub y02 y01 sub){M021}
\pnode(!x01 neg y01 neg){M001}
\pscircle[fillstyle=solid,fillcolor=gray!80](0,0){0.2}
\pscircle[fillstyle=solid,fillcolor=red!50](M001){0.25}
\pscircle[fillstyle=solid,fillcolor=blue](M021){0.07}
\rput(M021){\psline[unit=0.1,linecolor=blue,arrowinset=0.2,arrowsize=2]{->}(!v0x2 v0x1 sub v0y2 v0y1 sub)}
\rput(M001){\psline[unit=0.1,linecolor=red,arrowinset=0.2,arrowsize=2]{->}(!v0x1 neg v0y1 neg)}
\rput(-3,0){\psPrintValue[decimals=4]{periode1}\hphantom{000000}s}
\rput(-3,1){\psPrintValue[decimals=4]{periode2}\hphantom{000000}s}
\end{pspicture}
\end{center}
\subsection{Mouvements vus de $M_2$}
Les coordonnées de $M_1$ rapportées à $M_2$  sont :
\[
\left\{
\begin{array}[m]{l}
 x_{1/2}=x_1-x_2\\[1em]
 y_{1/2}=y_1-y_2\\[1em]
\end{array}
\right.
\]
Celles de $M$ sont $(-x_2,-y_2)$.
\begin{verbatim}
%   0    1    2    3    4    5    6    7
% y[0] y[1] y[2] y[3] y[4] y[5] y[6] y[7]
% x1    y1  x'1  y'1  x2   y2   x'2  y'2
\psequadiff[plotpoints=1000,algebraic,
            plotfuncx=y dup 4 get exch
                      0 get sub,
            plotfuncy=dup 5 get exch
                      1 get sub,
            tabname=XY12,
%            ,saveData,filename=IIIM12.dat
            ]{0}{periode2}{\InitCond}{\GravAlgIIIcorps}
\psequadiff[plotpoints=1000,algebraic,
            plotfuncx=y 4 get neg,
            plotfuncy=5 get neg,
            tabname=XYM02
%            ,saveData,filename=IIIXYM02.dat
            ]{0}{periode2}{\InitCond}{\GravAlgIIIcorps}
\end{verbatim}
\begin{center}
\def\InitCond{ x01  y01  v0x1  v0y1 x02 y02 v0x2  v0y2}
\psset{method=rk4}
\begin{pspicture}(-7,-6)(6,6)
\pstVerb{
        /G 1 def
        /x01 -1 def /y01 0 def
        /x02  5 def /y02 0 def
        /v0x1 0 def
        /v0y1 25 def
        /v0x2 0 def
        /v0y2 -10 def
        /M0 500 def
        /M1 10 def
        /M2 1 def
        /xG012 M1 x01 mul M2 x02 mul add M1 M2 add div def
        /yG012 M1 y01 mul M2 y02 mul add M1 M2 add div def
% calcul des périodes en supposant les interactions
% entre les 2 satellites négligeables
% satellite 1
        /Cste1 x01 v0y1 mul def % Cste des aires pour 1
        /Energie1 0.5 M1 mul v0y1 2 exp mul % 1/2 mv^2
                 G M1 M0 mul mul x01 abs div sub def
         /par1 Cste1 dup mul M0 div def %
         /exc1 2 Cste1 dup mul mul Energie1 mul M1 div M0 dup mul div 1 add sqrt def
         /a_1 par1 1 exc1 dup mul sub div def % demi-grand axe
         /periode1 6.283185 a_1 3 exp G div M0 div sqrt mul def
% satellite 2
        /Cste2 x02 v0y2 mul def % Cste des aires pour 2
        /Energie2 0.5 M2 mul v0y2 2 exp mul % 1/2 mv^2
                 G M2 M0 mul mul x02 abs div sub def
         /par2 Cste2 dup mul M0 div def %
         /exc2 2 Cste2 dup mul mul Energie2 mul M2 div M0 dup mul div 1 add sqrt def
         /a_2 par2 1 exc2 dup mul sub div def % demi-grand axe
         /periode2 6.283185 a_2 3 exp G div M0 div sqrt mul def
         }%
\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt](-7,-6)(6,6)%
\psequadiff[plotpoints=1000,algebraic,
            plotfuncx=y dup 0 get exch
                      4 get sub,
            plotfuncy=dup 1 get exch
                      5 get sub,
            tabname=XYM12
            ,saveData,filename=IIIXYM12.dat
            ]{0}{periode2}{\InitCond}{\GravAlgIIIcorps}
\listplot[unit=1,linecolor=blue]{XYM12 aload pop}
\psequadiff[plotpoints=1000,algebraic,
            plotfuncx=y 4 get neg,
            plotfuncy=5 get neg,
            tabname=XYM02
            ,saveData,filename=IIIXYM02.dat
            ]{0}{periode2}{\InitCond}{\GravAlgIIIcorps}
\listplot[unit=1,linecolor=red]{XYM02 aload pop}
\pnode(!x01 x02 sub y01 y02 sub){M012}
\pnode(!x02 neg y02 neg){M002}
\pscircle[fillstyle=solid,fillcolor=gray!80](0,0){0.2}
\pscircle[fillstyle=solid,fillcolor=red!50](M002){0.25}
\pscircle[fillstyle=solid,fillcolor=blue](M012){0.07}
\rput(M012){\psline[unit=0.1,linecolor=blue,arrowinset=0.2,arrowsize=2]{->}(!v0x1 v0x2 sub v0y1 v0y2 sub)}
\rput(M002){\psline[unit=0.1,linecolor=red,arrowinset=0.2,arrowsize=2]{->}(!v0x2 neg v0y2 neg)}
\rput(-3,0){\psPrintValue[decimals=4]{periode1}\hphantom{000000}s}
\rput(-3,1){\psPrintValue[decimals=4]{periode2}\hphantom{000000}s}
\end{pspicture}
\end{center}
\subsection{Animation}
\begin{center}
%
\def\nFrames{200}% 200 images
\begin{animateinline}[controls,timeline=pb2corps.dat,%loop,
                     begin={\begin{pspicture}(-7,-6)(7,6)},
                     end={\end{pspicture}}]{10}% 10 images/s
\pstVerb{/XY1 [(IIIXYM12.dat) run] def
         /XY2 [(IIIXYM02.dat) run] def}
\psframe*[linecolor={[cmyk]{1 1 0 0.7}}](-7,-6)(7,6)
\psgrid[subgriddiv=0,gridcolor=white,griddots=10,gridlabels=0pt]%
\listplot[linecolor=gray]{XY1 aload pop}
\listplot[linecolor=gray]{XY2 aload pop}
\pscircle[fillstyle=solid,fillcolor=gray!80](0,0){0.2}
\newframe
\multiframe{\nFrames}{i=0+10}{% 1 point sur 10
\pstVerb{/X1 XY1 \i\space get def
         /Y1 XY1 \i\space 1 add get def
         /X2 XY2 \i\space get def
         /Y2 XY2 \i\space 1 add get def
        }%
\pscircle[fillstyle=solid,fillcolor=yellow](!X1 Y1){0.1}
\pscircle[fillstyle=solid,fillcolor=red!50](!X2 Y2){0.4}
}
\end{animateinline}
\end{center}
\end{document} 