--- /dev/null
+\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}