
\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}
\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}}}
%timeline
\begin{filecontents}{pb2corps.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 : le problème des deux corps avec PSTricks\\ partie 2}
\date{5 juillet 2\,012}
\begin{document}
\maketitle
\section{Mouvements dans le repère inertiel}
On considère un système de deux corps en interaction gravitationnelle $M_1$ de masse $m_1$ et $M_2$ de masse $m_2$ dans le repère galiléen \textit{inertiel} $\mathcal{R}$. Ils sont supposés ponctuels.
\begin{center}
\begin{pspicture}(-3,-1)(5,6)
\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]%
\pstVerb{/x01 -2 def /y01 1 def
         /x02  5 def /y02 4 def
         /xr0 x02 x01 sub def
         /yr0 y02 y01 sub def
         /M1 3 def
         /M2 1 def
         /Mt M1 M2 add def
         /xG0 x01 M1 mul x02 M2 mul add Mt div def
         /yG0 y01 M1 mul y02 M2 mul add Mt div def}%
\pnode(0,0){O}
\pnode(!x01 y01){M1}
\pnode(!x02 y02){M2}
\pnode(!xG0 yG0){C}
\pscircle[fillstyle=solid,fillcolor=red!50](M1){0.21}
\pscircle[fillstyle=solid,fillcolor=blue!50](M2){0.07}
\psdot(C)
\psline[linestyle=dotted](M1)(M2)
\psline[linestyle=dashed](M2)(O)(M1)
\psline{<->}(5,0)(0,0)(0,6)
\uput[dl](O){$O$}
\uput[r](0,5.8){$y$}
\uput[u](4.8,0){$x$}
\uput[l](0,5.8){$\mathcal{R}$}
\uput{0.22}[l](M1){$M_1$}
\uput{0.1}[u](M2){$M_2$}
\uput[u](C){$C$}
\pcline[offset=5pt,linestyle=none](M1)(M2)
\ncput*[nrot=:U]{$r$}
\rput(M1){\psline[style=vecteurC]{->}(!xr0 5 div yr0 5 div)\uput[u](!xr0 10 div yr0 10 div){$\overrightarrow{F}_{2/1}$}}
\rput(M2){\psline[style=vecteurC]{->}(!xr0 5 div neg yr0 5 div neg)\uput[u](!xr0 10 div neg yr0 10 div neg){$\overrightarrow{F}_{1/2}$}}
\end{pspicture}
\end{center}
On note $\overrightarrow{r}=\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 :
\[
\overrightarrow{F}_{1/2}=-\mathcal{G}\frac{m_1m_2}{r^3}\overrightarrow{r}\quad \overrightarrow{F}_{2/1}=+\mathcal{G}\frac{m_1m_2}{r^3}\overrightarrow{r}
\]
La loi de Newton appliquée à chacune des particules s'écrit :
\[
m_1\frac{\mathrm{d}^2\overrightarrow{r}}{\mathrm{d}t^2}=\hphantom{-}\mathcal{G}\frac{m_1m_2}{r^3}\overrightarrow{r}
\]
\[
m_2\frac{\mathrm{d}^2\overrightarrow{r}}{\mathrm{d}t^2}=-\mathcal{G}\frac{m_1m_2}{r^3}\overrightarrow{r}
\]
Sachant que $r=\sqrt{(x_{2}-x_{1})^2+(y_{2}-y_{1})^2}$, nous obtenons un système de quatre équations différentielles :
\[
\left\{
\begin{array}[m]{l}
 \ddot{x_1}=\hphantom{-}\mathcal{G}\displaystyle\frac{m_2}{r^3}(x_2-x_1)\\[1em]
 \ddot{y_1}=\hphantom{-}\mathcal{G}\displaystyle\frac{m_2}{r^3}(y_2-y_1)\\[1em]
 \ddot{x_2}=-\mathcal{G}\displaystyle\frac{m_1}{r^3}(x_2-x_1)\\[1em]
 \ddot{y_2}=-\mathcal{G}\displaystyle\frac{m_1}{r^3}(y_2-y_1)
\end{array}
\right.
\]

%  x    y    x'   y'
% y[0] y[1] y[2] y[3]
\def\FictifAlg{%
  y[2]|y[3]|%
  -(M1+M2)*y[0]/(y[0]^2+y[1]^2)^1.5|%
  -(M1+M2)*y[1]/(y[0]^2+y[1]^2)^1.5}
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 :
\begin{verbatim}
%   0    1    2    3    4    5    6    7
% y[0] y[1] y[2] y[3] y[4] y[5] y[6] y[7]
% x1    y1  x'1  y'1  x2   y2   x'2  y'2
\def\GravAlgIIcorps{%
  y[2]|y[3]|%
  M2*(y[4]-y[0])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5|%
  M2*(y[5]-y[1])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5|%
  y[6]|y[7]|%
  M1*(y[0]-y[4])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5|%
  M1*(y[1]-y[5])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5}
\end{verbatim}
Les conditions initiales, positions et vitesses respectives des deux corps, déterminent l'évolution du système.
\begin{verbatim}
\def\InitCond{x01 y01 v0x1 v0y1 x02 y02 v0x2 v0y2}
\end{verbatim}
Elles sont définies dans l'environnement \verb+\begin{pspicture}+ :
\begin{verbatim}
\pstVerb{
        /x01 2 def
        /y01 2 def
        /v0x1 .2 def
        /v0y1 0.25 def
        /x02 -3 def
        /y02 0 def
        /v0x2 -0.25 def
        /v0y2 -0.5 def
        /M1 3 def
        /M2 1 def
        }%
\end{verbatim}
Dans cette même variable \verb+\pstVerb+, il est utile de faire d'autres calculs, comme la période du mouvement afin de bien choisir l'intervalle de temps sur lequel on dessine les trajectoires. Les calculs s'effectuent en \textsf{postscript} et l'affichage du résultat utilise, par exemple, le package \textsf{pst-tools}.
\begin{verbatim}
\rput(-2,0){\psPrintValue[decimals=4]{periode}\hphantom{00000}s}
\end{verbatim}
On désactivera cette ligne par \% au début, une fois la période connue. Dans l'exemple proposé $T=41,0688\,\mathrm{s}$. On étudie le mouvement sur trois périodes environ.
%%  0  1   2   3  4  5   6   7
%% x1 y1 x'1 y'1 x2 y2 x'2 y'2
\def\GravAlgIIcorps{%
  y[2]|y[3]|%
  M2*(y[4]-y[0])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5|%
  M2*(y[5]-y[1])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5|%
  y[6]|y[7]|%
  M1*(y[0]-y[4])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5|%
  M1*(y[1]-y[5])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5}
%\newpage
\begin{center}
% conditions initiales pour le point réduit
\def\InitCondred{ xr0 yr0 v0xr  v0yr}
\begin{pspicture}(-4,-5)(12,12)
\pstVerb{
/arccos {
   dup
   dup mul neg 1 add sqrt
   exch
   atan
} def
        /G 1 def
        /x01 2 def
        /y01 2 def
        /v0x1 .2 def
        /v0y1 0.25 def
        /x02 -3 def
        /y02 0 def
        /v0x2 -0.25 def
        /v0y2 -0.5 def
        /M1 3 def
        /M2 1 def
        /Mt M1 M2 add def
        /mu M1 M2 mul M1 M2 add div def % masse réduite
        /K G Mt mul def
        /xG0 M1 x01 mul M2 x02 mul add M1 M2 add div def
        /yG0 M1 y01 mul M2 y02 mul add M1 M2 add div def
        /vG0x M1 v0x1 mul M2 v0x2 mul add M1 M2 add div def
        /vG0y M1 v0y1 mul M2 v0y2 mul add M1 M2 add div def
        /K1 M2 Mt div neg def
        /K2 M1 Mt div def
% conditions initiales pour le point réduit
        /xr0 x02 x01 sub def
        /yr0 y02 y01 sub def
        /theta_0 yr0 xr0 atan def
        /r0 xr0 dup mul yr0 dup mul add sqrt def
        /v0xr v0x2 v0x1 sub def
        /v0yr v0y2 v0y1 sub def
        /v0r_2 v0xr dup mul v0yr dup mul add def
% constante des aires
        /Cste xr0 v0yr mul yr0 v0xr mul sub def % Cste des aires
        /Energie 0.5 mu mul v0r_2 mul % 1/2 mv^2
                 G M1 M2 mul mul r0 div sub
         def
         /par Cste dup mul K div def % p
         /exc 2 Cste dup mul mul Energie mul mu div K dup mul div 1 add sqrt def % e
         /a_2 par 1 exc dup mul sub div def % demi-grand axe
         /b_2 a_2 1 exc dup mul sub sqrt mul def % demi-petit axe
         /c_2 a_2 exc mul def % distance focale
         % phase
         /Phi theta_0 par r0 div 1 sub exc div arccos sub def
         /rP par 1 exc add div def
         /rA par 1 exc sub div def
% vitesses à l'apogée et au périgée
         /vA G Mt mul 2 rA div 1 a_2 div sub mul sqrt def
         /vP G Mt mul 2 rP div 1 a_2 div sub mul sqrt def
% positions du périgée et de l'apogée
         /xP rP Phi cos mul def
         /yP rP Phi sin mul def
         /xA rA Phi cos mul neg def
         /yA rA Phi sin mul neg def
         /xW xA xP add 2 div def
         /yW yA yP add 2 div def
% periode
       /periode 6.28 a_2 3 exp G div Mt div sqrt mul def
       /radius {par 1 exc t Phi sub cos mul add div} def
       /xE {radius t cos mul} def
       /yE {radius t sin mul}def
        }%
%%  0  1   2   3  4  5   6   7
%% x1 y1 x'1 y'1 x2 y2 x'2 y'2
\def\InitCondred{ x0  y0  v0x  v0y}
\def\InitCond{ x01  y01  v0x1  v0y1 x02 y02  v0x2   v0y2}
\psset{method=rk4}
\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=8pt](-4,-4)(12,12)%
\psequadiff[whichabs=0,whichord=1,
            plotpoints=1000,algebraic,
            tabname=X1Y1
%            ,saveData,filename=XYM1.dat
            ]{0}{120}{\InitCond}{\GravAlgIIcorps}
