gravitation/gravitation_01.tex

   1 \documentclass{article}
   2 \usepackage[a4paper,margin=2cm]{geometry}
   3 \usepackage[latin1]{inputenc}
   4 \usepackage{pstricks-add,pst-eqdf,pst-func}
   5 \usepackage{array,amsmath}
   6 \usepackage{listings}
   7 \newpsstyle{vecteurA}{arrowinset=0.05,arrowsize=0.125,linecolor={[rgb]{1 0.5 0}}}
   8 \newpsstyle{vecteurB}{arrowinset=0.05,arrowsize=0.1,linecolor={[rgb]{0 0.5 1}}}
   9 \newpsstyle{vecteurC}{arrowinset=0.1,arrowsize=0.2,linecolor={[rgb]{1 0.5 0}}}
  10 %%%%%%%%%%%%%%%%%%
  11 \title{Interaction gravitationnelle avec PSTricks}
  12 \date{19 juin 2\,012}
  13 \begin{document}
  14 \maketitle
  15 \def\eqsatellite{%
  16 y[2]|y[3]|-GM*y[0]/((sqrt(y[0]^2+y[1]^2))^3)|-GM*y[1]/((sqrt(y[0]^2+y[1]^2))^3)}%
  17 \section{Mise en orbite d'un satellite}
  18 \begin{figure}[htbp]
  19   \begin{center}
  20 \psset{unit=2}
  21 \begin{pspicture}(-6,-2)(2,5)
  22 \pstVerb{
  23     /GM 1 def
  24     /theta0 -35 def
  25     /r0 1 def
  26     /x0 r0 theta0 cos mul def
  27     /y0 r0 theta0 sin mul def
  28     /v0 1.3 def
  29     /v0x v0 theta0 sin mul neg def
  30     /v0y v0 theta0 cos mul def
  31     /Lc r0 v0 mul def % moment cinetique
  32     /par Lc dup mul GM div def % paramètre de l'ellipse
  33 % excentricité
  34     /exc 1 0.5 v0 4 exp mul r0 dup mul mul GM r0 mul v0 dup mul mul sub GM dup mul div 2 mul add sqrt def
  35 %%%%%%%%%%%%%%
  36     /a_2 par 1 exc dup mul sub div def % demi-grand axe
  37     /b_2 par 1 exc dup mul sub sqrt div def % demi-petit axe
  38     /periode 2 3.1416 mul a_2 3 exp GM div sqrt mul def
  39 % vitesse à l'apogée
  40     /vA GM par div sqrt 1 exc sub mul def
  41 % vitesse au périgée
  42     /vP GM par div sqrt 1 exc add mul def
  43 % coordonnées de vA
  44     /vAx vA theta0 90 add cos mul neg def
  45     /vAy vA theta0 90 add sin mul neg def
  46 % coordonnées de vP
  47     /vPx vP theta0 90 add cos mul def
  48     /vPy vP theta0 90 add sin mul def
  49 % Apogée
  50     /xA par 1 exc sub div theta0 cos mul neg def
  51     /yA par 1 exc sub div theta0 sin mul neg def
  52 % Périgée
  53    /xP par 1 exc add div theta0 cos mul def
  54    /yP par 1 exc add div theta0 sin mul def
  55 % Centre
  56    /xO xP xA add 2 div def
  57    /yO yP yA add 2 div def
  58 }%
  59 \pscircle[fillcolor=gray!70,fillstyle=solid](0,0){0.3}
  60 \psdot[dotstyle=+](0,0)
  61 \uput[d](0,0){O}
  62 \psgrid[subgriddiv=2,gridcolor=lightgray,gridlabels=8pt](-6,-2)(2,5)
  63 \rput(1,4){\psframebox[linestyle=none,fillstyle=solid,fillcolor=white]{T=\psPrintValue[decimals=2]{periode}\hphantom{0000}s}}
  64 \parametricplot[linecolor=red,plotpoints=360]{0}{360}{/radius par 1 exc t theta0 sub cos mul add div def
  65  radius t cos mul
  66  radius t sin mul}
  67 \pscircle*(!x0 y0){0.05}
  68 \pnode(!xP yP){P} % périgée
  69 \pnode(!xA yA){A} % Apogée
  70 \pnode(!xO yO){O} % centre
  71 \rput(!x0 y0){\psline[style=vecteurB]{->}(!v0x v0y)\uput[ur](!v0x v0y){$\overrightarrow{v_0}$}}
  72 \psline[linestyle=dashed](A)(P)
  73 \pscircle*(A){0.05}
  74 \uput[ul](A){$A$}
  75 \uput[dr](P){$P$}
  76 % position du satellite à un instant quelconque
  77 \pstVerb{/theta_i 170 def
  78  /radius par 1 exc theta_i theta0 sub cos mul add div def
  79  /xS radius theta_i cos mul def
  80  /yS radius theta_i sin mul def
  81  /ux theta_i cos 1 mul def
  82  /uy theta_i sin 1 mul def
  83  /xi xS ux sub def
  84  /yi yS uy sub def
  85  /xi2 xS ux 2 div sub def
  86  /yi2 yS uy 2 div sub def}%
  87 \pnode(!xi2 yi2){Mi2}
  88 \pnode(!xi yi){Mi}
  89 \pnode(!xS yS){S}
  90 \pscircle*(S){0.05}
  91 \psline[linestyle=dotted](S)(0,0)
  92 \psline[style=vecteurA]{->}(S)(Mi)
  93 \uput[l](S){S}
  94 \uput[u](Mi2){$\overrightarrow{F}$}
  95 \psarcn{->}(0,0){0.4}{0}{!theta0}
  96 \uput{0.5}[!theta0 2 div](0,0){$\theta_0$}
  97 \rput(A){\psline[style=vecteurB]{->}(!vAx vAy)}
  98 %\rput(P){\psline[style=vecteurB]{->}(!vPx vPy)}
  99 \psline[arrowinset=0.05,arrowsize=0.1]{<->}(2,0)(0,0)(0,5)
 100 \uput[r](0,4.9){$y$}\uput[u](1.9,0){$x$}
 101 \psdot(O)\uput[r](O){$\Omega$}
 102 \rput{!90 theta0 add}(O){\psline[linestyle=dashed](!b_2 neg 0)(!b_2 0)}
 103 \end{pspicture}
 104 \caption{Mouvement d'un satellite}
 105   \end{center}
 106 \end{figure}
 107 Soit (M) la masse de l'astre et (m) celle du satellite avec $m\ll M$. Le centre de masse du système \{M,m\} est confondu avec le centre de l'astre attracteur. La mise en orbite s'effectue à partir du point $M_0(x_0,y_0)$ avec une vitesse $\overrightarrow{v_0}(v_{0_x},v_{0_x})$. $\theta_0$ est l'angle que fait $\overrightarrow{Ox}$ avec $\overrightarrow{OM_0}$. Le satellite (S), supposé ponctuel, subit de la part de l'astre une force d'attraction gravitationnelle :
 108 \[
 109 \overrightarrow{F}=-\mathcal{G}\frac{Mm}{r^2}\overrightarrow{u}\qquad \text{avec}\quad \overrightarrow{u}=\frac{\overrightarrow{r}}{r}\quad \text{et}\quad \overrightarrow{r}=\overrightarrow{OS}
 110 \]
 111 Tous les \textit{bons} livres de mécanique\footnote{Comme celui, par exemple, de José-Philippe Pérez, aux éditions Masson.} établissent les relations suivantes :
 112 \[
 113 r=\frac{p}{1+\mathrm{e}\cos(\theta-\theta_0)}
 114 \]
 115 Paramètres et excentricité ont pour expressions respectives, avec les notations suivantes : $K=\mathcal{G}Mm$, $\mathcal{E}$ l'énergie du système et $L$ le moment cinétique.
 116 \[
 117 \mathrm{e}=\sqrt{1+\frac{2\mathcal{E}L^2}{mK^2}}\qquad p=\frac{L^2}{mK}
 118 \]
 119 On choisit une vitesse initiale $\overrightarrow{v_0}$ perpendiculaire à $\overrightarrow{OM_0}$, dans ces conditions le moment cinétique et l'énergie, qui restent constants, valent :
 120 \[
 121 L=mr_0v_0\qquad \mathcal{E}=-\frac{K}{r_0}+\frac{1}{2}mv_0^2
 122 \]
 123 En remplaçant $L$ et $\mathcal{E}$, on obtient pour l'excentricité et le paramètre les expressions suivantes :
 124 \[
 125 \mathrm{e}=\sqrt{1+\frac{1}{\mathcal{G}^2M^2}\Big(\frac{1}{2}v_0^4r_0^2-\mathcal{G}Mr_0v_0^2\Big)}
 126 \]
 127 \[
 128 p=\frac{v_0^2r_0^2}{\mathcal{G}M}
 129 \]
 130 On se limite au cas du mouvement elliptique, avec, en conséquence, la condition :
 131 \[
 132 \mathcal{E}=-\frac{K}{r_0}+\frac{1}{2}mv_0^2 < 0
 133 \]
 134 Le demi-grand axe $a$, le demi-petit axe $b$ sont :
 135 \[
 136 a=\frac{p}{1-\mathrm{e}^2}\qquad b=\frac{p}{\sqrt{1-\mathrm{e}^2}}
 137 \]
 138 La période $T$ qui obéit à la troisième loi de Képler :
 139 \[
 140 T^2=\frac{4\pi^2a^3}{\mathcal{G}M}
 141 \]
 142 La vitesse, en un point de l'ellipse, se calcule par :
 143 \[
 144 v^2=\mathcal{G}M\Big( \frac{2}{r}-\frac{1}{a}\Big)
 145 \]
 146 Sachant que $r_p=\dfrac{p}{1+\mathrm{e}}$ et $r_A=\dfrac{p}{1-\mathrm{e}}$, on en déduit les vitesses au périgée et à l'apogée :
 147 \[
 148 v_P=\sqrt{\frac{\mathcal{G}M}{p}}(1+\mathrm{e})\qquad v_A=\sqrt{\frac{\mathcal{G}M}{p}}(1-\mathrm{e})
 149 \]
 150
 151 \section{L'étude avec PSTricks}
 152 \subsection{La trajectoire}
 153 \def\parametres{
 154     /GM 1 def % 4e14 def % GxM
 155     /x0 6.5e6 def % position initiale
 156     /y0 0 def
 157     /vx0 0 def    % vitesse initiale
 158     /vy0 1e4 def
 159 }
 160 %  x0   y0   x'0   y'0
 161 % y[0] y[1] y[2]  y[3]
 162 \[
 163 \left\{
 164 \begin{array}{rcl}
 165 \ddot{x}&=&-\dfrac{GM}{r^3}x\\[1em]
 166 \ddot{y}&=&-\dfrac{GM}{r^3}y\\
 167 \end{array}
 168 \right.
 169 \label{eq1}
 170 \qquad r=\sqrt{x^2+y^2}
 171 \]
 172
 173 On peut dessiner la trajectoire du satellite et de ses caractéristiques de deux façons :
 174 \begin{itemize}
 175   \item par l'utilisation de \verb+\parametricplot+ ;
 176   \item ou celle de \verb+\psplotDiffEqn+.
 177 \end{itemize}
 178 \verb+\parametricplot+ utilise l'expression exacte de l'équation de la trajectoire en coordonnées polaires :
 179 \begin{verbatim}
 180    \parametricplot[linecolor=red,unit=2,plotpoints=360]{0}{360}{%
 181      /radius par 1 exc t theta0 sub cos mul add div def
 182       radius t cos mul
 183       radius t sin mul}
 184 \end{verbatim}
 185 L'excentricité, la période, demi-grand axe et demi-petit axe sont calculés par quelques lignes de code \textsf{postscript}. Il faut s'assurer que les conditions initiales choisies vérifient bien la condition d'une trajectoire elliptique, pour cela il faut que l'énergie initiale $\mathcal{E}_0<0$, sinon cela entraînera une erreur lors du passage à l'interpréteur \textsf{postscript}.
 186 \begin{verbatim}
 187 \pstVerb{
 188     /GM 1 def
 189     /theta0 -45 def
 190     /r0 0.5 def
 191     /x0 r0 theta0 cos mul def
 192     /y0 r0 theta0 sin mul def
 193     /v0 1.92 def
 194     /v0x v0 theta0 sin mul neg def
 195     /v0y v0 theta0 cos mul def
 196     /Lc r0 v0 mul def % moment cinetique
 197     /par Lc dup mul GM div def % paramètre de l'ellipse
 198 % excentricité
 199     /exc 1 0.5 v0 4 exp mul r0 dup mul mul GM r0 mul
 200          v0 dup mul mul sub GM dup mul div 2 mul add sqrt def
 201 %demi-grand axe
 202     /a_2 par 1 exc dup mul sub div def % demi-grand axe
 203 %demi-petit axe
 204     /b_2 par 1 exc dup mul sub sqrt div def % demi-petit axe
 205 % période
 206     /periode 2 3.1416 dup mul a_2 3 exp mul GM div sqrt mul def
 207 }%
 208 \end{verbatim}
 209 \verb+\psplotDiffEqn+ utilise les équations différentielles du mouvement, en notation algébrique :
 210 \begin{verbatim}
 211 %  x0   y0   x'0   y'0
 212 % y[0] y[1] y[2]  y[3]
 213 \def\eqsatellite{%
 214 y[2]|y[3]|-GM*y[0]/((sqrt(y[0]^2+y[1]^2))^3)|-GM*y[1]/((sqrt(y[0]^2+y[1]^2))^3)}
 215  \psplotDiffEqn[unit=2,whichabs=0,whichord=1,%
 216                 linecolor=blue,linewidth=0.1,%
 217                 method=rk4,plotpoints=1000,%
 218                 algebraic]{0}{37.8}{x0 y0 v0x v0y}{\eqsatellite}%
 219 \end{verbatim}
 220 Ce qui permet, par ailleurs de vérifier la qualité du tracé par la méthode numérique, en bleu, tandis que le tracé à partir de l'expression exacte est en trait fin en rouge.
 221 %
 222 \begin{center}
 223 \begin{pspicture}(-12,-2)(4,10)
 224 \pstVerb{
 225     /GM 1 def
 226     /theta0 -45 def
 227     /r0 0.5 def
 228     /x0 r0 theta0 cos mul def
 229     /y0 r0 theta0 sin mul def
 230     /v0 1.92 def
 231     /v0x v0 theta0 sin mul neg def
 232     /v0y v0 theta0 cos mul def
 233     /Lc r0 v0 mul def % moment cinetique
 234     /par Lc dup mul GM div def % paramètre de l'ellipse
 235 % excentricité
 236     /exc 1 0.5 v0 4 exp mul r0 dup mul mul GM r0 mul v0 dup mul mul sub GM dup mul div 2 mul add sqrt def
 237 %%%%%%%%%%%%%%
 238     /a_2 par 1 exc dup mul sub div def % demi-grand axe
 239     /b_2 par 1 exc dup mul sub sqrt div def % demi-petit axe
 240     /periode 2 3.1416 dup mul a_2 3 exp mul GM div sqrt mul def
 241 }%
 242 \psframe*[linecolor=white](-3,-0.2)(0,0.2)
 243 \pscircle[fillcolor=blue!50,fillstyle=solid](0,0){0.5}
 244 \psgrid[unit=2,subgriddiv=2,gridcolor=lightgray,gridlabels=8pt](-6,-1)(2,5)
 245 \rput(-2,0){T=\psPrintValue[decimals=2]{periode}\hphantom{00000}s}
 246   \psplotDiffEqn[unit=2,whichabs=0,whichord=1,linecolor=blue,linewidth=0.1,method=rk4,plotpoints=1000,algebraic]{0}{37.8}{x0 y0 v0x v0y}{\eqsatellite}%
 247 \parametricplot[linecolor=red,unit=2,plotpoints=360]{0}{360}{/radius par 1 exc t theta0 sub cos mul add div def
 248  radius t cos mul
 249  radius t sin mul}
 250 \psdot[unit=2,dotsize=0.12](!x0 y0)
 251 \rput(!x0 2 mul y0 2 mul){\psline[style=vecteurC]{->}(!v0x v0y)}
 252 \end{pspicture}
 253 \end{center}
 254 \subsection{La vitesse}
 255 \verb+\psplotDiffEqn+ permet de voir comment varie la vitesse sur l'ellipse :
 256 \begin{verbatim}
 257 %  x0   y0   x'0   y'0
 258 % y[0] y[1] y[2]  y[3]
 259  \psplotDiffEqn[xunit=0.2,yunit=5,%
 260                 plotfuncy=dup 2 get dup mul exch 3 get dup mul add sqrt,
 261                 linecolor=red,method=rk4,plotpoints=1000,
 262                 algebraic]{0}{50}{x0 y0 v0x v0y}{\eqsatellite}%
 263 \end{verbatim}
 264 \begin{center}
 265 \begin{pspicture}(0,-1)(10,7)
 266 \psgrid[subgriddiv=0,gridcolor=lightgray,griddots=10,gridlabels=0pt]
 267 \pstVerb{/GM 1 def
 268     /theta0 -35 def
 269     /r0 1 def
 270     /x0 r0 theta0 cos mul def
 271     /y0 r0 theta0 sin mul def
 272     /v0 1.3 def
 273     /v0x v0 theta0 sin mul neg def
 274     /v0y v0 theta0 cos mul def}%
 275 \psplotDiffEqn[xunit=0.2,yunit=5,
 276   plotfuncy=dup 2 get dup mul exch 3 get dup mul add sqrt,
 277   linecolor=red,method=rk4,plotpoints=1000,algebraic]{0}{50}{x0 y0 v0x v0y}{\eqsatellite}%
 278  \multido{\i=1+1,\I=5+5}{9}{\uput[u](\i,0){\I}}
 279 \pnode(! 36.4 5 div 0){P}
 280 \psdot(P)\uput[d](P){Périgée}
 281 \psline[linestyle=dashed](P)(! 36.4 5 div 7)
 282 \pnode(! 36.4 10 div 0){A}
 283 \psdot(A)\uput[d](A){Apogée}
 284 \psline[linestyle=dashed](A)(! 36.4 10 div 7)
 285 \uput[l](0,6.5){$v$}
 286 \psline[arrowinset=0.1,arrowsize=0.2]{<->}(10,0)(0,0)(0,7)
 287 \uput[u](10,0){$t$(s)}
 288 \uput[ul](0,0){0}
 289 \end{pspicture}
 290 \end{center}
 291 On peut obtenir les caractéristiques de la vitesse en un point quelconque, car le moment cinétique $L=mr\dot{\theta}$ étant constant, pour chaque valeur de $\theta$, on en déduit $r$ puis $\dot{\theta}$. En coordonnées polaires, la vitesse s'exprime par :
 292 \[
 293 \overrightarrow{v}=\dot{r}\overrightarrow{u_r}+r\dot{\theta}\overrightarrow{u_{\theta}}
 294 \]
 295 $\dot{\theta}$ et $\dot{r}$ s'obtiennent par les relations suivantes :
 296 \[
 297 \dot{\theta}=\frac{r_0v_0}{r^2}
 298 \]
 299 \[
 300 \dot{r}=-\frac{p(-\dot{\theta}\sin(\theta-\theta_0)}{(1+\mathrm{e}\cos(\theta-\theta_0))^2}=\frac{r_0v_0}{p}\sin(\theta-\theta_0)
 301 \]
 302 La chaîne de calculs est la suivante : $\theta\Longrightarrow r\Longrightarrow \dot{\theta}\Longrightarrow \dot{r}\Longrightarrow \overrightarrow{v}$.
 303
 304 Le package `\textsf{pst-eqdf}' comprend la commande \verb+\psequadiff+ qui est une version simplifiée de \verb+\psplotDiffEqn+, dont elle ne reprend que la méthode Runge-Kutta~4. Elle permet de sauvegarder sous forme de tableaux et/ou de fichiers toutes les variables et les dérivées de la fonction étudiée, cette possibilité est intéressante pour déterminer les caractéristiques de la vitesse au cours du temps et elle est particulièrement utile pour créer une animation. Nous allons l'utiliser pour dessiner le vecteur-vitesse à quelques instants.
 305
 306 On sauve successivement le tableau des positions et celui des vitesses.
 307 \begin{verbatim}
 308 \psequadiff[method=rk4,plotpoints=1000,
 309             algebraic,
 310             whichabs=0,whichord=1,
 311             tabname=XiYi]{0}{43}{x0 y0 v0x v0y}{\eqsatellite}%
 312 \end{verbatim}
 313 \begin{verbatim}
 314 \psequadiff[method=rk4,plotpoints=1000,
 315             algebraic,
 316             whichabs=2,whichord=3,
 317             tabname=vxvy]{0}{43}{x0 y0 v0x v0y}{\eqsatellite}%
 318 \end{verbatim}
 319
 320 Pour ensuite dessiner la trajectoire et les vecteurs-vitesse.
 321 \begin{verbatim}
 322 %\listplot[unit=1]{vxvy aload pop}
 323 % on dessine la vitesse un point sur 100
 324 \pscircle[fillcolor=gray!70,fillstyle=solid](0,0){0.3}
 325 \multido{\i=0+100}{20}{%
 326 \pstVerb{/vX vxvy \i\space get def
 327          /vY vxvy \i\space 1 add get def
 328          /xi XiYi \i\space get def
 329          /yi XiYi \i\space 1 add get def}%
 330 \rput(!xi yi){\psline[style=vecteurA]{->}(! vX 2 mul vY 2 mul)}}
 331 \end{verbatim}
 332 \begin{center}
 333 \begin{pspicture}(-10,-10)(6,7)
 334 \psset{unit=2}%
 335 \pstVerb{/GM 1 def
 336     /theta0 30 def
 337     /r0 2 def
 338     /x0 r0 theta0 cos mul def
 339     /y0 r0 theta0 sin mul def
 340     /v0 0.85 def
 341     /v0x v0 theta0 sin mul neg def
 342     /v0y v0 theta0 cos mul def
 343     /Lc r0 v0 mul def % moment cinetique
 344     /par Lc dup mul GM div def % paramètre de l'ellipse
 345 % excentricité
 346     /exc 1 0.5 v0 4 exp mul r0 dup mul mul GM r0 mul v0 dup mul mul sub GM dup mul div 2 mul add sqrt def
 347 %%%%%%%%%%%%%%
 348     /a_2 par 1 exc dup mul sub div def % demi-grand axe
 349     /b_2 par 1 exc dup mul sub sqrt div def % demi-petit axe
 350     /periode 2 3.1416 dup mul a_2 3 exp mul GM div sqrt mul def}%
 351 \rput(-2,0){T=\psPrintValue[decimals=2]{periode}\hphantom{00000}s}
 352 \psequadiff[method=rk4,
 353             plotpoints=1000,
 354             algebraic,
 355             whichabs=0,
 356             whichord=1,
 357             tabname=XiYi
 358 %           ,saveData,filename=XiYi.dat
 359 ]{0}{43}{x0 y0 v0x v0y}{\eqsatellite}%
 360 \listplot{XiYi aload pop}
 361 \psequadiff[method=rk4,
 362             plotpoints=1000,
 363             algebraic,
 364             whichabs=2,
 365             whichord=3,
 366             tabname=vxvy
 367 %           ,saveData,filename=vxvy.dat
 368 ]{0}{43}{x0 y0 v0x v0y}{\eqsatellite}%
 369 %\listplot[unit=1]{vxvy aload pop}
 370 % on dessine la vitesse un point sur 100
 371 \pscircle[fillcolor=gray!70,fillstyle=solid](0,0){0.3}
 372 \multido{\i=0+100}{20}{%
 373 \pstVerb{/vX vxvy \i\space get def
 374          /vY vxvy \i\space 1 add get def
 375          /xi XiYi \i\space get def
 376          /yi XiYi \i\space 1 add get def}%
 377 \rput(!xi yi){\psline[style=vecteurA]{->}(! vX 2 mul vY 2 mul)}}
 378 \psgrid[subgriddiv=2,gridcolor=lightgray,gridlabels=8pt](-5,-5)(3,3)
 379 \end{pspicture}
 380 \end{center}
 381 %
 382 \section{Mouvement circulaire}
 383  \begin{center}
 384 \psset{unit=1}
 385 \begin{pspicture}(-5,-5.5)(5,5.5)
 386 \pstVerb{
 387     /GM 1 def
 388     /theta0 60 def
 389     /r0 5 def
 390     /x0 r0 theta0 cos mul def
 391     /y0 r0 theta0 sin mul def
 392     /v0 GM r0 div sqrt def
 393     /v0x v0 theta0 sin mul neg def
 394     /v0y v0 theta0 cos mul def
 395     /Lc r0 v0 mul def % moment cinetique
 396     /par Lc dup mul GM div def % paramètre de l'ellipse
 397 % excentricité
 398     /exc 1 0.5 v0 4 exp mul r0 dup mul mul GM r0 mul v0 dup mul mul sub GM dup mul div 2 mul add sqrt def
 399 %%%%%%%%%%%%%%
 400     /a_2 par 1 exc dup mul sub div def % demi-grand axe
 401     /b_2 par 1 exc dup mul sub sqrt div def % demi-petit axe
 402     /periode 2 3.1416 mul a_2 3 exp GM div sqrt mul def
 403 % vitesse à l'apogée
 404     /vA GM par div sqrt 1 exc sub mul def
 405 % vitesse au périgée
 406     /vP GM par div sqrt 1 exc add mul def
 407 % coordonnées de vA
 408     /vAx vA theta0 90 add cos mul neg def
 409     /vAy vA theta0 90 add sin mul neg def
 410 % coordonnées de vP
 411     /vPx vP theta0 90 add cos mul def
 412     /vPy vP theta0 90 add sin mul def
 413 }%
 414 \pscircle[fillcolor=gray!70,fillstyle=solid](0,0){0.75}
 415 \psdot[dotstyle=+](0,0)
 416 \uput[d](0,0){O}
 417 \psgrid[subgriddiv=2,gridcolor=lightgray,gridlabels=8pt](-5,-5)(5,5)
 418 \rput(0,-2){\psframebox[linestyle=none,fillstyle=solid,fillcolor=white]{T=\psPrintValue[decimals=2]{periode}\hphantom{000000}s}}
 419 \parametricplot[linecolor=red,plotpoints=360]{0}{360}{/radius par 1 exc t theta0 sub cos mul add div def
 420  radius t cos mul
 421  radius t sin mul}
 422 \pscircle*(!x0 y0){0.1}
 423 %\pnode(!par 1 exc add div theta0 cos mul par 1 exc add div theta0 sin mul){P} % périgée
 424 %\pnode(!par 1 exc sub div theta0 cos mul neg par 1 exc sub div theta0 sin mul neg){A} % Apogée
 425 \rput(!x0 y0){\psline[arrowinset=0.1,arrowsize=0.2,linecolor={[rgb]{0 0.5 1}},unit=4]{->}(!v0x v0y)\uput[ur](!v0x 2 mul v0y 2 mul){$\overrightarrow{v_0}$}}
 426 \uput[ur](!x0 y0){$M_0$}
 427 \psline[linestyle=dashed](0,0)(!x0 y0)
 428 % position du satellite à un instant quelconque
 429 \pstVerb{/theta_i 170 def
 430  /radius par 1 exc theta_i theta0 sub cos mul add div def
 431  /xS radius theta_i cos mul def
 432  /yS radius theta_i sin mul def
 433  /ux theta_i cos 1 mul def
 434  /uy theta_i sin 1 mul def
 435  /xi xS ux sub def
 436  /yi yS uy sub def
 437  /xi2 xS ux 2 div sub def
 438  /yi2 yS uy 2 div sub def}%
 439 \pnode(!xi2 yi2){Mi2}
 440 \pnode(!xi yi){Mi}
 441 \pnode(!xS yS){S}
 442 \pscircle*(S){0.1}
 443 \psline[linestyle=dotted](S)(0,0)
 444 \psline[style=vecteurC]{->}(S)(Mi)
 445 \uput[l](S){S}
 446 \uput[u](Mi2){$\overrightarrow{F}$}
 447 \psarc{->}(0,0){1}{0}{!theta0}
 448 \uput{1.1}[!theta0 2 div](0,0){$\theta_0$}
 449 \psline[arrowinset=0.1,arrowsize=0.2]{<->}(5,0)(0,0)(0,5)
 450 \psline[arrowinset=0.05,arrowsize=0.1]{<->}(5,0)(0,0)(0,5)
 451 \uput[u](0,5){$y$}\uput[r](5,0){$x$}
 452 \end{pspicture}
 453  \end{center}
 454 Il s'obtient très facilement à partir de l'étude précédente si on sait que dans ce cas :
 455 \[
 456 v_0=\sqrt{\frac{\mathcal{G}M}{r_0}}
 457 \]
 458 \begin{verbatim}
 459 \pstVerb{
 460     /GM 1 def
 461     /theta0 60 def
 462     /r0 5 def
 463     /x0 r0 theta0 cos mul def
 464     /y0 r0 theta0 sin mul def
 465     /v0 GM r0 div sqrt def
 466     /v0x v0 theta0 sin mul neg def
 467     /v0y v0 theta0 cos mul def }%
 468 \end{verbatim}
 469 \section{L'hodographe du mouvement du satellite}
 470 On rappelle qu'en coordonnées polaires le vecteur-vitesse s'écrit :
 471 \[
 472 \overrightarrow{v}=\dot{r}\overrightarrow{u_r}+r\dot{\theta}\overrightarrow{u_{\theta}}
 473 \]
 474 \[
 475 \dot{\theta}=\frac{r_0v_0}{r^2}
 476 \]
 477 \[
 478 \dot{r}=-\frac{p(-\mathrm{e}\dot{\theta}\sin(\theta-\theta_0)}{(1+\mathrm{e}\cos(\theta-\theta_0))^2}=\frac{r_0v_0}{p}\mathrm{e}\sin(\theta-\theta_0)
 479 \]
 480 On peut exprimer $v$ uniquement en fonction de $\theta$ :
 481 \[
 482 \overrightarrow{v}=\frac{r_0v_0}{p}\mathrm{e}\sin(\theta-\theta_0)\overrightarrow{u_r}+\frac{r_0v_0}{p}\Big(1+\mathrm{e}\cos(\theta-\theta_0)\Big)\overrightarrow{u_\theta}
 483 \]
 484 Ses composantes dans la base $(\overrightarrow{u_r},u_\theta)$ sont :
 485 \[
 486 \left\{
 487 \begin{array}{rcl}
 488 v_r&=&\dfrac{r_0v_0}{p}\mathrm{e}\sin(\theta-\theta_0)\\[1em]
 489 v_\theta&=&\dfrac{r_0v_0}{p}\Big(1+\mathrm{e}\cos(\theta-\theta_0)\Big)
 490 \end{array}
 491 \right.
 492 \label{eqv1}
 493 \]
 494 On repasse aux coordonnées cartésiennes par une rotation d'angle $(-\theta)$.
 495 \[
 496 \left\{
 497 \begin{array}{rcl}
 498 \dot{x}&=&v_r\cos\theta-v_\theta\sin\theta\\[1em]
 499 \dot{y}&=&v_r\sin\theta+v_\theta\cos\theta
 500 \end{array}
 501 \right.
 502 \label{eqv2}
 503 \]
 504 \[
 505 \left\{
 506 \begin{array}{rcl}
 507 \dfrac{p}{r_0v_0}\dot{x}&=&\mathrm{e}\sin(\theta-\theta_0)\cos\theta-\Big(1+\mathrm{e}\cos(\theta-\theta_0)\Big)\sin\theta\\[1em]
 508 \dfrac{p}{r_0v_0}\dot{y}&=&\mathrm{e}\sin(\theta-\theta_0)\sin\theta+\Big(1+\mathrm{e}\cos(\theta-\theta_0)\Big)\cos\theta
 509 \end{array}
 510 \right.
 511 \label{eqv3}
 512 \]
 513 En utilisant les relations trigonométriques de soustraction :
 514 \[
 515 \sin\alpha\cos\beta-\cos\alpha\sin\beta=\sin(\alpha-\beta)\qquad \cos\alpha\cos\beta+\sin\alpha\sin\beta=\cos(\alpha-\beta)
 516 \]
 517 On obtient :
 518 \[
 519 \left\{
 520 \begin{array}{rcl}
 521 \dfrac{p}{r_0v_0}\dot{x}&=&-\mathrm{e}\sin\theta_0-\sin\theta\\[1em]
 522 \dfrac{p}{r_0v_0}\dot{y}&=&\hphantom{-}\mathrm{e}\cos\theta_0+\cos\theta
 523 \end{array}
 524 \right.
 525 \label{eqv4}
 526 \]
 527 \[
 528 \left\{
 529 \begin{array}{rcl}
 530 \dot{x}+\dfrac{r_0v_0}{p}\mathrm{e}\sin\theta_0&=&-\dfrac{r_0v_0}{p}\sin\theta\\[1em]
 531 \dot{y}-\dfrac{r_0v_0}{p}\mathrm{e}\cos\theta_0&=&\dfrac{r_0v_0}{p}\cos\theta
 532 \end{array}
 533 \right.
 534 \label{eqv5}
 535 \]
 536 L'équation de l'hodographe :
 537 \[
 538 \left(\dot{x}+\dfrac{r_0v_0}{p}\mathrm{e}\sin\theta_0\right)^2+\left(\dot{y}-\dfrac{r_0v_0}{p}\mathrm{e}\cos\theta_0\right)^2=\left(\dfrac{r_0v_0}{p}\right)^2
 539 \]
 540 est celle d'un cercle centré en $\Big(\dfrac{r_0v_0}{p}\mathrm{e}\sin\theta_0,\dfrac{r_0v_0}{p}\mathrm{e}\cos\theta_0\Big)$, de rayon $R=\dfrac{r_0v_0}{p}$. Ce qu'on vérifie sur le graphe suivant : en rouge l'hodographe a été obtenu à partir de l'équation différentielle et en noir avec un trait plus fin par l'expression exacte de l'équation du cercle.
 541 \begin{center}
 542 \begin{pspicture}(-12,-2)(4,8)
 543 \psset{unit=2}%
 544 \uput[r](0,4){$y$}\uput[u](2,0){$x$}
 545 \pstVerb{/GM 1 def
 546     /theta0 -35 def
 547     /r0 1 def
 548     /x0 r0 theta0 cos mul def
 549     /y0 r0 theta0 sin mul def
 550     /v0 1.3 def
 551     /v0x v0 theta0 sin mul neg def
 552     /v0y v0 theta0 cos mul def
 553     /Lc r0 v0 mul def % moment cinetique
 554     /par Lc dup mul GM div def % paramètre de l'ellipse
 555 % excentricité
 556     /exc 1 0.5 v0 4 exp mul r0 dup mul mul GM r0 mul v0 dup mul mul sub GM dup mul div 2 mul add sqrt def
 557 %%%%%%%%%%%%%%
 558     /a_2 par 1 exc dup mul sub div def % demi-grand axe
 559     /b_2 par 1 exc dup mul sub sqrt div def % demi-petit axe
 560     /periode 2 3.1416 dup mul a_2 3 exp mul GM div sqrt mul def
 561 % rayon de l'hodographe
 562    /rH r0 v0 mul par div def
 563 % les coordonnées du centre
 564    /xCH rH exc mul theta0 sin mul neg def
 565    /yCH rH exc mul theta0 cos mul def}%
 566 \psequadiff[method=rk4,
 567             plotpoints=1000,
 568             algebraic,
 569             whichabs=0,
 570             whichord=1,
 571             tabname=XiYi,
 572 %            saveData,filename=XiYi.dat
 573 ]{0}{36.5}{x0 y0 v0x v0y}{\eqsatellite}%
 574 \listplot{XiYi aload pop}
 575 \psequadiff[method=rk4,
 576             plotpoints=1000,
 577             algebraic,
 578             whichabs=2,
 579             whichord=3,
 580             tabname=vxvy,
 581 %            saveData,filename=vxvy.dat
 582 ]{0}{36.5}{x0 y0 v0x v0y}{\eqsatellite}%
 583 \listplot[unit=1,linecolor=red,linewidth=0.075]{vxvy aload pop}
 584 % on dessine la vitesse un point sur 50
 585 \multido{\i=0+50}{40}{%
 586 \pstVerb{ /vX vxvy \i\space get def
 587           /vY vxvy  \i\space 1 add get def
 588           /xi XiYi \i\space get def
 589           /yi XiYi \i\space 1 add get def}%
 590 \rput(!xi yi){\psline[unit=1,linecolor=blue]{->}(!vX vY)}}
 591 \pscircle[fillcolor=gray!70,fillstyle=solid](0,0){0.2}
 592 \psgrid[subgriddiv=2,gridcolor=lightgray,gridlabels=8pt](-6,-2)(2,4)
 593 \psline[arrowinset=0.05,arrowsize=0.1]{<->}(2,0)(0,0)(0,4)
 594 \pscircle(!xCH yCH){!rH}
 595 \psdot[dotstyle=+](!xCH yCH)
 596 \end{pspicture}
 597 \end{center}
 598 \end{document}