gravitation/pb_2corps_doc.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 \newpsstyle{vecteurA}{arrowinset=0.05,arrowsize=0.1,linecolor={[rgb]{1 0.5 0}}}
  10 \newpsstyle{vecteurB}{arrowinset=0.05,arrowsize=0.1,linecolor={[rgb]{0 0.5 1}}}
  11 \newpsstyle{vecteurC}{arrowinset=0.1,arrowsize=0.2,linecolor={[rgb]{1 0 0}}}
  12 %%%%%%%%%%%%%%%%%%
  13 \title{Gravitation : le problème des deux corps avec PSTricks\\ partie 1}
  14 \date{3 juillet 2\,012}
  15 \begin{document}
  16 \maketitle
  17 \section{Présentation}
  18  Cette première partie aborde uniquement le problème théorique et l'établissement des relations et formules indispensables à la réalisation des figures avec `PSTricks' et à l'animation avec le package `\textsf{animate}'. Elle pourra donc sembler superflue à tous ceux qui connaissent bien ce problème classique. De très bons livres traitent le \textit{problème des deux corps} de façon très claire comme celui de José-Philippe Pérez dans `\textit{Mécanique}' aux éditions Masson ou de façon très complète(et très claire aussi) comme celui publié par le \textsc{cnes} aux éditions \textsc{cepadues} : `\textit{Le mouvement du satellite}'.
  19
  20  La deuxième partie détaillera la procédure suivie pour réaliser les schémas et les animations. Toutefois, on peut, dès à présent, avoir accès au code des figures  dans le fichier source de ce document.
  21 \section{Étude théorique}
  22 \newcommand\GravAlg{%
  23   y[2]|y[3]|%
  24   M1*(y[4]-y[0])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5|%
  25   M1*(y[5]-y[1])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5|%
  26   y[6]|y[7]|%
  27   M2*(y[0]-y[4])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5|%
  28   M2*(y[1]-y[5])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5}
  29 %%  0  1   2   3  4  5   6   7
  30 %% x1 y1 x'1 y'1 x2 y2 x'2 y'2
  31 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.
  32 \begin{center}
  33 \begin{pspicture}(-3,-1)(5,6)
  34 \psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]%
  35 \pstVerb{/x01 -2 def /y01 1 def
  36          /x02  5 def /y02 4 def
  37          /xr0 x02 x01 sub def
  38          /yr0 y02 y01 sub def
  39          /M1 3 def
  40          /M2 1 def
  41          /Mt M1 M2 add def
  42          /xG0 x01 M1 mul x02 M2 mul add Mt div def
  43          /yG0 y01 M1 mul y02 M2 mul add Mt div def}%
  44 \pnode(0,0){O}
  45 \pnode(!x01 y01){M1}
  46 \pnode(!x02 y02){M2}
  47 \pnode(!xG0 yG0){C}
  48 \pscircle[fillstyle=solid,fillcolor=red!50](M1){0.21}
  49 \pscircle[fillstyle=solid,fillcolor=blue!50](M2){0.07}
  50 \psdot(C)
  51 \psline[linestyle=dotted](M1)(M2)
  52 \psline[linestyle=dashed](M2)(O)(M1)
  53 \psline{<->}(5,0)(0,0)(0,6)
  54 \uput[dl](O){$O$}
  55 \uput[r](0,5.8){$y$}
  56 \uput[u](4.8,0){$x$}
  57 \uput[l](0,5.8){$\mathcal{R}$}
  58 \uput{0.22}[l](M1){$M_1$}
  59 \uput{0.1}[u](M2){$M_2$}
  60 \uput[u](C){$C$}
  61 \pcline[offset=5pt,linestyle=none](M1)(M2)
  62 \ncput*[nrot=:U]{$r$}
  63 \rput(M1){\psline[style=vecteurC]{->}(!xr0 5 div yr0 5 div)\uput[u](!xr0 10 div yr0 10 div){$\overrightarrow{F}_{2/1}$}}
  64 \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}$}}
  65 \end{pspicture}
  66 \end{center}
  67 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 :
  68 \[
  69 \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}
  70 \]
  71 Bien sûr : $\overrightarrow{F}_{1/2}+\overrightarrow{F}_{2/1}=\overrightarrow{0}$. Par conséquent, le centre de masse $C$ du système $\{M_1,M_2\}$ est, dans le repère $\mathcal{R}$, soit immobile, soit en mouvement de translation rectiligne uniforme suivant les conditions initiales.
  72
  73 Dans $\mathcal{R}$ la position de $C$ est déterminée par la relation :
  74 \[
  75 \overrightarrow{OC}=\frac{m_1\overrightarrow{OM_1}+m_2\overrightarrow{OM_2}}{m_1+m_2}
  76 \]
  77 \begin{equation}
  78 \overrightarrow{v_{C/\mathcal{R}}}=\frac{m_1\overrightarrow{v_1}+m_2\overrightarrow{v_2}}{m_1+m_2}=\overrightarrow{v_0}
  79 \label{vc}
  80 \end{equation}
  81 $\overrightarrow{v_0}$ étant déterminé par les valeurs initiales de $\overrightarrow{v_1}$ et $\overrightarrow{v_2}$.
  82 L'étude peut donc se poursuivre dans le repère lié au centre de masse $\mathcal{R}^*$, lui-même galiléen.
  83
  84 On pose : $\overrightarrow{r_1}=\overrightarrow{CM_1}$ et $\overrightarrow{r_2}=\overrightarrow{CM_2}$. On a toujours : $\overrightarrow{r}=\overrightarrow{M_1M_2}=\overrightarrow{r_2}-\overrightarrow{r_1}$.
  85 \begin{center}
  86 \begin{pspicture}(-3.5,-1)(5.5,8.25)
  87 %\psframe(-3.5,-1)(5.5,9)
  88 \pstVerb{/x01 -2 def /y01 1 def
  89          /x02  4 def /y02 5 def
  90          /xr0 x02 x01 sub def
  91          /yr0 y02 y01 sub def
  92          /M1 3 def
  93          /M2 1 def
  94          /Mt M1 M2 add def
  95          /xG0 x01 M1 mul x02 M2 mul add Mt div def
  96          /yG0 y01 M1 mul y02 M2 mul add Mt div def}%
  97 \pnode(0,0){O}
  98 \pnode(!x01 y01){M1}
  99 \pnode(!x02 y02){M2}
 100 \pnode(!xG0 yG0){C}
 101 \pscircle[fillstyle=solid,fillcolor=red!50](M1){0.21}
 102 \pscircle[fillstyle=solid,fillcolor=blue!50](M2){0.07}
 103 \psdot(C)
 104 \psline[linestyle=dotted](M1)(M2)
 105 \psline[linestyle=dashed](M2)(O)(M1)
 106 \psline[linecolor=gray]{<->}(5,0)(0,0)(0,6)
 107 \rput(C){\psline[style=vecteurA]{<->}(5.5,0)(0,0)(0,6)%
 108          \uput[l](0,5.8){$\mathcal{R^*}$}
 109          \psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt](-3,-3)(6,6)%
 110          \uput[r](0,5.8){$y^*$}
 111          \uput[u](5.5,0){$x^*$}}
 112 \uput[dl](O){$O$}
 113 \uput[r](0,5.8){$y$}
 114 \uput[u](4.8,0){$x$}
 115 \uput[l](0,5.8){$\mathcal{R}$}
 116 \uput{0.22}[l](M1){$M_1$}
 117 \uput{0.1}[u](M2){$M_2$}
 118 \uput[ul](C){$C$}
 119 \pcline[offset=5pt,linestyle=none](M1)(M2)
 120 \ncput*[nrot=:U]{$r$}
 121 \rput(M1){\psline[style=vecteurC]{->}(!xr0 5 div yr0 5 div)\uput[u](!xr0 10 div yr0 10 div){$\overrightarrow{F}_{2/1}$}}
 122 \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}$}}
 123 \end{pspicture}
 124 \end{center}
 125
 126 Dans le repère $\mathcal{R}^*$, la quantité de mouvement du système $\{M_1,M_2\}$ est nulle : $\overrightarrow{p}^*=(m_1+m_2)\overrightarrow{v}_C^*=\overrightarrow{0}$.
 127 Pour chacun des deux corps, la quantité de mouvement dans $\mathcal{R}^*$ est :
 128 \[
 129 \overrightarrow{p}_1^*=m_1\overrightarrow{v}_1^*=m_1(\overrightarrow{v}_{1/\mathcal{R}}-\overrightarrow{v}_{C/\mathcal{R}})=-\overrightarrow{p}^*_2
 130 \]
 131 D'après (\ref{vc}) qui donne  $\overrightarrow{v}_{C/\mathcal{R}}$, on obtient :
 132 \[
 133 \overrightarrow{p}_1^*=m_1\left(\overrightarrow{v}_{1/\mathcal{R}}-\frac{m_1\overrightarrow{v}_{1/\mathcal{R}}+m_2\overrightarrow{v}_{2/\mathcal{R}}}{m_1+m_2}\right)=%
 134  \frac{m_1m_2}{m_1+m_2}(\overrightarrow{v}_{1/\mathcal{R}}-\overrightarrow{v}_{2/\mathcal{R}})=-\overrightarrow{p}^*_2
 135 \]
 136 Puisque, on a défini $\overrightarrow{r}$ par $\overrightarrow{r_2}-\overrightarrow{r_1}$, en posant $\mu\dfrac{m_1m_2}{m_1+m_2}$ :
 137 \[
 138 \overrightarrow{p}^*_1=-\mu\dfrac{\mathrm{d}\overrightarrow{r}}{\mathrm{d}t}\qquad \overrightarrow{p}^*_2=\mu\dfrac{\mathrm{d}\overrightarrow{r}}{\mathrm{d}t}
 139 \]
 140 La quantité de mouvement de $M_2$ est identique à celle d'une particule fictive de masse $\mu$, appelée \textit{masse réduite}, et de vitesse $\overrightarrow{v}_{2/\mathcal{R}}-\overrightarrow{v}_{1/\mathcal{R}}$ : vitesse relative de $M_2$ par rapport à $M_1$, idem pour $M_1$.
 141
 142 On considère cette particule fictive de masse $\mu$ située au point $M$ tel que : $\overrightarrow{CM}=\overrightarrow{r}$.
 143
 144 La relation fondamentale de la dynamique appliquée à chacun des deux corps donne :
 145 \[
 146 \left\{
 147 \begin{array}[m]{l}
 148  \dfrac{\mathrm{d}\overrightarrow{p}^*_1}{\mathrm{d}t}=\overrightarrow{F}_{2/1}\\[1em]
 149 \dfrac{\mathrm{d}\overrightarrow{p}^*_2}{\mathrm{d}t}=\overrightarrow{F}_{1/2}
 150 \end{array}
 151 \right.
 152 \]
 153 En remplaçant $\overrightarrow{p}^*_1$ (ou $\overrightarrow{p}^*_2$) par son expression en fonction de $\mu$, on a :
 154 \begin{equation}
 155  \mu\dfrac{\mathrm{d}^2\overrightarrow{r}}{\mathrm{d}t^2}=\overrightarrow{F}_{1/2}=-\mathcal{G}\frac{m_1m_2}{r^3}\overrightarrow{r}
 156 \label{RFD}
 157 \end{equation}
 158 Nous désignons cette force par $\overrightarrow{F}$. Remarquons qu'elle peut s'écrire, en notant la masse totale $(m_1+m_2)$  du système $\{M_1,M_2\}$ par $m_t$~:
 159 \begin{equation}
 160 \overrightarrow{F}=-\mathcal{G}\frac{\mu m_t}{r^3}\overrightarrow{r}
 161 \label{RF}
 162 \end{equation}
 163 En conclusion, cette particule fictive $M$, de masse $\mu$ positionnée en $\overrightarrow{CM}=\overrightarrow{r}$ est attirée par un corps fictif de masse égale à la masse totale du système $(m_1+m_2)$, placé en $C$. Sa position et sa vitesse initiales sont déterminées par les positions et vitesses initiales de $M_1$ et $M_2$.
 164 \[
 165 \overrightarrow{r}_0=\overrightarrow{r}_{2_0}-\overrightarrow{r}_{1_0} \quad \mbox{et} \quad \overrightarrow{v}^*_0=\overrightarrow{v}_{0/R}=\overrightarrow{v}_{2_0}^*-\overrightarrow{v}_{1_0}^*=
 166 \overrightarrow{v}_{{0_2}/R}-\overrightarrow{v}_{{0_1}/R}
 167 \]
 168 Si nous connaissons le mouvement de $M$ nous pourrons en déduire les mouvements respectifs de $M_1$ et $M_2$, puisque d'après la définition du centre de masse :
 169 \begin{equation}
 170 \overrightarrow{r_1}=-\frac{m_2}{m_1+m_2}\overrightarrow{r}\quad\mbox{et}\quad \overrightarrow{r_2}=\frac{m_1}{m_1+m_2}\overrightarrow{r}
 171 \label{r1r2}
 172 \end{equation}
 173 Les trajectoires de $M_1$ et $M_2$ se déduiront de celle de $M$ par une homothétie de centre C.
 174 \begin{center}
 175 \begin{pspicture}(-3.5,-1)(6.5,8.25)
 176 %\psframe(-3.5,-1)(5.5,9)
 177 \pstVerb{/x01 -2 def /y01 1 def
 178          /x02  4 def /y02 5 def
 179          /xr0 x02 x01 sub def
 180          /yr0 y02 y01 sub def
 181          /M1 3 def
 182          /M2 1 def
 183          /Mt M1 M2 add def
 184          /xG0 x01 M1 mul x02 M2 mul add Mt div def
 185          /yG0 y01 M1 mul y02 M2 mul add Mt div def}%
 186 \pnode(0,0){O}
 187 \pnode(!x01 y01){M1}
 188 \pnode(!x02 y02){M2}
 189 \pnode(!xG0 yG0){C}
 190 \pscircle[fillstyle=solid,fillcolor=red!50](M1){0.21}
 191 \pscircle[fillstyle=solid,fillcolor=blue!50](M2){0.07}
 192 \pscircle[hatchcolor=gray,fillstyle=vlines,hatchwidth=0.02,hatchsep=0.05](C){0.25}
 193 \psdot(C)
 194 \psline[linestyle=dotted](M1)(M2)
 195 \psline[linestyle=dashed](M2)(O)(M1)
 196 \psline[linecolor=gray]{<->}(5,0)(0,0)(0,6)
 197 \rput(C){\psline[style=vecteurA]{<->}(5.5,0)(0,0)(0,6)%
 198          \uput[l](0,5.8){$\mathcal{R^*}$}
 199          \psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt](-3,-3)(7,6)%
 200          \uput[r](0,5.8){$y^*$}
 201          \uput[u](5.5,0){$x^*$}
 202          \pnode(!xr0 yr0){M}
 203          \pscircle*(M){0.05}
 204          \pcline[offset=0.3,linewidth=0.01,arrowsize=0.1]{|->}(0,0)(M)
 205          \ncput*[nrot=:U]{$\overrightarrow{r}$}
 206          \psline[linestyle=dotted](M)\uput[r](M){$M(\mu)$}
 207          \rput(M){\psline[style=vecteurC]{->}(!xr0 5 div neg yr0 5 div neg)\uput[dr](!xr0 10 div neg yr0 10 div neg){$\overrightarrow{F}_{1/2}$}}}
 208 \uput[dl](O){$O$}
 209 \uput[r](0,5.8){$y$}
 210 \uput[u](4.8,0){$x$}
 211 \uput[l](0,5.8){$\mathcal{R}$}
 212 \uput{0.22}[l](M1){$M_1$}
 213 \uput{0.1}[dr](M2){$M_2$}
 214 \uput{0.3}[d](C){$C$}
 215 %\pcline[offset=10pt,linestyle=none](M1)(M2)
 216 %\ncput*[nrot=:U]{$r$}
 217 %\rput(M1){\psline[style=vecteurC]{->}(!xr0 5 div yr0 5 div)\uput[u](!xr0 10 div yr0 10 div){$\overrightarrow{F}_{2/1}$}}
 218 %\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}$}}
 219 \end{pspicture}
 220 \end{center}
 221 Pour étudier, dans $\mathcal{R}^*$, le mouvement de ce point matériel fictif $M$ de masse $(\mu)$ soumis à la force centrale $\overrightarrow{F}$, il est avantageux de passer en coordonnées polaires.
 222 \begin{center}
 223 \begin{pspicture}(-3.5,-1)(6.5,8.25)
 224 %\psframe(-3.5,-1)(5.5,9)
 225 \pstVerb{/x01 -2 def /y01 1 def
 226          /x02  4 def /y02 5 def
 227          /xr0 x02 x01 sub def
 228          /yr0 y02 y01 sub def
 229          /theta_0 yr0 xr0 atan def
 230          /M1 3 def
 231          /M2 1 def
 232          /Mt M1 M2 add def
 233          /xG0 x01 M1 mul x02 M2 mul add Mt div def
 234          /yG0 y01 M1 mul y02 M2 mul add Mt div def}%
 235 \pnode(0,0){O}
 236 \pnode(!x01 y01){M1}
 237 \pnode(!x02 y02){M2}
 238 \pnode(!xG0 yG0){C}
 239 \pscircle[fillstyle=solid,fillcolor=red!50](M1){0.21}
 240 \pscircle[fillstyle=solid,fillcolor=blue!50](M2){0.07}
 241 \pscircle[hatchcolor=gray,fillstyle=vlines,hatchwidth=0.02,hatchsep=0.05](C){0.25}
 242 \psdot(C)
 243 \psline[linestyle=dotted](M1)(M2)
 244 \psline[linestyle=dashed](M2)(O)(M1)
 245 \psline[linecolor=gray]{<->}(5,0)(0,0)(0,6)
 246 \rput(C){\psline[style=vecteurA]{<->}(5.5,0)(0,0)(0,6)%
 247          \uput[l](0,5.8){$\mathcal{R^*}$}
 248          \psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt](-3,-3)(7,6)%
 249          \uput[r](0,5.8){$y^*$}
 250          \uput[u](5.5,0){$x^*$}
 251          \pnode(!xr0 yr0){M}
 252          \pscircle*(M){0.05}
 253          \pcline[offset=0.3,linewidth=0.01,arrowsize=0.1]{|->}(0,0)(M)
 254          \ncput*[nrot=:U]{$\overrightarrow{r}$}
 255          \psline[linestyle=dotted](M)\uput[r](M){$M(\mu)$}
 256          \rput(M){\psline[style=vecteurC]{->}(!xr0 5 div neg yr0 5 div neg)\uput[dr](!xr0 10 div neg yr0 10 div neg){$\overrightarrow{F}$}}}
 257 \uput[dl](O){$O$}
 258 \uput[r](0,5.8){$y$}
 259 \uput[u](4.8,0){$x$}
 260 \uput[l](0,5.8){$\mathcal{R}$}
 261 \uput{0.22}[l](M1){$M_1$}
 262 \uput{0.1}[dr](M2){$M_2$}
 263 \uput{0.3}[d](C){$C$}
 264 \rput(C){\psline[arrowinset=0.1,arrowsize=0.15]{->}(!theta_0 90 add cos theta_0 90 add sin)
 265          \psline[arrowinset=0.1,arrowsize=0.15]{->}(!theta_0 cos theta_0 sin)
 266          \uput{.1}[u](!theta_0 90 add cos theta_0 90 add sin){$\overrightarrow{u}_\theta$}
 267          \uput[d](!theta_0 cos theta_0 sin){$\overrightarrow{u}_r$}
 268          \psarc{->}(0,0){1.5}{0}{!theta_0}
 269          \uput{1.6}[!theta_0 2 div](0,0){$\theta$}}
 270 \end{pspicture}
 271 \end{center}
 272 Pour la suite, retenons que dans $\mathcal{R}*$ :
 273 \[
 274 \overrightarrow{u}_r=
 275 \left|
 276 \begin{array}[m]{l}
 277 \cos\theta\\
 278 \sin\theta
 279 \end{array}
 280 \right.
 281 \quad
 282 \overrightarrow{u}_\theta=
 283 \left|
 284 \begin{array}[m]{l}
 285 -\sin\theta\\
 286 \cos\theta
 287 \end{array}
 288 \right.
 289 \]
 290 \[
 291 \frac{\mathrm{d}\overrightarrow{u}_r}{\mathrm{d}t}%=
 292 %\left|
 293 %\begin{array}[m]{l}
 294 %-\dot{\theta}\sin\theta\\
 295 %\dot{\theta}\cos\theta
 296 %\end{array}
 297 %\right.
 298 =\dot{\theta}\overrightarrow{u}_\theta
 299 \quad
 300 \frac{\mathrm{d}\overrightarrow{u}_\theta}{\mathrm{d}t}=
 301 %\left|
 302 %\begin{array}[m]{l}
 303 %-\dot{\theta}\cos\theta\\
 304 %-\dot{\theta}\sin\theta
 305 %\end{array}
 306 %\right.
 307 -\dot{\theta}\overrightarrow{u}_r
 308 \]
 309 Dans ce repère, la position de $M$ et sa vitesse ont pour expressions respectives :
 310 \[
 311 \overrightarrow{r}=r \overrightarrow{u}_r\quad ;\quad \overrightarrow{v}=\frac{\mathrm{d}\overrightarrow{r}}{\mathrm{d}t}=\frac{\mathrm{d}r}{\mathrm{d}t}\overrightarrow{u}_r + r\frac{\mathrm{d}\overrightarrow{u}_r}{\mathrm{d}t}=\frac{\mathrm{d}r}{\mathrm{d}t}\overrightarrow{u}_r + r\dot{\theta}\overrightarrow{u}_\theta
 312 \]
 313 l'accélération s'écrit :
 314 \[
 315 \overrightarrow{\gamma}=\frac{\mathrm{d}\overrightarrow{v}}{\mathrm{d}t}=\ddot{r}\overrightarrow{u}_r +
 316  \dot{r}\dot{\theta}\overrightarrow{u}_\theta+
 317  \dot{r}\dot{\theta}\overrightarrow{u}_\theta +
 318  r\ddot{\theta}\overrightarrow{u}_\theta -r\dot{\theta}^2\overrightarrow{u}_r
 319 \]
 320 \[
 321 =(\ddot{r} -r\dot{\theta}^2)\overrightarrow{u}_r +
 322  (2\dot{r}\dot{\theta}+r\ddot{\theta})\overrightarrow{u}_\theta
 323 \]
 324 En résumé :
 325 \[
 326 \overrightarrow{v}=
 327 \left|
 328 \begin{array}{l}
 329 \dot{r}\\
 330 r\dot{\theta}
 331 \end{array}
 332 \right.
 333 \quad  \quad
 334 \overrightarrow{\gamma}=
 335 \left|
 336 \begin{array}{l}
 337 \ddot{r} -r\dot{\theta}^2\\
 338 \dfrac{1}{r}\dfrac{\mathrm{d}}{\mathrm{d}t}(r^2\dot{\theta})
 339 \end{array}
 340 \right.
 341 \]
 342 L'équation (\ref{RFD}) nous donne, en divisant par $\mu$, l'accélération :
 343 \[
 344  \dfrac{\mathrm{d}^2\overrightarrow{r}}{\mathrm{d}t^2}=-\mathcal{G}\frac{m_t}{r^2}\overrightarrow{u}_r
 345 \]
 346 Et puisque l'accélération radiale est nulle :
 347 \[
 348 \dfrac{1}{r}\dfrac{\mathrm{d}}{\mathrm{d}t}(r^2\dot{\theta})=0
 349 \]
 350 D'où :
 351 \[
 352 r^2\dot{\theta}=\mathrm{C^{ste}}
 353 \]
 354 Il est intéressant ici de calculer le moment cinétique de $M$ et d'exprimer la loi des aires.
 355
 356 Le moment cinétique :
 357 \[
 358 \overrightarrow{\sigma}=\overrightarrow{r} \wedge \mu\overrightarrow{v}
 359 \]
 360 $\overrightarrow{F}$ étant dirigée vers $C$, son moment par rapport à $C$ est nul. Le moment cinétique est donc constant. Il en découle que le mouvement est plan, puisque le vecteur-position $\overrightarrow{r}$ reste toujours perpendiculaire à un vecteur constant. Ce plan est le plan perpendiculaire à $\sigma_0$ contenant $C$.
 361
 362 On peut calculer le moment cinétique $\overrightarrow{\sigma}$, à partir des coordonnées polaires :
 363 \[
 364 \overrightarrow{\sigma}=\mu
 365 \left|
 366 \begin{array}{l}
 367 r\\
 368 0\\
 369 0
 370 \end{array}
 371 \right.
 372 \wedge
 373 \left|
 374 \begin{array}{l}
 375 \dot{r}\\
 376 r\dot{\theta}\\
 377 0
 378 \end{array}
 379 \right.
 380 =
 381 \left|
 382 \begin{array}{l}
 383 0\\
 384 0\\
 385 r^2\dot{\theta}
 386 \end{array}
 387 \right.
 388 =\mu r^2\dot{\theta}\overrightarrow{u}_z
 389 \]
 390 La constante des aires(voir plus loin pourquoi elle s'appelle ainsi) est $\mathrm{C}=\dfrac{\sigma}{\mu}=r^2\dot{\theta}$, c'est aussi la valeur du moment cinétique par unité de masse.
 391
 392 On peut aussi calculer le moment cinétique en coordonnées cartésiennes. Les conditions initiales étant fixées, $\overrightarrow{v}_0$ est donné par ses coordonnées, ainsi que $\overrightarrow{r}_0$.
 393 \[
 394 \overrightarrow{\sigma}=\mu
 395 \left|
 396 \begin{array}{l}
 397 x_0\\
 398 y_0\\
 399 0
 400 \end{array}
 401 \right.
 402 \wedge
 403 \left|
 404 \begin{array}{l}
 405 v_{x_0}\\
 406 v_{y_0}\\
 407 0
 408 \end{array}
 409 \right.
 410 =
 411 \left|
 412 \begin{array}{l}
 413 0\\
 414 0\\
 415 x_0 v_{y_0} - y_0 v_{x_0}
 416 \end{array}
 417 \right.
 418 =\mu(x_0 v_{y_0} - y_0 v_{x_0})\overrightarrow{u}_z
 419 \]
 420 Les livres d'astronomie\footnote{page 43 dans \textit{Le mouvement du satellite} : conférences et exercices de mécanique spatiale. Cepadues-éditions 1983.} préfèrent définir la vitesse par un autre paramètre $\gamma$, \textit{angle de l'horizontale locale }$\overrightarrow{u}_z\wedge \overrightarrow{u}_r$ (c'est-à-dire $\overrightarrow{u}_\theta$) \textit{avec la vitesse}, en se plaçant toujours en coordonnées polaires.
 421 \begin{center}
 422 \begin{pspicture}(-3.5,-1)(6.5,8.25)
 423 %\psframe(-3.5,-1)(5.5,9)
 424 \pstVerb{/x01 -2 def /y01 1 def
 425          /x02  4 def /y02 5 def
 426          /xr0 x02 x01 sub def
 427          /yr0 y02 y01 sub def
 428          /theta_0 yr0 xr0 atan def
 429          /M1 3 def
 430          /M2 1 def
 431          /Mt M1 M2 add def
 432          /xG0 x01 M1 mul x02 M2 mul add Mt div def
 433          /yG0 y01 M1 mul y02 M2 mul add Mt div def}%
 434 \pnode(0,0){O}
 435 \pnode(!x01 y01){M1}
 436 \pnode(!x02 y02){M2}
 437 \pnode(!xG0 yG0){C}
 438 %\pscircle[fillstyle=solid,fillcolor=red!50](M1){0.21}
 439 %\pscircle[fillstyle=solid,fillcolor=blue!50](M2){0.07}
 440 \pscircle[hatchcolor=gray,fillstyle=vlines,hatchwidth=0.02,hatchsep=0.05](C){0.25}
 441 \psdot(C)
 442 %\psline[linestyle=dotted](M1)(M2)
 443 %\psline[linestyle=dashed](M2)(O)(M1)
 444 \psline[linecolor=gray]{<->}(5,0)(0,0)(0,6)
 445 \rput(C){\psline[style=vecteurA]{<->}(5.5,0)(0,0)(0,6)%
 446          \uput[l](0,5.8){$\mathcal{R^*}$}
 447          \psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt](-3,-3)(7,6)%
 448          \uput[r](0,5.8){$y^*$}
 449          \uput[u](5.5,0){$x^*$}
 450          \pnode(!xr0 yr0){M}
 451          \pscircle*(M){0.05}
 452 %         \pcline[offset=0.3,linewidth=0.01,arrowsize=0.1]{|->}(0,0)(M)
 453 %         \ncput*[nrot=:U]{$\overrightarrow{r}$}
 454          \psline[linestyle=dotted](M)\uput[r](M){$M(\mu)$}
 455          \rput(M){\psline[style=vecteurC]{->}(!xr0 5 div neg yr0 5 div neg)\uput[dr](!xr0 10 div neg yr0 10 div neg){$\overrightarrow{F}$}
 456          \psline[linestyle=dotted](!theta_0 90 add cos 2 mul theta_0 90 add sin 2 mul)
 457          \psline[arrowinset=0.1,arrowsize=0.15]{->}(!theta_0 90 add cos theta_0 90 add sin)
 458          \psline[arrowinset=0.1,arrowsize=0.15,linecolor=blue]{->}(!theta_0 45 add cos 2 mul theta_0 45 add sin 2 mul)
 459          \uput[r](!theta_0 45 add cos 1.75 mul theta_0 45 add sin 1.55 mul){$\overrightarrow{v}$}
 460          \psarc{->}(0,0){1.2}{!theta_0 45 add}{!theta_0 90 add}
 461          \uput{1.3}[!theta_0 67.5 add](0,0){$\gamma$}}
 462          }
 463 \uput[dl](O){$O$}
 464 \uput[r](0,5.8){$y$}
 465 \uput[u](4.8,0){$x$}
 466 \uput[l](0,5.8){$\mathcal{R}$}
 467 %\uput{0.22}[l](M1){$M_1$}
 468 %\uput{0.1}[dr](M2){$M_2$}
 469 \uput{0.3}[d](C){$C$}
 470 \rput(C){\psline[arrowinset=0.1,arrowsize=0.15]{->}(!theta_0 90 add cos theta_0 90 add sin)
 471          \psline[arrowinset=0.1,arrowsize=0.15]{->}(!theta_0 cos theta_0 sin)
 472          \uput{.1}[u](!theta_0 90 add cos theta_0 90 add sin){$\overrightarrow{u}_\theta$}
 473          \uput[d](!theta_0 cos theta_0 sin){$\overrightarrow{u}_r$}
 474          \psarc{->}(0,0){1.5}{0}{!theta_0}
 475          \uput{1.6}[!theta_0 2 div](0,0){$\theta$}}
 476 \end{pspicture}
 477 \end{center}
 478 Dans ce cas, on peut exprimer la constante des aires par :
 479 $
 480 \mathrm{C}=rv\cos\gamma
 481 $
 482 \begin{center}
 483 \begin{pspicture}(-2,-2)(7,5)
 484 \psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]
 485 \pstVerb{/par 2 def
 486          /exc 0.7 def
 487          /Phi 30 def
 488          /radius {par 1 exc t Phi sub cos mul add div} def
 489          /xE {radius t cos mul neg} def
 490          /yE {radius t sin mul neg} def}%
 491 \parametricplot[plotpoints=360]{0}{360}{
 492     xE
 493     yE}%
 494 \pscustom[fillstyle=vlines,hatchwidth=0.02,hatchcolor=gray]{%
 495     \psline(0,0)(! /t 200 def xE yE)
 496 \parametricplot[plotpoints=360]{200}{210}{
 497     xE
 498     yE}%
 499    \psline(! /t 210 def xE yE)(0,0)}
 500 \parametricplot[plotpoints=360,arrows=->]{200}{210}{
 501    /radius {2.25 1 exc t Phi sub cos mul add div} def
 502     xE
 503     yE}%
 504 \psarc{->}(0,0){3}{20}{30}
 505 \uput{3.2}[25](0,0){\psframebox[linestyle=none,fillstyle=solid,fillcolor=white,framesep=0pt]{d$\theta$}}
 506 \rput(5,2.3){\psframebox[linestyle=none,fillstyle=solid,fillcolor=white,framesep=0pt]{d$\mathcal{A}$}}
 507 %\psline[linearc=0.25]{->}(2.2,1.8)(3.5,1.8)(3.5,1.5)
 508 \pcline[offset=-0.3,linestyle=none](0,0)(! /t 200 def xE yE)
 509 \ncput*[nrot=:U]{$r$}
 510 \uput[l](0,4.8){$\mathcal{R^*}$}
 511 \uput[r](0,4.8){$y^*$}
 512 \uput[u](7,0){$x^*$}
 513 \uput[dl](0,0){$C$}
 514 \psline[style=vecteurA]{<->}(7,0)(0,0)(0,5)
 515 \end{pspicture}
 516 \end{center}
 517 Géométriquement lorsque $\theta$ varie de $\mathrm{d}\theta$, la surface balayée par le rayon vecteur vaut :
 518 \[
 519 \mathrm{d}\mathcal{A}=\frac{1}{2}r^2\mathrm{d}\theta
 520 \]
 521 On définit ainsi la vitesse \textit{aérolaire}, qui est l'aire balayée par unité de temps, en fonction de C :
 522 \[
 523 \frac{\mathrm{d}A}{\mathrm{d}t}=\frac{\mathrm{C}}{2}
 524 \]
 525 On retrouve ainsi la deuxième loi de Képler.
 526
 527 Pour déterminer l'équation de la trajectoire, on utilise la méthode de Binet, qui consiste en un changement de variable en posant $u=\dfrac{1}{r}$.
 528 \[
 529 \dfrac{\mathrm{d}r}{\mathrm{d}t}=-\dfrac{1}{u^2}\dfrac{\mathrm{d}u}{\mathrm{d}\theta}\dfrac{\mathrm{d}\theta}{\mathrm{d}t}
 530 \]
 531 Nous avons la constante des aires $\mathrm{C}=r^2\dot{\theta}=\dfrac{\dot{\theta}}{u^2}$, par conséquent :
 532 \[
 533 \dfrac{\mathrm{d}r}{\mathrm{d}t}=-\mathrm{C}\dfrac{\mathrm{d}u}{\mathrm{d}\theta}
 534 \]
 535 \[
 536 \dfrac{\mathrm{d}^2r}{\mathrm{d}t^2}=-\mathrm{C}\dfrac{\mathrm{d}^2u}{\mathrm{d}\theta^2}\dfrac{\mathrm{d}\theta}{\mathrm{d}t}
 537 \]
 538 On remplace $\dot{\theta}$ par $\mathrm{C}u^2$ :
 539 \[
 540 \dfrac{\mathrm{d}^2r}{\mathrm{d}t^2}=-\mathrm{C^2}u^2\dfrac{\mathrm{d}^2u}{\mathrm{d}\theta^2}
 541 \]
 542 Nous avons vu qu'en appliquant la loi de Newton, la composante de l'accélération $\gamma$ suivant $\sqrt[n]{u}_r$ vaut :
 543 \[
 544 \ddot{r} -r\dot{\theta}^2=-\mathcal{G}\frac{m_t}{r^2}
 545 \]
 546 Ce qui s'écrit avec la variable $u$ :
 547 \[
 548 -\mathrm{C^2}u^2\dfrac{\mathrm{d}^2u}{\mathrm{d}\theta^2}-\mathrm{C^2}u^3=-\mathcal{G}m_tu^2
 549 \]
 550 Ce qui, après simplification, donne :
 551 \begin{equation}
 552 \psframebox[linestyle=none,fillstyle=solid,fillcolor={[rgb]{1 1 0.5}}]{%
 553 \dfrac{\mathrm{d}^2u}{\mathrm{d}\theta^2}+u=\dfrac{\mathcal{G}m_t}{\mathrm{C^2}}}
 554 \label{equ}
 555 \end{equation}
 556 La solution générale de cette équation peut s'écrire :
 557 \[
 558 u=A\cos(\theta-\varphi)+\dfrac{\mathcal{G}m_t}{\mathrm{C^2}}
 559 \]
 560 $A$ et $\varphi$ sont fixés par les conditions initiales. Le rayon-vecteur a pour expression :
 561 \[
 562 r=\frac{1}{A\cos(\theta-\varphi)+\dfrac{\mathcal{G}m_t}{\mathrm{C^2}}}
 563 \]
 564 Il s'écrit encore ainsi :
 565 \[
 566 r=\dfrac{\dfrac{\mathrm{C^2}}{\mathcal{G}m_t}}{A\dfrac{\mathrm{C^2}}{\mathcal{G}m_t}\cos(\theta-\varphi)+1}
 567 \]
 568 Si on pose $p=\dfrac{\mathrm{C^2}}{\mathcal{G}m_t}$ et $\mathrm{e}=Ap$, on reconnaît l'équation d'une conique dont l'un des foyers est $C$, de paramètre $p$ et d'excentricité $\mathrm{e}$.
 569 \begin{equation}
 570 \psframebox[linestyle=none,fillstyle=solid,fillcolor={[rgb]{1 1 0.5}}]{%
 571 r=\dfrac{p}{1+e\cos(\theta-\varphi)}}
 572 \label{equr}
 573 \end{equation}
 574 Déterminons les constantes en fonction des conditions initiales.
 575
 576 Nous avons vu que la constante des aires C pouvait se calculer de deux façons suivant le choix fait pour définir la vitesse initiale, soit $
 577 \mathrm{C}=r_0v_0\cos\gamma_0$, soit $C=x_0 v_{y_0} - y_0 v_{x_0}$.  Par conséquent $p$ est bien déterminé.
 578
 579 Lorsque $\theta=\varphi$ alors $r=\dfrac{p}{1+\mathrm{e}}$, ce qui est la plus petite valeur possible de $r$. On est donc au périgée P de la conique, on note, en ce point, le rayon-vecteur par $\overrightarrow{r}_\mathrm{P}$. En ce point la vitesse $\overrightarrow{v}_\mathrm{P}$ est perpendiculaire à $\overrightarrow{r}_\mathrm{P}$ et la constante des aires peut s'écrire : $\mathrm{C}=r_\mathrm{P}v_\mathrm{P}$.
 580
 581 L'angle $\varphi$ est l'inclinaison du grand axe de la conique avec $Ox$, plus précisément l'angle du rayon-vecteur au périgée, $\overrightarrow{r}_P$ avec $Ox$.
 582
 583 L'énergie mécanique se conserve au cours du mouvement, elle est la somme de l'énergie cinétique et de l'énergie potentielle. En fonction des conditions initiales elle s'écrit(on rappelle que $\mu$ est la masse(\textit{réduite})) de la particule fictive $M$ étudiée :
 584 \[
 585 \mathcal{E}=\frac{1}{2}\mu v_0^2-\mathcal{G}\frac{\mu(m_1+m_2)}{r_0}\quad \text{ou}\quad \frac{1}{2}\mu v_0^2-\mathcal{G}\frac{m_1m_2}{r_0}
 586 \]
 587 Au périgée, elle s'exprimera par :
 588 \[
 589 \mathcal{E}=\frac{1}{2}\mu v_\mathrm{P}^2-\mathcal{G}\frac{m_1m_2}{r_\mathrm{P}}
 590 \]
 591 Pour exprimer l'énergie dans le cas général, calculons la vitesse en fonction de $u$. Rappelons que :
 592 \[
 593 \overrightarrow{v}=
 594 \left|
 595 \begin{array}{l}
 596 \dot{r}\\
 597 r\dot{\theta}
 598 \end{array}
 599 \right.
 600 \]
 601 Calculons séparément $\dot{r}$ et $r\dot{\theta}$ :
 602 \[
 603 \dot{r}=\dfrac{\mathrm{d}r}{\mathrm{d}u}\dfrac{\mathrm{d}u}{\mathrm{d}\theta}\dfrac{\mathrm{d}\theta}{\mathrm{d}t}=
 604    -\dfrac{1}{u^2}\dfrac{\mathrm{d}u}{\mathrm{d}\theta}\dot{\theta}=
 605    -\mathrm{C}\dfrac{\mathrm{d}u}{\mathrm{d}\theta}
 606 \]
 607 \[
 608 r\dot{\theta}=\dfrac{1}{u}\mathrm{C}u^2=\mathrm{C}u
 609 \]
 610 \[
 611 v^2=\mathrm{C^2}\left[\left(\frac{\mathrm{d}u}{\mathrm{d}\theta}\right)^2+u^2\right]
 612 \]
 613 \[
 614 v^2=\mathrm{C^2}\left[
 615  \frac{\mathrm{e^2}}{p^2}\sin^2(\theta-\varphi)+\frac{1}{p^2}\left[1+2\mathrm{e}\cos(\theta-\varphi)+\mathrm{e^2}\cos^2(\theta-\varphi)\right]
 616   \right]
 617 \]
 618 \[
 619 v^2=\frac{\mathrm{C^2}}{p^2}\left[1+2\mathrm{e}\cos(\theta-\varphi)+\mathrm{e^2} \right]
 620 \]
 621 L'énergie cinétique s'exprime par :
 622 \[
 623 E_C=\frac{1}{2}\mu\frac{\mathrm{C^2}}{p^2}\left[ 1+2\mathrm{e}\cos(\theta-\varphi)+\mathrm{e^2} \right]
 624 \]
 625 On se souvient que $p=\dfrac{\mathrm{C^2}}{\mathcal{G}m_t}$, il vient :
 626 \[
 627 E_C=\frac{1}{2}\frac{\mathcal{G}\mu m_t}{p}\left[ 1+2\mathrm{e}\cos(\theta-\varphi)+\mathrm{e^2} \right]
 628 \]
 629 Pour l'énergie potentielle nous avons :
 630 \[
 631 E_p=-\frac{\mathcal{G}\mu m_t}{r}=-\frac{\mathcal{G}\mu m_t}{p}\left[1+e\cos(\theta-\varphi)\right]
 632 \]
 633 et pour l'énergie mécanique, nous obtenons après simplification :
 634 \[
 635 \mathcal{E}=E_C+E_p=\frac{\mathcal{G}\mu m_t}{2p}(\mathrm{e^2}-1)
 636 \]
 637 On en déduit l'excentricité :
 638 \[
 639 \mathrm{e^2}=\frac{2p\mathcal{E}}{\mathcal{G}\mu m_t}+1\Longrightarrow \mathrm{e}=\sqrt{\frac{2p\mathcal{E}}{\mathcal{G}\mu m_t}+1}
 640 \]
 641 On peut l'exprimer en fonction de la constante des aires en posant $K=\mathcal{G}m_t$ et en remplaçant $p$ par $\dfrac{\mathrm{C^2}}{K}$ :
 642 \begin{equation}
 643 \psframebox[linestyle=none,fillstyle=solid,fillcolor={[rgb]{1 1 0.5}}]{%
 644 \mathrm{e}=\sqrt{\frac{2\mathrm{C^2}\mathcal{E}}{\mu K^2}+1}}
 645 \label{exc}
 646 \end{equation}
 647 Il reste à déterminer $\varphi$, c'est le problème le plus simple à résoudre.
 648 \[
 649 r_0=\frac{p}{1+\mathrm{e}\cos(\theta_0-\varphi)}
 650 \]
 651 On en déduit :
 652 \begin{equation}
 653 \psframebox[linestyle=none,fillstyle=solid,fillcolor={[rgb]{1 1 0.5}}]{%
 654 \varphi=\theta_0-\arccos\left[\left(\frac{p}{r_0}-1\right)\frac{1}{\mathrm{e}}\right]}
 655 \label{phi}
 656 \end{equation}
 657 On distingue trois cas :
 658 \begin{itemize}
 659 \item
 660   $\mathrm{e}<1$ ou $\mathcal{E}<0$ : ellipse
 661 \item
 662   $\mathrm{e}=1$ ou $\mathcal{E}=0$ : parabole
 663 \item
 664  $\mathrm{e}=>$ ou $\mathcal{E}>0$ : hyperbole
 665 \end{itemize}
 666 Et on s'intéresse maintenant au mouvement elliptique. Les conditions initiales $(\overrightarrow{r}_0,\overrightarrow{v}_0)$ doivent vérifier :
 667 \[
 668 \mathcal{E}=\frac{1}{2}\mu v_0^2-\mathcal{G} \frac{m_1 m_2}{r_0}<0
 669 \]
 670 \begin{center}
 671 \begin{pspicture}(-2,-3)(5,7)
 672 \psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]
 673 \pstVerb{
 674 /arccos {
 675    dup
 676    dup mul neg 1 add sqrt
 677    exch
 678    atan
 679 } def
 680         /G 1 def
 681         /x01 2 def
 682         /y01 2 def
 683         /v0x1 .05 def
 684         /v0y1 0.25 def
 685         /x02 -2 def
 686         /y02 -0.5 def
 687         /v0x2 -0.25 def
 688         /v0y2 -0.5 def
 689         /M1 3 def
 690         /M2 1 def
 691         /Mt M1 M2 add def
 692         /mu M1 M2 mul M1 M2 add div def % masse réduite
 693         /K G Mt mul def
 694         /xG0 M1 x01 mul M2 x02 mul add M1 M2 add div def
 695         /yG0 M1 y01 mul M2 y02 mul add M1 M2 add div def
 696 % conditions initiales pour le point réduit
 697         /xr0 x01 x02 sub def
 698         /yr0 y01 y02 sub def
 699         /theta_0 yr0 xr0 atan def
 700         /r0 xr0 dup mul yr0 dup mul add sqrt def
 701         /v0xr v0x1 v0x2 sub def
 702         /v0yr v0y1 v0y2 sub def
 703         /v0r_2 v0xr dup mul v0yr dup mul add def
 704 % constante des aires
 705         /Cste xr0 v0yr mul yr0 v0xr mul sub def % Cste des aires
 706         /Energie 0.5 mu mul v0r_2 mul % 1/2 mv^2
 707                  G M1 M2 mul mul r0 div sub
 708          def
 709          /par Cste dup mul K div def % p
 710          /exc 2 Cste dup mul mul Energie mul mu div K dup mul div 1 add sqrt def % e
 711          /a_2 par 1 exc dup mul sub div def % demi-grand axe
 712          /b_2 a_2 1 exc dup mul sub sqrt mul def % demi-petit axe
 713          /c_2 a_2 exc mul def % distance focale
 714          % phase
 715          /Phi theta_0 par r0 div 1 sub exc div arccos sub def
 716          /rP par 1 exc add div def
 717          /rA par 1 exc sub div def
 718 % vitesses à l'apogée et au périgée
 719          /vA G Mt mul 2 rA div 1 a_2 div sub mul sqrt def
 720          /vP G Mt mul 2 rP div 1 a_2 div sub mul sqrt def
 721 % positions du périgée et de l'apogée
 722          /xP rP Phi cos mul def
 723          /yP rP Phi sin mul def
 724          /xA rA Phi cos mul neg def
 725          /yA rA Phi sin mul neg def
 726          /xW xA xP add 2 div def
 727          /yW yA yP add 2 div def
 728 % periode
 729        /periode 6.28 a_2 3 exp G div Mt div sqrt mul def
 730       /radius {par 1 exc t Phi sub cos mul add div} def
 731       /xE {radius t cos mul} def
 732       /yE {radius t sin mul}def
 733         }%
 734 \parametricplot[plotpoints=360]{0}{360}{xE yE}%
 735 \uput[l](0,6.8){$\mathcal{R^*}$}
 736 \uput[r](0,6.8){$y^*$}
 737 \uput[u](5,0){$x^*$}
 738 \uput[l](0,0){$C$}
 739 \psdot(0,0)
 740 \psline[style=vecteurA]{<->}(5,0)(0,0)(0,7)
 741 \pnode(!xr0 yr0){M0}\psdot(M0)
 742 \rput(M0){\psline[unit=2,style=vecteurA]{->}(!v0xr v0yr)\uput[r](0,0){$M_0$}\uput[r](!v0xr 2 mul v0yr 2 mul){$\overrightarrow{v}_0$}}
 743 \psline(M0)
 744 \pcline[offset=0.25,linestyle=none](0,0)(M0)
 745 \ncput*[nrot=:U]{$r_0$}
 746 \psarc[style=vecteurB]{->}(0,0){1.5}{0}{!theta_0}
 747 \uput{1.6}[!theta_0 2 div](0,0){$\theta_0$}
 748 \pnode(!rA Phi cos mul neg rA Phi sin mul neg){A}
 749 \psline[linestyle=dotted](A)
 750 \uput[ur](A){$A$}
 751 \pnode(!rP Phi cos mul rP Phi sin mul){P}
 752 \psline[linestyle=dotted](P)
 753 \uput[dl](P){$P$}
 754 \rput(-2,0){\psPrintValue[decimals=4]{exc}}
 755 \pnode(!xW yW){W}\psdot(W)\uput[r](W){$\Omega$}
 756 \rput{!Phi}(W){\pnode(!0 b_2){B1}\pnode(!0 b_2 neg){B2}\psline[linestyle=dotted](B1)(B2)}
 757 \uput[ul](B2){$B$}
 758 \pcline[offset=-0.25,linestyle=none](W)(B2)
 759 \ncput*{$b$}
 760 \pcline[offset=0.25,linestyle=none](W)(A)
 761 \ncput*{$a$}
 762 \pcline[offset=0.25,linestyle=none](0,0)(W)
 763 \ncput*{$c$}
 764 \psline[linecolor=blue](W)(B2)
 765 \psline[linecolor=red](W)(A)
 766 \psline[linecolor=green](W)(0,0)
 767 \rput{! 180 Phi add}(W){\psdot(!c_2 0)}
 768 \psarcn{->}(0,0){0.5}{0}{!Phi}\uput{0.55}[!Phi 2 div](0,0){$\varphi$}
 769 \rput{!Phi}(P){\pnode(! 0 vP){vP}\psline[style=vecteurA]{->}(vP)}
 770 \uput[r](vP){$\overrightarrow{v}_P$}
 771 \rput{!Phi}(A){\pnode(! 0 vA neg){vA}\psline[style=vecteurA]{->}(vA)}
 772 \uput[u](vA){$\overrightarrow{v}_A$}
 773 \end{pspicture}
 774 \end{center}
 775 L'un des foyers de l'ellipse est le centre attracteur $C$. L'autre est le symétrique de $C$ par rapport à $\Omega$.
 776
 777 Au périgée, le rayon-vecteur vaut : $r_P=\dfrac{p}{1+\mathrm{e}}$, à l'apogée $r_A=\dfrac{p}{1-\mathrm{e}}$.
 778
 779 Le demi-grand axe est égal à : $a=\dfrac{r_P+r_A}{2}=\dfrac{p}{1-\mathrm{e^2}}$.
 780
 781 La distance focale est égale à : $c=a-r_P=\dfrac{p\mathrm{e}}{1-\mathrm{e^2}}=a\mathrm{e}$.
 782
 783 Pour calculer le demi-petit axe de l'ellipse $b$, on peut rappeler que l'ellipse est lieu des points dont la somme des distances aux foyers est égale à $2a$. Au point $B$ cette somme est égale à $2a$, par conséquent le théorème de Pythagore appliqué dans le triangle rectangle $C\Omega B$ donne la valeur du demi-axe :  $b=\sqrt{a^2-c^2}=a\sqrt{1-\mathrm{e^2}}$.
 784
 785 Rappelons que $\mathrm{e^2}=\dfrac{2p\mathcal{E}}{\mathcal{G}\mu m_t}+1$. On peut exprimer $\mathcal{E}$ en fonction de $\mathrm{e}$ :
 786
 787 \[
 788 \mathcal{E}=\frac{\mathcal{G}\mu m_t}{2p}(\mathrm{e^2}-1)=-\frac{\mathcal{G}\mu m_t}{2a}\quad\text{ou}\quad\mathcal{E}=-\frac{\mathcal{G}m_1m_2}{2a}
 789 \]
 790 Pour le système $\{M_1,M_2\}$,  son énergie s'exprime très simplement en fonction du grand axe de l'ellipse $2a$.
 791
 792 Calculons la vitesse tout au long de l'ellipse :
 793
 794 \[
 795 \mathcal{E}=\frac{1}{2}\mu v^2-\frac{\mathcal{G}\mu m_t}{r}=-\frac{\mathcal{G}\mu m_t}{2a}
 796 \]
 797 \[
 798 v^2=\mathcal{G}m_t\left(\frac{2}{r}-\frac{1}{a}\right)
 799 \]
 800 On en déduit les vitesses au périgée et à l'apogée :
 801 \[
 802 v_P=\sqrt{\frac{\mathcal{G}m_t}{p}}(1+\mathrm{e})\qquad v_A=\sqrt{\frac{\mathcal{G}m_t}{p}}(1-\mathrm{e})
 803 \]
 804 La vitesse est maximale au périgée et minimale à l'apogée.
 805
 806 La période de révolution se  détermine à partir de la vitesse \textit{aérolaire}. Cette vitesse vaut :
 807 \[
 808 \mathcal{V}=\frac{\mathrm{C}}{2}=\frac{\sqrt{p\mathcal{G}m_t}}{2}
 809 \]
 810 Sachant que surface de l'ellipse est $S=\pi a b$ :
 811 \[
 812 T=\frac{2\pi a b}{\sqrt{p\mathcal{G}m_t}}
 813 \]
 814 On remplace $b$ par $a\sqrt{(1-\mathrm{e^2})}$ et $p$ par $a(1-\mathrm{e^2})$ :
 815 \[
 816 T=\frac{2\pi a^2 \sqrt{(1-\mathrm{e^2})}}{\sqrt{a(1-\mathrm{e^2})\mathcal{G}m_t}}=\frac{2\pi a^2}{\sqrt{a\mathcal{G}m_t}}
 817 \]
 818 En élevant les deux membres au carré on retrouve la troisième loi de Képler :
 819 \begin{equation}
 820 \psframebox[linestyle=none,fillstyle=solid,fillcolor={[rgb]{1 1 0.5}}]{%
 821 \frac{T^2}{a^3}=\frac{4\pi^2}{\mathcal{G}m_t}}
 822 \label{Kepler3}
 823 \end{equation}
 824 \section{Le mouvement des 2 corps dans le repère du centre de masse}
 825 Nous avons vu (\ref{r1r2}), que les trajectoires de $M_1$ et $M_2$ se déduisent de celle de $M$ par une homothétie de centre~C.
 826 \[
 827 \overrightarrow{r_1}=-\frac{m_2}{m_1+m_2}\overrightarrow{r}\quad\mbox{et}\quad \overrightarrow{r_2}=\frac{m_1}{m_1+m_2}\overrightarrow{r}
 828 \]
 829 Les trajectoires sont donc deux ellipses dont les caractéristiques, paramètre, demi-grand axe et demi-petit axe, se déduisent de l'ellipse de la particule de masse réduite dans le même rapport.
 830
 831 \begin{itemize}
 832   \item Paramètre : $p_1=p\dfrac{m_2}{m_1+m_2}$.
 833   \item rayon-vecteur à l'apogée  : $r_{A_1}=r_A\dfrac{m_2}{m_1+m_2}$
 834   \item rayon-vecteur au périgée  : $r_{P_1}=r_P\dfrac{m_2}{m_1+m_2}$
 835   \item demi-grand axe : $a_1=\dfrac{r_{A_1}+r_{P_1}}{2}=\dfrac{m_2}{m_1+m_2}\dfrac{r_{A}+r_{P}}{2}=a\dfrac{m_2}{m_1+m_2}$
 836   \item excentricité, on vérifie facilement qu'elle est inchangée : $\mathrm{e_1}=\dfrac{r_{A_1}-r_{P_1}}{r_{A_1}+r_{P_1}}=\mathrm{e}$
 837 \end{itemize}
 838 pour la deuxième :
 839 \begin{itemize}
 840   \item Paramètre : $p_2=p\dfrac{m_1}{m_1+m_2}$.
 841   \item rayon-vecteur à l'apogée  : $r_{A_2}=r_A\dfrac{m_1}{m_1+m_2}$
 842   \item rayon-vecteur au périgée  : $r_{P_2}=r_P\dfrac{m_1}{m_1+m_2}$
 843   \item demi-grand axe : $a_1=\dfrac{r_{A_2}+r_{P_2}}{2}=\dfrac{m_1}{m_1+m_2}\dfrac{r_{A}+r_{P}}{2}=a\dfrac{m_1}{m_1+m_2}$
 844   \item excentricité : $\mathrm{e_2}=\dfrac{r_{A_2}-r_{P_2}}{r_{A_2}+r_{P_2}}=\mathrm{e}$
 845 \end{itemize}
 846 Leurs équations respectives s'écrivent :
 847 \[
 848 r_1=\dfrac{p_1}{1+\mathrm{e}\cos(\theta-\varphi)} \quad r_2=\dfrac{p_2}{1+\mathrm{e}\cos(\theta-\varphi)}
 849 \]
 850 Dans l'exemple suivant $\{m_1=3,m_2=1\}$. Les vitesses initiales dans $\mathcal{R}^*$ sont indiquées par une flèche. La trajectoire de $M_1$ est en bleu, celle de $M_2$ en rouge et celle du point fictif est en pointillés.
 851 \begin{center}
 852 \begin{pspicture}(-8,-10)(3,3)
 853 \psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]
 854 \pstVerb{
 855 /arccos {
 856    dup
 857    dup mul neg 1 add sqrt
 858    exch
 859    atan
 860 } def
 861         /G 1 def
 862         /x01 2 def
 863         /y01 2 def
 864         /v0x1 .2 def
 865         /v0y1 0.25 def
 866         /x02 -3 def
 867         /y02 0 def
 868         /v0x2 -0.25 def
 869         /v0y2 -0.5 def
 870         /M1 3 def
 871         /M2 1 def
 872         /Mt M1 M2 add def
 873         /mu M1 M2 mul M1 M2 add div def % masse réduite
 874         /K G Mt mul def
 875         /xG0 M1 x01 mul M2 x02 mul add M1 M2 add div def
 876         /yG0 M1 y01 mul M2 y02 mul add M1 M2 add div def
 877         /vG0x M1 v0x1 mul M2 v0x2 mul add M1 M2 add div def
 878         /vG0y M1 v0y1 mul M2 v0y2 mul add M1 M2 add div def
 879         /K1 M2 Mt div neg def
 880         /K2 M1 Mt div def
 881 % conditions initiales pour le point réduit
 882         /xr0 x02 x01 sub def
 883         /yr0 y02 y01 sub def
 884         /theta_0 yr0 xr0 atan def
 885         /r0 xr0 dup mul yr0 dup mul add sqrt def
 886         /v0xr v0x2 v0x1 sub def
 887         /v0yr v0y2 v0y1 sub def
 888         /v0r_2 v0xr dup mul v0yr dup mul add def
 889 % constante des aires
 890         /Cste xr0 v0yr mul yr0 v0xr mul sub def % Cste des aires
 891         /Energie 0.5 mu mul v0r_2 mul % 1/2 mv^2
 892                  G M1 M2 mul mul r0 div sub
 893          def
 894          /par Cste dup mul K div def % p
 895          /exc 2 Cste dup mul mul Energie mul mu div K dup mul div 1 add sqrt def % e
 896          /a_2 par 1 exc dup mul sub div def % demi-grand axe
 897          /b_2 a_2 1 exc dup mul sub sqrt mul def % demi-petit axe
 898          /c_2 a_2 exc mul def % distance focale
 899          % phase
 900          /Phi theta_0 par r0 div 1 sub exc div arccos sub def
 901          /rP par 1 exc add div def
 902          /rA par 1 exc sub div def
 903 % vitesses à l'apogée et au périgée
 904          /vA G Mt mul 2 rA div 1 a_2 div sub mul sqrt def
 905          /vP G Mt mul 2 rP div 1 a_2 div sub mul sqrt def
 906 % positions du périgée et de l'apogée
 907          /xP rP Phi cos mul def
 908          /yP rP Phi sin mul def
 909          /xA rA Phi cos mul neg def
 910          /yA rA Phi sin mul neg def
 911          /xW xA xP add 2 div def
 912          /yW yA yP add 2 div def
 913 % periode
 914        /periode 6.28 a_2 3 exp G div Mt div sqrt mul def
 915       /radius {par 1 exc t Phi sub cos mul add div} def
 916       /xE {radius t cos mul} def
 917       /yE {radius t sin mul}def
 918         }%
 919 \parametricplot[plotpoints=360,linestyle=dotted]{0}{360}{xE yE}%
 920 \parametricplot[plotpoints=360,linecolor=blue]{0}{360}{xE K1 mul yE K1 mul}%
 921 \parametricplot[plotpoints=360,linecolor=red]{0}{360}{xE K2 mul yE K2 mul}%
 922 \uput[l](0,2.8){$\mathcal{R^*}$}
 923 \uput[r](0,2.8){$y^*$}
 924 \uput[u](3,0){$x^*$}
 925 \uput[l](0,0){$C$}
 926 \psdot(0,0)
 927 \psline[style=vecteurA]{<->}(3,0)(0,0)(0,3)
 928 \pnode(!xr0 yr0){M0}\psdot(M0)
 929 \rput(M0){\psline[unit=2,style=vecteurA]{->}(!v0xr v0yr)\uput[l](0,0){$M_0$}\uput[l](!v0xr 2 mul v0yr 2 mul){$\overrightarrow{v}_0$}}
 930 \pcline[offset=0.25,linestyle=none](M0)(0,0)
 931 \ncput*[nrot=:U]{$r_0$}
 932 \psarcn[style=vecteurB]{->}(0,0){0.75}{0}{!theta_0}
 933 \uput{0.75}[!theta_0 2 div 180 add](0,0){$\theta_0$}
 934 \pnode(!rA Phi cos mul neg rA Phi sin mul neg){A}
 935 \psline[linestyle=dotted](A)
 936 %\uput[ur](A){$A$}
 937 \pnode(!rP Phi cos mul rP Phi sin mul){P}
 938 \psline[linestyle=dotted](P)
 939 %\uput[dl](P){$P$}
 940 %\rput(-2,0){\psPrintValue[decimals=4]{exc}}
 941 \pnode(!xW yW){W}\psdot(W)\uput[r](W){$\Omega$}
 942 %\psline[linecolor=blue](W)(B2)
 943 %\psline[linecolor=red](W)(A)
 944 %\psline[linecolor=green](W)(0,0)
 945 %\rput{! 180 Phi add}(W){\psdot(!c_2 0)}
 946 %\psarc{->}(0,0){0.5}{0}{!Phi}\uput{0.55}[!Phi 2 div](0,0){$\varphi$}
 947 %\rput{!Phi}(P){\pnode(! 0 vP){vP}\psline[style=vecteurA]{->}(vP)}
 948 %\uput[r](vP){$\overrightarrow{v}_P$}
 949 %\rput{!Phi}(A){\pnode(! 0 vA neg){vA}\psline[style=vecteurA]{->}(vA)}
 950 %\uput[u](vA){$\overrightarrow{v}_A$}
 951 \pnode(!x01 xG0 sub y01 yG0 sub){M01}\psdot(M01)\uput[r](M01){$M_{0_1}$}
 952 \pnode(!x02 xG0 sub y02 yG0 sub){M02}\psdot(M02)\uput[dr](M02){$M_{0_2}$}
 953 \psline[linestyle=dashed](M0)(M02)(M01)
 954 \rput(M01){\psline[unit=2,linecolor=blue]{->}(!v0x1 vG0x sub v0y1 vG0y sub)}
 955 \rput(M02){\psline[unit=2,linecolor=red]{->}(!v0x2 vG0x sub v0y2 vG0y sub)}
 956 \end{pspicture}
 957 \end{center}
 958 \section{Exemples à développer, réalisés avec pst-eqdf}
 959 \[
 960 \mu\frac{\mathrm{d}^2\overrightarrow{r}}{\mathrm{d}t^2}=-G\frac{m_1m_2}{r^3}\overrightarrow{r}
 961 \]
 962 \[
 963 \frac{\mathrm{d}^2\overrightarrow{r}}{\mathrm{d}t^2}=-G\frac{(m_1+m_2)}{r^3}\overrightarrow{r}
 964 \]
 965
 966 \[
 967 \left\{
 968 \begin{array}[m]{l}
 969  \ddot{x}=-G\displaystyle\frac{m_1+m_2}{r^3}x\\[1em]
 970  \ddot{y}=-G\displaystyle\frac{m_1+m_2}{r^3}y
 971 \end{array}
 972 \right.
 973 \]
 974 %  x    y    x'   y'
 975 % y[0] y[1] y[2] y[3]
 976 \def\FictifAlg{%
 977   y[2]|y[3]|%
 978   -(M1+M2)*y[0]/(y[0]^2+y[1]^2)^1.5|%
 979   -(M1+M2)*y[1]/(y[0]^2+y[1]^2)^1.5}
 980 \begin{center}
 981 \begin{pspicture}(-7,-4)(7,8)
 982 \pstVerb{/G 1 def
 983         /x01 1 def
 984         /y01 2 def
 985         /v0x1 .1 def
 986         /v0y1 0.1 def
 987         /x02 -1 def
 988         /y02 -2 def
 989         /v0x2 -2 def
 990         /v0y2 0 def
 991         /M1 2 def
 992         /M2 20 def
 993         /mu M1 M2 mul M1 M2 add div def % masse réduite
 994         /xG M1 x01 mul M2 x02 mul add M1 M2 add div def
 995         /yG M1 y01 mul M2 y02 mul add M1 M2 add div def
 996         /x0 x01 x02 sub def
 997         /y0 y01 y02 sub def
 998         /v0x v0x1 v0x2 sub def
 999         /v0y v0y1 v0y2 sub def
1000         }%
1001 %%  0  1   2   3  4  5   6   7
1002 %% x1 y1 x'1 y'1 x2 y2 x'2 y'2
1003 \def\InitCondred{ x0  y0  v0x  v0y}
1004 \def\InitCond{ x01  y01  v0x1  v0y1 x02 y02  v0x2   v0y2}
1005 \psset{method=rk4}
1006 \psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=8pt]%
1007 \psequadiff[whichabs=0,whichord=1,
1008             plotpoints=1000,algebraic,
1009             tabname=X1Y1]{0}{50}{\InitCond}{\GravAlg}
1010 \listplot[unit=1,linecolor=red]{X1Y1 aload pop}
1011 \psequadiff[whichabs=4,whichord=5,
1012             plotpoints=1000,algebraic,
1013             tabname=X2Y2]{0}{50}{\InitCond}{\GravAlg}
1014 \listplot[unit=1,linecolor=blue]{X2Y2 aload pop}
1015 \psequadiff[whichabs=0,whichord=1,
1016             plotpoints=1000,algebraic,
1017             tabname=XrYr]{0}{50}{\InitCondred}{\FictifAlg}
1018 \listplot[unit=1,linecolor=green]{XrYr aload pop}
1019 \psdots(!x01 y01)(!x02 y02)
1020 \psdot[linecolor=red](!xG yG)
1021 \psline{<->}(7,0)(0,0)(0,8)
1022 \end{pspicture}
1023 \end{center}
1024 %
1025 \newpage
1026 %  x    y    x'   y'
1027 % y[0] y[1] y[2] y[3]
1028 \def\FictifAlg{%
1029   y[2]|y[3]|%
1030   -(M1+M2)*y[0]/(y[0]^2+y[1]^2)^1.5|%
1031   -(M1+M2)*y[1]/(y[0]^2+y[1]^2)^1.5}
1032 % on se place au centre de masse
1033 \def\GravAlgIIcorps{%
1034   y[2]|y[3]|%
1035   M2*(y[4]-y[0])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5|%
1036   M2*(y[5]-y[1])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5|%
1037   y[6]|y[7]|%
1038   M1*(y[0]-y[4])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5|%
1039   M1*(y[1]-y[5])/((y[4]-y[0])^2+(y[5]-y[1])^2)^1.5}
1040 %%  0  1   2   3  4  5   6   7
1041 %% x1 y1 x'1 y'1 x2 y2 x'2 y'2
1042 \newpage
1043 \begin{center}
1044 \def\InitCondred{ xr0 yr0 v0xr  v0yr}
1045 \begin{pspicture}(-7,-4)(7,8)
1046 \pstVerb{/G 1 def
1047         /x01 2 def
1048         /y01 2 def
1049         /v0x1 .05 def
1050         /v0y1 0.25 def
1051         /x02 -2 def
1052         /y02 -0.35 def
1053         /v0x2 -0.25 def
1054         /v0y2 -0.5 def
1055         /M1 3 def
1056         /M2 1 def
1057         /Mt M1 M2 add def
1058         /mu M1 M2 mul M1 M2 add div def % masse réduite
1059         /xG M1 x01 mul M2 x02 mul add M1 M2 add div def
1060         /yG M1 y01 mul M2 y02 mul add M1 M2 add div def
1061         /x0 x01 x02 sub def
1062         /y0 y01 y02 sub def
1063         /v0x v0x1 v0x2 sub def
1064         /v0y v0y1 v0y2 sub def
1065         /mu M1 M2 mul M1 M2 add div def % masse réduite
1066         /xG M1 x01 mul M2 x02 mul add M1 M2 add div def
1067         /yG M1 y01 mul M2 y02 mul add M1 M2 add div def
1068         /vxG M1 v0x1 mul M2 v0x2 mul add M1 M2 add div def
1069         /vyG M1 v0y1 mul M2 v0y2 mul add M1 M2 add div def
1070 % conditions initiales pour le point réduit
1071         /xr0 x01 x02 sub def
1072         /yr0 y01 y02 sub def
1073         /r0 xr0 dup mul yr0 dup mul add sqrt def
1074         /v0xr v0x1 v0x2 sub def
1075         /v0yr v0y1 v0y2 sub def
1076         /Lc xr0 v0yr mul yr0 v0xr mul sub mu mul def % moment cinetique
1077         /K xr0 v0yr mul yr0 v0xr mul sub def
1078         /Energie 0.5 mu mul v0xr dup mul v0yr dup mul add mul % 1/2 mv^2
1079 %                 Lc dup mul 2 div mu div r0 dup mul div
1080                  G M1 M2 mul mul r0 div sub
1081          def
1082 %        /par xr0 v0yr mul yr0 v0xr mul sub dup mul G Mt mul div def % paramètre de l'ellipse
1083 %%%%%%%%%%%%%%
1084  %      /b_2 par 1 exc dup mul sub sqrt div def % demi-petit axe
1085         /a_2 G M1 M2 add mul mu mul Energie div 2 div neg def %
1086  % excentricité
1087         /exc 1 K 2 exp a_2 G mul Mt mul div sub sqrt def
1088 % paramètre
1089 %        /par2 a_2 1 exc dup mul sub mul def
1090 % periode
1091        /periode 6.28 a_2 3 exp G div Mt div sqrt mul def
1092         }%
1093 %%  0  1   2   3  4  5   6   7
1094 %% x1 y1 x'1 y'1 x2 y2 x'2 y'2
1095 \def\InitCondred{ x0  y0  v0x  v0y}
1096 \def\InitCond{ x01  y01  v0x1  v0y1 x02 y02  v0x2   v0y2}
1097 \psset{method=rk4}
1098 \psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=8pt]%
1099 \psequadiff[whichabs=0,whichord=1,
1100             plotpoints=1000,algebraic,
1101             tabname=X1Y1]{0}{100}{\InitCond}{\GravAlgIIcorps}
1102 \listplot[unit=1,linecolor=red]{X1Y1 aload pop}
1103 \psequadiff[whichabs=4,whichord=5,
1104             plotpoints=1000,algebraic,
1105             tabname=X2Y2]{0}{100}{\InitCond}{\GravAlgIIcorps}
1106 \listplot[unit=1,linecolor=blue]{X2Y2 aload pop}
1107 % mouvement de M2 par rapport à M1
1108 \psequadiff[plotpoints=1000,algebraic,
1109             plotfuncx=y dup 4 get exch 0 get sub ,
1110             plotfuncy=dup 5 get exch 1 get sub,
1111             tabname=XY1]{0}{22.7}{\InitCond}{\GravAlgIIcorps}
1112 \listplot[unit=1,linecolor=green]{XY1 aload pop}
1113 % mouvement de M1 par rapport à M2
1114 \psequadiff[plotpoints=1000,algebraic,
1115             plotfuncx=y dup 0 get exch 4 get sub ,
1116             plotfuncy=dup 1 get exch 5 get sub,
1117             tabname=XY2]{0}{22.7}{\InitCond}{\GravAlgIIcorps}
1118 \listplot[unit=1,linecolor=magenta]{XY2 aload pop}
1119 % mouvement de M1 par rapport à G
1120 %\psequadiff[plotpoints=1000,algebraic,
1121 %            plotfuncx=y 0 get
1122 %                      y 4 get M2 mul
1123 %                      y 0 get M1 mul add
1124 %                      Mt div sub ,
1125 %            plotfuncy=1 get
1126 %                      y 1 get
1127 %                      M1 mul
1128 %                      y 5 get M2 mul add
1129 %                      Mt div sub,
1130 %            tabname=XY3]{0}{22.7}{\InitCond}{\GravAlgIIcorps}
1131 %\listplot[unit=1]{XY3 aload pop}
1132 % mouvement de M2 par rapport à G
1133 %\psequadiff[plotpoints=1000,algebraic,
1134 %            plotfuncx=y 4 get
1135 %                      y 4 get M2 mul
1136 %                      y 0 get M1 mul add
1137 %                      Mt div sub ,
1138 %            plotfuncy=5 get
1139 %                      y 1 get M1 mul
1140 %                      y 5 get M2 mul add
1141 %                     Mt div sub,
1142 %            tabname=XY4]{0}{22.7}{\InitCond}{\GravAlgIIcorps}
1143 %\listplot[unit=1,linecolor=red]{XY4 aload pop}
1144 % centre de masse
1145 \psequadiff[plotpoints=1000,algebraic,
1146             plotfuncx=y dup 0 get M1 mul exch
1147                       4 get M2 mul add
1148                       Mt div,
1149             plotfuncy=dup 1 get M1 mul exch
1150                       5 get M2 mul add
1151                       Mt div,
1152             tabname=XGYG]{0}{100}{\InitCond}{\GravAlgIIcorps}
1153 \listplot[unit=1,linecolor=cyan]{XGYG aload pop}
1154 \psdots(!x01 y01)(!x02 y02)
1155 \psdot[linecolor=red](!xG yG)
1156 \psline{<->}(7,0)(0,0)(0,8)
1157 \rput(-2,0){\psPrintValue[decimals=4]{periode}\hphantom{00000}s}
1158 \end{pspicture}
1159 \end{center}
1160 %
1161 \newpage
1162 \begin{center}
1163 \def\InitCondred{ xr0 yr0 v0xr  v0yr}
1164 \begin{pspicture}(-7,-6)(7,8)
1165 \pstVerb{
1166 /arccos {
1167    dup
1168    dup mul neg 1 add sqrt
1169    exch
1170    atan
1171 } def
1172         /G 1 def
1173         /x01 2 def
1174         /y01 2 def
1175         /v0x1 .05 def
1176         /v0y1 0.25 def
1177         /x02 -2 def
1178         /y02 -0.5 def
1179         /v0x2 -0.25 def
1180         /v0y2 -0.5 def
1181         /M1 3 def
1182         /M2 1 def
1183         /Mt M1 M2 add def
1184         /mu M1 M2 mul M1 M2 add div def % masse réduite
1185         /xG0 M1 x01 mul M2 x02 mul add M1 M2 add div def
1186         /yG0 M1 y01 mul M2 y02 mul add M1 M2 add div def
1187         /x0 x01 x02 sub def
1188         /y0 y01 y02 sub def
1189         /theta_0 y0 x0 atan def
1190         /v0x v0x1 v0x2 sub def
1191         /v0y v0y1 v0y2 sub def
1192         /mu M1 M2 mul M1 M2 add div def % masse réduite
1193         /xG M1 x01 mul M2 x02 mul add M1 M2 add div def
1194         /yG M1 y01 mul M2 y02 mul add M1 M2 add div def
1195         /vxG M1 v0x1 mul M2 v0x2 mul add M1 M2 add div def
1196         /vyG M1 v0y1 mul M2 v0y2 mul add M1 M2 add div def
1197 % conditions initiales pour le point réduit
1198         /xr0 x01 x02 sub def
1199         /yr0 y01 y02 sub def
1200         /r0 xr0 dup mul yr0 dup mul add sqrt def
1201         /v0xr v0x1 v0x2 sub def
1202         /v0yr v0y1 v0y2 sub def
1203         /v0r_2 v0xr dup mul v0yr dup mul add def
1204 %         v0r 2 Mt mul r0 div G mul ge { a faire} if
1205         /Lc xr0 v0yr mul yr0 v0xr mul sub mu mul def % moment cinetique
1206         /K xr0 v0yr mul yr0 v0xr mul sub def
1207         /Energie 0.5 mu mul v0xr dup mul v0yr dup mul add mul % 1/2 mv^2
1208 %                 Lc dup mul 2 div mu div r0 dup mul div
1209                  G M1 M2 mul mul r0 div sub
1210          def
1211 %        /par xr0 v0yr mul yr0 v0xr mul sub dup mul G Mt mul div def % paramètre de l'ellipse
1212 %%%%%%%%%%%%%%
1213  %      /b_2 par 1 exc dup mul sub sqrt div def % demi-petit axe
1214         /a_2 G M1 M2 add mul mu mul Energie div 2 div neg def %
1215  % excentricité
1216         /exc 1 K 2 exp a_2 G mul Mt mul div sub sqrt def
1217 % paramètre
1218         /par2 a_2 1 exc dup mul sub mul def
1219 % periode
1220        /periode 6.28 a_2 3 exp G div Mt div sqrt mul def
1221 % phase
1222       /Phi theta_0 par2 r0 div 1 sub exc div arccos sub def
1223         }%
1224 %%  0  1   2   3  4  5   6   7
1225 %% x1 y1 x'1 y'1 x2 y2 x'2 y'2
1226 \def\InitCondred{ x0  y0  v0x  v0y}
1227 \def\InitCond{ x01  y01  v0x1  v0y1 x02 y02  v0x2   v0y2}
1228 \psset{method=rk4}
1229 %\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=8pt]%
1230 % mouvement de M1 par rapport à G
1231 \psequadiff[plotpoints=1000,algebraic,
1232             plotfuncx=y 0 get
1233                       y 4 get M2 mul
1234                       y 0 get M1 mul add
1235                       Mt div sub ,
1236             plotfuncy=1 get
1237                       y 1 get
1238                       M1 mul
1239                       y 5 get M2 mul add
1240                       Mt div sub,
1241             tabname=XY3]{0}{23.6}{\InitCond}{\GravAlgIIcorps}
1242 \listplot[unit=1,linecolor=red]{XY3 aload pop}
1243 % mouvement de M2 par rapport à G
1244 \psequadiff[plotpoints=1000,algebraic,
1245             plotfuncx=y 4 get
1246                       y 4 get M2 mul
1247                       y 0 get M1 mul add
1248                       Mt div sub ,
1249             plotfuncy=5 get
1250                       y 1 get M1 mul
1251                       y 5 get M2 mul add
1252                      Mt div sub,
1253             tabname=XY4]{0}{23.6}{\InitCond}{\GravAlgIIcorps}
1254 \listplot[unit=1,linecolor=blue]{XY4 aload pop}
1255 \psline{<->}(7,0)(0,0)(0,8)
1256 \rput(-2,2){\psPrintValue[decimals=4]{periode}\hphantom{0000000}s}
1257 %\parametricplot[plotpoints=360,linestyle=dotted]{0}{360}{%
1258 %    /radius par2 1 exc t theta_0 sub cos mul add div def
1259 %     radius t cos mul
1260 %     radius t sin mul}%
1261 \psdots(!x01 xG0 sub y01 yG0 sub)(!x02 xG0 sub y02 yG0 sub)
1262 \psline(!x01 xG0 sub y01 yG0 sub)(!x02 xG0 sub y02 yG0 sub)
1263 % mouvement de M (masse réduite) par rapport à G
1264 \psequadiff[plotpoints=1000,algebraic,
1265             plotfuncx=y 4 get
1266                       y 4 get M2 mul
1267                       y 0 get M1 mul add
1268                       Mt div sub
1269                       Mt mul M1 div,
1270             plotfuncy=5 get
1271                       y 1 get M1 mul
1272                       y 5 get M2 mul add
1273                       Mt div sub
1274                       Mt mul M1 div,
1275             tabname=XYM]{0}{23.6}{\InitCond}{\GravAlgIIcorps}%
1276 \listplot[unit=1,linecolor=gray!80,linewidth=0.1]{XYM aload pop}
1277 \rput(!x01 xG0 sub y01 yG0 sub){\psline[unit=5,linecolor=red!50]{->}(!v0x1 vxG sub v0y1 vyG sub)}
1278 \rput(!x02 xG0 sub y02 yG0 sub){\psline[unit=5,linecolor=blue!50]{->}(!v0x2 vxG sub v0y2 vyG sub)}
1279 % point "réduit"
1280 \pstVerb{/xr0 x01 xG0 sub Mt mul M2 div neg def
1281          /yr0 y01 yG0 sub Mt mul M2 div neg def}%
1282 \psdot[linecolor=gray](!xr0 yr0)
1283 \psline(!x01 xG0 sub y01 yG0 sub)(!x02 xG0 sub y02 yG0 sub)(!xr0 yr0)
1284 \uput[d](0,0){$C$}
1285 \uput[r](!x01 xG0 sub y01 yG0 sub){$M_1$}
1286 \uput[dr](!x02 xG0 sub y02 yG0 sub){$M_2$}
1287 \uput[l](!xr0 yr0){$M$}
1288 \rput(!xG0 neg yG0 neg){\psline{<->}(7,0)(0,0)(0,8)\psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=8pt]}
1289 \parametricplot[plotpoints=360,linecolor=white]{0}{360}{
1290     /radius par2 1 exc t Phi sub cos mul add div def
1291      radius t cos mul neg
1292      radius t sin mul neg}%
1293 \end{pspicture}
1294 \end{center}
1295 \end{document}