\listplot[unit=1,linecolor=red]{X1Y1 aload pop}
\psequadiff[whichabs=4,whichord=5,
            plotpoints=1000,algebraic,
            tabname=X2Y2
%            ,saveData,filename=XYM2.dat
            ]{0}{120}{\InitCond}{\GravAlgIIcorps}
\listplot[unit=1,linecolor=blue]{X2Y2 aload pop}
% mouvement de M2 par rapport à M1
%\psequadiff[plotpoints=1000,algebraic,
%            plotfuncx=y dup 4 get exch 0 get sub ,
%            plotfuncy=dup 5 get exch 1 get sub,
%            tabname=XY1]{0}{22.7}{\InitCond}{\GravAlgIIcorps}
%\listplot[unit=1,linecolor=green]{XY1 aload pop}
% mouvement de M1 par rapport à M2
%\psequadiff[plotpoints=1000,algebraic,
%            plotfuncx=y dup 0 get exch 4 get sub ,
%            plotfuncy=dup 1 get exch 5 get sub,
%            tabname=XY2]{0}{22.7}{\InitCond}{\GravAlgIIcorps}
%\listplot[unit=1,linecolor=magenta]{XY2 aload pop}
% mouvement de M1 par rapport à G
%\psequadiff[plotpoints=1000,algebraic,
%            plotfuncx=y 0 get
%                      y 4 get M2 mul
%                      y 0 get M1 mul add
%                      Mt div sub ,
%            plotfuncy=1 get
%                      y 1 get
%                      M1 mul
%                      y 5 get M2 mul add
%                      Mt div sub,
%            tabname=XY3]{0}{22.7}{\InitCond}{\GravAlgIIcorps}
%\listplot[unit=1]{XY3 aload pop}
% mouvement de M2 par rapport à G
%\psequadiff[plotpoints=1000,algebraic,
%            plotfuncx=y 4 get
%                      y 4 get M2 mul
%                      y 0 get M1 mul add
%                      Mt div sub ,
%            plotfuncy=5 get
%                      y 1 get M1 mul
%                      y 5 get M2 mul add
%                     Mt div sub,
%            tabname=XY4]{0}{22.7}{\InitCond}{\GravAlgIIcorps}
%\listplot[unit=1,linecolor=red]{XY4 aload pop}
% centre de masse
\psequadiff[plotpoints=1000,algebraic,
            plotfuncx=y dup 0 get M1 mul exch
                      4 get M2 mul add
                      Mt div,
            plotfuncy=dup 1 get M1 mul exch
                      5 get M2 mul add
                      Mt div,
            tabname=XGYG
%            ,saveData,filename=XYG.dat
            ]{0}{120}{\InitCond}{\GravAlgIIcorps}
