gravitation/Kepler16.tex

   1
   2 \documentclass{article}
   3 \usepackage[a4paper,margin=2cm]{geometry}
   4 \usepackage[T1]{fontenc}
   5 \usepackage[latin1]{inputenc}%
   6 \usepackage[garamond]{mathdesign}
   7 \usepackage{pst-eqdf,pst-node,pst-tools}
   8 \usepackage{array,amsmath}
   9 \usepackage{animate}
  10 \usepackage{wasysym}
  11 \newpsstyle{vecteurA}{arrowinset=0.05,arrowsize=0.1,linecolor={[rgb]{1 0.5 0}}}
  12 \newpsstyle{vecteurB}{arrowinset=0.05,arrowsize=0.1,linecolor={[rgb]{0 0.5 1}}}
  13 \newpsstyle{vecteurC}{arrowinset=0.1,arrowsize=0.2,linecolor={[rgb]{1 0 0}}}
  14 \makeatletter
  15 %% adapté de \psRandom du package pstricks-add
  16 %% pour rendre aléatoire la taille des étoiles
  17 %% Manuel Luque
  18 \newdimen\pssizeStar
  19 \def\psset@sizeStar#1{\pssetlength\pssizeStar{#1}}
  20 \psset@sizeStar{1pt}
  21 \define@key[psset]{pst-eqd}{randomPoints}[1000]{\def\psk@randomPoints{#1}}
  22 \psset[pst-eqd]{randomPoints=1000}
  23 \define@boolkey[psset]{pst-eqd}[Pst@]{color}[true]{}
  24 \psset[pst-eqd]{color=false}
  25 \def\psRandomStar{\pst@object{psRandomStar}}%
  26 \def\psRandomStar@i{\@ifnextchar({\psRandomStar@ii}{\psRandomStar@iii(0,0)(1,1)}}
  27 \def\psRandomStar@ii(#1){\@ifnextchar({\psRandomStar@iii(#1)}{\psRandomStar@iii(0,0)(#1)}}
  28 \def\psRandomStar@iii(#1)(#2)#3{%
  29   \def\pst@tempA{#3}%
  30   \ifx\pst@tempA\pst@empty\psclip{\psframe(#2)}\else\psclip{#3}\fi
  31   \pst@getcoor{#1}\pst@tempA
  32   \pst@getcoor{#2}\pst@tempB
  33   \begin@SpecialObj
  34   \addto@pscode{
  35     \pst@tempA\space /yMin exch def
  36     /xMin exch def
  37     \pst@tempB\space /yMax exch def
  38     /xMax exch def
  39     /dy yMax yMin sub def
  40     /dx xMax xMin sub def
  41     rrand srand                 % initializes the random generator
  42     /getRandReal { rand 2147483647 div } def
  43      \psk@randomPoints {
  44     /DS \pst@number\pssizeStar\space getRandReal mul def
  45     \@nameuse{psds@\psk@dotstyle}
  46      \ifPst@color getRandReal 1 1 sethsbcolor \fi
  47      getRandReal dx mul xMin add
  48      getRandReal dy mul yMin add
  49      Dot
  50      \ifx\psk@fillstyle\psfs@solid fill \fi stroke
  51     } repeat
  52   }%
  53   \end@SpecialObj
  54   \endpsclip
  55   \ignorespaces
  56 }
  57 \makeatother
  58 %%%%%%%%%%%%%%%%%%
  59 %timeline
  60 \begin{filecontents}{kepler16.dat}
  61 ::0x0
  62 ::1
  63 ::2
  64 ::3
  65 ::4
  66 ::5
  67 ::6
  68 ::7
  69 ::8
  70 ::9
  71 ::10
  72 ::11
  73 ::12
  74 ::13
  75 ::14
  76 ::15
  77 ::16
  78 ::17
  79 ::18
  80 ::19
  81 ::20
  82 ::21
  83 ::22
  84 ::23
  85 ::24
  86 ::25
  87 ::26
  88 ::27
  89 ::28
  90 ::29
  91 ::30
  92 ::31
  93 ::32
  94 ::33
  95 ::34
  96 ::35
  97 ::36
  98 ::37
  99 ::38
 100 ::39
 101 ::40
 102 ::41
 103 ::42
 104 ::43
 105 ::44
 106 ::45
 107 ::46
 108 ::47
 109 ::48
 110 ::49
 111 ::50
 112 ::51
 113 ::52
 114 ::53
 115 ::54
 116 ::55
 117 ::56
 118 ::57
 119 ::58
 120 ::59
 121 ::60
 122 ::61
 123 ::62
 124 ::63
 125 ::64
 126 ::65
 127 ::66
 128 ::67
 129 ::68
 130 ::69
 131 ::70
 132 ::71
 133 ::72
 134 ::73
 135 ::74
 136 ::75
 137 ::76
 138 ::77
 139 ::78
 140 ::79
 141 ::80
 142 ::81
 143 ::82
 144 ::83
 145 ::84
 146 ::85
 147 ::86
 148 ::87
 149 ::88
 150 ::89
 151 ::90
 152 ::91
 153 ::92
 154 ::93
 155 ::94
 156 ::95
 157 ::96
 158 ::97
 159 ::98
 160 ::99
 161 ::100
 162 ::101
 163 ::102
 164 ::103
 165 ::104
 166 ::105
 167 ::106
 168 ::107
 169 ::108
 170 ::109
 171 ::110
 172 ::111
 173 ::112
 174 ::113
 175 ::114
 176 ::115
 177 ::116
 178 ::117
 179 ::118
 180 ::119
 181 ::120
 182 ::121
 183 ::122
 184 ::123
 185 ::124
 186 ::125
 187 ::126
 188 ::127
 189 ::128
 190 ::129
 191 ::130
 192 ::131
 193 ::132
 194 ::133
 195 ::134
 196 ::135
 197 ::136
 198 ::137
 199 ::138
 200 ::139
 201 ::140
 202 ::141
 203 ::142
 204 ::143
 205 ::144
 206 ::145
 207 ::146
 208 ::147
 209 ::148
 210 ::149
 211 ::150
 212 ::151
 213 ::152
 214 ::153
 215 ::154
 216 ::155
 217 ::156
 218 ::157
 219 ::158
 220 ::159
 221 ::160
 222 ::161
 223 ::162
 224 ::163
 225 ::164
 226 ::165
 227 ::166
 228 ::167
 229 ::168
 230 ::169
 231 ::170
 232 ::171
 233 ::172
 234 ::173
 235 ::174
 236 ::175
 237 ::176
 238 ::177
 239 ::178
 240 ::179
 241 ::180
 242 ::181
 243 ::182
 244 ::183
 245 ::184
 246 ::185
 247 ::186
 248 ::187
 249 ::188
 250 ::189
 251 ::190
 252 ::191
 253 ::192
 254 ::193
 255 ::194
 256 ::195
 257 ::196
 258 ::197
 259 ::198
 260 ::199
 261 ::200
 262 \end{filecontents}
 263 %%%%%%%%%%%%%%%%%
 264 \title{Gravitation : une planète à deux soleils comme Tatooine celle de Luke Skywalker dans \textit{La Guerre des étoiles}}
 265 \date{10 juillet 2\,012}
 266 \begin{document}
 267 \maketitle
 268 \section{La représentation avec les données de la NASA}
 269 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}.
 270 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.
 271 \begin{center}
 272 \begin{tabular}{|ll|}
 273 \hline
 274 \textbf{ Paramètres}&\textbf{Valeurs}\\ \hline
 275   \textit{Étoile A}&\\ \hline
 276   Masse, $M_A(M_{\astrosun})$&0.6897\\
 277   Rayon, $R_A(R_{\astrosun})$&0.6489\\ \hline
 278   \textit{Étoile B}&\\ \hline
 279   Masse, $M_B(M_{\astrosun})$&0.20255\\
 280   Rayon, $R_B(R_{\astrosun})$&0.22623\\ \hline
 281   \textit{Planète b}&\\ \hline
 282   Masse, $M_b(M_{\jupiter})$&0.333\\
 283   Rayon, $R_b(R_{\jupiter})$&0.7538\\ \hline
 284   \textit{Orbite de l'étoile binaire}&\\ \hline
 285   Période(jours)&41.076\\
 286   Demi-grand axe(ua) $\mathrm{a_1}$&0.22431\\
 287   Excentricité $\mathrm{e_1}$&0.15944\\
 288   Argument du périastre(deg) $\omega_1$&263.464\\ \hline
 289   \textit{Orbite de la planète circumbinaire} &\\ \hline
 290    Période(jours)&228.776\\
 291   Demi-grand axe(ua)$\mathrm{a_2}$&0.7048\\
 292   Excentricité $\mathrm{e_2}$&0.0069\\
 293   Argument du périastre(deg) $\omega_2$&318\\ \hline
 294 \end{tabular}
 295 \end{center}
 296
 297 \begin{center}
 298 \begin{pspicture}(-5,-5)(5,5)
 299 \psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]%
 300 \pstVerb{% pour le système binaire d'étoiles
 301         /a1 0.22431 def % semi-major axis
 302         /e1 0.15944 def % eccentricity
 303         /w1 263.464 def % argument of periaspe
 304         /MA 0.6897 def  % masse of star A
 305         /MB 0.20255 def % masse of star B
 306         /Mt MA MB add def % masse totale
 307         /P1 41.076 def  % period (day)
 308         % pour la planète
 309         /Mb 0.3333 def  % masse of planet b (M_Jupiter)
 310         /a2 0.7048 def  % semi-major axis
 311         /e2 0.0069 def  % eccentricity
 312         /w2 318 def  % argument of periaspe
 313         /P2 228.776 def % period (day)
 314 % on en déduit le paramètre de chacune des ellipses
 315         /p1 a1 1 e1 dup mul add mul def
 316         /p2 a2 1 e2 dup mul add mul def
 317         }%
 318 \pnode(!/radius p1 1 e1 w1 w1 sub cos mul add div def
 319         /radius1 radius MB Mt div mul neg def
 320          radius1 w1 cos mul 5 mul
 321          radius1 w1 sin mul 5 mul){MA0}
 322 \parametricplot[plotpoints=360,unit=5,linecolor=red]{0}{360}{
 323     /radius p1 1 e1 t w1 sub cos mul add div def
 324     /radius1 radius MB Mt div mul neg def
 325      radius1 t cos mul
 326      radius1 t sin mul}
 327 \pscircle[fillstyle=solid,fillcolor=yellow!50](MA0){0.3}
 328 \parametricplot[plotpoints=360,unit=5,linecolor=green]{0}{360}{
 329     /radius p1 1 e1 t w1 sub cos mul add div def
 330     /radius2 radius MA Mt div mul def
 331      radius2 t cos mul
 332      radius2 t sin mul}
 333 \pnode(!/radius p1 1 e1 w1 w1 sub cos mul add div def
 334         /radius1 radius MA Mt div mul def
 335         radius1 w1 cos mul 5 mul
 336         radius1 w1 sin mul 5 mul){MB0}
 337 \pscircle[fillstyle=solid,fillcolor=red](MB0){0.15}
 338 \parametricplot[plotpoints=360,unit=5,linecolor=blue]{0}{360}{
 339     /radius p2 1 e2 t w2 sub cos mul add div def
 340      radius t cos mul
 341      radius t sin mul}
 342 \pnode(!/radius p2 1 e2 w2 w2 sub cos mul add div def
 343         radius w2 cos mul 5 mul
 344         radius w2 sin mul 5 mul){Mb0}
 345 \pscircle[fillstyle=solid,fillcolor=blue](Mb0){0.07}
 346 \end{pspicture}
 347 \end{center}
 348 \section{La simulation avec PSTricks}
 349 \begin{center}
 350 \begin{pspicture}(-3,-1)(6,7)
 351 \psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]%
 352 \pstVerb{/x01 -3 def /y01 0 def
 353          /x02  5 def /y02 0 def
 354          /xr012 x02 x01 sub def
 355          /yr012 y02 y01 sub def
 356          /x0 2 def /y0 5 def
 357          /xr002 x02 x0 sub def
 358          /yr002 y02 y0 sub def
 359          /xr001 x01 x0 sub def
 360          /yr001 y01 y0 sub def
 361          /M0 0.2 def
 362          /M1 0.07 def
 363          /M2 3.2e-4 def
 364          /Mt M1 M2 add def
 365          /xG012 x01 M1 mul x02 M2 mul add Mt div def
 366          /yG012 y01 M1 mul y02 M2 mul add Mt div def
 367          }%
 368 \pnode(0,0){O}
 369 \pnode(!x01 y01){M1}
 370 \pnode(!x02 y02){M2}
 371 \pnode(!x0 y0){M0}
 372 \pnode(!xG012 yG012){C12} % centre de masse de M1 et M2
 373 \pscircle[fillstyle=solid,fillcolor=yellow!50](M2){0.21}
 374 \pscircle[fillstyle=solid,fillcolor=red!50](M1){0.5}
 375 \pscircle[fillstyle=solid,fillcolor=blue!50](M0){0.07}
 376 \psline[linestyle=dotted](M1)(M2)
 377 \psline[linestyle=dashed](M2)(M0)(M1)
 378 \psline{<->}(6,0)(0,0)(0,7)
 379 \uput[u](M0){$M_0$}
 380 \uput[r](0,6.8){$y$}
 381 \uput[u](5.8,0){$x$}
 382 \uput[l](0,6.8){$\mathcal{R}$}
 383 \uput{0.6}[d](M1){$M_1$}
 384 \uput{0.25}[d](M2){$M_2$}
 385 \rput(M1){\psline[style=vecteurC]{->}(!xr012 3 div yr012 3 div)\uput[d](!xr012 3 div yr012 3 div){$\overrightarrow{F}_{2/1}$}
 386           \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}$}}
 387 \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}$}
 388              \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}$}}
 389 \rput(M0){\psline[style=vecteurC]{->}(!xr002 3 div yr002 3 div)\uput[r](!xr002 3 div yr002 3 div){$\overrightarrow{F}_{2/0}$}}
 390 \rput(M0){\psline[style=vecteurC]{->}(!xr001 5 div yr001 5 div)\uput[r](!xr001 5 div yr001 5 div){$\overrightarrow{F}_{1/0}$}}
 391 \end{pspicture}
 392 \end{center}
 393 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.
 394
 395 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 :
 396 \begin{itemize}
 397   \item $\overrightarrow{r}_{01}=\overrightarrow{r}_1-\overrightarrow{r}_0$
 398   \item $\overrightarrow{r}_{02}=\overrightarrow{r}_2-\overrightarrow{r}_0$
 399   \item $\overrightarrow{r}_{12}=\overrightarrow{r}_2-\overrightarrow{r}_1$
 400 \end{itemize}
 401 Chacune des forces s'exprime par :
 402 \[
 403 \left\{
 404 \begin{array}[m]{l}
 405 \overrightarrow{F}_{1/0}=\mathcal{G}\dfrac{m_1m_0}{r_{10}^3}\overrightarrow{r_{01}}=-\overrightarrow{F}_{0/1}\\[1em]
 406 \overrightarrow{F}_{2/0}=\mathcal{G}\dfrac{m_2m_0}{r_{20}^3}\overrightarrow{r_{02}}=-\overrightarrow{F}_{0/2}\\[1em]
 407 \overrightarrow{F}_{2/1}=\mathcal{G}\dfrac{m_1m_2}{r_{12}^3}\overrightarrow{r_{12}}=-\overrightarrow{F}_{1/2}
 408 \end{array}
 409 \right.
 410 \]
 411 L'application de la loi de Newton à chacun des corps donne :
 412 \[
 413 \left\{
 414 \begin{array}[m]{l}
 415 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}}+
 416                                           \mathcal{G}\dfrac{m_2m_0}{r_{20}^3}\overrightarrow{r_{02}}\\[1em]
 417 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}}+
 418                                           \mathcal{G}\dfrac{m_1m_2}{r_{12}^3}\overrightarrow{r_{12}}\\[1em]
 419 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}}-
 420                                           \mathcal{G}\dfrac{m_1m_2}{r_{12}^3}\overrightarrow{r_{12}}
 421 \end{array}
 422 \right.
 423 \]
 424 Ce qui conduit à un système de 6 équations différentielles :
 425 \[
 426 \left\{
 427 \begin{array}[m]{l}
 428 \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]
 429 \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]
 430 \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]
 431 \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]
 432 \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]
 433 \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)
 434 \end{array}
 435 \right.
 436 \]
 437
 438 \begin{verbatim}
 439 %   0    1    2    3    4    5    6    7    8   9     10    11
 440 % y[0] y[1] y[2] y[3] y[4] y[5] y[6] y[7] y[8] y[9] y[10] y[11]
 441 %  x0   y0  x'0  y'0   x1   y1  x'1  y'1   x2   y2   x'2   y'2
 442 \def\GravAlgIIIcorps{%
 443   y[2]|y[3]|%
 444   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|%
 445   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|%
 446   y[6]|y[7]|%
 447   -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|%
 448   -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|%
 449   y[10]|y[11]|%
 450   -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|%
 451   -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}
 452 \end{verbatim}
 453
 454
 455 \newpage
 456 % y[0] y[1] y[2] y[3] y[4] y[5] y[6] y[7] y[8] y[9] y[10] y[11]
 457 %  x0   y0  x'0  y'0   x1   y1  x'1  y'1   x2   y2   x'2   y'2
 458 \def\GravAlgIIIcorps{%
 459   y[2]|y[3]|%
 460   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|%
 461   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|%
 462   y[6]|y[7]|%
 463   -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|%
 464   -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|%
 465   y[10]|y[11]|%
 466   -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|%
 467   -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}
 468 %
 469
 470 \begin{center}
 471 \def\InitCond{x00 y00 v0x0 v0y0 x01 y01 v0x1 v0y1 x02 y02 v0x2 v0y2}
 472 \psset{method=rk4}
 473 \begin{pspicture}(-8,-8)(8,8)
 474 \psframe*[linecolor={[cmyk]{1 1 0 0.7}}](-8,-8)(8,8)
 475 \pstVerb{
 476 /arccos {
 477    dup
 478    dup mul neg 1 add sqrt
 479    exch
 480    atan
 481 } def
 482         /G 1 def
 483         /w2 55 def
 484         /x00 7 w2 cos mul def /y00 7 w2 sin mul def
 485         /x01 0.4962 def
 486         /y01 0 def
 487         /v0x1 0 def
 488         /v0y1 .135 def
 489         /x02 -1.69 def
 490         /y02 0 def
 491         /v0x2 0 def
 492         /v0y2 -.46 def
 493         /v0x0 0.36 w2 sin mul neg def
 494         /v0y0 0.36 w2 cos mul def
 495          /M1 0.7 def
 496          /M2 0.2 def
 497          /M0 3e-4 def
 498          /Mt M1 M2 add def
 499          /xG012 x01 M1 mul x02 M2 mul add Mt div def
 500          /yG012 y01 M1 mul y02 M2 mul add Mt div def
 501          }%
 502 \psgrid[subgriddiv=0,gridcolor=white,griddots=10,gridlabels=0pt](-8,-8)(8,8)%
 503 \psequadiff[plotpoints=1000,algebraic,
 504             plotfuncx=y 4 get
 505                       y 4 get M1 mul
 506                       y 8 get M2 mul add
 507                       M1 M2 add div
 508                       sub,
 509             plotfuncy=5 get
 510                       dup M1 mul
 511                       y 9 get M2 mul add
 512                       M1 M2 add div
 513                       sub,
 514             tabname=XAYA,
 515             saveData,filename=XAYA.dat
 516 %            ]{0}{18}{\InitCond}{\GravAlgIIIcorps}
 517             ]{0}{120}{\InitCond}{\GravAlgIIIcorps}
 518 \listplot[unit=1,linecolor=yellow]{XAYA aload pop}
 519 \psequadiff[plotpoints=1000,algebraic,
 520             plotfuncx=y 8 get
 521                       y 4 get M1 mul
 522                       y 8 get M2 mul add
 523                       M1 M2 add div
 524                       sub,
 525             plotfuncy=9 get
 526                       dup M2 mul
 527                       y 5 get M1 mul add
 528                       M1 M2 add div
 529                       sub,
 530             tabname=XBYB,
 531             saveData,filename=XBYB.dat
 532 %            ]{0}{18}{\InitCond}{\GravAlgIIIcorps}
 533             ]{0}{120}{\InitCond}{\GravAlgIIIcorps}
 534 \listplot[unit=1,linecolor=red]{XBYB aload pop}
 535 \psequadiff[plotpoints=1000,algebraic,
 536             plotfuncx=y 0 get
 537                       y 4 get M1 mul
 538                       y 8 get M2 mul add
 539                       M1 M2 add div
 540                       sub,
 541             plotfuncy=1 get
 542                       y 5 get M1 mul
 543                       y 9 get M2 mul add
 544                       M1 M2 add div
 545                       sub,
 546             tabname=XbYb,
 547             saveData,filename=XPYP.dat
 548             ]{0}{120}{\InitCond}{\GravAlgIIIcorps}
 549 \listplot[unit=1,linecolor=blue]{XbYb aload pop}
 550 \pnode(!x01 xG012 sub y01 yG012 sub){M01}
 551 \pnode(!x02 xG012 sub y02 yG012 sub){M02}
 552 \pnode(!x00  xG012 sub y00  yG012 sub){M00}
 553 \pscircle*[linecolor=yellow](M01){0.4}
 554 \pscircle*[linecolor=red](M02){0.16}
 555 \pscircle*[linecolor=white](M00){0.07}
 556 %\rput(M01){\psline[unit=2,linecolor=red]{->}(!v0x1 v0y1)}
 557 %\rput(M02){\psline[unit=2,linecolor=blue]{->}(!v0x2 v0y2)}
 558 \psdot(!xG012 yG012)
 559 \end{pspicture}
 560 \end{center}
 561 \section{Animation avec pst-eqdf et animate}
 562 \begin{center}
 563 \def\nFrames{200}% 200 images
 564 \begin{animateinline}[controls,timeline=kepler16.dat,loop,%
 565                      begin={\begin{pspicture}(-8,-8)(8,8)},
 566                      end={\end{pspicture}}]{5}% 5 images/s
 567 \pstVerb{/XY1 [(XAYA.dat) run] def
 568          /XY2 [(XBYB.dat) run] def
 569          /XY3 [(XPYP.dat) run] def
 570          }%
 571 \psframe*[linecolor={[cmyk]{1 1 0 0.7}}](-8,-8)(8,8)
 572 \psRandomStar[linecolor={[rgb]{1,1,0.5}},
 573     randomPoints=1000,sizeStar=1pt](-8,-8)(8,8){\psframe[linestyle=none](-8,-8)(8,8)}
 574 %\listplot[linecolor=gray]{XY1 aload pop}
 575 %\listplot[linecolor=gray]{XY2 aload pop}
 576 %\listplot[linecolor=gray]{XY3 aload pop}
 577 \newframe
 578 \multiframe{\nFrames}{i=0+10}{% 1 point sur 10
 579 \pstVerb{/X1 XY1 \i\space get def
 580          /Y1 XY1 \i\space 1 add get def
 581          /X2 XY2 \i\space get def
 582          /Y2 XY2 \i\space 1 add get def
 583          /X3 XY3 \i\space get def
 584          /Y3 XY3 \i\space 1 add get def
 585         }%
 586 %\psdot(!X1 Y1)
 587 \pscircle*[linecolor=yellow](!X1 Y1){0.5}
 588 \pscircle*[linecolor=red](!X2 Y2){0.2}
 589 \pscircle*[linecolor=white](!X3 Y3){0.07}
 590 }
 591 \end{animateinline}
 592 \end{center}
 593 \end{document}
 594