gravitation/pb_3corps_1.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 \newpsstyle{vecteurA}{arrowinset=0.05,arrowsize=0.1,linecolor={[rgb]{1 0.5 0}}}
  11 \newpsstyle{vecteurB}{arrowinset=0.05,arrowsize=0.1,linecolor={[rgb]{0 0.5 1}}}
  12 \newpsstyle{vecteurC}{arrowinset=0.1,arrowsize=0.2,linecolor={[rgb]{1 0 0}}}
  13 %timeline
  14 \begin{filecontents}{pb2corps.dat}
  15 ::0x0
  16 ::1
  17 ::2
  18 ::3
  19 ::4
  20 ::5
  21 ::6
  22 ::7
  23 ::8
  24 ::9
  25 ::10
  26 ::11
  27 ::12
  28 ::13
  29 ::14
  30 ::15
  31 ::16
  32 ::17
  33 ::18
  34 ::19
  35 ::20
  36 ::21
  37 ::22
  38 ::23
  39 ::24
  40 ::25
  41 ::26
  42 ::27
  43 ::28
  44 ::29
  45 ::30
  46 ::31
  47 ::32
  48 ::33
  49 ::34
  50 ::35
  51 ::36
  52 ::37
  53 ::38
  54 ::39
  55 ::40
  56 ::41
  57 ::42
  58 ::43
  59 ::44
  60 ::45
  61 ::46
  62 ::47
  63 ::48
  64 ::49
  65 ::50
  66 ::51
  67 ::52
  68 ::53
  69 ::54
  70 ::55
  71 ::56
  72 ::57
  73 ::58
  74 ::59
  75 ::60
  76 ::61
  77 ::62
  78 ::63
  79 ::64
  80 ::65
  81 ::66
  82 ::67
  83 ::68
  84 ::69
  85 ::70
  86 ::71
  87 ::72
  88 ::73
  89 ::74
  90 ::75
  91 ::76
  92 ::77
  93 ::78
  94 ::79
  95 ::80
  96 ::81
  97 ::82
  98 ::83
  99 ::84
 100 ::85
 101 ::86
 102 ::87
 103 ::88
 104 ::89
 105 ::90
 106 ::91
 107 ::92
 108 ::93
 109 ::94
 110 ::95
 111 ::96
 112 ::97
 113 ::98
 114 ::99
 115 ::100
 116 ::101
 117 ::102
 118 ::103
 119 ::104
 120 ::105
 121 ::106
 122 ::107
 123 ::108
 124 ::109
 125 ::110
 126 ::111
 127 ::112
 128 ::113
 129 ::114
 130 ::115
 131 ::116
 132 ::117
 133 ::118
 134 ::119
 135 ::120
 136 ::121
 137 ::122
 138 ::123
 139 ::124
 140 ::125
 141 ::126
 142 ::127
 143 ::128
 144 ::129
 145 ::130
 146 ::131
 147 ::132
 148 ::133
 149 ::134
 150 ::135
 151 ::136
 152 ::137
 153 ::138
 154 ::139
 155 ::140
 156 ::141
 157 ::142
 158 ::143
 159 ::144
 160 ::145
 161 ::146
 162 ::147
 163 ::148
 164 ::149
 165 ::150
 166 ::151
 167 ::152
 168 ::153
 169 ::154
 170 ::155
 171 ::156
 172 ::157
 173 ::158
 174 ::159
 175 ::160
 176 ::161
 177 ::162
 178 ::163
 179 ::164
 180 ::165
 181 ::166
 182 ::167
 183 ::168
 184 ::169
 185 ::170
 186 ::171
 187 ::172
 188 ::173
 189 ::174
 190 ::175
 191 ::176
 192 ::177
 193 ::178
 194 ::179
 195 ::180
 196 ::181
 197 ::182
 198 ::183
 199 ::184
 200 ::185
 201 ::186
 202 ::187
 203 ::188
 204 ::189
 205 ::190
 206 ::191
 207 ::192
 208 ::193
 209 ::194
 210 ::195
 211 ::196
 212 ::197
 213 ::198
 214 ::199
 215 ::200
 216 \end{filecontents}
 217 %%%%%%%%%%%%%%%%%%
 218 \title{Gravitation : le problème des trois corps avec PSTricks\\ partie 1}
 219 \date{7 juillet 2\,012}
 220 \begin{document}
 221 \maketitle
 222 \section{Système de trois corps, dont l'un est beaucoup plus massif que les deux autres}
 223 \subsection{Les équations}
 224 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.
 225 \begin{center}
 226 \begin{pspicture}(-3,-1)(6,7)
 227 \psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]%
 228 \pstVerb{/x01 -3 def /y01 4 def
 229          /x02  6 def /y02 6 def
 230          /xr012 x02 x01 sub def
 231          /yr012 y02 y01 sub def
 232          /M0 10 def
 233          /M1 3 def
 234          /M2 1 def
 235          /Mt M1 M2 add def
 236          /xG012 x01 M1 mul x02 M2 mul add Mt div def
 237          /yG012 y01 M1 mul y02 M2 mul add Mt div def
 238          }%
 239 \pnode(0,0){O}
 240 \pnode(!x01 y01){M1}
 241 \pnode(!x02 y02){M2}
 242 \pnode(!xG012 yG012){C12} % centre de masse de M1 et M2
 243 \pscircle[fillstyle=solid,fillcolor=yellow!50](0,0){0.5}
 244 \pscircle[fillstyle=solid,fillcolor=red!50](M1){0.21}
 245 \pscircle[fillstyle=solid,fillcolor=blue!50](M2){0.07}
 246 \psline[linestyle=dotted](M1)(M2)
 247 \psline[linestyle=dashed](M2)(O)(M1)
 248 \psline{<->}(6,0)(0,0)(0,7)
 249 \uput{0.6}[dl](O){$M$}
 250 \uput[r](0,6.8){$y$}
 251 \uput[u](5.8,0){$x$}
 252 \uput[l](0,6.8){$\mathcal{R}$}
 253 \uput{0.22}[l](M1){$M_1$}
 254 \uput{0.1}[u](M2){$M_2$}
 255 \rput(M1){\psline[style=vecteurC]{->}(!xr012 10 div yr012 10 div)\uput[u](!xr012 10 div yr012 10 div){$\overrightarrow{F}_{2/1}$}
 256           \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}$}}
 257 \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}$}
 258           \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}$}}
 259 \end{pspicture}
 260 \end{center}
 261 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 :
 262 \[
 263 \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}}
 264 \]
 265 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}$ :
 266 \[
 267 \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}
 268 \]
 269 La loi de Newton appliquée à chacune des particules $M_1$ et $M_2$ s'écrit :
 270 \[
 271 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}}-
 272  \mathcal{G}\frac{mm_1}{r^3_1}\overrightarrow{r_1}
 273 \]
 274 \[
 275 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}}-
 276  \mathcal{G}\frac{mm_2}{r^3_2}\overrightarrow{r_2}
 277 \]
 278 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 :
 279 \[
 280 \left\{
 281 \begin{array}[m]{l}
 282  \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]
 283  \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]
 284  \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]
 285  \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
 286 \end{array}
 287 \right.
 288 \]
 289 \subsection{La résolution numérique avec pst-eqdf}
 290 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 :
 291 \begin{verbatim}
 292 %   0    1    2    3    4    5    6    7
 293 % y[0] y[1] y[2] y[3] y[4] y[5] y[6] y[7]
 294 % x1    y1  x'1  y'1  x2   y2   x'2  y'2
 295 \def\GravAlgIIcorps{%
 296   y[2]|y[3]|%
 297   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|%
 298   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|%
 299   y[6]|y[7]|%
 300   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|%
 301   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}
 302 \end{verbatim}
 303 Les conditions initiales, positions et vitesses respectives de $M_1$ et $M_2$, déterminent l'évolution du système.
 304 \begin{verbatim}
 305 \def\InitCond{x01 y01 v0x1 v0y1 x02 y02 v0x2 v0y2}
 306 \end{verbatim}
 307 Elles sont définies dans l'environnement \verb+\begin{pspicture}+ :
 308 \begin{verbatim}
 309 \pstVerb{
 310         /x01 -1 def /y01 0 def
 311         /x02  5 def /y02 0 def
 312         /v0x1 0 def
 313         /v0y1 25 def
 314         /v0x2 0 def
 315         /v0y2 -10 def
 316         /M 500 def
 317         /M1 10 def
 318         /M2 1 def
 319         }%
 320 \end{verbatim}
 321 %   0    1    2    3    4    5    6    7
 322 % y[0] y[1] y[2] y[3] y[4] y[5] y[6] y[7]
 323 % x1    y1  x'1  y'1  x2   y2   x'2  y'2
 324 \def\GravAlgIIIcorps{%
 325   y[2]|y[3]|%
 326   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|%
 327   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|%
 328   y[6]|y[7]|%
 329   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|%
 330   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}
 331
 332 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.
 333 \begin{verbatim}
 334 \rput(-2,0){\psPrintValue[decimals=4]{periode1}\hphantom{000000}s}
 335 \rput(-2,1){\psPrintValue[decimals=4]{periode2}\hphantom{000000}s}
 336 \end{verbatim}
 337 La durée choisie qui doit être identique pour les deux, est la plus grande des périodes.
 338 \begin{verbatim}
 339 \psequadiff[whichabs=0,whichord=1,
 340             plotpoints=1000,algebraic,
 341             tabname=X1Y1
 342             ,saveData,filename=IIIXYM1.dat
 343             ]{0}{periode2}{\InitCond}{\GravAlgIIIcorps}
 344 \listplot[unit=1,linecolor=red]{X1Y1 aload pop}
 345 \psequadiff[whichabs=4,whichord=5,
 346             plotpoints=1000,algebraic,
 347             tabname=X2Y2
 348             ,saveData,filename=IIIXYM2.dat
 349             ]{0}{periode2}{\InitCond}{\GravAlgIIIcorps}
 350 \listplot[unit=1,linecolor=blue]{X2Y2 aload pop}
 351 \end{verbatim}
 352
 353 %\newpage
 354 \begin{center}
 355 \def\InitCond{ x01  y01  v0x1  v0y1 x02 y02 v0x2  v0y2}
 356 \psset{method=rk4}
 357 \begin{pspicture}(-6,-6)(6,6)
 358 \pstVerb{
 359 /arccos {
 360    dup
 361    dup mul neg 1 add sqrt
 362    exch
 363    atan
 364 } def
 365         /G 1 def
 366         /x01 -1 def /y01 0 def
 367         /x02  5 def /y02 0 def
 368         /v0x1 0 def
 369         /v0y1 25 def
 370         /v0x2 0 def
 371         /v0y2 -10 def
 372         /M0 500 def
 373         /M1 10 def
 374         /M2 1 def
 375         /xG012 M1 x01 mul M2 x02 mul add M1 M2 add div def
 376         /yG012 M1 y01 mul M2 y02 mul add M1 M2 add div def
 377 % calcul des périodes en supposant les interactions
 378 % entre les 2 satellites négligeables
 379 % satellite 1
 380         /Cste1 x01 v0y1 mul def % Cste des aires pour 1
 381         /Energie1 0.5 M1 mul v0y1 2 exp mul % 1/2 mv^2
 382                  G M1 M0 mul mul x01 abs div sub def
 383          /par1 Cste1 dup mul M0 div def %
 384          /exc1 2 Cste1 dup mul mul Energie1 mul M1 div M0 dup mul div 1 add sqrt def
 385          /a_1 par1 1 exc1 dup mul sub div def % demi-grand axe
 386          /periode1 6.283185 a_1 3 exp G div M0 div sqrt mul def
 387 % satellite 2
 388         /Cste2 x02 v0y2 mul def % Cste des aires pour 2
 389         /Energie2 0.5 M2 mul v0y2 2 exp mul % 1/2 mv^2
 390                  G M2 M0 mul mul x02 abs div sub def
 391          /par2 Cste2 dup mul M0 div def %
 392          /exc2 2 Cste2 dup mul mul Energie2 mul M2 div M0 dup mul div 1 add sqrt def
 393          /a_2 par2 1 exc2 dup mul sub div def % demi-grand axe
 394          /periode2 6.283185 a_2 3 exp G div M0 div sqrt mul def
 395          }%
 396 \psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=8pt](-6,-6)(6,6)%
 397 \psequadiff[whichabs=0,whichord=1,
 398             plotpoints=1000,algebraic,
 399             tabname=X1Y1
 400             ,saveData,filename=IIIXYM1.dat
 401             ]{0}{periode2}{\InitCond}{\GravAlgIIIcorps}
 402 \listplot[unit=1,linecolor=red]{X1Y1 aload pop}
 403 \psequadiff[whichabs=4,whichord=5,
 404             plotpoints=1000,algebraic,
 405             tabname=X2Y2
 406             ,saveData,filename=IIIXYM2.dat
 407             ]{0}{periode2}{\InitCond}{\GravAlgIIIcorps}
 408 \listplot[unit=1,linecolor=blue]{X2Y2 aload pop}
 409 \pnode(!x01 y01){M01}
 410 \pnode(!x02 y02){M02}
 411 \pscircle[fillstyle=solid,fillcolor=red!50](0,0){0.25}
 412 \pscircle[fillstyle=solid,fillcolor=gray!50](M01){0.15}
 413 \pscircle[fillstyle=solid,fillcolor=blue!50](M02){0.07}
 414 \rput(M01){\psline[unit=0.2,linecolor=red]{->}(!v0x1 v0y1)}
 415 \rput(M02){\psline[unit=0.2,linecolor=blue]{->}(!v0x2 v0y2)}
 416 \rput(-2,0){\psPrintValue[decimals=4]{periode1}\hphantom{000000}s}
 417 \rput(-2,1){\psPrintValue[decimals=4]{periode2}\hphantom{000000}s}
 418 \end{pspicture}
 419 \end{center}
 420 %\newpage
 421 \subsection{L'animation}
 422 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.
 423
 424 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.
 425
 426 \begin{center}
 427 %
 428 \def\nFrames{200}% 200 images
 429 \begin{animateinline}[controls,timeline=pb2corps.dat,%loop,
 430                      begin={\begin{pspicture}(-6,-6)(6,6)},
 431                      end={\end{pspicture}}]{10}% 10 images/s
 432 \pstVerb{/XY1 [(IIIXYM1.dat) run] def
 433          /XY2 [(IIIXYM2.dat) run] def}
 434 \psframe*[linecolor={[cmyk]{1 1 0 0.7}}](-6,-6)(6,6)
 435 \psgrid[subgriddiv=0,gridcolor=white,griddots=10,gridlabels=0pt]%
 436 \listplot[linecolor=gray]{XY1 aload pop}
 437 \listplot[linecolor=gray]{XY2 aload pop}
 438 \pscircle[fillstyle=solid,fillcolor=red!50](0,0){0.4}
 439 \newframe
 440 \multiframe{\nFrames}{i=0+10}{% 1 point sur 10
 441 \pstVerb{/X1 XY1 \i\space get def
 442          /Y1 XY1 \i\space 1 add get def
 443          /X2 XY2 \i\space get def
 444          /Y2 XY2 \i\space 1 add get def
 445         }%
 446 \pscircle[fillstyle=solid,fillcolor=gray!80](!X1 Y1){0.2}
 447 \pscircle[fillstyle=solid,fillcolor=yellow](!X2 Y2){0.1}
 448 }
 449 \end{animateinline}
 450 \end{center}
 451 \newpage
 452 \section{Mouvements relatifs}
 453 \subsection{Mouvements vus de $M_1$}
 454 Les coordonnées de $M_2$ rapportées à $M_1$  sont :
 455 \[
 456 \left\{
 457 \begin{array}[m]{l}
 458  x_{2/1}=x_2-x_1\\[1em]
 459  y_{2/1}=y_2-y_1\\[1em]
 460 \end{array}
 461 \right.
 462 \]
 463 Celles de $M$ sont $(-x_1,-y_1)$.
 464 \begin{verbatim}
 465 %   0    1    2    3    4    5    6    7
 466 % y[0] y[1] y[2] y[3] y[4] y[5] y[6] y[7]
 467 % x1    y1  x'1  y'1  x2   y2   x'2  y'2
 468 \psequadiff[plotpoints=1000,algebraic,
 469             plotfuncx=y dup 0 get exch
 470                       4 get sub,
 471             plotfuncy=dup 1 get exch
 472                       5 get sub,
 473             tabname=XY21,
 474 %            ,saveData,filename=IIIM21.dat
 475             ]{0}{periode2}{\InitCond}{\GravAlgIIIcorps}
 476 \psequadiff[plotpoints=1000,algebraic,
 477             plotfuncx=y 0 get neg,
 478             plotfuncy=1 get neg,
 479             tabname=XY01
 480 %            ,saveData,filename=IIIXY01.dat
 481             ]{0}{periode2}{\InitCond}{\GravAlgIIIcorps}
 482 \end{verbatim}
 483 \begin{center}
 484 \def\InitCond{ x01  y01  v0x1  v0y1 x02 y02 v0x2  v0y2}
 485 \psset{method=rk4}
 486 \begin{pspicture}(-7,-6)(6,6)
 487 \pstVerb{
 488         /G 1 def
 489         /x01 -1 def /y01 0 def
 490         /x02  5 def /y02 0 def
 491         /v0x1 0 def
 492         /v0y1 25 def
 493         /v0x2 0 def
 494         /v0y2 -10 def
 495         /M0 500 def
 496         /M1 10 def
 497         /M2 1 def
 498         /xG012 M1 x01 mul M2 x02 mul add M1 M2 add div def
 499         /yG012 M1 y01 mul M2 y02 mul add M1 M2 add div def
 500 % calcul des périodes en supposant les interactions
 501 % entre les 2 satellites négligeables
 502 % satellite 1
 503         /Cste1 x01 v0y1 mul def % Cste des aires pour 1
 504         /Energie1 0.5 M1 mul v0y1 2 exp mul % 1/2 mv^2
 505                  G M1 M0 mul mul x01 abs div sub def
 506          /par1 Cste1 dup mul M0 div def %
 507          /exc1 2 Cste1 dup mul mul Energie1 mul M1 div M0 dup mul div 1 add sqrt def
 508          /a_1 par1 1 exc1 dup mul sub div def % demi-grand axe
 509          /periode1 6.283185 a_1 3 exp G div M0 div sqrt mul def
 510 % satellite 2
 511         /Cste2 x02 v0y2 mul def % Cste des aires pour 2
 512         /Energie2 0.5 M2 mul v0y2 2 exp mul % 1/2 mv^2
 513                  G M2 M0 mul mul x02 abs div sub def
 514          /par2 Cste2 dup mul M0 div def %
 515          /exc2 2 Cste2 dup mul mul Energie2 mul M2 div M0 dup mul div 1 add sqrt def
 516          /a_2 par2 1 exc2 dup mul sub div def % demi-grand axe
 517          /periode2 6.283185 a_2 3 exp G div M0 div sqrt mul def
 518          }%
 519 \psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt](-7,-6)(6,6)%
 520 \psequadiff[plotpoints=1000,algebraic,
 521             plotfuncx=y dup 4 get exch
 522                       0 get sub,
 523             plotfuncy=dup 5 get exch
 524                       1 get sub,
 525             tabname=XY21
 526 %            ,saveData,filename=IIIM21.dat
 527             ]{0}{periode2}{\InitCond}{\GravAlgIIIcorps}
 528 \listplot[unit=1,linecolor=blue]{XY21 aload pop}
 529 \psequadiff[plotpoints=1000,algebraic,
 530             plotfuncx=y 0 get neg,
 531             plotfuncy=1 get neg,
 532             tabname=XYM01
 533 %            ,saveData,filename=IIIXYM01.dat
 534             ]{0}{periode2}{\InitCond}{\GravAlgIIIcorps}
 535 \listplot[unit=1,linecolor=red]{XYM01 aload pop}
 536 \pnode(!x02 x01 sub y02 y01 sub){M021}
 537 \pnode(!x01 neg y01 neg){M001}
 538 \pscircle[fillstyle=solid,fillcolor=gray!80](0,0){0.2}
 539 \pscircle[fillstyle=solid,fillcolor=red!50](M001){0.25}
 540 \pscircle[fillstyle=solid,fillcolor=blue](M021){0.07}
 541 \rput(M021){\psline[unit=0.1,linecolor=blue,arrowinset=0.2,arrowsize=2]{->}(!v0x2 v0x1 sub v0y2 v0y1 sub)}
 542 \rput(M001){\psline[unit=0.1,linecolor=red,arrowinset=0.2,arrowsize=2]{->}(!v0x1 neg v0y1 neg)}
 543 \rput(-3,0){\psPrintValue[decimals=4]{periode1}\hphantom{000000}s}
 544 \rput(-3,1){\psPrintValue[decimals=4]{periode2}\hphantom{000000}s}
 545 \end{pspicture}
 546 \end{center}
 547 \subsection{Mouvements vus de $M_2$}
 548 Les coordonnées de $M_1$ rapportées à $M_2$  sont :
 549 \[
 550 \left\{
 551 \begin{array}[m]{l}
 552  x_{1/2}=x_1-x_2\\[1em]
 553  y_{1/2}=y_1-y_2\\[1em]
 554 \end{array}
 555 \right.
 556 \]
 557 Celles de $M$ sont $(-x_2,-y_2)$.
 558 \begin{verbatim}
 559 %   0    1    2    3    4    5    6    7
 560 % y[0] y[1] y[2] y[3] y[4] y[5] y[6] y[7]
 561 % x1    y1  x'1  y'1  x2   y2   x'2  y'2
 562 \psequadiff[plotpoints=1000,algebraic,
 563             plotfuncx=y dup 4 get exch
 564                       0 get sub,
 565             plotfuncy=dup 5 get exch
 566                       1 get sub,
 567             tabname=XY12,
 568 %            ,saveData,filename=IIIM12.dat
 569             ]{0}{periode2}{\InitCond}{\GravAlgIIIcorps}
 570 \psequadiff[plotpoints=1000,algebraic,
 571             plotfuncx=y 4 get neg,
 572             plotfuncy=5 get neg,
 573             tabname=XYM02
 574 %            ,saveData,filename=IIIXYM02.dat
 575             ]{0}{periode2}{\InitCond}{\GravAlgIIIcorps}
 576 \end{verbatim}
 577 \begin{center}
 578 \def\InitCond{ x01  y01  v0x1  v0y1 x02 y02 v0x2  v0y2}
 579 \psset{method=rk4}
 580 \begin{pspicture}(-7,-6)(6,6)
 581 \pstVerb{
 582         /G 1 def
 583         /x01 -1 def /y01 0 def
 584         /x02  5 def /y02 0 def
 585         /v0x1 0 def
 586         /v0y1 25 def
 587         /v0x2 0 def
 588         /v0y2 -10 def
 589         /M0 500 def
 590         /M1 10 def
 591         /M2 1 def
 592         /xG012 M1 x01 mul M2 x02 mul add M1 M2 add div def
 593         /yG012 M1 y01 mul M2 y02 mul add M1 M2 add div def
 594 % calcul des périodes en supposant les interactions
 595 % entre les 2 satellites négligeables
 596 % satellite 1
 597         /Cste1 x01 v0y1 mul def % Cste des aires pour 1
 598         /Energie1 0.5 M1 mul v0y1 2 exp mul % 1/2 mv^2
 599                  G M1 M0 mul mul x01 abs div sub def
 600          /par1 Cste1 dup mul M0 div def %
 601          /exc1 2 Cste1 dup mul mul Energie1 mul M1 div M0 dup mul div 1 add sqrt def
 602          /a_1 par1 1 exc1 dup mul sub div def % demi-grand axe
 603          /periode1 6.283185 a_1 3 exp G div M0 div sqrt mul def
 604 % satellite 2
 605         /Cste2 x02 v0y2 mul def % Cste des aires pour 2
 606         /Energie2 0.5 M2 mul v0y2 2 exp mul % 1/2 mv^2
 607                  G M2 M0 mul mul x02 abs div sub def
 608          /par2 Cste2 dup mul M0 div def %
 609          /exc2 2 Cste2 dup mul mul Energie2 mul M2 div M0 dup mul div 1 add sqrt def
 610          /a_2 par2 1 exc2 dup mul sub div def % demi-grand axe
 611          /periode2 6.283185 a_2 3 exp G div M0 div sqrt mul def
 612          }%
 613 \psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt](-7,-6)(6,6)%
 614 \psequadiff[plotpoints=1000,algebraic,
 615             plotfuncx=y dup 0 get exch
 616                       4 get sub,
 617             plotfuncy=dup 1 get exch
 618                       5 get sub,
 619             tabname=XYM12
 620             ,saveData,filename=IIIXYM12.dat
 621             ]{0}{periode2}{\InitCond}{\GravAlgIIIcorps}
 622 \listplot[unit=1,linecolor=blue]{XYM12 aload pop}
 623 \psequadiff[plotpoints=1000,algebraic,
 624             plotfuncx=y 4 get neg,
 625             plotfuncy=5 get neg,
 626             tabname=XYM02
 627             ,saveData,filename=IIIXYM02.dat
 628             ]{0}{periode2}{\InitCond}{\GravAlgIIIcorps}
 629 \listplot[unit=1,linecolor=red]{XYM02 aload pop}
 630 \pnode(!x01 x02 sub y01 y02 sub){M012}
 631 \pnode(!x02 neg y02 neg){M002}
 632 \pscircle[fillstyle=solid,fillcolor=gray!80](0,0){0.2}
 633 \pscircle[fillstyle=solid,fillcolor=red!50](M002){0.25}
 634 \pscircle[fillstyle=solid,fillcolor=blue](M012){0.07}
 635 \rput(M012){\psline[unit=0.1,linecolor=blue,arrowinset=0.2,arrowsize=2]{->}(!v0x1 v0x2 sub v0y1 v0y2 sub)}
 636 \rput(M002){\psline[unit=0.1,linecolor=red,arrowinset=0.2,arrowsize=2]{->}(!v0x2 neg v0y2 neg)}
 637 \rput(-3,0){\psPrintValue[decimals=4]{periode1}\hphantom{000000}s}
 638 \rput(-3,1){\psPrintValue[decimals=4]{periode2}\hphantom{000000}s}
 639 \end{pspicture}
 640 \end{center}
 641 \subsection{Animation}
 642 \begin{center}
 643 %
 644 \def\nFrames{200}% 200 images
 645 \begin{animateinline}[controls,timeline=pb2corps.dat,%loop,
 646                      begin={\begin{pspicture}(-7,-6)(7,6)},
 647                      end={\end{pspicture}}]{10}% 10 images/s
 648 \pstVerb{/XY1 [(IIIXYM12.dat) run] def
 649          /XY2 [(IIIXYM02.dat) run] def}
 650 \psframe*[linecolor={[cmyk]{1 1 0 0.7}}](-7,-6)(7,6)
 651 \psgrid[subgriddiv=0,gridcolor=white,griddots=10,gridlabels=0pt]%
 652 \listplot[linecolor=gray]{XY1 aload pop}
 653 \listplot[linecolor=gray]{XY2 aload pop}
 654 \pscircle[fillstyle=solid,fillcolor=gray!80](0,0){0.2}
 655 \newframe
 656 \multiframe{\nFrames}{i=0+10}{% 1 point sur 10
 657 \pstVerb{/X1 XY1 \i\space get def
 658          /Y1 XY1 \i\space 1 add get def
 659          /X2 XY2 \i\space get def
 660          /Y2 XY2 \i\space 1 add get def
 661         }%
 662 \pscircle[fillstyle=solid,fillcolor=yellow](!X1 Y1){0.1}
 663 \pscircle[fillstyle=solid,fillcolor=red!50](!X2 Y2){0.4}
 664 }
 665 \end{animateinline}
 666 \end{center}
 667 \end{document}