\listplot[unit=1,linecolor=cyan]{XGYG aload pop}
\pnode(!x01 y01){M01}
\pnode(!x02 y02){M02}
\pnode(!xG0 yG0){G0}
\psdots(M01)(M02)
\psdot[linecolor=cyan](G0)
\rput(M01){\psline[unit=5,linecolor=red]{->}(!v0x1 v0y1)}
\rput(M02){\psline[unit=5,linecolor=blue]{->}(!v0x2 v0y2)}
\rput(G0){\psline[unit=5]{->}(!vG0x vG0y)}
\psline{<->}(12,0)(0,0)(0,12)
\rput(-2,0){\psPrintValue[decimals=4]{periode}\hphantom{00000}s}
\end{pspicture}
\end{center}
La partie suivante est une animation en \textit{temps réel}, ce qui permet de visualiser l'effet de fronde (accélération), lorsque les deux corps passent au voisinage l'un de l'autre.

Pour faciliter les calculs les positions successives des deux corps et du centre de masse ont été enregistrées dans trois fichiers. Voici pour le premier :
\begin{verbatim}
\psequadiff[whichabs=0,whichord=1,
            plotpoints=1000,algebraic,
            tabname=X1Y1,
            saveData,filename=XYM1.dat]{0}{120}{\InitCond}{\GravAlgIIcorps}
