\documentclass{article}
\usepackage[a4paper,margin=2cm]{geometry}
\usepackage[latin1]{inputenc}
\usepackage{pstricks,pst-eqdf}
\usepackage{animate}
\newpsstyle{vecteurA}{arrowinset=0.05,arrowsize=0.125,linecolor={[rgb]{1 0.5 0}}}
\newpsstyle{vecteurB}{arrowinset=0.05,arrowsize=0.1,linecolor={[rgb]{0 0.7 1}}}
\newpsstyle{vecteurC}{arrowinset=0.1,arrowsize=0.2,linecolor={[rgb]{1 0.5 0}}}
%%%%%%%%%%%%%%%%%%
\title{Animation du mouvement d'un satellite}
\date{23 juin 2\,012}
%timeline
\begin{filecontents}{satellite.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
\end{filecontents}
\newcommand{\Jupiter}{%
  \psclip{\pscircle[fillstyle=solid,fillcolor=yellow]{1.5}}
  \psset{fillstyle=solid,fillcolor={[cmyk]{0 0.2 0.4 0}},linestyle=none}
    \psframe(-2,.6)(2,1)
    \psframe(-2,0.45)(2,0.52)
    \psframe(-2,.05)(2,.3)
    \psframe(-2,-0.35)(2,-0.2)
    \psframe(-2,-0.9)(2,-.5)
  \endpsclip
  \psellipse[fillstyle=solid,fillcolor={[cmyk]{0 0.4 0.6 0}},linestyle=none](-0.3,-0.6)(0.35,0.2)
  \pscircle{1.5}}
\begin{document}
\def\eqsatellite{%
y[2]|y[3]|-GM*y[0]/((sqrt(y[0]^2+y[1]^2))^3)|-GM*y[1]/((sqrt(y[0]^2+y[1]^2))^3)}%
\begin{center}
%\begin{pspicture}(-10,-10)(6,6)
\def\nFrames{50}% 50 images
\begin{animateinline}[controls,loop,timeline=satellite.dat,%
                     begin={\begin{pspicture}(-10,-10)(6,6)},
                     end={\end{pspicture}}]{5}% 5 images/s
%\uput[r](0,3){$y$}\uput[u](3,0){$x$}
%\psset{unit=2}
\psframe*[linecolor={[rgb]{0 0 0.5}}](-10,-10)(6,6)
\pstVerb{/tabXiYi [(XiYi.dat) run] def}%
\pstVerb{/GM 1 def
    /theta0 30 def
    /r0 2 def
    /x0 r0 theta0 cos mul def
    /y0 r0 theta0 sin mul def
    /v0 0.85 def
    /v0x v0 theta0 sin mul neg def
    /v0y v0 theta0 cos mul def
    /Lc r0 v0 mul def % moment cinetique
    /par Lc dup mul GM div def % paramètre de l'ellipse
% excentricité
    /exc 1 0.5 v0 4 exp mul r0 dup mul mul GM r0 mul v0 dup mul mul sub GM dup mul div 2 mul add sqrt def
%%%%%%%%%%%%%%
    /a_2 par 1 exc dup mul sub div def % demi-grand axe
    /b_2 par 1 exc dup mul sub sqrt div def % demi-petit axe
    /periode 2 3.1416 dup mul a_2 3 exp mul GM div sqrt mul def}%
\psequadiff[method=rk4,
            plotpoints=1000,
            algebraic,
            whichabs=0,
            whichord=1,
            tabname=XiYi,
%            saveData,filename=XiYi.dat
]{0}{43}{x0 y0 v0x v0y}{\eqsatellite}%
\listplot[linecolor=gray,unit=2]{XiYi aload pop}
\psequadiff[method=rk4,
            plotpoints=1000,
            algebraic,
            whichabs=2,
            whichord=3,
            tabname=vxvy,
%            saveData,filename=vxvy.dat
]{0}{43}{x0 y0 v0x v0y}{\eqsatellite}%
%\listplot[unit=2,linecolor=red]{vxvy aload pop}
% on dessine la vitesse un point sur 50
%\pscircle[fillcolor=gray!70,fillstyle=solid,unit=2](0,0){0.4}
\Jupiter%
%\psgrid[subgriddiv=2,gridcolor=lightgray,gridlabels=8pt,unit=2](-5,-5)(3,3)
%\psline[arrowinset=0.05,arrowsize=0.1,unit=2]{<->}(3,0)(0,0)(0,3)
\newframe
\multiframe{\nFrames}{i=0+40}{
\psset{unit=2}%
\pstVerb{/tabXYpartiel {
          2 2 \i\space {/I exch def
          tabXiYi I get
          tabXiYi I 1 add get
          } for
          } def
         /vX vxvy \i\space get def
         /vY vxvy \i\space 1 add get def
         /xi tabXiYi \i\space get def
         /yi tabXiYi \i\space 1 add get def
         /ri3 xi dup mul yi dup mul add sqrt 3 exp def
         /Fx xi ri3 div neg def
         /Fy yi ri3 div neg def}%
\listplot[linecolor=red,linewidth=0.05]{tabXYpartiel}
\rput(!xi yi){\pscircle[fillstyle=solid,fillcolor=gray!50]{0.075}\psline[style=vecteurA]{->}(! vX 2 mul vY 2 mul)\psline[style=vecteurB]{->}(! Fx 5 mul Fy 5 mul)}}
\newframe
\listplot[linecolor=gray,unit=2]{XiYi aload pop}
\listplot[linecolor=red,unit=2,linewidth=0.05]{XiYi aload pop}
\psset{unit=2}%
\rput(!x0 y0){\pscircle[fillstyle=solid,fillcolor=gray!50]{0.075}\psline[style=vecteurA]{->}(! v0x 2 mul v0y 2 mul)\psline[style=vecteurB]{->}(! x0 r0 3 exp div 5 mul neg y0 r0 3 exp div 5 mul neg)}
\end{animateinline}
\end{center}
\end{document}
