
\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}
\usepackage{wasysym}
\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}}}
\makeatletter
%% adapté de \psRandom du package pstricks-add
%% pour rendre aléatoire la taille des étoiles
%% Manuel Luque
\newdimen\pssizeStar
\def\psset@sizeStar#1{\pssetlength\pssizeStar{#1}}
\psset@sizeStar{1pt}
\define@key[psset]{pst-eqd}{randomPoints}[1000]{\def\psk@randomPoints{#1}}
\psset[pst-eqd]{randomPoints=1000}
\define@boolkey[psset]{pst-eqd}[Pst@]{color}[true]{}
\psset[pst-eqd]{color=false}
\def\psRandomStar{\pst@object{psRandomStar}}%
\def\psRandomStar@i{\@ifnextchar({\psRandomStar@ii}{\psRandomStar@iii(0,0)(1,1)}}
\def\psRandomStar@ii(#1){\@ifnextchar({\psRandomStar@iii(#1)}{\psRandomStar@iii(0,0)(#1)}}
\def\psRandomStar@iii(#1)(#2)#3{%
  \def\pst@tempA{#3}%
  \ifx\pst@tempA\pst@empty\psclip{\psframe(#2)}\else\psclip{#3}\fi
  \pst@getcoor{#1}\pst@tempA
  \pst@getcoor{#2}\pst@tempB
  \begin@SpecialObj
  \addto@pscode{
    \pst@tempA\space /yMin exch def
    /xMin exch def
    \pst@tempB\space /yMax exch def
    /xMax exch def
    /dy yMax yMin sub def
    /dx xMax xMin sub def
    rrand srand                 % initializes the random generator
    /getRandReal { rand 2147483647 div } def
     \psk@randomPoints {
    /DS \pst@number\pssizeStar\space getRandReal mul def
    \@nameuse{psds@\psk@dotstyle}
     \ifPst@color getRandReal 1 1 sethsbcolor \fi
     getRandReal dx mul xMin add
     getRandReal dy mul yMin add
     Dot
     \ifx\psk@fillstyle\psfs@solid fill \fi stroke
    } repeat
  }%
  \end@SpecialObj
  \endpsclip
  \ignorespaces
}
\makeatother
%%%%%%%%%%%%%%%%%%
%timeline
\begin{filecontents}{kepler16.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 : une planète à deux soleils comme Tatooine celle de Luke Skywalker dans \textit{La Guerre des étoiles}}
\date{10 juillet 2\,012}
\begin{document}
\maketitle
\section{La représentation avec les données de la NASA}
Ceci est une tentative pour essayer de schématiser une planète orbitant autour d'une étoile binaire, comme Kepler-16b. On parle dans ce cas-là d'une planète \textit{circumbinaire}.
Le site de la NASA dédié à cette planète\footnote{http://kepler.nasa.gov/Mission/discoveries/kepler16b/}, fournit un grand nombre de renseignements sur les étoiles et la planète, ce qui permet de reconstituer leurs trajectoires respectives. Voici les caractéristiques utiles pour le schéma et l'animation.
\begin{center}
\begin{tabular}{|ll|}
\hline
\textbf{ Paramètres}&\textbf{Valeurs}\\ \hline
  \textit{Étoile A}&\\ \hline
  Masse, $M_A(M_{\astrosun})$&0.6897\\
  Rayon, $R_A(R_{\astrosun})$&0.6489\\ \hline
  \textit{Étoile B}&\\ \hline
  Masse, $M_B(M_{\astrosun})$&0.20255\\
  Rayon, $R_B(R_{\astrosun})$&0.22623\\ \hline
  \textit{Planète b}&\\ \hline
  Masse, $M_b(M_{\jupiter})$&0.333\\
  Rayon, $R_b(R_{\jupiter})$&0.7538\\ \hline
  \textit{Orbite de l'étoile binaire}&\\ \hline
  Période(jours)&41.076\\
  Demi-grand axe(ua) $\mathrm{a_1}$&0.22431\\
  Excentricité $\mathrm{e_1}$&0.15944\\
  Argument du périastre(deg) $\omega_1$&263.464\\ \hline
  \textit{Orbite de la planète circumbinaire} &\\ \hline
   Période(jours)&228.776\\
  Demi-grand axe(ua)$\mathrm{a_2}$&0.7048\\
  Excentricité $\mathrm{e_2}$&0.0069\\
  Argument du périastre(deg) $\omega_2$&318\\ \hline
\end{tabular}
\end{center}

\begin{center}
\begin{pspicture}(-5,-5)(5,5)
\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]%
\pstVerb{% pour le système binaire d'étoiles
        /a1 0.22431 def % semi-major axis
        /e1 0.15944 def % eccentricity
        /w1 263.464 def % argument of periaspe
        /MA 0.6897 def  % masse of star A
        /MB 0.20255 def % masse of star B
        /Mt MA MB add def % masse totale
        /P1 41.076 def  % period (day)
        % pour la planète
        /Mb 0.3333 def  % masse of planet b (M_Jupiter)
        /a2 0.7048 def  % semi-major axis
        /e2 0.0069 def  % eccentricity
        /w2 318 def  % argument of periaspe
        /P2 228.776 def % period (day)
% on en déduit le paramètre de chacune des ellipses
        /p1 a1 1 e1 dup mul add mul def
        /p2 a2 1 e2 dup mul add mul def
        }%
\pnode(!/radius p1 1 e1 w1 w1 sub cos mul add div def
        /radius1 radius MB Mt div mul neg def
         radius1 w1 cos mul 5 mul
         radius1 w1 sin mul 5 mul){MA0}
\parametricplot[plotpoints=360,unit=5,linecolor=red]{0}{360}{
    /radius p1 1 e1 t w1 sub cos mul add div def
    /radius1 radius MB Mt div mul neg def
     radius1 t cos mul
     radius1 t sin mul}
\pscircle[fillstyle=solid,fillcolor=yellow!50](MA0){0.3}
\parametricplot[plotpoints=360,unit=5,linecolor=green]{0}{360}{
    /radius p1 1 e1 t w1 sub cos mul add div def
    /radius2 radius MA Mt div mul def
     radius2 t cos mul
     radius2 t sin mul}
\pnode(!/radius p1 1 e1 w1 w1 sub cos mul add div def
        /radius1 radius MA Mt div mul def
        radius1 w1 cos mul 5 mul
        radius1 w1 sin mul 5 mul){MB0}
\pscircle[fillstyle=solid,fillcolor=red](MB0){0.15}
\parametricplot[plotpoints=360,unit=5,linecolor=blue]{0}{360}{
    /radius p2 1 e2 t w2 sub cos mul add div def
     radius t cos mul
     radius t sin mul}
\pnode(!/radius p2 1 e2 w2 w2 sub cos mul add div def
        radius w2 cos mul 5 mul
        radius w2 sin mul 5 mul){Mb0}
\pscircle[fillstyle=solid,fillcolor=blue](Mb0){0.07}
\end{pspicture}
\end{center}
\section{La simulation avec PSTricks}
\begin{center}
\begin{pspicture}(-3,-1)(6,7)
\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]%
\pstVerb{/x01 -3 def /y01 0 def
         /x02  5 def /y02 0 def
         /xr012 x02 x01 sub def
         /yr012 y02 y01 sub def
         /x0 2 def /y0 5 def
         /xr002 x02 x0 sub def
         /yr002 y02 y0 sub def
         /xr001 x01 x0 sub def
         /yr001 y01 y0 sub def
         /M0 0.2 def
         /M1 0.07 def
         /M2 3.2e-4 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(!x0 y0){M0}
\pnode(!xG012 yG012){C12} % centre de masse de M1 et M2
\pscircle[fillstyle=solid,fillcolor=yellow!50](M2){0.21}
\pscircle[fillstyle=solid,fillcolor=red!50](M1){0.5}
\pscircle[fillstyle=solid,fillcolor=blue!50](M0){0.07}
\psline[linestyle=dotted](M1)(M2)
\psline[linestyle=dashed](M2)(M0)(M1)
\psline{<->}(6,0)(0,0)(0,7)
\uput[u](M0){$M_0$}
\uput[r](0,6.8){$y$}
\uput[u](5.8,0){$x$}
\uput[l](0,6.8){$\mathcal{R}$}
\uput{0.6}[d](M1){$M_1$}
\uput{0.25}[d](M2){$M_2$}
\rput(M1){\psline[style=vecteurC]{->}(!xr012 3 div yr012 3 div)\uput[d](!xr012 3 div yr012 3 div){$\overrightarrow{F}_{2/1}$}
          \psline[style=vecteurC]{->}(!xr001 5 div neg yr001 5 div neg)\uput[r](!xr001 5 div neg yr001 5 div neg){$\overrightarrow{F}_{0/1}$}}
\rput(M2){\psline[style=vecteurC]{->}(!xr012 3 div neg yr012 3 div neg)\uput[d](!xr012 3 div neg yr012 3 div neg){$\overrightarrow{F}_{1/2}$}
             \psline[style=vecteurC]{->}(!xr002 3 div neg yr002 3 div neg)\uput[r](!xr002 3 div neg yr002 3 div neg){$\overrightarrow{F}_{0/2}$}}
\rput(M0){\psline[style=vecteurC]{->}(!xr002 3 div yr002 3 div)\uput[r](!xr002 3 div yr002 3 div){$\overrightarrow{F}_{2/0}$}}
\rput(M0){\psline[style=vecteurC]{->}(!xr001 5 div yr001 5 div)\uput[r](!xr001 5 div yr001 5 div){$\overrightarrow{F}_{1/0}$}}
\end{pspicture}
\end{center}
Pour cela, on considère un système de trois corps en interaction gravitationnelle : l'étoile $M_1$ de masse $m_1$, l'étoile $M_2$ de masse $m_2$ constituant l'étoile binaire et une planète $P$ notée $M_0$, de masse $m$ orbitant autour de cette étoile double.

Avec : $\overrightarrow{r_{01}}=\overrightarrow{M_0M_1}$, $\overrightarrow{r_{02}}=\overrightarrow{M_0M_2}$ et $\overrightarrow{r_{12}}=\overrightarrow{M_1M_2}$ on pose :
\begin{itemize}
  \item $\overrightarrow{r}_{01}=\overrightarrow{r}_1-\overrightarrow{r}_0$
  \item $\overrightarrow{r}_{02}=\overrightarrow{r}_2-\overrightarrow{r}_0$
  \item $\overrightarrow{r}_{12}=\overrightarrow{r}_2-\overrightarrow{r}_1$
\end{itemize}
Chacune des forces s'exprime par :
\[
\left\{
\begin{array}[m]{l}
\overrightarrow{F}_{1/0}=\mathcal{G}\dfrac{m_1m_0}{r_{10}^3}\overrightarrow{r_{01}}=-\overrightarrow{F}_{0/1}\\[1em]
\overrightarrow{F}_{2/0}=\mathcal{G}\dfrac{m_2m_0}{r_{20}^3}\overrightarrow{r_{02}}=-\overrightarrow{F}_{0/2}\\[1em]
\overrightarrow{F}_{2/1}=\mathcal{G}\dfrac{m_1m_2}{r_{12}^3}\overrightarrow{r_{12}}=-\overrightarrow{F}_{1/2}
\end{array}
\right.
\]
L'application de la loi de Newton à chacun des corps donne :
\[
\left\{
\begin{array}[m]{l}
m_0\dfrac{\mathrm{d^2}\overrightarrow{r_0}}{\mathrm{d}t^2}=\hphantom{-}\mathcal{G}\dfrac{m_1m_0}{r_{10}^3}\overrightarrow{r_{01}}+
                                          \mathcal{G}\dfrac{m_2m_0}{r_{20}^3}\overrightarrow{r_{02}}\\[1em]
m_1\dfrac{\mathrm{d^2}\overrightarrow{r_1}}{\mathrm{d}t^2}=-\mathcal{G}\dfrac{m_0m_1}{r_{10}^3}\overrightarrow{r_{01}}+
                                          \mathcal{G}\dfrac{m_1m_2}{r_{12}^3}\overrightarrow{r_{12}}\\[1em]
m_2\dfrac{\mathrm{d^2}\overrightarrow{r_2}}{\mathrm{d}t^2}=-\mathcal{G}\dfrac{m_0m_2}{r_{20}^3}\overrightarrow{r_{02}}-
                                          \mathcal{G}\dfrac{m_1m_2}{r_{12}^3}\overrightarrow{r_{12}}
\end{array}
\right.
\]
Ce qui conduit à un système de 6 équations différentielles :
\[
\left\{
\begin{array}[m]{l}
\ddot{x_0}=\hphantom{-}\mathcal{G}\dfrac{m_1}{r_{10}^3}(x_1-x_0)+\mathcal{G}\dfrac{m_2}{r_{20}^3}(x_2-x_0)\\[1em]
\ddot{y_0}=\hphantom{-}\mathcal{G}\dfrac{m_1}{r_{10}^3}(y_1-y_0)+\mathcal{G}\dfrac{m_2}{r_{20}^3}(y_2-y_0)\\[1em]
\ddot{x_1}=-\mathcal{G}\dfrac{m_0}{r_{10}^3}(x_1-x_0)+\mathcal{G}\dfrac{m_2}{r_{12}^3}(x_2-x_1)\\[1em]
\ddot{y_1}=-\mathcal{G}\dfrac{m_0}{r_{10}^3}(y_1-y_0)+\mathcal{G}\dfrac{m_2}{r_{12}^3}(y_2-y_1)\\[1em]
\ddot{x_2}=-\mathcal{G}\dfrac{m_0}{r_{10}^3}(x_2-x_0)-\mathcal{G}\dfrac{m_1}{r_{12}^3}(x_2-x_1)\\[1em]
\ddot{y_2}=-\mathcal{G}\dfrac{m_0}{r_{02}^3}(y_2-y_0)-\mathcal{G}\dfrac{m_1}{r_{12}^3}(y_2-y_1)
\end{array}
\right.
\]

\begin{verbatim}
%   0    1    2    3    4    5    6    7    8   9     10    11
% y[0] y[1] y[2] y[3] y[4] y[5] y[6] y[7] y[8] y[9] y[10] y[11]
%  x0   y0  x'0  y'0   x1   y1  x'1  y'1   x2   y2   x'2   y'2
\def\GravAlgIIIcorps{%
  y[2]|y[3]|%
  M1*(y[4]-y[0])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5+M2*(y[8]-y[0])/((y[8]-y[0])^2+(y[9]-y[1])^2)^1.5|%
  M1*(y[5]-y[1])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5+M2*(y[9]-y[1])/((y[8]-y[0])^2+(y[9]-y[1])^2)^1.5|%
  y[6]|y[7]|%
  -M0*(y[4]-y[0])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5+M2*(y[8]-y[4])/((y[8]-y[4])^2+(y[9]-y[5])^2)^1.5|%
  -M0*(y[5]-y[1])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5+M2*(y[9]-y[5])/((y[8]-y[4])^2+(y[9]-y[5])^2)^1.5|%
  y[10]|y[11]|%
  -M0*(y[8]-y[0])/((y[8]-y[0])^2+(y[9]-y[1])^2)^1.5-M1*(y[8]-y[4])/((y[8]-y[4])^2+(y[9]-y[5])^2)^1.5|%
  -M0*(y[9]-y[1])/((y[8]-y[0])^2+(y[9]-y[1])^2)^1.5-M1*(y[9]-y[5])/((y[8]-y[4])^2+(y[9]-y[5])^2)^1.5}
\end{verbatim}


\newpage
% y[0] y[1] y[2] y[3] y[4] y[5] y[6] y[7] y[8] y[9] y[10] y[11]
%  x0   y0  x'0  y'0   x1   y1  x'1  y'1   x2   y2   x'2   y'2
\def\GravAlgIIIcorps{%
  y[2]|y[3]|%
  M1*(y[4]-y[0])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5+M2*(y[8]-y[0])/((y[8]-y[0])^2+(y[9]-y[1])^2)^1.5|%
  M1*(y[5]-y[1])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5+M2*(y[9]-y[1])/((y[8]-y[0])^2+(y[9]-y[1])^2)^1.5|%
  y[6]|y[7]|%
  -M0*(y[4]-y[0])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5+M2*(y[8]-y[4])/((y[8]-y[4])^2+(y[9]-y[5])^2)^1.5|%
  -M0*(y[5]-y[1])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5+M2*(y[9]-y[5])/((y[8]-y[4])^2+(y[9]-y[5])^2)^1.5|%
  y[10]|y[11]|%
  -M0*(y[8]-y[0])/((y[8]-y[0])^2+(y[9]-y[1])^2)^1.5-M1*(y[8]-y[4])/((y[8]-y[4])^2+(y[9]-y[5])^2)^1.5|%
  -M0*(y[9]-y[1])/((y[8]-y[0])^2+(y[9]-y[1])^2)^1.5-M1*(y[9]-y[5])/((y[8]-y[4])^2+(y[9]-y[5])^2)^1.5}
%

\begin{center}
\def\InitCond{x00 y00 v0x0 v0y0 x01 y01 v0x1 v0y1 x02 y02 v0x2 v0y2}
\psset{method=rk4}
\begin{pspicture}(-8,-8)(8,8)
\psframe*[linecolor={[cmyk]{1 1 0 0.7}}](-8,-8)(8,8)
\pstVerb{
/arccos {
   dup
   dup mul neg 1 add sqrt
   exch
   atan
} def
        /G 1 def
        /w2 55 def
        /x00 7 w2 cos mul def /y00 7 w2 sin mul def
        /x01 0.4962 def
        /y01 0 def
        /v0x1 0 def
        /v0y1 .135 def
        /x02 -1.69 def
        /y02 0 def
        /v0x2 0 def
        /v0y2 -.46 def
        /v0x0 0.36 w2 sin mul neg def
        /v0y0 0.36 w2 cos mul def
         /M1 0.7 def
         /M2 0.2 def
         /M0 3e-4 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
         }%
\psgrid[subgriddiv=0,gridcolor=white,griddots=10,gridlabels=0pt](-8,-8)(8,8)%
\psequadiff[plotpoints=1000,algebraic,
            plotfuncx=y 4 get
                      y 4 get M1 mul
                      y 8 get M2 mul add
                      M1 M2 add div
                      sub,
            plotfuncy=5 get
                      dup M1 mul
                      y 9 get M2 mul add
                      M1 M2 add div
                      sub,
            tabname=XAYA,
            saveData,filename=XAYA.dat
%            ]{0}{18}{\InitCond}{\GravAlgIIIcorps}
            ]{0}{120}{\InitCond}{\GravAlgIIIcorps}
\listplot[unit=1,linecolor=yellow]{XAYA aload pop}
\psequadiff[plotpoints=1000,algebraic,
            plotfuncx=y 8 get
                      y 4 get M1 mul
                      y 8 get M2 mul add
                      M1 M2 add div
                      sub,
            plotfuncy=9 get
                      dup M2 mul
                      y 5 get M1 mul add
                      M1 M2 add div
                      sub,
            tabname=XBYB,
            saveData,filename=XBYB.dat
%            ]{0}{18}{\InitCond}{\GravAlgIIIcorps}
            ]{0}{120}{\InitCond}{\GravAlgIIIcorps}
\listplot[unit=1,linecolor=red]{XBYB aload pop}
\psequadiff[plotpoints=1000,algebraic,
            plotfuncx=y 0 get
                      y 4 get M1 mul
                      y 8 get M2 mul add
                      M1 M2 add div
                      sub,
            plotfuncy=1 get
                      y 5 get M1 mul
                      y 9 get M2 mul add
                      M1 M2 add div
                      sub,
            tabname=XbYb,
            saveData,filename=XPYP.dat
            ]{0}{120}{\InitCond}{\GravAlgIIIcorps}
\listplot[unit=1,linecolor=blue]{XbYb aload pop}
\pnode(!x01 xG012 sub y01 yG012 sub){M01}
\pnode(!x02 xG012 sub y02 yG012 sub){M02}
\pnode(!x00  xG012 sub y00  yG012 sub){M00}
\pscircle*[linecolor=yellow](M01){0.4}
\pscircle*[linecolor=red](M02){0.16}
\pscircle*[linecolor=white](M00){0.07}
%\rput(M01){\psline[unit=2,linecolor=red]{->}(!v0x1 v0y1)}
%\rput(M02){\psline[unit=2,linecolor=blue]{->}(!v0x2 v0y2)}
\psdot(!xG012 yG012)
\end{pspicture}
\end{center}
\section{Animation avec pst-eqdf et animate}
\begin{center}
\def\nFrames{200}% 200 images
\begin{animateinline}[controls,timeline=kepler16.dat,loop,%
                     begin={\begin{pspicture}(-8,-8)(8,8)},
                     end={\end{pspicture}}]{5}% 5 images/s
\pstVerb{/XY1 [(XAYA.dat) run] def
         /XY2 [(XBYB.dat) run] def
         /XY3 [(XPYP.dat) run] def
         }%
\psframe*[linecolor={[cmyk]{1 1 0 0.7}}](-8,-8)(8,8)
\psRandomStar[linecolor={[rgb]{1,1,0.5}},
    randomPoints=1000,sizeStar=1pt](-8,-8)(8,8){\psframe[linestyle=none](-8,-8)(8,8)}
%\listplot[linecolor=gray]{XY1 aload pop}
%\listplot[linecolor=gray]{XY2 aload pop}
%\listplot[linecolor=gray]{XY3 aload pop}
\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
         /X3 XY3 \i\space get def
         /Y3 XY3 \i\space 1 add get def
        }%
%\psdot(!X1 Y1)
\pscircle*[linecolor=yellow](!X1 Y1){0.5}
\pscircle*[linecolor=red](!X2 Y2){0.2}
\pscircle*[linecolor=white](!X3 Y3){0.07}
}
\end{animateinline}
\end{center}
\end{document}