\end{verbatim}
Cette option est ensuite désactivée. Pour l'animation, les fichiers sont lus et les données placées dans un tableau.
\begin{verbatim}
\pstVerb{/XY1 [(XYM1.dat) run] def
         /XY2 [(XYM2.dat) run] def
         /XYG [(XYG.dat) run] def}%
\end{verbatim}
On pourra tracer un segment joignant les deux corps en activant la ligne :
\begin{verbatim}
    %\psline[linestyle=dashed,linecolor=white](!X1 Y1)(!X2 Y2)
\end{verbatim}

\newpage
\section{L'animation : version 1}
\begin{center}
\pstVerb{/XY1 [(XYM1.dat) run] def
         /XY2 [(XYM2.dat) run] def
         /XYG [(XYG.dat) run] def}%
\def\nFrames{200}% 200 images
\begin{animateinline}[controls,loop,timeline=pb2corps.dat,%
                     begin={\begin{pspicture}(-4,-5)(13,12)},
                     end={\end{pspicture}}]{10}% 10 images/s
\psframe*(-4,-5)(13,12)
\psgrid[subgriddiv=0,gridcolor=white,griddots=10,gridlabels=0pt]%
\listplot[linecolor=gray]{XY1 aload pop}
\listplot[linecolor=gray]{XY2 aload pop}
\listplot[linecolor=gray]{XYG aload pop}
\newframe
\multiframe{\nFrames}{i=0+10,I=2+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
         /XG XYG \i\space get def
         /YG XYG \i\space 1 add get def
        }%
%\psdot(!X1 Y1)
\pscircle*[linecolor=red](!X1 Y1){0.3}
\pscircle*[linecolor=blue!50](!X2 Y2){0.1}
\pscircle*[linecolor=yellow](!XG YG){0.05}
%\psline[linestyle=dashed,linecolor=white](!X1 Y1)(!X2 Y2)
}
\end{animateinline}
\end{center}
\newpage
\section{L'animation : version 2}
\begin{center}
\pstVerb{/XY1 [(XYM1.dat) run] def
         /XY2 [(XYM2.dat) run] def
         /XYG [(XYG.dat) run] def}%
\def\nFrames{200}% 200 images
\begin{animateinline}[controls,loop,timeline=pb2corps.dat,%
                     begin={\begin{pspicture}(-4,-5)(13,12)},
                     end={\end{pspicture}}]{10}% 10 images/s
\psframe*(-4,-5)(13,12)
\psgrid[subgriddiv=0,gridcolor=white,griddots=10,gridlabels=0pt]%
\listplot[linecolor=gray]{XY1 aload pop}
\listplot[linecolor=gray]{XY2 aload pop}
\listplot[linecolor=gray]{XYG aload pop}
\newframe
\multiframe{\nFrames}{i=0+10,I=2+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
         /XG XYG \i\space get def
         /YG XYG \i\space 1 add get def
     /xy1 [
         0 2 \I\space {/i exch def
         XY1 i get
          XY1 i 1 add get
        } for
         ] def
     /xy2 [
         0 2 \I\space {/i exch def
         XY2 i get
         XY2 i 1 add get
        } for
         ] def
     /xyG [
         0 2 \I\space {/i exch def
         XYG i get
         XYG i 1 add get
        } for
         ] def
}%
\psgrid[subgriddiv=0,gridcolor=white,griddots=10,gridlabels=0pt]%
\listplot[linecolor=red]{xy1 aload pop}
\listplot[linecolor=blue!50]{xy2 aload pop}
\listplot[linecolor=yellow]{xyG aload pop}
\pscircle*[linecolor=red](!X1 Y1){0.3}
\pscircle*[linecolor=blue!50](!X2 Y2){0.1}
\pscircle*[linecolor=yellow](!XG YG){0.05}
%\psline[linestyle=dashed,linecolor=white](!X1 Y1)(!X2 Y2)
}
\end{animateinline}
\end{center}
\end{document} 